1 /**********************************************************************
2  *
3  * rttopo - topology library
4  * http://git.osgeo.org/gitea/rttopo/librttopo
5  *
6  * rttopo 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  * rttopo 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 rttopo.  If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright (C) 2001-2006 Refractions Research Inc.
22  *
23  **********************************************************************/
24 
25 
26 
27 #include "rttopo_config.h"
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <stdarg.h>
31 
32 #include "librttopo_geom_internal.h"
33 #include "rtgeom_log.h"
34 
35 
36 /** Force Right-hand-rule on RTGEOM polygons **/
37 void
rtgeom_force_clockwise(const RTCTX * ctx,RTGEOM * rtgeom)38 rtgeom_force_clockwise(const RTCTX *ctx, RTGEOM *rtgeom)
39 {
40   RTCOLLECTION *coll;
41   int i;
42 
43   switch (rtgeom->type)
44   {
45   case RTPOLYGONTYPE:
46     rtpoly_force_clockwise(ctx, (RTPOLY *)rtgeom);
47     return;
48 
49   case RTTRIANGLETYPE:
50     rttriangle_force_clockwise(ctx, (RTTRIANGLE *)rtgeom);
51     return;
52 
53     /* Not handle POLYHEDRALSURFACE and TIN
54        as they are supposed to be well oriented */
55   case RTMULTIPOLYGONTYPE:
56   case RTCOLLECTIONTYPE:
57     coll = (RTCOLLECTION *)rtgeom;
58     for (i=0; i<coll->ngeoms; i++)
59       rtgeom_force_clockwise(ctx, coll->geoms[i]);
60     return;
61   }
62 }
63 
64 /** Reverse vertex order of RTGEOM **/
65 void
rtgeom_reverse(const RTCTX * ctx,RTGEOM * rtgeom)66 rtgeom_reverse(const RTCTX *ctx, RTGEOM *rtgeom)
67 {
68   int i;
69   RTCOLLECTION *col;
70 
71   switch (rtgeom->type)
72   {
73   case RTLINETYPE:
74     rtline_reverse(ctx, (RTLINE *)rtgeom);
75     return;
76   case RTPOLYGONTYPE:
77     rtpoly_reverse(ctx, (RTPOLY *)rtgeom);
78     return;
79   case RTTRIANGLETYPE:
80     rttriangle_reverse(ctx, (RTTRIANGLE *)rtgeom);
81     return;
82   case RTMULTILINETYPE:
83   case RTMULTIPOLYGONTYPE:
84   case RTPOLYHEDRALSURFACETYPE:
85   case RTTINTYPE:
86   case RTCOLLECTIONTYPE:
87     col = (RTCOLLECTION *)rtgeom;
88     for (i=0; i<col->ngeoms; i++)
89       rtgeom_reverse(ctx, col->geoms[i]);
90     return;
91   }
92 }
93 
94 RTPOINT *
rtgeom_as_rtpoint(const RTCTX * ctx,const RTGEOM * rtgeom)95 rtgeom_as_rtpoint(const RTCTX *ctx, const RTGEOM *rtgeom)
96 {
97   if ( rtgeom == NULL ) return NULL;
98   if ( rtgeom->type == RTPOINTTYPE )
99     return (RTPOINT *)rtgeom;
100   else return NULL;
101 }
102 
103 RTLINE *
rtgeom_as_rtline(const RTCTX * ctx,const RTGEOM * rtgeom)104 rtgeom_as_rtline(const RTCTX *ctx, const RTGEOM *rtgeom)
105 {
106   if ( rtgeom == NULL ) return NULL;
107   if ( rtgeom->type == RTLINETYPE )
108     return (RTLINE *)rtgeom;
109   else return NULL;
110 }
111 
112 RTCIRCSTRING *
rtgeom_as_rtcircstring(const RTCTX * ctx,const RTGEOM * rtgeom)113 rtgeom_as_rtcircstring(const RTCTX *ctx, const RTGEOM *rtgeom)
114 {
115   if ( rtgeom == NULL ) return NULL;
116   if ( rtgeom->type == RTCIRCSTRINGTYPE )
117     return (RTCIRCSTRING *)rtgeom;
118   else return NULL;
119 }
120 
121 RTCOMPOUND *
rtgeom_as_rtcompound(const RTCTX * ctx,const RTGEOM * rtgeom)122 rtgeom_as_rtcompound(const RTCTX *ctx, const RTGEOM *rtgeom)
123 {
124   if ( rtgeom == NULL ) return NULL;
125   if ( rtgeom->type == RTCOMPOUNDTYPE )
126     return (RTCOMPOUND *)rtgeom;
127   else return NULL;
128 }
129 
130 RTCURVEPOLY *
rtgeom_as_rtcurvepoly(const RTCTX * ctx,const RTGEOM * rtgeom)131 rtgeom_as_rtcurvepoly(const RTCTX *ctx, const RTGEOM *rtgeom)
132 {
133   if ( rtgeom == NULL ) return NULL;
134   if ( rtgeom->type == RTCURVEPOLYTYPE )
135     return (RTCURVEPOLY *)rtgeom;
136   else return NULL;
137 }
138 
139 RTPOLY *
rtgeom_as_rtpoly(const RTCTX * ctx,const RTGEOM * rtgeom)140 rtgeom_as_rtpoly(const RTCTX *ctx, const RTGEOM *rtgeom)
141 {
142   if ( rtgeom == NULL ) return NULL;
143   if ( rtgeom->type == RTPOLYGONTYPE )
144     return (RTPOLY *)rtgeom;
145   else return NULL;
146 }
147 
148 RTTRIANGLE *
rtgeom_as_rttriangle(const RTCTX * ctx,const RTGEOM * rtgeom)149 rtgeom_as_rttriangle(const RTCTX *ctx, const RTGEOM *rtgeom)
150 {
151   if ( rtgeom == NULL ) return NULL;
152   if ( rtgeom->type == RTTRIANGLETYPE )
153     return (RTTRIANGLE *)rtgeom;
154   else return NULL;
155 }
156 
157 RTCOLLECTION *
rtgeom_as_rtcollection(const RTCTX * ctx,const RTGEOM * rtgeom)158 rtgeom_as_rtcollection(const RTCTX *ctx, const RTGEOM *rtgeom)
159 {
160   if ( rtgeom == NULL ) return NULL;
161   if ( rtgeom_is_collection(ctx, rtgeom) )
162     return (RTCOLLECTION*)rtgeom;
163   else return NULL;
164 }
165 
166 RTMPOINT *
rtgeom_as_rtmpoint(const RTCTX * ctx,const RTGEOM * rtgeom)167 rtgeom_as_rtmpoint(const RTCTX *ctx, const RTGEOM *rtgeom)
168 {
169   if ( rtgeom == NULL ) return NULL;
170   if ( rtgeom->type == RTMULTIPOINTTYPE )
171     return (RTMPOINT *)rtgeom;
172   else return NULL;
173 }
174 
175 RTMLINE *
rtgeom_as_rtmline(const RTCTX * ctx,const RTGEOM * rtgeom)176 rtgeom_as_rtmline(const RTCTX *ctx, const RTGEOM *rtgeom)
177 {
178   if ( rtgeom == NULL ) return NULL;
179   if ( rtgeom->type == RTMULTILINETYPE )
180     return (RTMLINE *)rtgeom;
181   else return NULL;
182 }
183 
184 RTMPOLY *
rtgeom_as_rtmpoly(const RTCTX * ctx,const RTGEOM * rtgeom)185 rtgeom_as_rtmpoly(const RTCTX *ctx, const RTGEOM *rtgeom)
186 {
187   if ( rtgeom == NULL ) return NULL;
188   if ( rtgeom->type == RTMULTIPOLYGONTYPE )
189     return (RTMPOLY *)rtgeom;
190   else return NULL;
191 }
192 
193 RTPSURFACE *
rtgeom_as_rtpsurface(const RTCTX * ctx,const RTGEOM * rtgeom)194 rtgeom_as_rtpsurface(const RTCTX *ctx, const RTGEOM *rtgeom)
195 {
196   if ( rtgeom->type == RTPOLYHEDRALSURFACETYPE )
197     return (RTPSURFACE *)rtgeom;
198   else return NULL;
199 }
200 
201 RTTIN *
rtgeom_as_rttin(const RTCTX * ctx,const RTGEOM * rtgeom)202 rtgeom_as_rttin(const RTCTX *ctx, const RTGEOM *rtgeom)
203 {
204   if ( rtgeom->type == RTTINTYPE )
205     return (RTTIN *)rtgeom;
206   else return NULL;
207 }
208 
rttin_as_rtgeom(const RTCTX * ctx,const RTTIN * obj)209 RTGEOM * rttin_as_rtgeom(const RTCTX *ctx, const RTTIN *obj)
210 {
211   return (RTGEOM *)obj;
212 }
213 
rtpsurface_as_rtgeom(const RTCTX * ctx,const RTPSURFACE * obj)214 RTGEOM * rtpsurface_as_rtgeom(const RTCTX *ctx, const RTPSURFACE *obj)
215 {
216   return (RTGEOM *)obj;
217 }
218 
rtmpoly_as_rtgeom(const RTCTX * ctx,const RTMPOLY * obj)219 RTGEOM * rtmpoly_as_rtgeom(const RTCTX *ctx, const RTMPOLY *obj)
220 {
221   if ( obj == NULL ) return NULL;
222   return (RTGEOM *)obj;
223 }
rtmline_as_rtgeom(const RTCTX * ctx,const RTMLINE * obj)224 RTGEOM * rtmline_as_rtgeom(const RTCTX *ctx, const RTMLINE *obj)
225 {
226   if ( obj == NULL ) return NULL;
227   return (RTGEOM *)obj;
228 }
rtmpoint_as_rtgeom(const RTCTX * ctx,const RTMPOINT * obj)229 RTGEOM * rtmpoint_as_rtgeom(const RTCTX *ctx, const RTMPOINT *obj)
230 {
231   if ( obj == NULL ) return NULL;
232   return (RTGEOM *)obj;
233 }
rtcollection_as_rtgeom(const RTCTX * ctx,const RTCOLLECTION * obj)234 RTGEOM * rtcollection_as_rtgeom(const RTCTX *ctx, const RTCOLLECTION *obj)
235 {
236   if ( obj == NULL ) return NULL;
237   return (RTGEOM *)obj;
238 }
rtcircstring_as_rtgeom(const RTCTX * ctx,const RTCIRCSTRING * obj)239 RTGEOM * rtcircstring_as_rtgeom(const RTCTX *ctx, const RTCIRCSTRING *obj)
240 {
241   if ( obj == NULL ) return NULL;
242   return (RTGEOM *)obj;
243 }
rtcurvepoly_as_rtgeom(const RTCTX * ctx,const RTCURVEPOLY * obj)244 RTGEOM * rtcurvepoly_as_rtgeom(const RTCTX *ctx, const RTCURVEPOLY *obj)
245 {
246   if ( obj == NULL ) return NULL;
247   return (RTGEOM *)obj;
248 }
rtcompound_as_rtgeom(const RTCTX * ctx,const RTCOMPOUND * obj)249 RTGEOM * rtcompound_as_rtgeom(const RTCTX *ctx, const RTCOMPOUND *obj)
250 {
251   if ( obj == NULL ) return NULL;
252   return (RTGEOM *)obj;
253 }
rtpoly_as_rtgeom(const RTCTX * ctx,const RTPOLY * obj)254 RTGEOM * rtpoly_as_rtgeom(const RTCTX *ctx, const RTPOLY *obj)
255 {
256   if ( obj == NULL ) return NULL;
257   return (RTGEOM *)obj;
258 }
rttriangle_as_rtgeom(const RTCTX * ctx,const RTTRIANGLE * obj)259 RTGEOM * rttriangle_as_rtgeom(const RTCTX *ctx, const RTTRIANGLE *obj)
260 {
261   if ( obj == NULL ) return NULL;
262   return (RTGEOM *)obj;
263 }
rtline_as_rtgeom(const RTCTX * ctx,const RTLINE * obj)264 RTGEOM * rtline_as_rtgeom(const RTCTX *ctx, const RTLINE *obj)
265 {
266   if ( obj == NULL ) return NULL;
267   return (RTGEOM *)obj;
268 }
rtpoint_as_rtgeom(const RTCTX * ctx,const RTPOINT * obj)269 RTGEOM * rtpoint_as_rtgeom(const RTCTX *ctx, const RTPOINT *obj)
270 {
271   if ( obj == NULL ) return NULL;
272   return (RTGEOM *)obj;
273 }
274 
275 
276 /**
277 ** Look-up for the correct MULTI* type promotion for singleton types.
278 */
279 uint8_t RTMULTITYPE[RTNUMTYPES] =
280 {
281   0,
282   RTMULTIPOINTTYPE,        /*  1 */
283   RTMULTILINETYPE,         /*  2 */
284   RTMULTIPOLYGONTYPE,      /*  3 */
285   0,0,0,0,
286   RTMULTICURVETYPE,        /*  8 */
287   RTMULTICURVETYPE,        /*  9 */
288   RTMULTISURFACETYPE,      /* 10 */
289   RTPOLYHEDRALSURFACETYPE, /* 11 */
290   0, 0,
291   RTTINTYPE,               /* 14 */
292   0
293 };
294 
295 /**
296 * Create a new RTGEOM of the appropriate MULTI* type.
297 */
298 RTGEOM *
rtgeom_as_multi(const RTCTX * ctx,const RTGEOM * rtgeom)299 rtgeom_as_multi(const RTCTX *ctx, const RTGEOM *rtgeom)
300 {
301   RTGEOM **ogeoms;
302   RTGEOM *ogeom = NULL;
303   RTGBOX *box = NULL;
304   int type;
305 
306   type = rtgeom->type;
307 
308   if ( ! RTMULTITYPE[type] ) return rtgeom_clone(ctx, rtgeom);
309 
310   if( rtgeom_is_empty(ctx, rtgeom) )
311   {
312     ogeom = (RTGEOM *)rtcollection_construct_empty(ctx,
313       RTMULTITYPE[type],
314       rtgeom->srid,
315       RTFLAGS_GET_Z(rtgeom->flags),
316       RTFLAGS_GET_M(rtgeom->flags)
317     );
318   }
319   else
320   {
321     ogeoms = rtalloc(ctx, sizeof(RTGEOM*));
322     ogeoms[0] = rtgeom_clone(ctx, rtgeom);
323 
324     /* Sub-geometries are not allowed to have bboxes or SRIDs, move the bbox to the collection */
325     box = ogeoms[0]->bbox;
326     ogeoms[0]->bbox = NULL;
327     ogeoms[0]->srid = SRID_UNKNOWN;
328 
329     ogeom = (RTGEOM *)rtcollection_construct(ctx, RTMULTITYPE[type], rtgeom->srid, box, 1, ogeoms);
330   }
331 
332   return ogeom;
333 }
334 
335 /**
336 * Create a new RTGEOM of the appropriate CURVE* type.
337 */
338 RTGEOM *
rtgeom_as_curve(const RTCTX * ctx,const RTGEOM * rtgeom)339 rtgeom_as_curve(const RTCTX *ctx, const RTGEOM *rtgeom)
340 {
341   RTGEOM *ogeom;
342   int type = rtgeom->type;
343   /*
344   int hasz = RTFLAGS_GET_Z(rtgeom->flags);
345   int hasm = RTFLAGS_GET_M(rtgeom->flags);
346   int srid = rtgeom->srid;
347   */
348 
349   switch(type)
350   {
351     case RTLINETYPE:
352       /* turn to COMPOUNDCURVE */
353       ogeom = (RTGEOM*)rtcompound_construct_from_rtline(ctx, (RTLINE*)rtgeom);
354       break;
355     case RTPOLYGONTYPE:
356       ogeom = (RTGEOM*)rtcurvepoly_construct_from_rtpoly(ctx, rtgeom_as_rtpoly(ctx, rtgeom));
357       break;
358     case RTMULTILINETYPE:
359       /* turn to MULTICURVE */
360       ogeom = rtgeom_clone(ctx, rtgeom);
361       ogeom->type = RTMULTICURVETYPE;
362       break;
363     case RTMULTIPOLYGONTYPE:
364       /* turn to MULTISURFACE */
365       ogeom = rtgeom_clone(ctx, rtgeom);
366       ogeom->type = RTMULTISURFACETYPE;
367       break;
368     case RTCOLLECTIONTYPE:
369     default:
370       ogeom = rtgeom_clone(ctx, rtgeom);
371       break;
372   }
373 
374   /* TODO: copy bbox from input geom ? */
375 
376   return ogeom;
377 }
378 
379 
380 /**
381 * Free the containing RTGEOM and the associated BOX. Leave the underlying
382 * geoms/points/point objects intact. Useful for functions that are stripping
383 * out subcomponents of complex objects, or building up new temporary objects
384 * on top of subcomponents.
385 */
386 void
rtgeom_release(const RTCTX * ctx,RTGEOM * rtgeom)387 rtgeom_release(const RTCTX *ctx, RTGEOM *rtgeom)
388 {
389   if ( ! rtgeom )
390     rterror(ctx, "rtgeom_release: someone called on 0x0");
391 
392   RTDEBUGF(ctx, 3, "releasing type %s", rttype_name(ctx, rtgeom->type));
393 
394   /* Drop bounding box (artays a copy) */
395   if ( rtgeom->bbox )
396   {
397     RTDEBUGF(ctx, 3, "rtgeom_release: releasing bbox. %p", rtgeom->bbox);
398     rtfree(ctx, rtgeom->bbox);
399   }
400   rtfree(ctx, rtgeom);
401 
402 }
403 
404 
405 /* @brief Clone RTGEOM object. Serialized point lists are not copied.
406  *
407  * @see ptarray_clone
408  */
409 RTGEOM *
rtgeom_clone(const RTCTX * ctx,const RTGEOM * rtgeom)410 rtgeom_clone(const RTCTX *ctx, const RTGEOM *rtgeom)
411 {
412   RTDEBUGF(ctx, 2, "rtgeom_clone called with %p, %s",
413            rtgeom, rttype_name(ctx, rtgeom->type));
414 
415   switch (rtgeom->type)
416   {
417   case RTPOINTTYPE:
418     return (RTGEOM *)rtpoint_clone(ctx, (RTPOINT *)rtgeom);
419   case RTLINETYPE:
420     return (RTGEOM *)rtline_clone(ctx, (RTLINE *)rtgeom);
421   case RTCIRCSTRINGTYPE:
422     return (RTGEOM *)rtcircstring_clone(ctx, (RTCIRCSTRING *)rtgeom);
423   case RTPOLYGONTYPE:
424     return (RTGEOM *)rtpoly_clone(ctx, (RTPOLY *)rtgeom);
425   case RTTRIANGLETYPE:
426     return (RTGEOM *)rttriangle_clone(ctx, (RTTRIANGLE *)rtgeom);
427   case RTCOMPOUNDTYPE:
428   case RTCURVEPOLYTYPE:
429   case RTMULTICURVETYPE:
430   case RTMULTISURFACETYPE:
431   case RTMULTIPOINTTYPE:
432   case RTMULTILINETYPE:
433   case RTMULTIPOLYGONTYPE:
434   case RTPOLYHEDRALSURFACETYPE:
435   case RTTINTYPE:
436   case RTCOLLECTIONTYPE:
437     return (RTGEOM *)rtcollection_clone(ctx, (RTCOLLECTION *)rtgeom);
438   default:
439     rterror(ctx, "rtgeom_clone: Unknown geometry type: %s", rttype_name(ctx, rtgeom->type));
440     return NULL;
441   }
442 }
443 
444 /**
445 * Deep-clone an #RTGEOM object. #RTPOINTARRAY <em>are</em> copied.
446 */
447 RTGEOM *
rtgeom_clone_deep(const RTCTX * ctx,const RTGEOM * rtgeom)448 rtgeom_clone_deep(const RTCTX *ctx, const RTGEOM *rtgeom)
449 {
450   RTDEBUGF(ctx, 2, "rtgeom_clone called with %p, %s",
451            rtgeom, rttype_name(ctx, rtgeom->type));
452 
453   switch (rtgeom->type)
454   {
455   case RTPOINTTYPE:
456   case RTLINETYPE:
457   case RTCIRCSTRINGTYPE:
458   case RTTRIANGLETYPE:
459     return (RTGEOM *)rtline_clone_deep(ctx, (RTLINE *)rtgeom);
460   case RTPOLYGONTYPE:
461     return (RTGEOM *)rtpoly_clone_deep(ctx, (RTPOLY *)rtgeom);
462   case RTCOMPOUNDTYPE:
463   case RTCURVEPOLYTYPE:
464   case RTMULTICURVETYPE:
465   case RTMULTISURFACETYPE:
466   case RTMULTIPOINTTYPE:
467   case RTMULTILINETYPE:
468   case RTMULTIPOLYGONTYPE:
469   case RTPOLYHEDRALSURFACETYPE:
470   case RTTINTYPE:
471   case RTCOLLECTIONTYPE:
472     return (RTGEOM *)rtcollection_clone_deep(ctx, (RTCOLLECTION *)rtgeom);
473   default:
474     rterror(ctx, "rtgeom_clone_deep: Unknown geometry type: %s", rttype_name(ctx, rtgeom->type));
475     return NULL;
476   }
477 }
478 
479 
480 /**
481  * Return an alloced string
482  */
483 char*
rtgeom_to_ewkt(const RTCTX * ctx,const RTGEOM * rtgeom)484 rtgeom_to_ewkt(const RTCTX *ctx, const RTGEOM *rtgeom)
485 {
486   char* wkt = NULL;
487   size_t wkt_size = 0;
488 
489   wkt = rtgeom_to_wkt(ctx, rtgeom, RTWKT_EXTENDED, 12, &wkt_size);
490 
491   if ( ! wkt )
492   {
493     rterror(ctx, "Error writing geom %p to RTWKT", rtgeom);
494   }
495 
496   return wkt;
497 }
498 
499 /**
500  * @brief geom1 same as geom2
501  *    iff
502  *        + have same type
503  *      + have same # objects
504  *        + have same bvol
505  *        + each object in geom1 has a corresponding object in geom2 (see above)
506  *  @param rtgeom1
507  *  @param rtgeom2
508  */
509 char
rtgeom_same(const RTCTX * ctx,const RTGEOM * rtgeom1,const RTGEOM * rtgeom2)510 rtgeom_same(const RTCTX *ctx, const RTGEOM *rtgeom1, const RTGEOM *rtgeom2)
511 {
512   RTDEBUGF(ctx, 2, "rtgeom_same(ctx, %s, %s) called",
513            rttype_name(ctx, rtgeom1->type),
514            rttype_name(ctx, rtgeom2->type));
515 
516   if ( rtgeom1->type != rtgeom2->type )
517   {
518     RTDEBUG(ctx, 3, " type differ");
519 
520     return RT_FALSE;
521   }
522 
523   if ( RTFLAGS_GET_ZM(rtgeom1->flags) != RTFLAGS_GET_ZM(rtgeom2->flags) )
524   {
525     RTDEBUG(ctx, 3, " ZM flags differ");
526 
527     return RT_FALSE;
528   }
529 
530   /* Check boxes if both already computed  */
531   if ( rtgeom1->bbox && rtgeom2->bbox )
532   {
533     /*rtnotice(ctx, "bbox1:%p, bbox2:%p", rtgeom1->bbox, rtgeom2->bbox);*/
534     if ( ! gbox_same(ctx, rtgeom1->bbox, rtgeom2->bbox) )
535     {
536       RTDEBUG(ctx, 3, " bounding boxes differ");
537 
538       return RT_FALSE;
539     }
540   }
541 
542   /* geoms have same type, invoke type-specific function */
543   switch (rtgeom1->type)
544   {
545   case RTPOINTTYPE:
546     return rtpoint_same(ctx, (RTPOINT *)rtgeom1,
547                         (RTPOINT *)rtgeom2);
548   case RTLINETYPE:
549     return rtline_same(ctx, (RTLINE *)rtgeom1,
550                        (RTLINE *)rtgeom2);
551   case RTPOLYGONTYPE:
552     return rtpoly_same(ctx, (RTPOLY *)rtgeom1,
553                        (RTPOLY *)rtgeom2);
554   case RTTRIANGLETYPE:
555     return rttriangle_same(ctx, (RTTRIANGLE *)rtgeom1,
556                            (RTTRIANGLE *)rtgeom2);
557   case RTCIRCSTRINGTYPE:
558     return rtcircstring_same(ctx, (RTCIRCSTRING *)rtgeom1,
559            (RTCIRCSTRING *)rtgeom2);
560   case RTMULTIPOINTTYPE:
561   case RTMULTILINETYPE:
562   case RTMULTIPOLYGONTYPE:
563   case RTMULTICURVETYPE:
564   case RTMULTISURFACETYPE:
565   case RTCOMPOUNDTYPE:
566   case RTCURVEPOLYTYPE:
567   case RTPOLYHEDRALSURFACETYPE:
568   case RTTINTYPE:
569   case RTCOLLECTIONTYPE:
570     return rtcollection_same(ctx, (RTCOLLECTION *)rtgeom1,
571                              (RTCOLLECTION *)rtgeom2);
572   default:
573     rterror(ctx, "rtgeom_same: unsupported geometry type: %s",
574             rttype_name(ctx, rtgeom1->type));
575     return RT_FALSE;
576   }
577 
578 }
579 
580 int
rtpoint_inside_circle(const RTCTX * ctx,const RTPOINT * p,double cx,double cy,double rad)581 rtpoint_inside_circle(const RTCTX *ctx, const RTPOINT *p, double cx, double cy, double rad)
582 {
583   const RTPOINT2D *pt;
584   RTPOINT2D center;
585 
586   if ( ! p || ! p->point )
587     return RT_FALSE;
588 
589   pt = rt_getPoint2d_cp(ctx, p->point, 0);
590 
591   center.x = cx;
592   center.y = cy;
593 
594   if ( distance2d_pt_pt(ctx, pt, &center) < rad )
595     return RT_TRUE;
596 
597   return RT_FALSE;
598 }
599 
600 void
rtgeom_drop_bbox(const RTCTX * ctx,RTGEOM * rtgeom)601 rtgeom_drop_bbox(const RTCTX *ctx, RTGEOM *rtgeom)
602 {
603   if ( rtgeom->bbox ) rtfree(ctx, rtgeom->bbox);
604   rtgeom->bbox = NULL;
605   RTFLAGS_SET_BBOX(rtgeom->flags, 0);
606 }
607 
608 /**
609  * Ensure there's a box in the RTGEOM.
610  * If the box is already there just return,
611  * else compute it.
612  */
613 void
rtgeom_add_bbox(const RTCTX * ctx,RTGEOM * rtgeom)614 rtgeom_add_bbox(const RTCTX *ctx, RTGEOM *rtgeom)
615 {
616   /* an empty RTGEOM has no bbox */
617   if ( rtgeom_is_empty(ctx, rtgeom) ) return;
618 
619   if ( rtgeom->bbox ) return;
620   RTFLAGS_SET_BBOX(rtgeom->flags, 1);
621   rtgeom->bbox = gbox_new(ctx, rtgeom->flags);
622   rtgeom_calculate_gbox(ctx, rtgeom, rtgeom->bbox);
623 }
624 
625 void
rtgeom_add_bbox_deep(const RTCTX * ctx,RTGEOM * rtgeom,RTGBOX * gbox)626 rtgeom_add_bbox_deep(const RTCTX *ctx, RTGEOM *rtgeom, RTGBOX *gbox)
627 {
628   if ( rtgeom_is_empty(ctx, rtgeom) ) return;
629 
630   RTFLAGS_SET_BBOX(rtgeom->flags, 1);
631 
632   if ( ! ( gbox || rtgeom->bbox ) )
633   {
634     rtgeom->bbox = gbox_new(ctx, rtgeom->flags);
635     rtgeom_calculate_gbox(ctx, rtgeom, rtgeom->bbox);
636   }
637   else if ( gbox && ! rtgeom->bbox )
638   {
639     rtgeom->bbox = gbox_clone(ctx, gbox);
640   }
641 
642   if ( rtgeom_is_collection(ctx, rtgeom) )
643   {
644     int i;
645     RTCOLLECTION *rtcol = (RTCOLLECTION*)rtgeom;
646 
647     for ( i = 0; i < rtcol->ngeoms; i++ )
648     {
649       rtgeom_add_bbox_deep(ctx, rtcol->geoms[i], rtgeom->bbox);
650     }
651   }
652 }
653 
654 const RTGBOX *
rtgeom_get_bbox(const RTCTX * ctx,const RTGEOM * rtg)655 rtgeom_get_bbox(const RTCTX *ctx, const RTGEOM *rtg)
656 {
657   /* add it if not already there */
658   rtgeom_add_bbox(ctx, (RTGEOM *)rtg);
659   return rtg->bbox;
660 }
661 
662 
663 /**
664 * Calculate the gbox for this goemetry, a cartesian box or
665 * geodetic box, depending on how it is flagged.
666 */
rtgeom_calculate_gbox(const RTCTX * ctx,const RTGEOM * rtgeom,RTGBOX * gbox)667 int rtgeom_calculate_gbox(const RTCTX *ctx, const RTGEOM *rtgeom, RTGBOX *gbox)
668 {
669   gbox->flags = rtgeom->flags;
670   if( RTFLAGS_GET_GEODETIC(rtgeom->flags) )
671     return rtgeom_calculate_gbox_geodetic(ctx, rtgeom, gbox);
672   else
673     return rtgeom_calculate_gbox_cartesian(ctx, rtgeom, gbox);
674 }
675 
676 void
rtgeom_drop_srid(const RTCTX * ctx,RTGEOM * rtgeom)677 rtgeom_drop_srid(const RTCTX *ctx, RTGEOM *rtgeom)
678 {
679   rtgeom->srid = SRID_UNKNOWN;  /* TODO: To be changed to SRID_UNKNOWN */
680 }
681 
682 RTGEOM *
rtgeom_segmentize2d(const RTCTX * ctx,RTGEOM * rtgeom,double dist)683 rtgeom_segmentize2d(const RTCTX *ctx, RTGEOM *rtgeom, double dist)
684 {
685   switch (rtgeom->type)
686   {
687   case RTLINETYPE:
688     return (RTGEOM *)rtline_segmentize2d(ctx, (RTLINE *)rtgeom,
689                                          dist);
690   case RTPOLYGONTYPE:
691     return (RTGEOM *)rtpoly_segmentize2d(ctx, (RTPOLY *)rtgeom,
692                                          dist);
693   case RTMULTILINETYPE:
694   case RTMULTIPOLYGONTYPE:
695   case RTCOLLECTIONTYPE:
696     return (RTGEOM *)rtcollection_segmentize2d(ctx,
697                (RTCOLLECTION *)rtgeom, dist);
698 
699   default:
700     return rtgeom_clone(ctx, rtgeom);
701   }
702 }
703 
704 RTGEOM*
rtgeom_force_2d(const RTCTX * ctx,const RTGEOM * geom)705 rtgeom_force_2d(const RTCTX *ctx, const RTGEOM *geom)
706 {
707   return rtgeom_force_dims(ctx, geom, 0, 0);
708 }
709 
710 RTGEOM*
rtgeom_force_3dz(const RTCTX * ctx,const RTGEOM * geom)711 rtgeom_force_3dz(const RTCTX *ctx, const RTGEOM *geom)
712 {
713   return rtgeom_force_dims(ctx, geom, 1, 0);
714 }
715 
716 RTGEOM*
rtgeom_force_3dm(const RTCTX * ctx,const RTGEOM * geom)717 rtgeom_force_3dm(const RTCTX *ctx, const RTGEOM *geom)
718 {
719   return rtgeom_force_dims(ctx, geom, 0, 1);
720 }
721 
722 RTGEOM*
rtgeom_force_4d(const RTCTX * ctx,const RTGEOM * geom)723 rtgeom_force_4d(const RTCTX *ctx, const RTGEOM *geom)
724 {
725   return rtgeom_force_dims(ctx, geom, 1, 1);
726 }
727 
728 RTGEOM*
rtgeom_force_dims(const RTCTX * ctx,const RTGEOM * geom,int hasz,int hasm)729 rtgeom_force_dims(const RTCTX *ctx, const RTGEOM *geom, int hasz, int hasm)
730 {
731   switch(geom->type)
732   {
733     case RTPOINTTYPE:
734       return rtpoint_as_rtgeom(ctx, rtpoint_force_dims(ctx, (RTPOINT*)geom, hasz, hasm));
735     case RTCIRCSTRINGTYPE:
736     case RTLINETYPE:
737     case RTTRIANGLETYPE:
738       return rtline_as_rtgeom(ctx, rtline_force_dims(ctx, (RTLINE*)geom, hasz, hasm));
739     case RTPOLYGONTYPE:
740       return rtpoly_as_rtgeom(ctx, rtpoly_force_dims(ctx, (RTPOLY*)geom, hasz, hasm));
741     case RTCOMPOUNDTYPE:
742     case RTCURVEPOLYTYPE:
743     case RTMULTICURVETYPE:
744     case RTMULTISURFACETYPE:
745     case RTMULTIPOINTTYPE:
746     case RTMULTILINETYPE:
747     case RTMULTIPOLYGONTYPE:
748     case RTPOLYHEDRALSURFACETYPE:
749     case RTTINTYPE:
750     case RTCOLLECTIONTYPE:
751       return rtcollection_as_rtgeom(ctx, rtcollection_force_dims(ctx, (RTCOLLECTION*)geom, hasz, hasm));
752     default:
753       rterror(ctx, "rtgeom_force_2d: unsupported geom type: %s", rttype_name(ctx, geom->type));
754       return NULL;
755   }
756 }
757 
758 RTGEOM*
rtgeom_force_sfs(const RTCTX * ctx,RTGEOM * geom,int version)759 rtgeom_force_sfs(const RTCTX *ctx, RTGEOM *geom, int version)
760 {
761   RTCOLLECTION *col;
762   int i;
763   RTGEOM *g;
764 
765   /* SFS 1.2 version */
766   if (version == 120)
767   {
768     switch(geom->type)
769     {
770       /* SQL/MM types */
771       case RTCIRCSTRINGTYPE:
772       case RTCOMPOUNDTYPE:
773       case RTCURVEPOLYTYPE:
774       case RTMULTICURVETYPE:
775       case RTMULTISURFACETYPE:
776         return rtgeom_stroke(ctx, geom, 32);
777 
778       case RTCOLLECTIONTYPE:
779         col = (RTCOLLECTION*)geom;
780         for ( i = 0; i < col->ngeoms; i++ )
781           col->geoms[i] = rtgeom_force_sfs(ctx, (RTGEOM*)col->geoms[i], version);
782 
783         return rtcollection_as_rtgeom(ctx, (RTCOLLECTION*)geom);
784 
785       default:
786         return (RTGEOM *)geom;
787     }
788   }
789 
790 
791   /* SFS 1.1 version */
792   switch(geom->type)
793   {
794     /* SQL/MM types */
795     case RTCIRCSTRINGTYPE:
796     case RTCOMPOUNDTYPE:
797     case RTCURVEPOLYTYPE:
798     case RTMULTICURVETYPE:
799     case RTMULTISURFACETYPE:
800       return rtgeom_stroke(ctx, geom, 32);
801 
802     /* SFS 1.2 types */
803     case RTTRIANGLETYPE:
804       g = rtpoly_as_rtgeom(ctx, rtpoly_from_rtlines(ctx, (RTLINE*)geom, 0, NULL));
805       rtgeom_free(ctx, geom);
806       return g;
807 
808     case RTTINTYPE:
809       col = (RTCOLLECTION*) geom;
810       for ( i = 0; i < col->ngeoms; i++ )
811       {
812         g = rtpoly_as_rtgeom(ctx, rtpoly_from_rtlines(ctx, (RTLINE*)col->geoms[i], 0, NULL));
813         rtgeom_free(ctx, col->geoms[i]);
814         col->geoms[i] = g;
815       }
816       col->type = RTCOLLECTIONTYPE;
817       return rtmpoly_as_rtgeom(ctx, (RTMPOLY*)geom);
818 
819     case RTPOLYHEDRALSURFACETYPE:
820       geom->type = RTCOLLECTIONTYPE;
821       return (RTGEOM *)geom;
822 
823     /* Collection */
824     case RTCOLLECTIONTYPE:
825       col = (RTCOLLECTION*)geom;
826       for ( i = 0; i < col->ngeoms; i++ )
827         col->geoms[i] = rtgeom_force_sfs(ctx, (RTGEOM*)col->geoms[i], version);
828 
829       return rtcollection_as_rtgeom(ctx, (RTCOLLECTION*)geom);
830 
831     default:
832       return (RTGEOM *)geom;
833   }
834 }
835 
836 int32_t
rtgeom_get_srid(const RTCTX * ctx,const RTGEOM * geom)837 rtgeom_get_srid(const RTCTX *ctx, const RTGEOM *geom)
838 {
839   if ( ! geom ) return SRID_UNKNOWN;
840   return geom->srid;
841 }
842 
843 uint32_t
rtgeom_get_type(const RTCTX * ctx,const RTGEOM * geom)844 rtgeom_get_type(const RTCTX *ctx, const RTGEOM *geom)
845 {
846   if ( ! geom ) return 0;
847   return geom->type;
848 }
849 
850 int
rtgeom_has_z(const RTCTX * ctx,const RTGEOM * geom)851 rtgeom_has_z(const RTCTX *ctx, const RTGEOM *geom)
852 {
853   if ( ! geom ) return RT_FALSE;
854   return RTFLAGS_GET_Z(geom->flags);
855 }
856 
857 int
rtgeom_has_m(const RTCTX * ctx,const RTGEOM * geom)858 rtgeom_has_m(const RTCTX *ctx, const RTGEOM *geom)
859 {
860   if ( ! geom ) return RT_FALSE;
861   return RTFLAGS_GET_M(geom->flags);
862 }
863 
864 int
rtgeom_ndims(const RTCTX * ctx,const RTGEOM * geom)865 rtgeom_ndims(const RTCTX *ctx, const RTGEOM *geom)
866 {
867   if ( ! geom ) return 0;
868   return RTFLAGS_NDIMS(geom->flags);
869 }
870 
871 
872 void
rtgeom_set_geodetic(const RTCTX * ctx,RTGEOM * geom,int value)873 rtgeom_set_geodetic(const RTCTX *ctx, RTGEOM *geom, int value)
874 {
875   RTPOINT *pt;
876   RTLINE *ln;
877   RTPOLY *ply;
878   RTCOLLECTION *col;
879   int i;
880 
881   RTFLAGS_SET_GEODETIC(geom->flags, value);
882   if ( geom->bbox )
883     RTFLAGS_SET_GEODETIC(geom->bbox->flags, value);
884 
885   switch(geom->type)
886   {
887     case RTPOINTTYPE:
888       pt = (RTPOINT*)geom;
889       if ( pt->point )
890         RTFLAGS_SET_GEODETIC(pt->point->flags, value);
891       break;
892     case RTLINETYPE:
893       ln = (RTLINE*)geom;
894       if ( ln->points )
895         RTFLAGS_SET_GEODETIC(ln->points->flags, value);
896       break;
897     case RTPOLYGONTYPE:
898       ply = (RTPOLY*)geom;
899       for ( i = 0; i < ply->nrings; i++ )
900         RTFLAGS_SET_GEODETIC(ply->rings[i]->flags, value);
901       break;
902     case RTMULTIPOINTTYPE:
903     case RTMULTILINETYPE:
904     case RTMULTIPOLYGONTYPE:
905     case RTCOLLECTIONTYPE:
906       col = (RTCOLLECTION*)geom;
907       for ( i = 0; i < col->ngeoms; i++ )
908         rtgeom_set_geodetic(ctx, col->geoms[i], value);
909       break;
910     default:
911       rterror(ctx, "rtgeom_set_geodetic: unsupported geom type: %s", rttype_name(ctx, geom->type));
912       return;
913   }
914 }
915 
916 void
rtgeom_longitude_shift(const RTCTX * ctx,RTGEOM * rtgeom)917 rtgeom_longitude_shift(const RTCTX *ctx, RTGEOM *rtgeom)
918 {
919   int i;
920   switch (rtgeom->type)
921   {
922     RTPOINT *point;
923     RTLINE *line;
924     RTPOLY *poly;
925     RTTRIANGLE *triangle;
926     RTCOLLECTION *coll;
927 
928   case RTPOINTTYPE:
929     point = (RTPOINT *)rtgeom;
930     ptarray_longitude_shift(ctx, point->point);
931     return;
932   case RTLINETYPE:
933     line = (RTLINE *)rtgeom;
934     ptarray_longitude_shift(ctx, line->points);
935     return;
936   case RTPOLYGONTYPE:
937     poly = (RTPOLY *)rtgeom;
938     for (i=0; i<poly->nrings; i++)
939       ptarray_longitude_shift(ctx, poly->rings[i]);
940     return;
941   case RTTRIANGLETYPE:
942     triangle = (RTTRIANGLE *)rtgeom;
943     ptarray_longitude_shift(ctx, triangle->points);
944     return;
945   case RTMULTIPOINTTYPE:
946   case RTMULTILINETYPE:
947   case RTMULTIPOLYGONTYPE:
948   case RTPOLYHEDRALSURFACETYPE:
949   case RTTINTYPE:
950   case RTCOLLECTIONTYPE:
951     coll = (RTCOLLECTION *)rtgeom;
952     for (i=0; i<coll->ngeoms; i++)
953       rtgeom_longitude_shift(ctx, coll->geoms[i]);
954     return;
955   default:
956     rterror(ctx, "rtgeom_longitude_shift: unsupported geom type: %s",
957             rttype_name(ctx, rtgeom->type));
958   }
959 }
960 
961 int
rtgeom_is_closed(const RTCTX * ctx,const RTGEOM * geom)962 rtgeom_is_closed(const RTCTX *ctx, const RTGEOM *geom)
963 {
964   int type = geom->type;
965 
966   if( rtgeom_is_empty(ctx, geom) )
967     return RT_FALSE;
968 
969   /* Test linear types for closure */
970   switch (type)
971   {
972   case RTLINETYPE:
973     return rtline_is_closed(ctx, (RTLINE*)geom);
974   case RTPOLYGONTYPE:
975     return rtpoly_is_closed(ctx, (RTPOLY*)geom);
976   case RTCIRCSTRINGTYPE:
977     return rtcircstring_is_closed(ctx, (RTCIRCSTRING*)geom);
978   case RTCOMPOUNDTYPE:
979     return rtcompound_is_closed(ctx, (RTCOMPOUND*)geom);
980   case RTTINTYPE:
981     return rttin_is_closed(ctx, (RTTIN*)geom);
982   case RTPOLYHEDRALSURFACETYPE:
983     return rtpsurface_is_closed(ctx, (RTPSURFACE*)geom);
984   }
985 
986   /* Recurse into collections and see if anything is not closed */
987   if ( rtgeom_is_collection(ctx, geom) )
988   {
989     RTCOLLECTION *col = rtgeom_as_rtcollection(ctx, geom);
990     int i;
991     int closed;
992     for ( i = 0; i < col->ngeoms; i++ )
993     {
994       closed = rtgeom_is_closed(ctx, col->geoms[i]);
995       if ( ! closed )
996         return RT_FALSE;
997     }
998     return RT_TRUE;
999   }
1000 
1001   /* All non-linear non-collection types we will call closed */
1002   return RT_TRUE;
1003 }
1004 
1005 int
rtgeom_is_collection(const RTCTX * ctx,const RTGEOM * geom)1006 rtgeom_is_collection(const RTCTX *ctx, const RTGEOM *geom)
1007 {
1008   if( ! geom ) return RT_FALSE;
1009   return rttype_is_collection(ctx, geom->type);
1010 }
1011 
1012 /** Return TRUE if the geometry may contain sub-geometries, i.e. it is a MULTI* or COMPOUNDCURVE */
1013 int
rttype_is_collection(const RTCTX * ctx,uint8_t type)1014 rttype_is_collection(const RTCTX *ctx, uint8_t type)
1015 {
1016 
1017   switch (type)
1018   {
1019   case RTMULTIPOINTTYPE:
1020   case RTMULTILINETYPE:
1021   case RTMULTIPOLYGONTYPE:
1022   case RTCOLLECTIONTYPE:
1023   case RTCURVEPOLYTYPE:
1024   case RTCOMPOUNDTYPE:
1025   case RTMULTICURVETYPE:
1026   case RTMULTISURFACETYPE:
1027   case RTPOLYHEDRALSURFACETYPE:
1028   case RTTINTYPE:
1029     return RT_TRUE;
1030     break;
1031 
1032   default:
1033     return RT_FALSE;
1034   }
1035 }
1036 
1037 /**
1038 * Given an rttype number, what homogeneous collection can hold it?
1039 */
1040 int
rttype_get_collectiontype(const RTCTX * ctx,uint8_t type)1041 rttype_get_collectiontype(const RTCTX *ctx, uint8_t type)
1042 {
1043   switch (type)
1044   {
1045     case RTPOINTTYPE:
1046       return RTMULTIPOINTTYPE;
1047     case RTLINETYPE:
1048       return RTMULTILINETYPE;
1049     case RTPOLYGONTYPE:
1050       return RTMULTIPOLYGONTYPE;
1051     case RTCIRCSTRINGTYPE:
1052       return RTMULTICURVETYPE;
1053     case RTCOMPOUNDTYPE:
1054       return RTMULTICURVETYPE;
1055     case RTCURVEPOLYTYPE:
1056       return RTMULTISURFACETYPE;
1057     case RTTRIANGLETYPE:
1058       return RTTINTYPE;
1059     default:
1060       return RTCOLLECTIONTYPE;
1061   }
1062 }
1063 
1064 
rtgeom_free(const RTCTX * ctx,RTGEOM * rtgeom)1065 void rtgeom_free(const RTCTX *ctx, RTGEOM *rtgeom)
1066 {
1067 
1068   /* There's nothing here to free... */
1069   if( ! rtgeom ) return;
1070 
1071   RTDEBUGF(ctx, 5,"freeing a %s",rttype_name(ctx, rtgeom->type));
1072 
1073   switch (rtgeom->type)
1074   {
1075   case RTPOINTTYPE:
1076     rtpoint_free(ctx, (RTPOINT *)rtgeom);
1077     break;
1078   case RTLINETYPE:
1079     rtline_free(ctx, (RTLINE *)rtgeom);
1080     break;
1081   case RTPOLYGONTYPE:
1082     rtpoly_free(ctx, (RTPOLY *)rtgeom);
1083     break;
1084   case RTCIRCSTRINGTYPE:
1085     rtcircstring_free(ctx, (RTCIRCSTRING *)rtgeom);
1086     break;
1087   case RTTRIANGLETYPE:
1088     rttriangle_free(ctx, (RTTRIANGLE *)rtgeom);
1089     break;
1090   case RTMULTIPOINTTYPE:
1091     rtmpoint_free(ctx, (RTMPOINT *)rtgeom);
1092     break;
1093   case RTMULTILINETYPE:
1094     rtmline_free(ctx, (RTMLINE *)rtgeom);
1095     break;
1096   case RTMULTIPOLYGONTYPE:
1097     rtmpoly_free(ctx, (RTMPOLY *)rtgeom);
1098     break;
1099   case RTPOLYHEDRALSURFACETYPE:
1100     rtpsurface_free(ctx, (RTPSURFACE *)rtgeom);
1101     break;
1102   case RTTINTYPE:
1103     rttin_free(ctx, (RTTIN *)rtgeom);
1104     break;
1105   case RTCURVEPOLYTYPE:
1106   case RTCOMPOUNDTYPE:
1107   case RTMULTICURVETYPE:
1108   case RTMULTISURFACETYPE:
1109   case RTCOLLECTIONTYPE:
1110     rtcollection_free(ctx, (RTCOLLECTION *)rtgeom);
1111     break;
1112   default:
1113     rterror(ctx, "rtgeom_free called with unknown type (%d) %s", rtgeom->type, rttype_name(ctx, rtgeom->type));
1114   }
1115   return;
1116 }
1117 
rtgeom_needs_bbox(const RTCTX * ctx,const RTGEOM * geom)1118 int rtgeom_needs_bbox(const RTCTX *ctx, const RTGEOM *geom)
1119 {
1120   assert(geom);
1121   if ( geom->type == RTPOINTTYPE )
1122   {
1123     return RT_FALSE;
1124   }
1125   else if ( geom->type == RTLINETYPE )
1126   {
1127     if ( rtgeom_count_vertices(ctx, geom) <= 2 )
1128       return RT_FALSE;
1129     else
1130       return RT_TRUE;
1131   }
1132   else if ( geom->type == RTMULTIPOINTTYPE )
1133   {
1134     if ( ((RTCOLLECTION*)geom)->ngeoms == 1 )
1135       return RT_FALSE;
1136     else
1137       return RT_TRUE;
1138   }
1139   else if ( geom->type == RTMULTILINETYPE )
1140   {
1141     if ( ((RTCOLLECTION*)geom)->ngeoms == 1 && rtgeom_count_vertices(ctx, geom) <= 2 )
1142       return RT_FALSE;
1143     else
1144       return RT_TRUE;
1145   }
1146   else
1147   {
1148     return RT_TRUE;
1149   }
1150 }
1151 
1152 /**
1153 * Count points in an #RTGEOM.
1154 */
rtgeom_count_vertices(const RTCTX * ctx,const RTGEOM * geom)1155 int rtgeom_count_vertices(const RTCTX *ctx, const RTGEOM *geom)
1156 {
1157   int result = 0;
1158 
1159   /* Null? Zero. */
1160   if( ! geom ) return 0;
1161 
1162   RTDEBUGF(ctx, 4, "rtgeom_count_vertices got type %s",
1163            rttype_name(ctx, geom->type));
1164 
1165   /* Empty? Zero. */
1166   if( rtgeom_is_empty(ctx, geom) ) return 0;
1167 
1168   switch (geom->type)
1169   {
1170   case RTPOINTTYPE:
1171     result = 1;
1172     break;
1173   case RTTRIANGLETYPE:
1174   case RTCIRCSTRINGTYPE:
1175   case RTLINETYPE:
1176     result = rtline_count_vertices(ctx, (RTLINE *)geom);
1177     break;
1178   case RTPOLYGONTYPE:
1179     result = rtpoly_count_vertices(ctx, (RTPOLY *)geom);
1180     break;
1181   case RTCOMPOUNDTYPE:
1182   case RTCURVEPOLYTYPE:
1183   case RTMULTICURVETYPE:
1184   case RTMULTISURFACETYPE:
1185   case RTMULTIPOINTTYPE:
1186   case RTMULTILINETYPE:
1187   case RTMULTIPOLYGONTYPE:
1188   case RTPOLYHEDRALSURFACETYPE:
1189   case RTTINTYPE:
1190   case RTCOLLECTIONTYPE:
1191     result = rtcollection_count_vertices(ctx, (RTCOLLECTION *)geom);
1192     break;
1193   default:
1194     rterror(ctx, "%s: unsupported input geometry type: %s",
1195             __func__, rttype_name(ctx, geom->type));
1196     break;
1197   }
1198   RTDEBUGF(ctx, 3, "counted %d vertices", result);
1199   return result;
1200 }
1201 
1202 /**
1203 * For an #RTGEOM, returns 0 for points, 1 for lines,
1204 * 2 for polygons, 3 for volume, and the max dimension
1205 * of a collection.
1206 */
rtgeom_dimension(const RTCTX * ctx,const RTGEOM * geom)1207 int rtgeom_dimension(const RTCTX *ctx, const RTGEOM *geom)
1208 {
1209 
1210   /* Null? Zero. */
1211   if( ! geom ) return -1;
1212 
1213   RTDEBUGF(ctx, 4, "rtgeom_dimension got type %s",
1214            rttype_name(ctx, geom->type));
1215 
1216   /* Empty? Zero. */
1217   /* if( rtgeom_is_empty(ctx, geom) ) return 0; */
1218 
1219   switch (geom->type)
1220   {
1221   case RTPOINTTYPE:
1222   case RTMULTIPOINTTYPE:
1223     return 0;
1224   case RTCIRCSTRINGTYPE:
1225   case RTLINETYPE:
1226   case RTCOMPOUNDTYPE:
1227   case RTMULTICURVETYPE:
1228   case RTMULTILINETYPE:
1229     return 1;
1230   case RTTRIANGLETYPE:
1231   case RTPOLYGONTYPE:
1232   case RTCURVEPOLYTYPE:
1233   case RTMULTISURFACETYPE:
1234   case RTMULTIPOLYGONTYPE:
1235   case RTTINTYPE:
1236     return 2;
1237   case RTPOLYHEDRALSURFACETYPE:
1238   {
1239     /* A closed polyhedral surface contains a volume. */
1240     int closed = rtpsurface_is_closed(ctx, (RTPSURFACE*)geom);
1241     return ( closed ? 3 : 2 );
1242   }
1243   case RTCOLLECTIONTYPE:
1244   {
1245     int maxdim = 0, i;
1246     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1247     for( i = 0; i < col->ngeoms; i++ )
1248     {
1249       int dim = rtgeom_dimension(ctx, col->geoms[i]);
1250       maxdim = ( dim > maxdim ? dim : maxdim );
1251     }
1252     return maxdim;
1253   }
1254   default:
1255     rterror(ctx, "%s: unsupported input geometry type: %s",
1256             __func__, rttype_name(ctx, geom->type));
1257   }
1258   return -1;
1259 }
1260 
1261 /**
1262 * Count rings in an #RTGEOM.
1263 */
rtgeom_count_rings(const RTCTX * ctx,const RTGEOM * geom)1264 int rtgeom_count_rings(const RTCTX *ctx, const RTGEOM *geom)
1265 {
1266   int result = 0;
1267 
1268   /* Null? Empty? Zero. */
1269   if( ! geom || rtgeom_is_empty(ctx, geom) )
1270     return 0;
1271 
1272   switch (geom->type)
1273   {
1274   case RTPOINTTYPE:
1275   case RTCIRCSTRINGTYPE:
1276   case RTCOMPOUNDTYPE:
1277   case RTMULTICURVETYPE:
1278   case RTMULTIPOINTTYPE:
1279   case RTMULTILINETYPE:
1280   case RTLINETYPE:
1281     result = 0;
1282     break;
1283   case RTTRIANGLETYPE:
1284     result = 1;
1285     break;
1286   case RTPOLYGONTYPE:
1287     result = ((RTPOLY *)geom)->nrings;
1288     break;
1289   case RTCURVEPOLYTYPE:
1290     result = ((RTCURVEPOLY *)geom)->nrings;
1291     break;
1292   case RTMULTISURFACETYPE:
1293   case RTMULTIPOLYGONTYPE:
1294   case RTPOLYHEDRALSURFACETYPE:
1295   case RTTINTYPE:
1296   case RTCOLLECTIONTYPE:
1297   {
1298     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1299     int i = 0;
1300     for( i = 0; i < col->ngeoms; i++ )
1301       result += rtgeom_count_rings(ctx, col->geoms[i]);
1302     break;
1303   }
1304   default:
1305     rterror(ctx, "rtgeom_count_rings: unsupported input geometry type: %s", rttype_name(ctx, geom->type));
1306     break;
1307   }
1308   RTDEBUGF(ctx, 3, "counted %d rings", result);
1309   return result;
1310 }
1311 
rtgeom_is_empty(const RTCTX * ctx,const RTGEOM * geom)1312 int rtgeom_is_empty(const RTCTX *ctx, const RTGEOM *geom)
1313 {
1314   int result = RT_FALSE;
1315   RTDEBUGF(ctx, 4, "rtgeom_is_empty: got type %s",
1316            rttype_name(ctx, geom->type));
1317 
1318   switch (geom->type)
1319   {
1320   case RTPOINTTYPE:
1321     return rtpoint_is_empty(ctx, (RTPOINT*)geom);
1322     break;
1323   case RTLINETYPE:
1324     return rtline_is_empty(ctx, (RTLINE*)geom);
1325     break;
1326   case RTCIRCSTRINGTYPE:
1327     return rtcircstring_is_empty(ctx, (RTCIRCSTRING*)geom);
1328     break;
1329   case RTPOLYGONTYPE:
1330     return rtpoly_is_empty(ctx, (RTPOLY*)geom);
1331     break;
1332   case RTTRIANGLETYPE:
1333     return rttriangle_is_empty(ctx, (RTTRIANGLE*)geom);
1334     break;
1335   case RTMULTIPOINTTYPE:
1336   case RTMULTILINETYPE:
1337   case RTMULTIPOLYGONTYPE:
1338   case RTCOMPOUNDTYPE:
1339   case RTCURVEPOLYTYPE:
1340   case RTMULTICURVETYPE:
1341   case RTMULTISURFACETYPE:
1342   case RTPOLYHEDRALSURFACETYPE:
1343   case RTTINTYPE:
1344   case RTCOLLECTIONTYPE:
1345     return rtcollection_is_empty(ctx, (RTCOLLECTION *)geom);
1346     break;
1347   default:
1348     rterror(ctx, "rtgeom_is_empty: unsupported input geometry type: %s",
1349             rttype_name(ctx, geom->type));
1350     break;
1351   }
1352   return result;
1353 }
1354 
rtgeom_has_srid(const RTCTX * ctx,const RTGEOM * geom)1355 int rtgeom_has_srid(const RTCTX *ctx, const RTGEOM *geom)
1356 {
1357   if ( geom->srid != SRID_UNKNOWN )
1358     return RT_TRUE;
1359 
1360   return RT_FALSE;
1361 }
1362 
1363 
rtcollection_dimensionality(const RTCTX * ctx,RTCOLLECTION * col)1364 static int rtcollection_dimensionality(const RTCTX *ctx, RTCOLLECTION *col)
1365 {
1366   int i;
1367   int dimensionality = 0;
1368   for ( i = 0; i < col->ngeoms; i++ )
1369   {
1370     int d = rtgeom_dimensionality(ctx, col->geoms[i]);
1371     if ( d > dimensionality )
1372       dimensionality = d;
1373   }
1374   return dimensionality;
1375 }
1376 
rtgeom_dimensionality(const RTCTX * ctx,RTGEOM * geom)1377 extern int rtgeom_dimensionality(const RTCTX *ctx, RTGEOM *geom)
1378 {
1379   int dim;
1380 
1381   RTDEBUGF(ctx, 3, "rtgeom_dimensionality got type %s",
1382            rttype_name(ctx, geom->type));
1383 
1384   switch (geom->type)
1385   {
1386   case RTPOINTTYPE:
1387   case RTMULTIPOINTTYPE:
1388     return 0;
1389     break;
1390   case RTLINETYPE:
1391   case RTCIRCSTRINGTYPE:
1392   case RTMULTILINETYPE:
1393   case RTCOMPOUNDTYPE:
1394   case RTMULTICURVETYPE:
1395     return 1;
1396     break;
1397   case RTPOLYGONTYPE:
1398   case RTTRIANGLETYPE:
1399   case RTCURVEPOLYTYPE:
1400   case RTMULTIPOLYGONTYPE:
1401   case RTMULTISURFACETYPE:
1402     return 2;
1403     break;
1404 
1405   case RTPOLYHEDRALSURFACETYPE:
1406   case RTTINTYPE:
1407     dim = rtgeom_is_closed(ctx, geom)?3:2;
1408     return dim;
1409     break;
1410 
1411   case RTCOLLECTIONTYPE:
1412     return rtcollection_dimensionality(ctx, (RTCOLLECTION *)geom);
1413     break;
1414   default:
1415     rterror(ctx, "rtgeom_dimensionality: unsupported input geometry type: %s",
1416             rttype_name(ctx, geom->type));
1417     break;
1418   }
1419   return 0;
1420 }
1421 
rtgeom_remove_repeated_points(const RTCTX * ctx,const RTGEOM * in,double tolerance)1422 extern RTGEOM* rtgeom_remove_repeated_points(const RTCTX *ctx, const RTGEOM *in, double tolerance)
1423 {
1424   RTDEBUGF(ctx, 4, "rtgeom_remove_repeated_points got type %s",
1425            rttype_name(ctx, in->type));
1426 
1427   if(rtgeom_is_empty(ctx, in))
1428   {
1429     return rtgeom_clone_deep(ctx, in);
1430   }
1431 
1432   switch (in->type)
1433   {
1434   case RTMULTIPOINTTYPE:
1435     return rtmpoint_remove_repeated_points(ctx, (RTMPOINT*)in, tolerance);
1436     break;
1437   case RTLINETYPE:
1438     return rtline_remove_repeated_points(ctx, (RTLINE*)in, tolerance);
1439 
1440   case RTMULTILINETYPE:
1441   case RTCOLLECTIONTYPE:
1442   case RTMULTIPOLYGONTYPE:
1443   case RTPOLYHEDRALSURFACETYPE:
1444     return rtcollection_remove_repeated_points(ctx, (RTCOLLECTION *)in, tolerance);
1445 
1446   case RTPOLYGONTYPE:
1447     return rtpoly_remove_repeated_points(ctx, (RTPOLY *)in, tolerance);
1448     break;
1449 
1450   case RTPOINTTYPE:
1451   case RTTRIANGLETYPE:
1452   case RTTINTYPE:
1453     /* No point is repeated for a single point, or for Triangle or TIN */
1454     return rtgeom_clone_deep(ctx, in);
1455 
1456   case RTCIRCSTRINGTYPE:
1457   case RTCOMPOUNDTYPE:
1458   case RTMULTICURVETYPE:
1459   case RTCURVEPOLYTYPE:
1460   case RTMULTISURFACETYPE:
1461     /* Dunno how to handle these, will return untouched */
1462     return rtgeom_clone_deep(ctx, in);
1463 
1464   default:
1465     rtnotice(ctx, "%s: unsupported geometry type: %s",
1466              __func__, rttype_name(ctx, in->type));
1467     return rtgeom_clone_deep(ctx, in);
1468     break;
1469   }
1470   return 0;
1471 }
1472 
rtgeom_flip_coordinates(const RTCTX * ctx,RTGEOM * in)1473 RTGEOM* rtgeom_flip_coordinates(const RTCTX *ctx, RTGEOM *in)
1474 {
1475   rtgeom_swap_ordinates(ctx, in,RTORD_X,RTORD_Y);
1476   return in;
1477 }
1478 
rtgeom_swap_ordinates(const RTCTX * ctx,RTGEOM * in,RTORD o1,RTORD o2)1479 void rtgeom_swap_ordinates(const RTCTX *ctx, RTGEOM *in, RTORD o1, RTORD o2)
1480 {
1481   RTCOLLECTION *col;
1482   RTPOLY *poly;
1483   int i;
1484 
1485 #if PARANOIA_LEVEL > 0
1486   assert(o1 < 4);
1487   assert(o2 < 4);
1488 #endif
1489 
1490   if ( (!in) || rtgeom_is_empty(ctx, in) ) return;
1491 
1492   /* TODO: check for rtgeom NOT having the specified dimension ? */
1493 
1494   RTDEBUGF(ctx, 4, "rtgeom_flip_coordinates, got type: %s",
1495            rttype_name(ctx, in->type));
1496 
1497   switch (in->type)
1498   {
1499   case RTPOINTTYPE:
1500     ptarray_swap_ordinates(ctx, rtgeom_as_rtpoint(ctx, in)->point, o1, o2);
1501     break;
1502 
1503   case RTLINETYPE:
1504     ptarray_swap_ordinates(ctx, rtgeom_as_rtline(ctx, in)->points, o1, o2);
1505     break;
1506 
1507   case RTCIRCSTRINGTYPE:
1508     ptarray_swap_ordinates(ctx, rtgeom_as_rtcircstring(ctx, in)->points, o1, o2);
1509     break;
1510 
1511   case RTPOLYGONTYPE:
1512     poly = (RTPOLY *) in;
1513     for (i=0; i<poly->nrings; i++)
1514     {
1515       ptarray_swap_ordinates(ctx, poly->rings[i], o1, o2);
1516     }
1517     break;
1518 
1519   case RTTRIANGLETYPE:
1520     ptarray_swap_ordinates(ctx, rtgeom_as_rttriangle(ctx, in)->points, o1, o2);
1521     break;
1522 
1523   case RTMULTIPOINTTYPE:
1524   case RTMULTILINETYPE:
1525   case RTMULTIPOLYGONTYPE:
1526   case RTCOLLECTIONTYPE:
1527   case RTCOMPOUNDTYPE:
1528   case RTCURVEPOLYTYPE:
1529   case RTMULTISURFACETYPE:
1530   case RTMULTICURVETYPE:
1531   case RTPOLYHEDRALSURFACETYPE:
1532   case RTTINTYPE:
1533     col = (RTCOLLECTION *) in;
1534     for (i=0; i<col->ngeoms; i++)
1535     {
1536       rtgeom_swap_ordinates(ctx, col->geoms[i], o1, o2);
1537     }
1538     break;
1539 
1540   default:
1541     rterror(ctx, "rtgeom_swap_ordinates: unsupported geometry type: %s",
1542             rttype_name(ctx, in->type));
1543     return;
1544   }
1545 
1546   /* only refresh bbox if X or Y changed */
1547   if ( in->bbox && (o1 < 2 || o2 < 2) )
1548   {
1549     rtgeom_drop_bbox(ctx, in);
1550     rtgeom_add_bbox(ctx, in);
1551   }
1552 }
1553 
rtgeom_set_srid(const RTCTX * ctx,RTGEOM * geom,int32_t srid)1554 void rtgeom_set_srid(const RTCTX *ctx, RTGEOM *geom, int32_t srid)
1555 {
1556   int i;
1557 
1558   RTDEBUGF(ctx, 4,"entered with srid=%d",srid);
1559 
1560   geom->srid = srid;
1561 
1562   if ( rtgeom_is_collection(ctx, geom) )
1563   {
1564     /* All the children are set to the same SRID value */
1565     RTCOLLECTION *col = rtgeom_as_rtcollection(ctx, geom);
1566     for ( i = 0; i < col->ngeoms; i++ )
1567     {
1568       rtgeom_set_srid(ctx, col->geoms[i], srid);
1569     }
1570   }
1571 }
1572 
rtgeom_simplify(const RTCTX * ctx,const RTGEOM * igeom,double dist,int preserve_collapsed)1573 RTGEOM* rtgeom_simplify(const RTCTX *ctx, const RTGEOM *igeom, double dist, int preserve_collapsed)
1574 {
1575   switch (igeom->type)
1576   {
1577   case RTPOINTTYPE:
1578   case RTMULTIPOINTTYPE:
1579     return rtgeom_clone(ctx, igeom);
1580   case RTLINETYPE:
1581     return (RTGEOM*)rtline_simplify(ctx, (RTLINE*)igeom, dist, preserve_collapsed);
1582   case RTPOLYGONTYPE:
1583     return (RTGEOM*)rtpoly_simplify(ctx, (RTPOLY*)igeom, dist, preserve_collapsed);
1584   case RTMULTILINETYPE:
1585   case RTMULTIPOLYGONTYPE:
1586   case RTCOLLECTIONTYPE:
1587     return (RTGEOM*)rtcollection_simplify(ctx, (RTCOLLECTION *)igeom, dist, preserve_collapsed);
1588   default:
1589     rterror(ctx, "%s: unsupported geometry type: %s", __func__, rttype_name(ctx, igeom->type));
1590   }
1591   return NULL;
1592 }
1593 
rtgeom_area(const RTCTX * ctx,const RTGEOM * geom)1594 double rtgeom_area(const RTCTX *ctx, const RTGEOM *geom)
1595 {
1596   int type = geom->type;
1597 
1598   if ( type == RTPOLYGONTYPE )
1599     return rtpoly_area(ctx, (RTPOLY*)geom);
1600   else if ( type == RTCURVEPOLYTYPE )
1601     return rtcurvepoly_area(ctx, (RTCURVEPOLY*)geom);
1602   else if (type ==  RTTRIANGLETYPE )
1603     return rttriangle_area(ctx, (RTTRIANGLE*)geom);
1604   else if ( rtgeom_is_collection(ctx, geom) )
1605   {
1606     double area = 0.0;
1607     int i;
1608     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1609     for ( i = 0; i < col->ngeoms; i++ )
1610       area += rtgeom_area(ctx, col->geoms[i]);
1611     return area;
1612   }
1613   else
1614     return 0.0;
1615 }
1616 
rtgeom_perimeter(const RTCTX * ctx,const RTGEOM * geom)1617 double rtgeom_perimeter(const RTCTX *ctx, const RTGEOM *geom)
1618 {
1619   int type = geom->type;
1620   if ( type == RTPOLYGONTYPE )
1621     return rtpoly_perimeter(ctx, (RTPOLY*)geom);
1622   else if ( type == RTCURVEPOLYTYPE )
1623     return rtcurvepoly_perimeter(ctx, (RTCURVEPOLY*)geom);
1624   else if ( type == RTTRIANGLETYPE )
1625     return rttriangle_perimeter(ctx, (RTTRIANGLE*)geom);
1626   else if ( rtgeom_is_collection(ctx, geom) )
1627   {
1628     double perimeter = 0.0;
1629     int i;
1630     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1631     for ( i = 0; i < col->ngeoms; i++ )
1632       perimeter += rtgeom_perimeter(ctx, col->geoms[i]);
1633     return perimeter;
1634   }
1635   else
1636     return 0.0;
1637 }
1638 
rtgeom_perimeter_2d(const RTCTX * ctx,const RTGEOM * geom)1639 double rtgeom_perimeter_2d(const RTCTX *ctx, const RTGEOM *geom)
1640 {
1641   int type = geom->type;
1642   if ( type == RTPOLYGONTYPE )
1643     return rtpoly_perimeter_2d(ctx, (RTPOLY*)geom);
1644   else if ( type == RTCURVEPOLYTYPE )
1645     return rtcurvepoly_perimeter_2d(ctx, (RTCURVEPOLY*)geom);
1646   else if ( type == RTTRIANGLETYPE )
1647     return rttriangle_perimeter_2d(ctx, (RTTRIANGLE*)geom);
1648   else if ( rtgeom_is_collection(ctx, geom) )
1649   {
1650     double perimeter = 0.0;
1651     int i;
1652     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1653     for ( i = 0; i < col->ngeoms; i++ )
1654       perimeter += rtgeom_perimeter_2d(ctx, col->geoms[i]);
1655     return perimeter;
1656   }
1657   else
1658     return 0.0;
1659 }
1660 
rtgeom_length(const RTCTX * ctx,const RTGEOM * geom)1661 double rtgeom_length(const RTCTX *ctx, const RTGEOM *geom)
1662 {
1663   int type = geom->type;
1664   if ( type == RTLINETYPE )
1665     return rtline_length(ctx, (RTLINE*)geom);
1666   else if ( type == RTCIRCSTRINGTYPE )
1667     return rtcircstring_length(ctx, (RTCIRCSTRING*)geom);
1668   else if ( type == RTCOMPOUNDTYPE )
1669     return rtcompound_length(ctx, (RTCOMPOUND*)geom);
1670   else if ( rtgeom_is_collection(ctx, geom) )
1671   {
1672     double length = 0.0;
1673     int i;
1674     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1675     for ( i = 0; i < col->ngeoms; i++ )
1676       length += rtgeom_length(ctx, col->geoms[i]);
1677     return length;
1678   }
1679   else
1680     return 0.0;
1681 }
1682 
rtgeom_length_2d(const RTCTX * ctx,const RTGEOM * geom)1683 double rtgeom_length_2d(const RTCTX *ctx, const RTGEOM *geom)
1684 {
1685   int type = geom->type;
1686   if ( type == RTLINETYPE )
1687     return rtline_length_2d(ctx, (RTLINE*)geom);
1688   else if ( type == RTCIRCSTRINGTYPE )
1689     return rtcircstring_length_2d(ctx, (RTCIRCSTRING*)geom);
1690   else if ( type == RTCOMPOUNDTYPE )
1691     return rtcompound_length_2d(ctx, (RTCOMPOUND*)geom);
1692   else if ( rtgeom_is_collection(ctx, geom) )
1693   {
1694     double length = 0.0;
1695     int i;
1696     RTCOLLECTION *col = (RTCOLLECTION*)geom;
1697     for ( i = 0; i < col->ngeoms; i++ )
1698       length += rtgeom_length_2d(ctx, col->geoms[i]);
1699     return length;
1700   }
1701   else
1702     return 0.0;
1703 }
1704 
1705 void
rtgeom_affine(const RTCTX * ctx,RTGEOM * geom,const RTAFFINE * affine)1706 rtgeom_affine(const RTCTX *ctx, RTGEOM *geom, const RTAFFINE *affine)
1707 {
1708   int type = geom->type;
1709   int i;
1710 
1711   switch(type)
1712   {
1713     /* Take advantage of fact tht pt/ln/circ/tri have same memory structure */
1714     case RTPOINTTYPE:
1715     case RTLINETYPE:
1716     case RTCIRCSTRINGTYPE:
1717     case RTTRIANGLETYPE:
1718     {
1719       RTLINE *l = (RTLINE*)geom;
1720       ptarray_affine(ctx, l->points, affine);
1721       break;
1722     }
1723     case RTPOLYGONTYPE:
1724     {
1725       RTPOLY *p = (RTPOLY*)geom;
1726       for( i = 0; i < p->nrings; i++ )
1727         ptarray_affine(ctx, p->rings[i], affine);
1728       break;
1729     }
1730     case RTCURVEPOLYTYPE:
1731     {
1732       RTCURVEPOLY *c = (RTCURVEPOLY*)geom;
1733       for( i = 0; i < c->nrings; i++ )
1734         rtgeom_affine(ctx, c->rings[i], affine);
1735       break;
1736     }
1737     default:
1738     {
1739       if( rtgeom_is_collection(ctx, geom) )
1740       {
1741         RTCOLLECTION *c = (RTCOLLECTION*)geom;
1742         for( i = 0; i < c->ngeoms; i++ )
1743         {
1744           rtgeom_affine(ctx, c->geoms[i], affine);
1745         }
1746       }
1747       else
1748       {
1749         rterror(ctx, "rtgeom_affine: unable to handle type '%s'", rttype_name(ctx, type));
1750       }
1751     }
1752   }
1753 
1754 }
1755 
1756 void
rtgeom_scale(const RTCTX * ctx,RTGEOM * geom,const RTPOINT4D * factor)1757 rtgeom_scale(const RTCTX *ctx, RTGEOM *geom, const RTPOINT4D *factor)
1758 {
1759   int type = geom->type;
1760   int i;
1761 
1762   switch(type)
1763   {
1764     /* Take advantage of fact tht pt/ln/circ/tri have same memory structure */
1765     case RTPOINTTYPE:
1766     case RTLINETYPE:
1767     case RTCIRCSTRINGTYPE:
1768     case RTTRIANGLETYPE:
1769     {
1770       RTLINE *l = (RTLINE*)geom;
1771       ptarray_scale(ctx, l->points, factor);
1772       break;
1773     }
1774     case RTPOLYGONTYPE:
1775     {
1776       RTPOLY *p = (RTPOLY*)geom;
1777       for( i = 0; i < p->nrings; i++ )
1778         ptarray_scale(ctx, p->rings[i], factor);
1779       break;
1780     }
1781     case RTCURVEPOLYTYPE:
1782     {
1783       RTCURVEPOLY *c = (RTCURVEPOLY*)geom;
1784       for( i = 0; i < c->nrings; i++ )
1785         rtgeom_scale(ctx, c->rings[i], factor);
1786       break;
1787     }
1788     default:
1789     {
1790       if( rtgeom_is_collection(ctx, geom) )
1791       {
1792         RTCOLLECTION *c = (RTCOLLECTION*)geom;
1793         for( i = 0; i < c->ngeoms; i++ )
1794         {
1795           rtgeom_scale(ctx, c->geoms[i], factor);
1796         }
1797       }
1798       else
1799       {
1800         rterror(ctx, "rtgeom_scale: unable to handle type '%s'", rttype_name(ctx, type));
1801       }
1802     }
1803   }
1804 
1805   /* Recompute bbox if needed */
1806 
1807   if ( geom->bbox )
1808   {
1809     /* TODO: expose a gbox_scale function */
1810     geom->bbox->xmin *= factor->x;
1811     geom->bbox->xmax *= factor->x;
1812     geom->bbox->ymin *= factor->y;
1813     geom->bbox->ymax *= factor->y;
1814     geom->bbox->zmin *= factor->z;
1815     geom->bbox->zmax *= factor->z;
1816     geom->bbox->mmin *= factor->m;
1817     geom->bbox->mmax *= factor->m;
1818   }
1819 }
1820 
1821 RTGEOM*
rtgeom_construct_empty(const RTCTX * ctx,uint8_t type,int srid,char hasz,char hasm)1822 rtgeom_construct_empty(const RTCTX *ctx, uint8_t type, int srid, char hasz, char hasm)
1823 {
1824   switch(type)
1825   {
1826     case RTPOINTTYPE:
1827       return rtpoint_as_rtgeom(ctx, rtpoint_construct_empty(ctx, srid, hasz, hasm));
1828     case RTLINETYPE:
1829       return rtline_as_rtgeom(ctx, rtline_construct_empty(ctx, srid, hasz, hasm));
1830     case RTPOLYGONTYPE:
1831       return rtpoly_as_rtgeom(ctx, rtpoly_construct_empty(ctx, srid, hasz, hasm));
1832     case RTCURVEPOLYTYPE:
1833       return rtcurvepoly_as_rtgeom(ctx, rtcurvepoly_construct_empty(ctx, srid, hasz, hasm));
1834     case RTCIRCSTRINGTYPE:
1835       return rtcircstring_as_rtgeom(ctx, rtcircstring_construct_empty(ctx, srid, hasz, hasm));
1836     case RTTRIANGLETYPE:
1837       return rttriangle_as_rtgeom(ctx, rttriangle_construct_empty(ctx, srid, hasz, hasm));
1838     case RTCOMPOUNDTYPE:
1839     case RTMULTIPOINTTYPE:
1840     case RTMULTILINETYPE:
1841     case RTMULTIPOLYGONTYPE:
1842     case RTCOLLECTIONTYPE:
1843       return rtcollection_as_rtgeom(ctx, rtcollection_construct_empty(ctx, type, srid, hasz, hasm));
1844     default:
1845       rterror(ctx, "rtgeom_construct_empty: unsupported geometry type: %s",
1846               rttype_name(ctx, type));
1847       return NULL;
1848   }
1849 }
1850 
1851 int
rtgeom_startpoint(const RTCTX * ctx,const RTGEOM * rtgeom,RTPOINT4D * pt)1852 rtgeom_startpoint(const RTCTX *ctx, const RTGEOM* rtgeom, RTPOINT4D* pt)
1853 {
1854   if ( ! rtgeom )
1855     return RT_FAILURE;
1856 
1857   switch( rtgeom->type )
1858   {
1859     case RTPOINTTYPE:
1860       return ptarray_startpoint(ctx, ((RTPOINT*)rtgeom)->point, pt);
1861     case RTTRIANGLETYPE:
1862     case RTCIRCSTRINGTYPE:
1863     case RTLINETYPE:
1864       return ptarray_startpoint(ctx, ((RTLINE*)rtgeom)->points, pt);
1865     case RTPOLYGONTYPE:
1866       return rtpoly_startpoint(ctx, (RTPOLY*)rtgeom, pt);
1867     case RTCURVEPOLYTYPE:
1868     case RTCOMPOUNDTYPE:
1869     case RTMULTIPOINTTYPE:
1870     case RTMULTILINETYPE:
1871     case RTMULTIPOLYGONTYPE:
1872     case RTCOLLECTIONTYPE:
1873       return rtcollection_startpoint(ctx, (RTCOLLECTION*)rtgeom, pt);
1874     default:
1875       rterror(ctx, "int: unsupported geometry type: %s",
1876               rttype_name(ctx, rtgeom->type));
1877       return RT_FAILURE;
1878   }
1879 }
1880 
1881 
1882 RTGEOM *
rtgeom_grid(const RTCTX * ctx,const RTGEOM * rtgeom,const gridspec * grid)1883 rtgeom_grid(const RTCTX *ctx, const RTGEOM *rtgeom, const gridspec *grid)
1884 {
1885   switch ( rtgeom->type )
1886   {
1887     case RTPOINTTYPE:
1888       return (RTGEOM *)rtpoint_grid(ctx, (RTPOINT *)rtgeom, grid);
1889     case RTLINETYPE:
1890       return (RTGEOM *)rtline_grid(ctx, (RTLINE *)rtgeom, grid);
1891     case RTPOLYGONTYPE:
1892       return (RTGEOM *)rtpoly_grid(ctx, (RTPOLY *)rtgeom, grid);
1893     case RTMULTIPOINTTYPE:
1894     case RTMULTILINETYPE:
1895     case RTMULTIPOLYGONTYPE:
1896     case RTCOLLECTIONTYPE:
1897     case RTCOMPOUNDTYPE:
1898       return (RTGEOM *)rtcollection_grid(ctx, (RTCOLLECTION *)rtgeom, grid);
1899     case RTCIRCSTRINGTYPE:
1900       return (RTGEOM *)rtcircstring_grid(ctx, (RTCIRCSTRING *)rtgeom, grid);
1901     default:
1902       rterror(ctx, "rtgeom_grid: Unsupported geometry type: %s",
1903               rttype_name(ctx, rtgeom->type));
1904       return NULL;
1905   }
1906 }
1907 
1908 
1909 /* Prototype for recursion */
1910 static int
1911 rtgeom_subdivide_recursive(const RTCTX *ctx, const RTGEOM *geom, int maxvertices, int depth, RTCOLLECTION *col, const RTGBOX *clip);
1912 
1913 static int
rtgeom_subdivide_recursive(const RTCTX * ctx,const RTGEOM * geom,int maxvertices,int depth,RTCOLLECTION * col,const RTGBOX * clip)1914 rtgeom_subdivide_recursive(const RTCTX *ctx, const RTGEOM *geom, int maxvertices, int depth, RTCOLLECTION *col, const RTGBOX *clip)
1915 {
1916   const int maxdepth = 50;
1917   int nvertices = 0;
1918   int i, n = 0;
1919   double width = clip->xmax - clip->xmin;
1920   double height = clip->ymax - clip->ymin;
1921   RTGBOX subbox1, subbox2;
1922   RTGEOM *clipped1, *clipped2;
1923 
1924   if ( geom->type == RTPOLYHEDRALSURFACETYPE || geom->type == RTTINTYPE )
1925   {
1926     rterror(ctx, "%s: unsupported geometry type '%s'", __func__, rttype_name(ctx, geom->type));
1927   }
1928 
1929   if ( width == 0.0 && height == 0.0 )
1930     return 0;
1931 
1932   /* Artays just recurse into collections */
1933   if ( rtgeom_is_collection(ctx, geom) )
1934   {
1935     RTCOLLECTION *incol = (RTCOLLECTION*)geom;
1936     int n = 0;
1937     for ( i = 0; i < incol->ngeoms; i++ )
1938     {
1939       /* Don't increment depth yet, since we aren't actually subdividing geomtries yet */
1940       n += rtgeom_subdivide_recursive(ctx, incol->geoms[i], maxvertices, depth, col, clip);
1941     }
1942     return n;
1943   }
1944 
1945   /* But don't go too far. 2^25 = 33M, that's enough subdivision */
1946   /* Signal callers above that we depth'ed out with a negative */
1947   /* return value */
1948   if ( depth > maxdepth )
1949   {
1950     return 0;
1951   }
1952 
1953   nvertices = rtgeom_count_vertices(ctx, geom);
1954   /* Skip empties entirely */
1955   if ( nvertices == 0 )
1956   {
1957     return 0;
1958   }
1959 
1960   /* If it is under the vertex tolerance, just add it, we're done */
1961   if ( nvertices < maxvertices )
1962   {
1963     rtcollection_add_rtgeom(ctx, col, rtgeom_clone_deep(ctx, geom));
1964     return 1;
1965   }
1966 
1967   subbox1 = subbox2 = *clip;
1968   if ( width > height )
1969   {
1970     subbox1.xmax = subbox2.xmin = (clip->xmin + clip->xmax)/2;
1971   }
1972   else
1973   {
1974     subbox1.ymax = subbox2.ymin = (clip->ymin + clip->ymax)/2;
1975   }
1976 
1977   if ( height == 0 )
1978   {
1979     subbox1.ymax += FP_TOLERANCE;
1980     subbox2.ymax += FP_TOLERANCE;
1981     subbox1.ymin -= FP_TOLERANCE;
1982     subbox2.ymin -= FP_TOLERANCE;
1983   }
1984 
1985   if ( width == 0 )
1986   {
1987     subbox1.xmax += FP_TOLERANCE;
1988     subbox2.xmax += FP_TOLERANCE;
1989     subbox1.xmin -= FP_TOLERANCE;
1990     subbox2.xmin -= FP_TOLERANCE;
1991   }
1992 
1993   clipped1 = rtgeom_clip_by_rect(ctx, geom, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
1994   clipped2 = rtgeom_clip_by_rect(ctx, geom, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
1995 
1996   if ( clipped1 )
1997   {
1998     n += rtgeom_subdivide_recursive(ctx, clipped1, maxvertices, ++depth, col, &subbox1);
1999     rtgeom_free(ctx, clipped1);
2000   }
2001 
2002   if ( clipped2 )
2003   {
2004     n += rtgeom_subdivide_recursive(ctx, clipped2, maxvertices, ++depth, col, &subbox2);
2005     rtgeom_free(ctx, clipped2);
2006   }
2007 
2008   return n;
2009 
2010 }
2011 
2012 RTCOLLECTION *
rtgeom_subdivide(const RTCTX * ctx,const RTGEOM * geom,int maxvertices)2013 rtgeom_subdivide(const RTCTX *ctx, const RTGEOM *geom, int maxvertices)
2014 {
2015   static int startdepth = 0;
2016   static int minmaxvertices = 8;
2017   RTCOLLECTION *col;
2018   RTGBOX clip;
2019 
2020   col = rtcollection_construct_empty(ctx, RTCOLLECTIONTYPE, geom->srid, rtgeom_has_z(ctx, geom), rtgeom_has_m(ctx, geom));
2021 
2022   if ( rtgeom_is_empty(ctx, geom) )
2023     return col;
2024 
2025   if ( maxvertices < minmaxvertices )
2026   {
2027     rtcollection_free(ctx, col);
2028     rterror(ctx, "%s: cannot subdivide to fewer than %d vertices per output", __func__, minmaxvertices);
2029   }
2030 
2031   clip = *(rtgeom_get_bbox(ctx, geom));
2032   rtgeom_subdivide_recursive(ctx, geom, maxvertices, startdepth, col, &clip);
2033   rtgeom_set_srid(ctx, (RTGEOM*)col, geom->srid);
2034   return col;
2035 }
2036 
2037 
2038 int
rtgeom_is_trajectory(const RTCTX * ctx,const RTGEOM * geom)2039 rtgeom_is_trajectory(const RTCTX *ctx, const RTGEOM *geom)
2040 {
2041   int type = geom->type;
2042 
2043   if( type != RTLINETYPE )
2044   {
2045     rtnotice(ctx, "Geometry is not a LINESTRING");
2046     return RT_FALSE;
2047   }
2048   return rtline_is_trajectory(ctx, (RTLINE*)geom);
2049 }
2050 
2051