1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2012-2013 Oslandia <infos@oslandia.com>
22 *
23 **********************************************************************/
24
25 #include "lwgeom_sfcgal.h"
26
27 static int SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_t type);
28 static POINTARRAY *ptarray_from_SFCGAL(const sfcgal_geometry_t *geom, int force3D);
29 static sfcgal_geometry_t *ptarray_to_SFCGAL(const POINTARRAY *pa, int type);
30
31 /* Return SFCGAL version string */
32 const char *
lwgeom_sfcgal_version()33 lwgeom_sfcgal_version()
34 {
35 const char *version = sfcgal_version();
36 return version;
37 }
38
39 /*
40 * Mapping between SFCGAL and PostGIS types
41 *
42 * Throw an error if type is unsupported
43 */
44 static int
SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_t type)45 SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_t type)
46 {
47 switch (type)
48 {
49 case SFCGAL_TYPE_POINT:
50 return POINTTYPE;
51
52 case SFCGAL_TYPE_LINESTRING:
53 return LINETYPE;
54
55 case SFCGAL_TYPE_POLYGON:
56 return POLYGONTYPE;
57
58 case SFCGAL_TYPE_MULTIPOINT:
59 return MULTIPOINTTYPE;
60
61 case SFCGAL_TYPE_MULTILINESTRING:
62 return MULTILINETYPE;
63
64 case SFCGAL_TYPE_MULTIPOLYGON:
65 return MULTIPOLYGONTYPE;
66
67 case SFCGAL_TYPE_MULTISOLID:
68 return COLLECTIONTYPE; /* Nota: PolyhedralSurface closed inside
69 aim is to use true solid type as soon
70 as available in OGC SFS */
71
72 case SFCGAL_TYPE_GEOMETRYCOLLECTION:
73 return COLLECTIONTYPE;
74
75 #if 0
76 case SFCGAL_TYPE_CIRCULARSTRING:
77 return CIRCSTRINGTYPE;
78
79 case SFCGAL_TYPE_COMPOUNDCURVE:
80 return COMPOUNDTYPE;
81
82 case SFCGAL_TYPE_CURVEPOLYGON:
83 return CURVEPOLYTYPE;
84
85 case SFCGAL_TYPE_MULTICURVE:
86 return MULTICURVETYPE;
87
88 case SFCGAL_TYPE_MULTISURFACE:
89 return MULTISURFACETYPE;
90 #endif
91
92 case SFCGAL_TYPE_POLYHEDRALSURFACE:
93 return POLYHEDRALSURFACETYPE;
94
95 case SFCGAL_TYPE_TRIANGULATEDSURFACE:
96 return TINTYPE;
97
98 case SFCGAL_TYPE_TRIANGLE:
99 return TRIANGLETYPE;
100
101 default:
102 lwerror("SFCGAL_type_to_lwgeom_type: Unknown Type");
103 return 0;
104 }
105 }
106
107 /*
108 * Return a PostGIS pointarray from a simple SFCGAL geometry:
109 * POINT, LINESTRING or TRIANGLE
110 *
111 * Throw an error on others types
112 */
113 static POINTARRAY *
ptarray_from_SFCGAL(const sfcgal_geometry_t * geom,int want3d)114 ptarray_from_SFCGAL(const sfcgal_geometry_t *geom, int want3d)
115 {
116 POINT4D point;
117 uint32_t i, npoints;
118 POINTARRAY *pa = NULL;
119
120 assert(geom);
121
122 switch (sfcgal_geometry_type_id(geom))
123 {
124 case SFCGAL_TYPE_POINT:
125 {
126 pa = ptarray_construct(want3d, 0, 1);
127 point.x = sfcgal_point_x(geom);
128 point.y = sfcgal_point_y(geom);
129
130 if (sfcgal_geometry_is_3d(geom))
131 point.z = sfcgal_point_z(geom);
132 else if (want3d)
133 point.z = 0.0;
134
135 ptarray_set_point4d(pa, 0, &point);
136 }
137 break;
138
139 case SFCGAL_TYPE_LINESTRING:
140 {
141 npoints = sfcgal_linestring_num_points(geom);
142 pa = ptarray_construct(want3d, 0, npoints);
143
144 for (i = 0; i < npoints; i++)
145 {
146 const sfcgal_geometry_t *pt = sfcgal_linestring_point_n(geom, i);
147 point.x = sfcgal_point_x(pt);
148 point.y = sfcgal_point_y(pt);
149
150 if (sfcgal_geometry_is_3d(geom))
151 point.z = sfcgal_point_z(pt);
152 else if (want3d)
153 point.z = 0.0;
154
155 ptarray_set_point4d(pa, i, &point);
156 }
157 }
158 break;
159
160 case SFCGAL_TYPE_TRIANGLE:
161 {
162 pa = ptarray_construct(want3d, 0, 4);
163
164 for (i = 0; i < 4; i++)
165 {
166 const sfcgal_geometry_t *pt = sfcgal_triangle_vertex(geom, (i % 3));
167 point.x = sfcgal_point_x(pt);
168 point.y = sfcgal_point_y(pt);
169
170 if (sfcgal_geometry_is_3d(geom))
171 point.z = sfcgal_point_z(pt);
172 else if (want3d)
173 point.z = 0.0;
174
175 ptarray_set_point4d(pa, i, &point);
176 }
177 }
178 break;
179
180 /* Other types should not be called directly ... */
181 default:
182 lwerror("ptarray_from_SFCGAL: Unknown Type");
183 break;
184 }
185 return pa;
186 }
187
188 /*
189 * Convert a PostGIS pointarray to SFCGAL structure
190 *
191 * Used for simple LWGEOM geometry POINT, LINESTRING, TRIANGLE
192 * and POLYGON rings
193 */
194 static sfcgal_geometry_t *
ptarray_to_SFCGAL(const POINTARRAY * pa,int type)195 ptarray_to_SFCGAL(const POINTARRAY *pa, int type)
196 {
197 POINT3DZ point;
198 int is_3d;
199 uint32_t i;
200
201 assert(pa);
202
203 is_3d = FLAGS_GET_Z(pa->flags) != 0;
204
205 switch (type)
206 {
207 case POINTTYPE:
208 {
209 getPoint3dz_p(pa, 0, &point);
210 if (is_3d)
211 return sfcgal_point_create_from_xyz(point.x, point.y, point.z);
212 else
213 return sfcgal_point_create_from_xy(point.x, point.y);
214 }
215 break;
216
217 case LINETYPE:
218 {
219 sfcgal_geometry_t *line = sfcgal_linestring_create();
220
221 for (i = 0; i < pa->npoints; i++)
222 {
223 getPoint3dz_p(pa, i, &point);
224 if (is_3d)
225 {
226 sfcgal_linestring_add_point(line,
227 sfcgal_point_create_from_xyz(point.x, point.y, point.z));
228 }
229 else
230 {
231 sfcgal_linestring_add_point(line, sfcgal_point_create_from_xy(point.x, point.y));
232 }
233 }
234
235 return line;
236 }
237 break;
238
239 case TRIANGLETYPE:
240 {
241 sfcgal_geometry_t *triangle = sfcgal_triangle_create();
242
243 getPoint3dz_p(pa, 0, &point);
244 if (is_3d)
245 sfcgal_triangle_set_vertex_from_xyz(triangle, 0, point.x, point.y, point.z);
246 else
247 sfcgal_triangle_set_vertex_from_xy(triangle, 0, point.x, point.y);
248
249 getPoint3dz_p(pa, 1, &point);
250 if (is_3d)
251 sfcgal_triangle_set_vertex_from_xyz(triangle, 1, point.x, point.y, point.z);
252 else
253 sfcgal_triangle_set_vertex_from_xy(triangle, 1, point.x, point.y);
254
255 getPoint3dz_p(pa, 2, &point);
256 if (is_3d)
257 sfcgal_triangle_set_vertex_from_xyz(triangle, 2, point.x, point.y, point.z);
258 else
259 sfcgal_triangle_set_vertex_from_xy(triangle, 2, point.x, point.y);
260
261 return triangle;
262 }
263 break;
264
265 /* Other SFCGAL types should not be called directly ... */
266 default:
267 lwerror("ptarray_from_SFCGAL: Unknown Type");
268 return NULL;
269 }
270 }
271
272 /*
273 * Convert a SFCGAL structure to PostGIS LWGEOM
274 *
275 * Throws an error on unsupported type
276 */
277 LWGEOM *
SFCGAL2LWGEOM(const sfcgal_geometry_t * geom,int force3D,int32_t srid)278 SFCGAL2LWGEOM(const sfcgal_geometry_t *geom, int force3D, int32_t srid)
279 {
280 uint32_t ngeoms, nshells, nsolids;
281 uint32_t i, j, k;
282 int want3d;
283
284 assert(geom);
285
286 want3d = force3D || sfcgal_geometry_is_3d(geom);
287
288 switch (sfcgal_geometry_type_id(geom))
289 {
290 case SFCGAL_TYPE_POINT:
291 {
292 if (sfcgal_geometry_is_empty(geom))
293 return (LWGEOM *)lwpoint_construct_empty(srid, want3d, 0);
294
295 POINTARRAY *pa = ptarray_from_SFCGAL(geom, want3d);
296 return (LWGEOM *)lwpoint_construct(srid, NULL, pa);
297 }
298
299 case SFCGAL_TYPE_LINESTRING:
300 {
301 if (sfcgal_geometry_is_empty(geom))
302 return (LWGEOM *)lwline_construct_empty(srid, want3d, 0);
303
304 POINTARRAY *pa = ptarray_from_SFCGAL(geom, want3d);
305 return (LWGEOM *)lwline_construct(srid, NULL, pa);
306 }
307
308 case SFCGAL_TYPE_TRIANGLE:
309 {
310 if (sfcgal_geometry_is_empty(geom))
311 return (LWGEOM *)lwtriangle_construct_empty(srid, want3d, 0);
312
313 POINTARRAY *pa = ptarray_from_SFCGAL(geom, want3d);
314 return (LWGEOM *)lwtriangle_construct(srid, NULL, pa);
315 }
316
317 case SFCGAL_TYPE_POLYGON:
318 {
319 if (sfcgal_geometry_is_empty(geom))
320 return (LWGEOM *)lwpoly_construct_empty(srid, want3d, 0);
321
322 uint32_t nrings = sfcgal_polygon_num_interior_rings(geom) + 1;
323 POINTARRAY **pa = (POINTARRAY **)lwalloc(sizeof(POINTARRAY *) * nrings);
324
325 pa[0] = ptarray_from_SFCGAL(sfcgal_polygon_exterior_ring(geom), want3d);
326 for (i = 1; i < nrings; i++)
327 pa[i] = ptarray_from_SFCGAL(sfcgal_polygon_interior_ring_n(geom, i - 1), want3d);
328
329 return (LWGEOM *)lwpoly_construct(srid, NULL, nrings, pa);
330 }
331
332 case SFCGAL_TYPE_MULTIPOINT:
333 case SFCGAL_TYPE_MULTILINESTRING:
334 case SFCGAL_TYPE_MULTIPOLYGON:
335 case SFCGAL_TYPE_MULTISOLID:
336 case SFCGAL_TYPE_GEOMETRYCOLLECTION:
337 {
338 ngeoms = sfcgal_geometry_collection_num_geometries(geom);
339 LWGEOM **geoms = NULL;
340 if (ngeoms)
341 {
342 nsolids = 0;
343 geoms = (LWGEOM **)lwalloc(sizeof(LWGEOM *) * ngeoms);
344 for (i = 0; i < ngeoms; i++)
345 {
346 const sfcgal_geometry_t *g = sfcgal_geometry_collection_geometry_n(geom, i);
347 geoms[i] = SFCGAL2LWGEOM(g, 0, srid);
348 if (FLAGS_GET_SOLID(geoms[i]->flags))
349 ++nsolids;
350 }
351 geoms = (LWGEOM **)lwrealloc(geoms, sizeof(LWGEOM *) * ngeoms);
352 }
353 LWGEOM *rgeom = (LWGEOM *)lwcollection_construct(
354 SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_id(geom)), srid, NULL, ngeoms, geoms);
355 if (ngeoms)
356 {
357 if (ngeoms == nsolids)
358 FLAGS_SET_SOLID(rgeom->flags, 1);
359 else if (nsolids)
360 lwnotice(
361 "SFCGAL2LWGEOM: SOLID in heterogeneous collection will be treated as a POLYHEDRALSURFACETYPE");
362 }
363 return rgeom;
364 }
365
366 #if 0
367 case SFCGAL_TYPE_CIRCULARSTRING:
368 case SFCGAL_TYPE_COMPOUNDCURVE:
369 case SFCGAL_TYPE_CURVEPOLYGON:
370 case SFCGAL_TYPE_MULTICURVE:
371 case SFCGAL_TYPE_MULTISURFACE:
372 case SFCGAL_TYPE_CURVE:
373 case SFCGAL_TYPE_SURFACE:
374
375 /* TODO curve types handling */
376 #endif
377
378 case SFCGAL_TYPE_POLYHEDRALSURFACE:
379 {
380 ngeoms = sfcgal_polyhedral_surface_num_polygons(geom);
381
382 LWGEOM **geoms = NULL;
383 if (ngeoms)
384 {
385 geoms = (LWGEOM **)lwalloc(sizeof(LWGEOM *) * ngeoms);
386 for (i = 0; i < ngeoms; i++)
387 {
388 const sfcgal_geometry_t *g = sfcgal_polyhedral_surface_polygon_n(geom, i);
389 geoms[i] = SFCGAL2LWGEOM(g, 0, srid);
390 }
391 }
392 return (LWGEOM *)lwcollection_construct(POLYHEDRALSURFACETYPE, srid, NULL, ngeoms, geoms);
393 }
394
395 /* Solid is map as a closed PolyhedralSurface (for now) */
396 case SFCGAL_TYPE_SOLID:
397 {
398 nshells = sfcgal_solid_num_shells(geom);
399
400 for (ngeoms = 0, i = 0; i < nshells; i++)
401 ngeoms += sfcgal_polyhedral_surface_num_polygons(sfcgal_solid_shell_n(geom, i));
402
403 LWGEOM **geoms = 0;
404 if (ngeoms)
405 {
406 geoms = (LWGEOM **)lwalloc(sizeof(LWGEOM *) * ngeoms);
407 for (i = 0, k = 0; i < nshells; i++)
408 {
409 const sfcgal_geometry_t *shell = sfcgal_solid_shell_n(geom, i);
410 ngeoms = sfcgal_polyhedral_surface_num_polygons(shell);
411
412 for (j = 0; j < ngeoms; j++)
413 {
414 const sfcgal_geometry_t *g = sfcgal_polyhedral_surface_polygon_n(shell, j);
415 geoms[k] = SFCGAL2LWGEOM(g, 1, srid);
416 k++;
417 }
418 }
419 }
420 LWGEOM *rgeom = (LWGEOM *)lwcollection_construct(POLYHEDRALSURFACETYPE, srid, NULL, ngeoms, geoms);
421 if (ngeoms)
422 FLAGS_SET_SOLID(rgeom->flags, 1);
423 return rgeom;
424 }
425
426 case SFCGAL_TYPE_TRIANGULATEDSURFACE:
427 {
428 ngeoms = sfcgal_triangulated_surface_num_triangles(geom);
429 LWGEOM **geoms = NULL;
430 if (ngeoms)
431 {
432 geoms = (LWGEOM **)lwalloc(sizeof(LWGEOM *) * ngeoms);
433 for (i = 0; i < ngeoms; i++)
434 {
435 const sfcgal_geometry_t *g = sfcgal_triangulated_surface_triangle_n(geom, i);
436 geoms[i] = SFCGAL2LWGEOM(g, 0, srid);
437 }
438 }
439 return (LWGEOM *)lwcollection_construct(TINTYPE, srid, NULL, ngeoms, geoms);
440 }
441
442 default:
443 lwerror("SFCGAL2LWGEOM: Unknown Type");
444 return NULL;
445 }
446 }
447
448 sfcgal_geometry_t *
LWGEOM2SFCGAL(const LWGEOM * geom)449 LWGEOM2SFCGAL(const LWGEOM *geom)
450 {
451 uint32_t i;
452 sfcgal_geometry_t *ret_geom = NULL;
453
454 assert(geom);
455
456 switch (geom->type)
457 {
458 case POINTTYPE:
459 {
460 const LWPOINT *lwp = (const LWPOINT *)geom;
461 if (lwgeom_is_empty(geom))
462 return sfcgal_point_create();
463
464 return ptarray_to_SFCGAL(lwp->point, POINTTYPE);
465 }
466 break;
467
468 case LINETYPE:
469 {
470 const LWLINE *line = (const LWLINE *)geom;
471 if (lwgeom_is_empty(geom))
472 return sfcgal_linestring_create();
473
474 return ptarray_to_SFCGAL(line->points, LINETYPE);
475 }
476 break;
477
478 case TRIANGLETYPE:
479 {
480 const LWTRIANGLE *triangle = (const LWTRIANGLE *)geom;
481 if (lwgeom_is_empty(geom))
482 return sfcgal_triangle_create();
483 return ptarray_to_SFCGAL(triangle->points, TRIANGLETYPE);
484 }
485 break;
486
487 case POLYGONTYPE:
488 {
489 const LWPOLY *poly = (const LWPOLY *)geom;
490 uint32_t nrings = poly->nrings - 1;
491
492 if (lwgeom_is_empty(geom))
493 return sfcgal_polygon_create();
494
495 sfcgal_geometry_t *exterior_ring = ptarray_to_SFCGAL(poly->rings[0], LINETYPE);
496 ret_geom = sfcgal_polygon_create_from_exterior_ring(exterior_ring);
497
498 for (i = 0; i < nrings; i++)
499 {
500 sfcgal_geometry_t *ring = ptarray_to_SFCGAL(poly->rings[i + 1], LINETYPE);
501 sfcgal_polygon_add_interior_ring(ret_geom, ring);
502 }
503 return ret_geom;
504 }
505 break;
506
507 case MULTIPOINTTYPE:
508 case MULTILINETYPE:
509 case MULTIPOLYGONTYPE:
510 case COLLECTIONTYPE:
511 {
512 if (geom->type == MULTIPOINTTYPE)
513 ret_geom = sfcgal_multi_point_create();
514 else if (geom->type == MULTILINETYPE)
515 ret_geom = sfcgal_multi_linestring_create();
516 else if (geom->type == MULTIPOLYGONTYPE)
517 ret_geom = sfcgal_multi_polygon_create();
518 else
519 ret_geom = sfcgal_geometry_collection_create();
520
521 const LWCOLLECTION *lwc = (const LWCOLLECTION *)geom;
522 for (i = 0; i < lwc->ngeoms; i++)
523 {
524 sfcgal_geometry_t *g = LWGEOM2SFCGAL(lwc->geoms[i]);
525 sfcgal_geometry_collection_add_geometry(ret_geom, g);
526 }
527
528 return ret_geom;
529 }
530 break;
531
532 case POLYHEDRALSURFACETYPE:
533 {
534 const LWPSURFACE *lwp = (const LWPSURFACE *)geom;
535 ret_geom = sfcgal_polyhedral_surface_create();
536
537 for (i = 0; i < lwp->ngeoms; i++)
538 {
539 sfcgal_geometry_t *g = LWGEOM2SFCGAL((const LWGEOM *)lwp->geoms[i]);
540 sfcgal_polyhedral_surface_add_polygon(ret_geom, g);
541 }
542 /* We treat polyhedral surface as the only exterior shell,
543 since we can't distinguish exterior from interior shells ... */
544 if (FLAGS_GET_SOLID(lwp->flags))
545 {
546 return sfcgal_solid_create_from_exterior_shell(ret_geom);
547 }
548
549 return ret_geom;
550 }
551 break;
552
553 case TINTYPE:
554 {
555 const LWTIN *lwp = (const LWTIN *)geom;
556 ret_geom = sfcgal_triangulated_surface_create();
557
558 for (i = 0; i < lwp->ngeoms; i++)
559 {
560 sfcgal_geometry_t *g = LWGEOM2SFCGAL((const LWGEOM *)lwp->geoms[i]);
561 sfcgal_triangulated_surface_add_triangle(ret_geom, g);
562 }
563
564 return ret_geom;
565 }
566 break;
567
568 default:
569 lwerror("LWGEOM2SFCGAL: Unknown geometry type !");
570 return NULL;
571 }
572 }
573
574 /*
575 * No Operation SFCGAL function, used (only) for cunit tests
576 * Take a PostGIS geometry, send it to SFCGAL and return it unchanged (in theory)
577 */
578 LWGEOM *
lwgeom_sfcgal_noop(const LWGEOM * geom_in)579 lwgeom_sfcgal_noop(const LWGEOM *geom_in)
580 {
581 sfcgal_geometry_t *converted;
582
583 assert(geom_in);
584
585 converted = LWGEOM2SFCGAL(geom_in);
586 assert(converted);
587
588 LWGEOM *geom_out = SFCGAL2LWGEOM(converted, 0, SRID_UNKNOWN);
589 sfcgal_geometry_delete(converted);
590
591 /* copy SRID (SFCGAL does not store the SRID) */
592 geom_out->srid = geom_in->srid;
593 return geom_out;
594 }
595