1 /******************************************************************************
2  * $Id$
3  *
4  * Project:  MapServer
5  * Purpose:  MapServer-GEOS integration.
6  * Author:   Steve Lime and the MapServer team.
7  *
8  ******************************************************************************
9  * Copyright (c) 1996-2005 Regents of the University of Minnesota.
10  *
11  * Permission is hereby granted, free of charge, to any person obtaining a
12  * copy of this software and associated documentation files (the "Software"),
13  * to deal in the Software without restriction, including without limitation
14  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15  * and/or sell copies of the Software, and to permit persons to whom the
16  * Software is furnished to do so, subject to the following conditions:
17  *
18  * The above copyright notice and this permission notice shall be included in
19  * all copies of this Software or works derived from this Software.
20  *
21  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27  * DEALINGS IN THE SOFTWARE.
28  *****************************************************************************/
29 
30 #include "mapserver.h"
31 
32 #ifdef USE_GEOS
33 
34 // To avoid accidental use of non reentrant GEOS API.
35 // (check only effective in GEOS >= 3.5)
36 #define GEOS_USE_ONLY_R_API
37 
38 #include <geos_c.h>
39 
40 /*
41 ** Error handling...
42 */
msGEOSError(const char * format,...)43 static void msGEOSError(const char *format, ...)
44 {
45   va_list args;
46 
47   va_start (args, format);
48   msSetError(MS_GEOSERR, format, "msGEOSError()", args); /* just pass along to MapServer error handling */
49   va_end(args);
50 
51   return;
52 }
53 
msGEOSNotice(const char * fmt,...)54 static void msGEOSNotice(const char *fmt, ...)
55 {
56   return; /* do nothing with notices at this point */
57 }
58 
59 #ifndef USE_THREAD
60 static GEOSContextHandle_t geos_handle;
61 
msGetGeosContextHandle()62 static inline GEOSContextHandle_t msGetGeosContextHandle()
63 {
64   return geos_handle;
65 }
66 
67 #else
68 #include "mapthread.h"
69 typedef struct geos_thread_info {
70   struct geos_thread_info *next;
71   void*             thread_id;
72   GEOSContextHandle_t        geos_handle;
73 } geos_thread_info_t;
74 
75 static geos_thread_info_t *geos_list = NULL;
76 
msGetGeosContextHandle()77 static GEOSContextHandle_t msGetGeosContextHandle()
78 {
79   geos_thread_info_t *link;
80   GEOSContextHandle_t ret_obj;
81   void*        thread_id;
82 
83   msAcquireLock( TLOCK_GEOS );
84 
85   thread_id = msGetThreadId();
86 
87   /* find link for this thread */
88 
89   for( link = geos_list;
90        link != NULL && link->thread_id != thread_id
91        && link->next != NULL && link->next->thread_id != thread_id;
92        link = link->next ) {}
93 
94   /* If the target thread link is already at the head of the list were ok */
95   if( geos_list != NULL && geos_list->thread_id == thread_id ) {
96   }
97 
98   /* We don't have one ... initialize one. */
99   else if( link == NULL || link->next == NULL ) {
100     geos_thread_info_t *new_link;
101     new_link = (geos_thread_info_t *) malloc(sizeof(geos_thread_info_t));
102     new_link->next = geos_list;
103     new_link->thread_id = thread_id;
104     new_link->geos_handle = initGEOS_r(msGEOSNotice, msGEOSError);
105 
106     geos_list = new_link;
107   }
108 
109   /* If the link is not already at the head of the list, promote it */
110   else if( link != NULL && link->next != NULL ) {
111     geos_thread_info_t *target = link->next;
112 
113     link->next = link->next->next;
114     target->next = geos_list;
115     geos_list = target;
116   }
117 
118   ret_obj = geos_list->geos_handle;
119 
120   msReleaseLock( TLOCK_GEOS );
121 
122   return ret_obj;
123 }
124 #endif
125 
126 /*
127 ** Setup/Cleanup wrappers
128 */
msGEOSSetup()129 void msGEOSSetup()
130 {
131 #ifndef USE_THREAD
132   geos_handle = initGEOS_r(msGEOSNotice, msGEOSError);
133 #endif
134 }
135 
msGEOSCleanup()136 void msGEOSCleanup()
137 {
138 #ifndef USE_THREAD
139   finishGEOS_r(geos_handle);
140   geos_handle = NULL;
141 #else
142   geos_thread_info_t *link;
143   msAcquireLock( TLOCK_GEOS );
144   for( link = geos_list; link != NULL;) {
145     geos_thread_info_t *cur = link;
146     link = link->next;
147     finishGEOS_r(cur->geos_handle);
148     free(cur);
149   }
150   geos_list = NULL;
151   msReleaseLock( TLOCK_GEOS );
152 #endif
153 }
154 
155 /*
156 ** Translation functions
157 */
msGEOSShape2Geometry_point(pointObj * point)158 static GEOSGeom msGEOSShape2Geometry_point(pointObj *point)
159 {
160   GEOSCoordSeq coords;
161   GEOSGeom g;
162   GEOSContextHandle_t handle = msGetGeosContextHandle();
163 
164   if(!point) return NULL;
165 
166   coords = GEOSCoordSeq_create_r(handle,1, 2); /* todo handle z's */
167   if(!coords) return NULL;
168 
169   GEOSCoordSeq_setX_r(handle,coords, 0, point->x);
170   GEOSCoordSeq_setY_r(handle,coords, 0, point->y);
171   /* GEOSCoordSeq_setY(coords, 0, point->z); */
172 
173   g = GEOSGeom_createPoint_r(handle,coords); /* g owns the coordinate in coords */
174 
175   return g;
176 }
177 
msGEOSShape2Geometry_multipoint(lineObj * multipoint)178 static GEOSGeom msGEOSShape2Geometry_multipoint(lineObj *multipoint)
179 {
180   int i;
181   GEOSGeom g;
182   GEOSGeom *points;
183   GEOSContextHandle_t handle = msGetGeosContextHandle();
184 
185   if(!multipoint) return NULL;
186 
187   points = malloc(multipoint->numpoints*sizeof(GEOSGeom));
188   if(!points) return NULL;
189 
190   for(i=0; i<multipoint->numpoints; i++)
191     points[i] = msGEOSShape2Geometry_point(&(multipoint->point[i]));
192 
193   g = GEOSGeom_createCollection_r(handle,GEOS_MULTIPOINT, points, multipoint->numpoints);
194 
195   free(points);
196 
197   return g;
198 }
199 
msGEOSShape2Geometry_line(lineObj * line)200 static GEOSGeom msGEOSShape2Geometry_line(lineObj *line)
201 {
202   int i;
203   GEOSGeom g;
204   GEOSCoordSeq coords;
205   GEOSContextHandle_t handle = msGetGeosContextHandle();
206 
207   if(!line) return NULL;
208 
209   coords = GEOSCoordSeq_create_r(handle,line->numpoints, 2); /* todo handle z's */
210   if(!coords) return NULL;
211 
212   for(i=0; i<line->numpoints; i++) {
213     GEOSCoordSeq_setX_r(handle,coords, i, line->point[i].x);
214     GEOSCoordSeq_setY_r(handle,coords, i, line->point[i].y);
215     /* GEOSCoordSeq_setZ(coords, i, line->point[i].z); */
216   }
217 
218   g = GEOSGeom_createLineString_r(handle,coords); /* g owns the coordinates in coords */
219 
220   return g;
221 }
222 
msGEOSShape2Geometry_multiline(shapeObj * multiline)223 static GEOSGeom msGEOSShape2Geometry_multiline(shapeObj *multiline)
224 {
225   int i;
226   GEOSGeom g;
227   GEOSGeom *lines;
228   GEOSContextHandle_t handle = msGetGeosContextHandle();
229 
230   if(!multiline) return NULL;
231 
232   lines = malloc(multiline->numlines*sizeof(GEOSGeom));
233   if(!lines) return NULL;
234 
235   for(i=0; i<multiline->numlines; i++)
236     lines[i] = msGEOSShape2Geometry_line(&(multiline->line[i]));
237 
238   g = GEOSGeom_createCollection_r(handle,GEOS_MULTILINESTRING, lines, multiline->numlines);
239 
240   free(lines);
241 
242   return g;
243 }
244 
msGEOSShape2Geometry_simplepolygon(shapeObj * shape,int r,int * outerList)245 static GEOSGeom msGEOSShape2Geometry_simplepolygon(shapeObj *shape, int r, int *outerList)
246 {
247   int i, j, k;
248   GEOSCoordSeq coords;
249   GEOSGeom g;
250   GEOSGeom outerRing;
251   GEOSGeom *innerRings=NULL;
252   int numInnerRings=0, *innerList;
253   GEOSContextHandle_t handle = msGetGeosContextHandle();
254 
255   if(!shape || !outerList) return NULL;
256 
257   /* build the outer shell */
258   coords = GEOSCoordSeq_create_r(handle,shape->line[r].numpoints, 2); /* todo handle z's */
259   if(!coords) return NULL;
260 
261   for(i=0; i<shape->line[r].numpoints; i++) {
262     GEOSCoordSeq_setX_r(handle,coords, i, shape->line[r].point[i].x);
263     GEOSCoordSeq_setY_r(handle,coords, i, shape->line[r].point[i].y);
264     /* GEOSCoordSeq_setZ(coords, i, shape->line[r].point[i].z); */
265   }
266 
267   outerRing = GEOSGeom_createLinearRing_r(handle,coords); /* outerRing owns the coordinates in coords */
268 
269   /* build the holes */
270   innerList = msGetInnerList(shape, r, outerList);
271   for(j=0; j<shape->numlines; j++)
272     if(innerList[j] == MS_TRUE) numInnerRings++;
273 
274   if(numInnerRings > 0) {
275     k = 0; /* inner ring counter */
276 
277     innerRings = msSmallMalloc(numInnerRings*sizeof(GEOSGeom));
278 
279     for(j=0; j<shape->numlines; j++) {
280       if(innerList[j] == MS_FALSE) continue;
281 
282       coords = GEOSCoordSeq_create_r(handle,shape->line[j].numpoints, 2); /* todo handle z's */
283       if(!coords) {
284         free(innerRings);
285         free(innerList);
286         return NULL; /* todo, this will leak memory (shell + allocated holes) */
287       }
288 
289       for(i=0; i<shape->line[j].numpoints; i++) {
290         GEOSCoordSeq_setX_r(handle,coords, i, shape->line[j].point[i].x);
291         GEOSCoordSeq_setY_r(handle,coords, i, shape->line[j].point[i].y);
292         /* GEOSCoordSeq_setZ(coords, i, shape->line[j].point[i].z); */
293       }
294 
295       innerRings[k] = GEOSGeom_createLinearRing_r(handle,coords); /* innerRings[k] owns the coordinates in coords */
296       k++;
297     }
298   }
299 
300   g = GEOSGeom_createPolygon_r(handle,outerRing, innerRings, numInnerRings);
301 
302   free(innerList); /* clean up */
303   free(innerRings); /* clean up */
304 
305   return g;
306 }
307 
msGEOSShape2Geometry_polygon(shapeObj * shape)308 static GEOSGeom msGEOSShape2Geometry_polygon(shapeObj *shape)
309 {
310   int i, j;
311   GEOSGeom *polygons;
312   int *outerList, numOuterRings=0, lastOuterRing=0;
313   GEOSGeom g;
314   GEOSContextHandle_t handle = msGetGeosContextHandle();
315 
316   outerList = msGetOuterList(shape);
317   for(i=0; i<shape->numlines; i++) {
318     if(outerList[i] == MS_TRUE) {
319       numOuterRings++;
320       lastOuterRing = i; /* save for the simple case */
321     }
322   }
323 
324   if(numOuterRings == 1) {
325     g = msGEOSShape2Geometry_simplepolygon(shape, lastOuterRing, outerList);
326   } else { /* a true multipolygon */
327     polygons = msSmallMalloc(numOuterRings*sizeof(GEOSGeom));
328 
329     j = 0; /* part counter */
330     for(i=0; i<shape->numlines; i++) {
331       if(outerList[i] == MS_FALSE) continue;
332       polygons[j] = msGEOSShape2Geometry_simplepolygon(shape, i, outerList); /* TODO: account for NULL return values */
333       j++;
334     }
335 
336     g = GEOSGeom_createCollection_r(handle,GEOS_MULTIPOLYGON, polygons, numOuterRings);
337     free(polygons);
338   }
339 
340   free(outerList);
341   return g;
342 }
343 
msGEOSShape2Geometry(shapeObj * shape)344 GEOSGeom msGEOSShape2Geometry(shapeObj *shape)
345 {
346   if(!shape)
347     return NULL; /* a NULL shape generates a NULL geometry */
348 
349   switch(shape->type) {
350     case MS_SHAPE_POINT:
351       if(shape->numlines == 0 || shape->line[0].numpoints == 0) /* not enough info for a point */
352         return NULL;
353 
354       if(shape->line[0].numpoints == 1) /* simple point */
355         return msGEOSShape2Geometry_point(&(shape->line[0].point[0]));
356       else /* multi-point */
357         return msGEOSShape2Geometry_multipoint(&(shape->line[0]));
358       break;
359     case MS_SHAPE_LINE:
360       if(shape->numlines == 0 || shape->line[0].numpoints < 2) /* not enough info for a line */
361         return NULL;
362 
363       if(shape->numlines == 1) /* simple line */
364         return msGEOSShape2Geometry_line(&(shape->line[0]));
365       else /* multi-line */
366         return msGEOSShape2Geometry_multiline(shape);
367       break;
368     case MS_SHAPE_POLYGON:
369       if(shape->numlines == 0 || shape->line[0].numpoints < 4) /* not enough info for a polygon (first=last) */
370         return NULL;
371 
372       return msGEOSShape2Geometry_polygon(shape); /* simple and multipolygon cases are addressed */
373       break;
374     default:
375       break;
376   }
377 
378   return NULL; /* should not get here */
379 }
380 
msGEOSGeometry2Shape_point(GEOSGeom g)381 static shapeObj *msGEOSGeometry2Shape_point(GEOSGeom g)
382 {
383   GEOSCoordSeq coords;
384   shapeObj *shape=NULL;
385   GEOSContextHandle_t handle = msGetGeosContextHandle();
386 
387   if(!g) return NULL;
388 
389   shape = (shapeObj *) malloc(sizeof(shapeObj));
390   msInitShape(shape);
391 
392   shape->type = MS_SHAPE_POINT;
393   shape->line = (lineObj *) malloc(sizeof(lineObj));
394   shape->numlines = 1;
395   shape->line[0].point = (pointObj *) malloc(sizeof(pointObj));
396   shape->line[0].numpoints = 1;
397   shape->geometry = (GEOSGeom) g;
398 
399   coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,g);
400 
401   GEOSCoordSeq_getX_r(handle,coords, 0, &(shape->line[0].point[0].x));
402   GEOSCoordSeq_getY_r(handle,coords, 0, &(shape->line[0].point[0].y));
403   /* GEOSCoordSeq_getZ(coords, 0, &(shape->line[0].point[0].z)); */
404 
405   shape->bounds.minx = shape->bounds.maxx = shape->line[0].point[0].x;
406   shape->bounds.miny = shape->bounds.maxy = shape->line[0].point[0].y;
407 
408   return shape;
409 }
410 
msGEOSGeometry2Shape_multipoint(GEOSGeom g)411 static shapeObj *msGEOSGeometry2Shape_multipoint(GEOSGeom g)
412 {
413   int i;
414   int numPoints;
415   GEOSCoordSeq coords;
416   GEOSGeom point;
417 
418   shapeObj *shape=NULL;
419   GEOSContextHandle_t handle = msGetGeosContextHandle();
420 
421   if(!g) return NULL;
422   numPoints = GEOSGetNumGeometries_r(handle,g); /* each geometry has 1 point */
423 
424   shape = (shapeObj *) malloc(sizeof(shapeObj));
425   msInitShape(shape);
426 
427   shape->type = MS_SHAPE_POINT;
428   shape->line = (lineObj *) malloc(sizeof(lineObj));
429   shape->numlines = 1;
430   shape->line[0].point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
431   shape->line[0].numpoints = numPoints;
432   shape->geometry = (GEOSGeom) g;
433 
434   for(i=0; i<numPoints; i++) {
435     point = (GEOSGeom) GEOSGetGeometryN_r(handle,g, i);
436     coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,point);
437 
438     GEOSCoordSeq_getX_r(handle,coords, 0, &(shape->line[0].point[i].x));
439     GEOSCoordSeq_getY_r(handle,coords, 0, &(shape->line[0].point[i].y));
440     /* GEOSCoordSeq_getZ(coords, 0, &(shape->line[0].point[i].z)); */
441   }
442 
443   msComputeBounds(shape);
444 
445   return shape;
446 }
447 
msGEOSGeometry2Shape_line(GEOSGeom g)448 static shapeObj *msGEOSGeometry2Shape_line(GEOSGeom g)
449 {
450   shapeObj *shape=NULL;
451   GEOSContextHandle_t handle = msGetGeosContextHandle();
452 
453   int i;
454   int numPoints;
455   GEOSCoordSeq coords;
456 
457   if(!g) return NULL;
458   numPoints = GEOSGetNumCoordinates_r(handle,g);
459   coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,g);
460 
461   shape = (shapeObj *) malloc(sizeof(shapeObj));
462   msInitShape(shape);
463 
464   shape->type = MS_SHAPE_LINE;
465   shape->line = (lineObj *) malloc(sizeof(lineObj));
466   shape->numlines = 1;
467   shape->line[0].point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
468   shape->line[0].numpoints = numPoints;
469   shape->geometry = (GEOSGeom) g;
470 
471   for(i=0; i<numPoints; i++) {
472     GEOSCoordSeq_getX_r(handle,coords, i, &(shape->line[0].point[i].x));
473     GEOSCoordSeq_getY_r(handle,coords, i, &(shape->line[0].point[i].y));
474     /* GEOSCoordSeq_getZ(coords, i, &(shape->line[0].point[i].z)); */
475   }
476 
477   msComputeBounds(shape);
478 
479   return shape;
480 }
481 
msGEOSGeometry2Shape_multiline(GEOSGeom g)482 static shapeObj *msGEOSGeometry2Shape_multiline(GEOSGeom g)
483 {
484   int i, j;
485   int numPoints, numLines;
486   GEOSCoordSeq coords;
487   GEOSGeom lineString;
488 
489   shapeObj *shape=NULL;
490   lineObj line;
491   GEOSContextHandle_t handle = msGetGeosContextHandle();
492 
493   if(!g) return NULL;
494   numLines = GEOSGetNumGeometries_r(handle,g);
495 
496   shape = (shapeObj *) malloc(sizeof(shapeObj));
497   msInitShape(shape);
498 
499   shape->type = MS_SHAPE_LINE;
500   shape->geometry = (GEOSGeom) g;
501 
502   for(j=0; j<numLines; j++) {
503     lineString = (GEOSGeom) GEOSGetGeometryN_r(handle,g, j);
504     numPoints = GEOSGetNumCoordinates_r(handle,lineString);
505     coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,lineString);
506 
507     line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
508     line.numpoints = numPoints;
509 
510     for(i=0; i<numPoints; i++) {
511       GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
512       GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
513       /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
514     }
515 
516     msAddLineDirectly(shape, &line);
517   }
518 
519   msComputeBounds(shape);
520 
521   return shape;
522 }
523 
msGEOSGeometry2Shape_polygon(GEOSGeom g)524 static shapeObj *msGEOSGeometry2Shape_polygon(GEOSGeom g)
525 {
526   shapeObj *shape=NULL;
527   lineObj line;
528   int numPoints, numRings;
529   int i, j;
530 
531   GEOSCoordSeq coords;
532   GEOSGeom ring;
533   GEOSContextHandle_t handle = msGetGeosContextHandle();
534 
535   if(!g) return NULL;
536 
537   shape = (shapeObj *) malloc(sizeof(shapeObj));
538   msInitShape(shape);
539   shape->type = MS_SHAPE_POLYGON;
540   shape->geometry = (GEOSGeom) g;
541 
542   /* exterior ring */
543   ring = (GEOSGeom) GEOSGetExteriorRing_r(handle,g);
544   numPoints = GEOSGetNumCoordinates_r(handle,ring);
545   coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
546 
547   line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
548   line.numpoints = numPoints;
549 
550   for(i=0; i<numPoints; i++) {
551     GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
552     GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
553     /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
554   }
555   msAddLineDirectly(shape, &line);
556 
557   /* interior rings */
558   numRings = GEOSGetNumInteriorRings_r(handle,g);
559   for(j=0; j<numRings; j++) {
560     ring = (GEOSGeom) GEOSGetInteriorRingN_r(handle,g, j);
561     if(GEOSisRing_r(handle,ring) != 1) continue; /* skip it */
562 
563     numPoints = GEOSGetNumCoordinates_r(handle,ring);
564     coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
565 
566     line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
567     line.numpoints = numPoints;
568 
569     for(i=0; i<numPoints; i++) {
570       GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
571       GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
572       /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
573     }
574     msAddLineDirectly(shape, &line);
575   }
576 
577   msComputeBounds(shape);
578 
579   return shape;
580 }
581 
msGEOSGeometry2Shape_multipolygon(GEOSGeom g)582 static shapeObj *msGEOSGeometry2Shape_multipolygon(GEOSGeom g)
583 {
584   int i, j, k;
585   shapeObj *shape=NULL;
586   lineObj line;
587   int numPoints, numRings, numPolygons;
588 
589   GEOSCoordSeq coords;
590   GEOSGeom polygon, ring;
591   GEOSContextHandle_t handle = msGetGeosContextHandle();
592 
593   if(!g) return NULL;
594   numPolygons = GEOSGetNumGeometries_r(handle,g);
595 
596   shape = (shapeObj *) malloc(sizeof(shapeObj));
597   msInitShape(shape);
598   shape->type = MS_SHAPE_POLYGON;
599   shape->geometry = (GEOSGeom) g;
600 
601   for(k=0; k<numPolygons; k++) { /* for each polygon */
602     polygon = (GEOSGeom) GEOSGetGeometryN_r(handle,g, k);
603 
604     /* exterior ring */
605     ring = (GEOSGeom) GEOSGetExteriorRing_r(handle,polygon);
606     numPoints = GEOSGetNumCoordinates_r(handle,ring);
607     coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
608 
609     line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
610     line.numpoints = numPoints;
611 
612     for(i=0; i<numPoints; i++) {
613       GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
614       GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
615       /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
616     }
617     msAddLineDirectly(shape, &line);
618 
619     /* interior rings */
620     numRings = GEOSGetNumInteriorRings_r(handle,polygon);
621 
622     for(j=0; j<numRings; j++) {
623       ring = (GEOSGeom) GEOSGetInteriorRingN_r(handle,polygon, j);
624       if(GEOSisRing_r(handle,ring) != 1) continue; /* skip it */
625 
626       numPoints = GEOSGetNumCoordinates_r(handle,ring);
627       coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,ring);
628 
629       line.point = (pointObj *) malloc(sizeof(pointObj)*numPoints);
630       line.numpoints = numPoints;
631 
632       for(i=0; i<numPoints; i++) {
633         GEOSCoordSeq_getX_r(handle,coords, i, &(line.point[i].x));
634         GEOSCoordSeq_getY_r(handle,coords, i, &(line.point[i].y));
635         /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
636       }
637       msAddLineDirectly(shape, &line);
638     }
639   } /* next polygon */
640 
641   msComputeBounds(shape);
642 
643   return shape;
644 }
645 
msGEOSGeometry2Shape(GEOSGeom g)646 shapeObj *msGEOSGeometry2Shape(GEOSGeom g)
647 {
648   int type;
649   GEOSContextHandle_t handle = msGetGeosContextHandle();
650 
651   if(!g)
652     return NULL; /* a NULL geometry generates a NULL shape */
653 
654   type = GEOSGeomTypeId_r(handle,g);
655   switch(type) {
656     case GEOS_POINT:
657       return msGEOSGeometry2Shape_point(g);
658       break;
659     case GEOS_MULTIPOINT:
660       return msGEOSGeometry2Shape_multipoint(g);
661       break;
662     case GEOS_LINESTRING:
663       return msGEOSGeometry2Shape_line(g);
664       break;
665     case GEOS_MULTILINESTRING:
666       return msGEOSGeometry2Shape_multiline(g);
667       break;
668     case GEOS_POLYGON:
669       return msGEOSGeometry2Shape_polygon(g);
670       break;
671     case GEOS_MULTIPOLYGON:
672       return msGEOSGeometry2Shape_multipolygon(g);
673       break;
674     case GEOS_GEOMETRYCOLLECTION:
675       if (!GEOSisEmpty_r(handle,g))
676       {
677         int i, j, numGeoms;
678         shapeObj* shape;
679 
680         numGeoms = GEOSGetNumGeometries_r(handle,g);
681 
682         shape = (shapeObj *) malloc(sizeof(shapeObj));
683         msInitShape(shape);
684         shape->type = MS_SHAPE_LINE;
685         shape->geometry = (GEOSGeom) g;
686 
687         numGeoms = GEOSGetNumGeometries_r(handle,g);
688         for(i = 0; i < numGeoms; i++) { /* for each geometry */
689            shapeObj* shape2 = msGEOSGeometry2Shape((GEOSGeom)GEOSGetGeometryN_r(handle,g, i));
690            if (shape2) {
691               for (j = 0; j < shape2->numlines; j++)
692                  msAddLineDirectly(shape, &shape2->line[j]);
693               shape2->numlines = 0;
694               shape2->geometry = NULL; /* not owned */
695               msFreeShape(shape2);
696            }
697         }
698         msComputeBounds(shape);
699         return shape;
700       }
701       break;
702     default:
703       msSetError(MS_GEOSERR, "Unsupported GEOS geometry type (%d).", "msGEOSGeometry2Shape()", type);
704   }
705   return NULL;
706 }
707 #endif
708 
709 /*
710 ** Maintenence functions exposed to MapServer/MapScript.
711 */
712 
msGEOSFreeGeometry(shapeObj * shape)713 void msGEOSFreeGeometry(shapeObj *shape)
714 {
715 #ifdef USE_GEOS
716   GEOSGeom g=NULL;
717   GEOSContextHandle_t handle = msGetGeosContextHandle();
718 
719   if(!shape || !shape->geometry)
720     return;
721 
722   g = (GEOSGeom) shape->geometry;
723   GEOSGeom_destroy_r(handle,g);
724   shape->geometry = NULL;
725 #else
726   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSFreeGEOSGeom()");
727   return;
728 #endif
729 }
730 
731 /*
732 ** WKT input and output functions
733 */
msGEOSShapeFromWKT(const char * wkt)734 shapeObj *msGEOSShapeFromWKT(const char *wkt)
735 {
736 #ifdef USE_GEOS
737   GEOSGeom g;
738   GEOSContextHandle_t handle = msGetGeosContextHandle();
739 
740   if(!wkt)
741     return NULL;
742 
743   g = GEOSGeomFromWKT_r(handle,wkt);
744   if(!g) {
745     msSetError(MS_GEOSERR, "Error reading WKT geometry \"%s\".", "msGEOSShapeFromWKT()", wkt);
746     return NULL;
747   } else {
748     return msGEOSGeometry2Shape(g);
749   }
750 
751 #else
752   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSShapeFromWKT()");
753   return NULL;
754 #endif
755 }
756 
757 /* Return should be freed with msGEOSFreeWKT */
msGEOSShapeToWKT(shapeObj * shape)758 char *msGEOSShapeToWKT(shapeObj *shape)
759 {
760 #ifdef USE_GEOS
761   GEOSGeom g;
762   GEOSContextHandle_t handle = msGetGeosContextHandle();
763 
764   if(!shape)
765     return NULL;
766 
767   /* if we have a geometry, we should update it*/
768   msGEOSFreeGeometry(shape);
769 
770   shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
771   g = (GEOSGeom) shape->geometry;
772   if(!g) return NULL;
773 
774   return GEOSGeomToWKT_r(handle,g);
775 #else
776   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSShapeToWKT()");
777   return NULL;
778 #endif
779 }
780 
msGEOSFreeWKT(char * pszGEOSWKT)781 void msGEOSFreeWKT(char* pszGEOSWKT)
782 {
783 #ifdef USE_GEOS
784   GEOSContextHandle_t handle = msGetGeosContextHandle();
785 #if GEOS_VERSION_MAJOR > 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 2)
786   GEOSFree_r(handle,pszGEOSWKT);
787 #endif
788 #else
789   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSFreeWKT()");
790 #endif
791 }
792 
msGEOSOffsetCurve(shapeObj * p,double offset)793 shapeObj *msGEOSOffsetCurve(shapeObj *p, double offset) {
794 #if defined USE_GEOS && (GEOS_VERSION_MAJOR > 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 3))
795   GEOSGeom g1, g2;
796   GEOSContextHandle_t handle = msGetGeosContextHandle();
797 
798   if(!p)
799     return NULL;
800 
801   if(!p->geometry) /* if no geometry for the shape then build one */
802     p->geometry = (GEOSGeom) msGEOSShape2Geometry(p);
803 
804   g1 = (GEOSGeom) p->geometry;
805   if(!g1) return NULL;
806 
807   if (GEOSGeomTypeId_r(handle,g1) == GEOS_MULTILINESTRING)
808   {
809     GEOSGeom *lines = malloc(p->numlines*sizeof(GEOSGeom));
810     if (!lines) return NULL;
811     for(int i=0; i<p->numlines; i++)
812     {
813       lines[i] = GEOSOffsetCurve_r(handle, GEOSGetGeometryN_r(handle,g1,i),
814                                    offset, 4, GEOSBUF_JOIN_MITRE, fabs(offset*1.5));
815     }
816     g2 = GEOSGeom_createCollection_r(handle,GEOS_MULTILINESTRING, lines, p->numlines);
817     free(lines);
818   }
819   else
820   {
821     g2 = GEOSOffsetCurve_r(handle,g1, offset, 4, GEOSBUF_JOIN_MITRE, fabs(offset*1.5));
822   }
823   return msGEOSGeometry2Shape(g2);
824 #else
825   msSetError(MS_GEOSERR, "GEOS Offset Curve support is not available.", "msGEOSingleSidedBuffer()");
826   return NULL;
827 #endif
828 }
829 
830 /*
831 ** Analytical functions exposed to MapServer/MapScript.
832 */
833 
msGEOSBuffer(shapeObj * shape,double width)834 shapeObj *msGEOSBuffer(shapeObj *shape, double width)
835 {
836 #ifdef USE_GEOS
837   GEOSGeom g1, g2;
838   GEOSContextHandle_t handle = msGetGeosContextHandle();
839 
840   if(!shape)
841     return NULL;
842 
843   if(!shape->geometry) /* if no geometry for the shape then build one */
844     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
845 
846   g1 = (GEOSGeom) shape->geometry;
847   if(!g1) return NULL;
848 
849   g2 = GEOSBuffer_r(handle,g1, width, 30);
850   return msGEOSGeometry2Shape(g2);
851 #else
852   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBuffer()");
853   return NULL;
854 #endif
855 }
856 
msGEOSSimplify(shapeObj * shape,double tolerance)857 shapeObj *msGEOSSimplify(shapeObj *shape, double tolerance)
858 {
859 #ifdef USE_GEOS
860   GEOSGeom g1, g2;
861   GEOSContextHandle_t handle = msGetGeosContextHandle();
862 
863   if(!shape)
864     return NULL;
865 
866   if(!shape->geometry) /* if no geometry for the shape then build one */
867     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
868 
869   g1 = (GEOSGeom) shape->geometry;
870   if(!g1) return NULL;
871 
872   g2 = GEOSSimplify_r(handle,g1, tolerance);
873   return msGEOSGeometry2Shape(g2);
874 #else
875   msSetError(MS_GEOSERR, "GEOS Simplify support is not available.", "msGEOSSimplify()");
876   return NULL;
877 #endif
878 }
879 
msGEOSTopologyPreservingSimplify(shapeObj * shape,double tolerance)880 shapeObj *msGEOSTopologyPreservingSimplify(shapeObj *shape, double tolerance)
881 {
882 #ifdef USE_GEOS
883   GEOSGeom g1, g2;
884   GEOSContextHandle_t handle = msGetGeosContextHandle();
885 
886   if(!shape)
887     return NULL;
888 
889   if(!shape->geometry) /* if no geometry for the shape then build one */
890     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
891 
892   g1 = (GEOSGeom) shape->geometry;
893   if(!g1) return NULL;
894 
895   g2 = GEOSTopologyPreserveSimplify_r(handle,g1, tolerance);
896   return msGEOSGeometry2Shape(g2);
897 #else
898   msSetError(MS_GEOSERR, "GEOS Simplify support is not available.", "msGEOSTopologyPreservingSimplify()");
899   return NULL;
900 #endif
901 }
902 
msGEOSConvexHull(shapeObj * shape)903 shapeObj *msGEOSConvexHull(shapeObj *shape)
904 {
905 #ifdef USE_GEOS
906   GEOSGeom g1, g2;
907   GEOSContextHandle_t handle = msGetGeosContextHandle();
908 
909   if(!shape) return NULL;
910 
911   if(!shape->geometry) /* if no geometry for the shape then build one */
912     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
913   g1 = (GEOSGeom) shape->geometry;
914   if(!g1) return NULL;
915 
916   g2 = GEOSConvexHull_r(handle,g1);
917   return msGEOSGeometry2Shape(g2);
918 #else
919   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSConvexHull()");
920   return NULL;
921 #endif
922 }
923 
msGEOSBoundary(shapeObj * shape)924 shapeObj *msGEOSBoundary(shapeObj *shape)
925 {
926 #ifdef USE_GEOS
927   GEOSGeom g1, g2;
928   GEOSContextHandle_t handle = msGetGeosContextHandle();
929 
930   if(!shape) return NULL;
931 
932   if(!shape->geometry) /* if no geometry for the shape then build one */
933     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
934   g1 = (GEOSGeom) shape->geometry;
935   if(!g1) return NULL;
936 
937   g2 = GEOSBoundary_r(handle,g1);
938   return msGEOSGeometry2Shape(g2);
939 #else
940   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBoundary()");
941   return NULL;
942 #endif
943 }
944 
msGEOSGetCentroid(shapeObj * shape)945 pointObj *msGEOSGetCentroid(shapeObj *shape)
946 {
947 #ifdef USE_GEOS
948   GEOSGeom g1, g2;
949   GEOSCoordSeq coords;
950   pointObj *point;
951   GEOSContextHandle_t handle = msGetGeosContextHandle();
952 
953   if(!shape) return NULL;
954 
955   if(!shape->geometry) /* if no geometry for the shape then build one */
956     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
957   g1 = (GEOSGeom) shape->geometry;
958   if(!g1) return NULL;
959 
960   g2 = GEOSGetCentroid_r(handle,g1);
961   if (!g2) return NULL;
962 
963   point = (pointObj *) malloc(sizeof(pointObj));
964 
965   coords = (GEOSCoordSeq) GEOSGeom_getCoordSeq_r(handle,g2);
966 
967   GEOSCoordSeq_getX_r(handle,coords, 0, &(point->x));
968   GEOSCoordSeq_getY_r(handle,coords, 0, &(point->y));
969   /* GEOSCoordSeq_getZ(coords, 0, &(point->z)); */
970 
971   GEOSGeom_destroy_r(handle, g2);
972 
973   return point;
974 #else
975   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSGetCentroid()");
976   return NULL;
977 #endif
978 }
979 
msGEOSUnion(shapeObj * shape1,shapeObj * shape2)980 shapeObj *msGEOSUnion(shapeObj *shape1, shapeObj *shape2)
981 {
982 #ifdef USE_GEOS
983   GEOSGeom g1, g2, g3;
984   GEOSContextHandle_t handle = msGetGeosContextHandle();
985 
986   if(!shape1 || !shape2)
987     return NULL;
988 
989   if(!shape1->geometry) /* if no geometry for the shape then build one */
990     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
991   g1 = (GEOSGeom) shape1->geometry;
992   if(!g1) return NULL;
993 
994   if(!shape2->geometry) /* if no geometry for the shape then build one */
995     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
996   g2 = (GEOSGeom) shape2->geometry;
997   if(!g2) return NULL;
998 
999   g3 = GEOSUnion_r(handle,g1, g2);
1000   return msGEOSGeometry2Shape(g3);
1001 #else
1002   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSUnion()");
1003   return NULL;
1004 #endif
1005 }
1006 
msGEOSIntersection(shapeObj * shape1,shapeObj * shape2)1007 shapeObj *msGEOSIntersection(shapeObj *shape1, shapeObj *shape2)
1008 {
1009 #ifdef USE_GEOS
1010   GEOSGeom g1, g2, g3;
1011   GEOSContextHandle_t handle = msGetGeosContextHandle();
1012 
1013   if(!shape1 || !shape2)
1014     return NULL;
1015 
1016   if(!shape1->geometry) /* if no geometry for the shape then build one */
1017     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1018   g1 = (GEOSGeom) shape1->geometry;
1019   if(!g1) return NULL;
1020 
1021   if(!shape2->geometry) /* if no geometry for the shape then build one */
1022     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1023   g2 = (GEOSGeom) shape2->geometry;
1024   if(!g2) return NULL;
1025 
1026   g3 = GEOSIntersection_r(handle,g1, g2);
1027   return msGEOSGeometry2Shape(g3);
1028 #else
1029   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSIntersection()");
1030   return NULL;
1031 #endif
1032 }
1033 
msGEOSDifference(shapeObj * shape1,shapeObj * shape2)1034 shapeObj *msGEOSDifference(shapeObj *shape1, shapeObj *shape2)
1035 {
1036 #ifdef USE_GEOS
1037   GEOSGeom g1, g2, g3;
1038   GEOSContextHandle_t handle = msGetGeosContextHandle();
1039 
1040   if(!shape1 || !shape2)
1041     return NULL;
1042 
1043   if(!shape1->geometry) /* if no geometry for the shape then build one */
1044     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1045   g1 = (GEOSGeom) shape1->geometry;
1046   if(!g1) return NULL;
1047 
1048   if(!shape2->geometry) /* if no geometry for the shape then build one */
1049     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1050   g2 = (GEOSGeom) shape2->geometry;
1051   if(!g2) return NULL;
1052 
1053   g3 = GEOSDifference_r(handle,g1, g2);
1054   return msGEOSGeometry2Shape(g3);
1055 #else
1056   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSDifference()");
1057   return NULL;
1058 #endif
1059 }
1060 
msGEOSSymDifference(shapeObj * shape1,shapeObj * shape2)1061 shapeObj *msGEOSSymDifference(shapeObj *shape1, shapeObj *shape2)
1062 {
1063 #ifdef USE_GEOS
1064   GEOSGeom g1, g2, g3;
1065   GEOSContextHandle_t handle = msGetGeosContextHandle();
1066 
1067   if(!shape1 || !shape2)
1068     return NULL;
1069 
1070   if(!shape1->geometry) /* if no geometry for the shape then build one */
1071     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1072   g1 = (GEOSGeom) shape1->geometry;
1073   if(!g1) return NULL;
1074 
1075   if(!shape2->geometry) /* if no geometry for the shape then build one */
1076     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1077   g2 = (GEOSGeom) shape2->geometry;
1078   if(!g2) return NULL;
1079 
1080   g3 = GEOSSymDifference_r(handle,g1, g2);
1081   return msGEOSGeometry2Shape(g3);
1082 #else
1083   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSSymDifference()");
1084   return NULL;
1085 #endif
1086 }
1087 
1088 /*
1089 ** Binary predicates exposed to MapServer/MapScript
1090 */
1091 
1092 /*
1093 ** Does shape1 contain shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1094 */
msGEOSContains(shapeObj * shape1,shapeObj * shape2)1095 int msGEOSContains(shapeObj *shape1, shapeObj *shape2)
1096 {
1097 #ifdef USE_GEOS
1098   GEOSGeom g1, g2;
1099   int result;
1100   GEOSContextHandle_t handle = msGetGeosContextHandle();
1101 
1102   if(!shape1 || !shape2)
1103     return -1;
1104 
1105   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1106     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1107   g1 = shape1->geometry;
1108   if(!g1) return -1;
1109 
1110   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1111     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1112   g2 = shape2->geometry;
1113   if(!g2) return -1;
1114 
1115   result = GEOSContains_r(handle,g1, g2);
1116   return ((result==2) ? -1 : result);
1117 #else
1118   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSContains()");
1119   return -1;
1120 #endif
1121 }
1122 
1123 /*
1124 ** Does shape1 overlap shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1125 */
msGEOSOverlaps(shapeObj * shape1,shapeObj * shape2)1126 int msGEOSOverlaps(shapeObj *shape1, shapeObj *shape2)
1127 {
1128 #ifdef USE_GEOS
1129   GEOSGeom g1, g2;
1130   int result;
1131   GEOSContextHandle_t handle = msGetGeosContextHandle();
1132 
1133   if(!shape1 || !shape2)
1134     return -1;
1135 
1136   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1137     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1138   g1 = shape1->geometry;
1139   if(!g1) return -1;
1140 
1141   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1142     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1143   g2 = shape2->geometry;
1144   if(!g2) return -1;
1145 
1146   result = GEOSOverlaps_r(handle,g1, g2);
1147   return ((result==2) ? -1 : result);
1148 #else
1149   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSOverlaps()");
1150   return -1;
1151 #endif
1152 }
1153 
1154 /*
1155 ** Is shape1 within shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1156 */
msGEOSWithin(shapeObj * shape1,shapeObj * shape2)1157 int msGEOSWithin(shapeObj *shape1, shapeObj *shape2)
1158 {
1159 #ifdef USE_GEOS
1160   GEOSGeom g1, g2;
1161   int result;
1162   GEOSContextHandle_t handle = msGetGeosContextHandle();
1163 
1164   if(!shape1 || !shape2)
1165     return -1;
1166 
1167   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1168     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1169   g1 = shape1->geometry;
1170   if(!g1) return -1;
1171 
1172   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1173     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1174   g2 = shape2->geometry;
1175   if(!g2) return -1;
1176 
1177   result = GEOSWithin_r(handle,g1, g2);
1178   return ((result==2) ? -1 : result);
1179 #else
1180   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSWithin()");
1181   return -1;
1182 #endif
1183 }
1184 
1185 /*
1186 ** Does shape1 cross shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1187 */
msGEOSCrosses(shapeObj * shape1,shapeObj * shape2)1188 int msGEOSCrosses(shapeObj *shape1, shapeObj *shape2)
1189 {
1190 #ifdef USE_GEOS
1191   GEOSGeom g1, g2;
1192   int result;
1193   GEOSContextHandle_t handle = msGetGeosContextHandle();
1194 
1195   if(!shape1 || !shape2)
1196     return -1;
1197 
1198   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1199     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1200   g1 = shape1->geometry;
1201   if(!g1) return -1;
1202 
1203   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1204     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1205   g2 = shape2->geometry;
1206   if(!g2) return -1;
1207 
1208   result = GEOSCrosses_r(handle,g1, g2);
1209   return ((result==2) ? -1 : result);
1210 #else
1211   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSCrosses()");
1212   return -1;
1213 #endif
1214 }
1215 
1216 /*
1217 ** Does shape1 intersect shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1218 */
msGEOSIntersects(shapeObj * shape1,shapeObj * shape2)1219 int msGEOSIntersects(shapeObj *shape1, shapeObj *shape2)
1220 {
1221 #ifdef USE_GEOS
1222   GEOSGeom g1, g2;
1223   int result;
1224   GEOSContextHandle_t handle = msGetGeosContextHandle();
1225 
1226   if(!shape1 || !shape2)
1227     return -1;
1228 
1229   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1230     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1231   g1 = (GEOSGeom) shape1->geometry;
1232   if(!g1) return -1;
1233 
1234   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1235     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1236   g2 = (GEOSGeom) shape2->geometry;
1237   if(!g2) return -1;
1238 
1239   result = GEOSIntersects_r(handle,g1, g2);
1240   return ((result==2) ? -1 : result);
1241 #else
1242   if(!shape1 || !shape2)
1243     return -1;
1244 
1245   switch(shape1->type) { /* todo: deal with point shapes */
1246     case(MS_SHAPE_LINE):
1247       switch(shape2->type) {
1248         case(MS_SHAPE_LINE):
1249           return msIntersectPolylines(shape1, shape2);
1250         case(MS_SHAPE_POLYGON):
1251           return msIntersectPolylinePolygon(shape1, shape2);
1252       }
1253       break;
1254     case(MS_SHAPE_POLYGON):
1255       switch(shape2->type) {
1256         case(MS_SHAPE_LINE):
1257           return msIntersectPolylinePolygon(shape2, shape1);
1258         case(MS_SHAPE_POLYGON):
1259           return msIntersectPolygons(shape1, shape2);
1260       }
1261       break;
1262   }
1263 
1264   return -1;
1265 #endif
1266 }
1267 
1268 /*
1269 ** Does shape1 touch shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1270 */
msGEOSTouches(shapeObj * shape1,shapeObj * shape2)1271 int msGEOSTouches(shapeObj *shape1, shapeObj *shape2)
1272 {
1273 #ifdef USE_GEOS
1274   GEOSGeom g1, g2;
1275   int result;
1276   GEOSContextHandle_t handle = msGetGeosContextHandle();
1277 
1278   if(!shape1 || !shape2)
1279     return -1;
1280 
1281   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1282     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1283   g1 = (GEOSGeom) shape1->geometry;
1284   if(!g1) return -1;
1285 
1286   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1287     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1288   g2 = (GEOSGeom) shape2->geometry;
1289   if(!g2) return -1;
1290 
1291   result = GEOSTouches_r(handle,g1, g2);
1292   return ((result==2) ? -1 : result);
1293 #else
1294   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSTouches()");
1295   return -1;
1296 #endif
1297 }
1298 
1299 /*
1300 ** Does shape1 equal shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
1301 */
msGEOSEquals(shapeObj * shape1,shapeObj * shape2)1302 int msGEOSEquals(shapeObj *shape1, shapeObj *shape2)
1303 {
1304 #ifdef USE_GEOS
1305   GEOSGeom g1, g2;
1306   int result;
1307   GEOSContextHandle_t handle = msGetGeosContextHandle();
1308 
1309   if(!shape1 || !shape2)
1310     return -1;
1311 
1312   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1313     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1314   g1 = (GEOSGeom) shape1->geometry;
1315   if(!g1) return -1;
1316 
1317   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1318     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1319   g2 = (GEOSGeom) shape2->geometry;
1320   if(!g2) return -1;
1321 
1322   result = GEOSEquals_r(handle,g1, g2);
1323   return ((result==2) ? -1 : result);
1324 #else
1325   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSEquals()");
1326   return -1;
1327 #endif
1328 }
1329 
1330 /*
1331 ** Are shape1 and shape2 disjoint, returns MS_TRUE/MS_FALSE or -1 for an error.
1332 */
msGEOSDisjoint(shapeObj * shape1,shapeObj * shape2)1333 int msGEOSDisjoint(shapeObj *shape1, shapeObj *shape2)
1334 {
1335 #ifdef USE_GEOS
1336   GEOSGeom g1, g2;
1337   int result;
1338   GEOSContextHandle_t handle = msGetGeosContextHandle();
1339 
1340   if(!shape1 || !shape2)
1341     return -1;
1342 
1343   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1344     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1345   g1 = (GEOSGeom) shape1->geometry;
1346   if(!g1) return -1;
1347 
1348   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1349     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1350   g2 = (GEOSGeom) shape2->geometry;
1351   if(!g2) return -1;
1352 
1353   result = GEOSDisjoint_r(handle,g1, g2);
1354   return ((result==2) ? -1 : result);
1355 #else
1356   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSDisjoint()");
1357   return -1;
1358 #endif
1359 }
1360 
1361 /*
1362 ** Useful misc. functions that return -1 on failure.
1363 */
msGEOSArea(shapeObj * shape)1364 double msGEOSArea(shapeObj *shape)
1365 {
1366 #if defined(USE_GEOS) && defined(GEOS_CAPI_VERSION_MAJOR) && defined(GEOS_CAPI_VERSION_MINOR) && (GEOS_CAPI_VERSION_MAJOR > 1 || GEOS_CAPI_VERSION_MINOR >= 1)
1367   GEOSGeom g;
1368   double area;
1369   int result;
1370   GEOSContextHandle_t handle = msGetGeosContextHandle();
1371 
1372   if(!shape) return -1;
1373 
1374   if(!shape->geometry) /* if no geometry for the shape then build one */
1375     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
1376   g = (GEOSGeom) shape->geometry;
1377   if(!g) return -1;
1378 
1379   result = GEOSArea_r(handle,g, &area);
1380   return  ((result==0) ? -1 : area);
1381 #elif defined(USE_GEOS)
1382   msSetError(MS_GEOSERR, "GEOS support enabled, but old version lacks GEOSArea().", "msGEOSArea()");
1383   return -1;
1384 #else
1385   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSArea()");
1386   return -1;
1387 #endif
1388 }
1389 
msGEOSLength(shapeObj * shape)1390 double msGEOSLength(shapeObj *shape)
1391 {
1392 #if defined(USE_GEOS) && defined(GEOS_CAPI_VERSION_MAJOR) && defined(GEOS_CAPI_VERSION_MINOR) && (GEOS_CAPI_VERSION_MAJOR > 1 || GEOS_CAPI_VERSION_MINOR >= 1)
1393 
1394   GEOSGeom g;
1395   double length;
1396   int result;
1397   GEOSContextHandle_t handle = msGetGeosContextHandle();
1398 
1399   if(!shape) return -1;
1400 
1401   if(!shape->geometry) /* if no geometry for the shape then build one */
1402     shape->geometry = (GEOSGeom) msGEOSShape2Geometry(shape);
1403   g = (GEOSGeom) shape->geometry;
1404   if(!g) return -1;
1405 
1406   result = GEOSLength_r(handle,g, &length);
1407   return  ((result==0) ? -1 : length);
1408 #elif defined(USE_GEOS)
1409   msSetError(MS_GEOSERR, "GEOS support enabled, but old version lacks GEOSLength().", "msGEOSLength()");
1410   return -1;
1411 #else
1412   msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSLength()");
1413   return -1;
1414 #endif
1415 }
1416 
msGEOSDistance(shapeObj * shape1,shapeObj * shape2)1417 double msGEOSDistance(shapeObj *shape1, shapeObj *shape2)
1418 {
1419 #ifdef USE_GEOS
1420   GEOSGeom g1, g2;
1421   double distance;
1422   int result;
1423   GEOSContextHandle_t handle = msGetGeosContextHandle();
1424 
1425   if(!shape1 || !shape2)
1426     return -1;
1427 
1428   if(!shape1->geometry) /* if no geometry for shape1 then build one */
1429     shape1->geometry = (GEOSGeom) msGEOSShape2Geometry(shape1);
1430   g1 = (GEOSGeom) shape1->geometry;
1431   if(!g1) return -1;
1432 
1433   if(!shape2->geometry) /* if no geometry for shape2 then build one */
1434     shape2->geometry = (GEOSGeom) msGEOSShape2Geometry(shape2);
1435   g2 = (GEOSGeom) shape2->geometry;
1436   if(!g2) return -1;
1437 
1438   result = GEOSDistance_r(handle,g1, g2, &distance);
1439   return ((result==0) ? -1 : distance);
1440 #else
1441   return msDistanceShapeToShape(shape1, shape2); /* fall back on brute force method (for MapScript) */
1442 #endif
1443 }
1444