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
26 #include "lwgeom_sfcgal.h"
27
28 static int SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_t type);
29 static POINTARRAY* ptarray_from_SFCGAL(const sfcgal_geometry_t* geom, int force3D);
30 static sfcgal_geometry_t* ptarray_to_SFCGAL(const POINTARRAY* pa, int type);
31
32
33
34 /* Return SFCGAL version string */
35 const char*
lwgeom_sfcgal_version()36 lwgeom_sfcgal_version()
37 {
38 const char *version = sfcgal_version();
39 return version;
40 }
41
42
43 /*
44 * Mapping between SFCGAL and PostGIS types
45 *
46 * Throw an error if type is unsupported
47 */
48 static int
SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_t type)49 SFCGAL_type_to_lwgeom_type(sfcgal_geometry_type_t type)
50 {
51 switch (type)
52 {
53 case SFCGAL_TYPE_POINT:
54 return POINTTYPE;
55
56 case SFCGAL_TYPE_LINESTRING:
57 return LINETYPE;
58
59 case SFCGAL_TYPE_POLYGON:
60 return POLYGONTYPE;
61
62 case SFCGAL_TYPE_MULTIPOINT:
63 return MULTIPOINTTYPE;
64
65 case SFCGAL_TYPE_MULTILINESTRING:
66 return MULTILINETYPE;
67
68 case SFCGAL_TYPE_MULTIPOLYGON:
69 return MULTIPOLYGONTYPE;
70
71 case SFCGAL_TYPE_MULTISOLID:
72 return COLLECTIONTYPE; /* Nota: PolyhedralSurface closed inside
73 aim is to use true solid type as soon
74 as available in OGC SFS */
75
76 case SFCGAL_TYPE_GEOMETRYCOLLECTION:
77 return COLLECTIONTYPE;
78
79 #if 0
80 case SFCGAL_TYPE_CIRCULARSTRING:
81 return CIRCSTRINGTYPE;
82
83 case SFCGAL_TYPE_COMPOUNDCURVE:
84 return COMPOUNDTYPE;
85
86 case SFCGAL_TYPE_CURVEPOLYGON:
87 return CURVEPOLYTYPE;
88
89 case SFCGAL_TYPE_MULTICURVE:
90 return MULTICURVETYPE;
91
92 case SFCGAL_TYPE_MULTISURFACE:
93 return MULTISURFACETYPE;
94 #endif
95
96 case SFCGAL_TYPE_POLYHEDRALSURFACE:
97 return POLYHEDRALSURFACETYPE;
98
99 case SFCGAL_TYPE_TRIANGULATEDSURFACE:
100 return TINTYPE;
101
102 case SFCGAL_TYPE_TRIANGLE:
103 return TRIANGLETYPE;
104
105 default:
106 lwerror("SFCGAL_type_to_lwgeom_type: Unknown Type");
107 return 0;
108 }
109 }
110
111
112 /*
113 * Return a PostGIS pointarray from a simple SFCGAL geometry:
114 * POINT, LINESTRING or TRIANGLE
115 *
116 * Throw an error on others types
117 */
118 static POINTARRAY*
ptarray_from_SFCGAL(const sfcgal_geometry_t * geom,int want3d)119 ptarray_from_SFCGAL(const sfcgal_geometry_t* geom, int want3d)
120 {
121 POINT4D point;
122 uint32_t i, npoints;
123 POINTARRAY* pa = NULL;
124
125 assert(geom);
126
127 switch (sfcgal_geometry_type_id(geom))
128 {
129 case SFCGAL_TYPE_POINT:
130 {
131 pa = ptarray_construct(want3d, 0, 1);
132 point.x = sfcgal_point_x(geom);
133 point.y = sfcgal_point_y(geom);
134
135 if (sfcgal_geometry_is_3d(geom))
136 point.z = sfcgal_point_z(geom);
137 else if (want3d)
138 point.z = 0.0;
139
140 ptarray_set_point4d(pa, 0, &point);
141 }
142 break;
143
144 case SFCGAL_TYPE_LINESTRING:
145 {
146 npoints = sfcgal_linestring_num_points(geom);
147 pa = ptarray_construct(want3d, 0, npoints);
148
149 for (i = 0; i < npoints; i++)
150 {
151 const sfcgal_geometry_t* pt = sfcgal_linestring_point_n(geom, i);
152 point.x = sfcgal_point_x(pt);
153 point.y = sfcgal_point_y(pt);
154
155 if (sfcgal_geometry_is_3d(geom))
156 point.z = sfcgal_point_z(pt);
157 else if (want3d)
158 point.z = 0.0;
159
160 ptarray_set_point4d(pa, i, &point);
161 }
162 }
163 break;
164
165 case SFCGAL_TYPE_TRIANGLE:
166 {
167 pa = ptarray_construct(want3d, 0, 4);
168
169 for (i = 0; i < 4; i++)
170 {
171 const sfcgal_geometry_t* pt = sfcgal_triangle_vertex(geom, (i%3));
172 point.x = sfcgal_point_x(pt);
173 point.y = sfcgal_point_y(pt);
174
175 if ( sfcgal_geometry_is_3d(geom))
176 point.z = sfcgal_point_z(pt);
177 else if (want3d)
178 point.z = 0.0;
179
180 ptarray_set_point4d(pa, i, &point);
181 }
182 }
183 break;
184
185 /* Other types should not be called directly ... */
186 default:
187 lwerror("ptarray_from_SFCGAL: Unknown Type");
188 break;
189 }
190 return pa;
191 }
192
193
194 /*
195 * Convert a PostGIS pointarray to SFCGAL structure
196 *
197 * Used for simple LWGEOM geometry POINT, LINESTRING, TRIANGLE
198 * and POLYGON rings
199 */
200 static sfcgal_geometry_t*
ptarray_to_SFCGAL(const POINTARRAY * pa,int type)201 ptarray_to_SFCGAL(const POINTARRAY* pa, int type)
202 {
203 POINT3DZ point;
204 int is_3d;
205 uint32_t i;
206
207 assert(pa);
208
209 is_3d = FLAGS_GET_Z(pa->flags) != 0;
210
211 switch (type)
212 {
213 case POINTTYPE:
214 {
215 getPoint3dz_p(pa, 0, &point);
216 if (is_3d) return sfcgal_point_create_from_xyz(point.x, point.y, point.z);
217 else return sfcgal_point_create_from_xy(point.x, point.y);
218 }
219 break;
220
221 case LINETYPE:
222 {
223 sfcgal_geometry_t* line = sfcgal_linestring_create();
224
225 for (i = 0; i < pa->npoints; i++)
226 {
227 getPoint3dz_p(pa, i, &point);
228 if (is_3d)
229 {
230 sfcgal_linestring_add_point(line,
231 sfcgal_point_create_from_xyz(point.x, point.y, point.z));
232 }
233 else
234 {
235 sfcgal_linestring_add_point(line,
236 sfcgal_point_create_from_xy(point.x, point.y));
237 }
238 }
239
240 return line;
241 }
242 break;
243
244 case TRIANGLETYPE:
245 {
246 sfcgal_geometry_t* triangle = sfcgal_triangle_create();
247
248 getPoint3dz_p(pa, 0, &point);
249 if (is_3d) sfcgal_triangle_set_vertex_from_xyz(triangle, 0, point.x, point.y, point.z);
250 else sfcgal_triangle_set_vertex_from_xy (triangle, 0, point.x, point.y);
251
252 getPoint3dz_p(pa, 1, &point);
253 if (is_3d) sfcgal_triangle_set_vertex_from_xyz(triangle, 1, point.x, point.y, point.z);
254 else sfcgal_triangle_set_vertex_from_xy (triangle, 1, point.x, point.y);
255
256 getPoint3dz_p(pa, 2, &point);
257 if (is_3d) sfcgal_triangle_set_vertex_from_xyz(triangle, 2, point.x, point.y, point.z);
258 else sfcgal_triangle_set_vertex_from_xy (triangle, 2, point.x, point.y);
259
260 return triangle;
261 }
262 break;
263
264 /* Other SFCGAL types should not be called directly ... */
265 default:
266 lwerror("ptarray_from_SFCGAL: Unknown Type");
267 return NULL;
268 }
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,int srid)278 SFCGAL2LWGEOM(const sfcgal_geometry_t* geom, int force3D, int 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 ) ) ++nsolids;
349 }
350 geoms = (LWGEOM**) lwrealloc(geoms, sizeof(LWGEOM*) * ngeoms);
351 }
352 LWGEOM* rgeom = (LWGEOM*) lwcollection_construct(SFCGAL_type_to_lwgeom_type(
353 sfcgal_geometry_type_id(geom)), srid, NULL, ngeoms, geoms);
354 if ( ngeoms )
355 {
356 if ( ngeoms == nsolids ) FLAGS_SET_SOLID( rgeom->flags, 1);
357 else if ( nsolids ) lwnotice("SFCGAL2LWGEOM: SOLID in heterogeneous collection will be treated as a POLYHEDRALSURFACETYPE");
358 }
359 return rgeom;
360 }
361
362 #if 0
363 case SFCGAL_TYPE_CIRCULARSTRING:
364 case SFCGAL_TYPE_COMPOUNDCURVE:
365 case SFCGAL_TYPE_CURVEPOLYGON:
366 case SFCGAL_TYPE_MULTICURVE:
367 case SFCGAL_TYPE_MULTISURFACE:
368 case SFCGAL_TYPE_CURVE:
369 case SFCGAL_TYPE_SURFACE:
370
371 /* TODO curve types handling */
372 #endif
373
374 case SFCGAL_TYPE_POLYHEDRALSURFACE:
375 {
376 ngeoms = sfcgal_polyhedral_surface_num_polygons(geom);
377
378 LWGEOM** geoms = NULL;
379 if (ngeoms)
380 {
381 geoms = (LWGEOM**) lwalloc(sizeof(LWGEOM*) * ngeoms);
382 for (i = 0; i < ngeoms; i++)
383 {
384 const sfcgal_geometry_t* g = sfcgal_polyhedral_surface_polygon_n( geom, i );
385 geoms[i] = SFCGAL2LWGEOM(g, 0, srid);
386 }
387 }
388 return (LWGEOM*)lwcollection_construct(POLYHEDRALSURFACETYPE, srid, NULL, ngeoms, geoms);
389 }
390
391 /* Solid is map as a closed PolyhedralSurface (for now) */
392 case SFCGAL_TYPE_SOLID:
393 {
394 nshells = sfcgal_solid_num_shells(geom);
395
396 for (ngeoms = 0, i = 0; i < nshells; i++)
397 ngeoms += sfcgal_polyhedral_surface_num_polygons(sfcgal_solid_shell_n(geom, i));
398
399 LWGEOM** geoms = 0;
400 if (ngeoms)
401 {
402 geoms = (LWGEOM**) lwalloc( sizeof(LWGEOM*) * ngeoms);
403 for (i = 0, k =0 ; i < nshells; i++)
404 {
405 const sfcgal_geometry_t* shell = sfcgal_solid_shell_n(geom, i);
406 ngeoms = sfcgal_polyhedral_surface_num_polygons(shell);
407
408 for (j = 0; j < ngeoms; j++)
409 {
410 const sfcgal_geometry_t* g = sfcgal_polyhedral_surface_polygon_n(shell, j);
411 geoms[k] = SFCGAL2LWGEOM(g, 1, srid);
412 k++;
413 }
414 }
415 }
416 LWGEOM* rgeom = (LWGEOM*) lwcollection_construct(POLYHEDRALSURFACETYPE, srid, NULL, ngeoms, geoms);
417 if (ngeoms) FLAGS_SET_SOLID( rgeom->flags, 1);
418 return rgeom;
419 }
420
421 case SFCGAL_TYPE_TRIANGULATEDSURFACE:
422 {
423 ngeoms = sfcgal_triangulated_surface_num_triangles(geom);
424 LWGEOM** geoms = NULL;
425 if (ngeoms)
426 {
427 geoms = (LWGEOM**) lwalloc(sizeof(LWGEOM*) * ngeoms);
428 for (i = 0; i < ngeoms; i++)
429 {
430 const sfcgal_geometry_t* g = sfcgal_triangulated_surface_triangle_n(geom, i);
431 geoms[i] = SFCGAL2LWGEOM(g, 0, srid);
432 }
433 }
434 return (LWGEOM*) lwcollection_construct(TINTYPE, srid, NULL, ngeoms, geoms);
435 }
436
437 default:
438 lwerror("SFCGAL2LWGEOM: Unknown Type");
439 return NULL;
440 }
441 }
442
443
444 sfcgal_geometry_t*
LWGEOM2SFCGAL(const LWGEOM * geom)445 LWGEOM2SFCGAL(const LWGEOM* geom)
446 {
447 uint32_t i;
448 sfcgal_geometry_t* ret_geom = NULL;
449
450 assert(geom);
451
452 switch (geom->type)
453 {
454 case POINTTYPE:
455 {
456 const LWPOINT* lwp = (const LWPOINT*) geom;
457 if (lwgeom_is_empty(geom)) return sfcgal_point_create();
458
459 return ptarray_to_SFCGAL(lwp->point, POINTTYPE);
460 }
461 break;
462
463 case LINETYPE:
464 {
465 const LWLINE* line = (const LWLINE*) geom;
466 if (lwgeom_is_empty(geom)) return sfcgal_linestring_create();
467
468 return ptarray_to_SFCGAL(line->points, LINETYPE);
469 }
470 break;
471
472 case TRIANGLETYPE:
473 {
474 const LWTRIANGLE* triangle = (const LWTRIANGLE*) geom;
475 if (lwgeom_is_empty(geom)) return sfcgal_triangle_create();
476 return ptarray_to_SFCGAL(triangle->points, TRIANGLETYPE);
477 }
478 break;
479
480 case POLYGONTYPE:
481 {
482 const LWPOLY* poly = (const LWPOLY*) geom;
483 uint32_t nrings = poly->nrings - 1;
484
485 if (lwgeom_is_empty(geom)) return sfcgal_polygon_create();
486
487 sfcgal_geometry_t* exterior_ring = ptarray_to_SFCGAL(poly->rings[0], LINETYPE);
488 ret_geom = sfcgal_polygon_create_from_exterior_ring(exterior_ring);
489
490 for (i = 0; i < nrings; i++)
491 {
492 sfcgal_geometry_t* ring = ptarray_to_SFCGAL(poly->rings[i + 1], LINETYPE);
493 sfcgal_polygon_add_interior_ring(ret_geom, ring);
494 }
495 return ret_geom;
496 }
497 break;
498
499 case MULTIPOINTTYPE:
500 case MULTILINETYPE:
501 case MULTIPOLYGONTYPE:
502 case COLLECTIONTYPE:
503 {
504 if (geom->type == MULTIPOINTTYPE) ret_geom = sfcgal_multi_point_create();
505 else if (geom->type == MULTILINETYPE) ret_geom = sfcgal_multi_linestring_create();
506 else if (geom->type == MULTIPOLYGONTYPE) ret_geom = sfcgal_multi_polygon_create();
507 else ret_geom = sfcgal_geometry_collection_create();
508
509 const LWCOLLECTION* lwc = (const LWCOLLECTION*)geom;
510 for (i = 0; i < lwc->ngeoms; i++)
511 {
512 sfcgal_geometry_t* g = LWGEOM2SFCGAL(lwc->geoms[i]);
513 sfcgal_geometry_collection_add_geometry(ret_geom, g);
514 }
515
516 return ret_geom;
517 }
518 break;
519
520 case POLYHEDRALSURFACETYPE:
521 {
522 const LWPSURFACE* lwp = (const LWPSURFACE*) geom;
523 ret_geom = sfcgal_polyhedral_surface_create();
524
525 for (i = 0; i < lwp->ngeoms; i++)
526 {
527 sfcgal_geometry_t* g = LWGEOM2SFCGAL((const LWGEOM*) lwp->geoms[i]);
528 sfcgal_polyhedral_surface_add_polygon(ret_geom, g);
529 }
530 /* We treat polyhedral surface as the only exterior shell,
531 since we can't distinguish exterior from interior shells ... */
532 if (FLAGS_GET_SOLID(lwp->flags))
533 {
534 return sfcgal_solid_create_from_exterior_shell(ret_geom);
535 }
536
537 return ret_geom;
538 }
539 break;
540
541 case TINTYPE:
542 {
543 const LWTIN* lwp = (const LWTIN*) geom;
544 ret_geom = sfcgal_triangulated_surface_create();
545
546 for (i = 0; i < lwp->ngeoms; i++)
547 {
548 sfcgal_geometry_t* g = LWGEOM2SFCGAL((const LWGEOM*) lwp->geoms[i]);
549 sfcgal_triangulated_surface_add_triangle(ret_geom, g);
550 }
551
552 return ret_geom;
553 }
554 break;
555
556 default:
557 lwerror("LWGEOM2SFCGAL: Unknown geometry type !");
558 return NULL;
559 }
560 }
561
562
563 /*
564 * No Operation SFCGAL function, used (only) for cunit tests
565 * Take a PostGIS geometry, send it to SFCGAL and return it unchanged (in theory)
566 */
lwgeom_sfcgal_noop(const LWGEOM * geom_in)567 LWGEOM* lwgeom_sfcgal_noop(const LWGEOM* geom_in)
568 {
569 sfcgal_geometry_t* converted;
570
571 assert(geom_in);
572
573 converted = LWGEOM2SFCGAL(geom_in);
574 assert(converted);
575
576 LWGEOM* geom_out = SFCGAL2LWGEOM(converted, 0, SRID_UNKNOWN);
577 sfcgal_geometry_delete(converted);
578
579 /* copy SRID (SFCGAL does not store the SRID) */
580 geom_out->srid = geom_in->srid;
581 return geom_out;
582 }
583