source: trunk/liblwgeom/lwgeom_api.c @ 5274

Last change on this file since 5274 was 5274, checked in by colivier, 6 years ago

make astyle session

  • Property svn:keywords set to Author Date Id Revision
File size: 43.9 KB
Line 
1
2#include <math.h>
3#include <float.h>
4#include <string.h>
5#include <stdio.h>
6#include <errno.h>
7
8#include "liblwgeom.h"
9#include "wktparse.h"
10
11/*
12 * Lower this to reduce integrity checks
13 */
14#define PARANOIA_LEVEL 1
15
16
17
18/**********************************************************************
19 * BOX routines
20 *
21 * returns the float thats very close to the input, but <=
22 *  handles the funny differences in float4 and float8 reps.
23 **********************************************************************/
24
25
26/*
27 * These are taken from glibc
28 * some machines do *not* have these functions defined, so we give
29 *  an implementation of them here.
30 */
31typedef int int32_tt;
32typedef unsigned int u_int32_tt;
33
34typedef union
35{
36        float value;
37        u_int32_tt word;
38} ieee_float_shape_type;
39
40#define GET_FLOAT_WORD(i,d)                     \
41        do {                                    \
42                ieee_float_shape_type gf_u;     \
43                gf_u.value = (d);               \
44                (i) = gf_u.word;                \
45        } while (0)
46
47
48#define SET_FLOAT_WORD(d,i)                     \
49        do {                                    \
50                ieee_float_shape_type sf_u;     \
51                sf_u.word = (i);                \
52                (d) = sf_u.value;               \
53        } while (0)
54
55
56/*
57 * Returns the next smaller or next larger float
58 * from x (in direction of y).
59 */
60float
61nextafterf_custom(float x, float y)
62{
63        int32_tt hx,hy,ix,iy;
64
65        GET_FLOAT_WORD(hx,x);
66        GET_FLOAT_WORD(hy,y);
67        ix = hx&0x7fffffff;             /* |x| */
68        iy = hy&0x7fffffff;             /* |y| */
69
70        if ((ix>0x7f800000) ||   /* x is nan */
71                (iy>0x7f800000))     /* y is nan */
72                return x+y;
73        if (x==y) return y;              /* x=y, return y */
74        if (ix==0)
75        {
76                /* x == 0 */
77                SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
78                y = x*x;
79                if (y==x) return y;
80                else return x;   /* raise underflow flag */
81        }
82        if (hx>=0)
83        {
84                /* x > 0 */
85                if (hx>hy)
86                {
87                        /* x > y, x -= ulp */
88                        hx -= 1;
89                }
90                else
91                {
92                        /* x < y, x += ulp */
93                        hx += 1;
94                }
95        }
96        else
97        {
98                /* x < 0 */
99                if (hy>=0||hx>hy)
100                {
101                        /* x < y, x -= ulp */
102                        hx -= 1;
103                }
104                else
105                {
106                        /* x > y, x += ulp */
107                        hx += 1;
108                }
109        }
110        hy = hx&0x7f800000;
111        if (hy>=0x7f800000) return x+x;  /* overflow  */
112        if (hy<0x00800000)
113        {
114                /* underflow */
115                y = x*x;
116                if (y!=x)
117                {
118                        /* raise underflow flag */
119                        SET_FLOAT_WORD(y,hx);
120                        return y;
121                }
122        }
123        SET_FLOAT_WORD(x,hx);
124        return x;
125}
126
127
128float nextDown_f(double d)
129{
130        float result  = d;
131
132        if ( ((double) result) <=d)
133                return result;
134
135        return nextafterf_custom(result, result - 1000000);
136
137}
138
139/*
140 * Returns the float thats very close to the input, but >=.
141 * handles the funny differences in float4 and float8 reps.
142 */
143float
144nextUp_f(double d)
145{
146        float result  = d;
147
148        if ( ((double) result) >=d)
149                return result;
150
151        return nextafterf_custom(result, result + 1000000);
152}
153
154
155/*
156 * Returns the double thats very close to the input, but <.
157 * handles the funny differences in float4 and float8 reps.
158 */
159double
160nextDown_d(float d)
161{
162        double result  = d;
163
164        if ( result < d)
165                return result;
166
167        return nextafterf_custom(result, result - 1000000);
168}
169
170/*
171 * Returns the double thats very close to the input, but >
172 * handles the funny differences in float4 and float8 reps.
173 */
174double
175nextUp_d(float d)
176{
177        double result  = d;
178
179        if ( result > d)
180                return result;
181
182        return nextafterf_custom(result, result + 1000000);
183}
184
185
186
187/*
188 * Convert BOX3D to BOX2D
189 * returned box2d is allocated with 'lwalloc'
190 */
191BOX2DFLOAT4 *
192box3d_to_box2df(BOX3D *box)
193{
194        BOX2DFLOAT4 *result = (BOX2DFLOAT4*) lwalloc(sizeof(BOX2DFLOAT4));
195
196#if PARANOIA_LEVEL > 0
197        if (box == NULL)
198        {
199                lwerror("box3d_to_box2df got NUL box");
200                return NULL;
201        }
202#endif
203
204        result->xmin = nextDown_f(box->xmin);
205        result->ymin = nextDown_f(box->ymin);
206
207        result->xmax = nextUp_f(box->xmax);
208        result->ymax = nextUp_f(box->ymax);
209
210        return result;
211}
212
213/*
214 * Convert BOX3D to BOX2D using pre-allocated BOX2D
215 * returned box2d is allocated with 'lwalloc'
216 * return 0 on error (NULL input box)
217 */
218int
219box3d_to_box2df_p(BOX3D *box, BOX2DFLOAT4 *result)
220{
221#if PARANOIA_LEVEL > 0
222        if (box == NULL)
223        {
224                lwerror("box3d_to_box2df got NUL box");
225                return 0;
226        }
227#endif
228
229        result->xmin = nextDown_f(box->xmin);
230        result->ymin = nextDown_f(box->ymin);
231
232        result->xmax = nextUp_f(box->xmax);
233        result->ymax = nextUp_f(box->ymax);
234
235        return 1;
236}
237
238
239
240/*
241 * Convert BOX2D to BOX3D
242 * zmin and zmax are set to NO_Z_VALUE
243 */
244BOX3D
245box2df_to_box3d(BOX2DFLOAT4 *box)
246{
247        BOX3D result;
248
249#if PARANOIA_LEVEL > 0
250        if (box == NULL)
251                lwerror("box2df_to_box3d got NULL box");
252#endif
253
254        result.xmin = box->xmin;
255        result.ymin = box->ymin;
256
257        result.xmax = box->xmax;
258        result.ymax = box->ymax;
259
260        result.zmin = result.zmax = NO_Z_VALUE;
261
262        return result;
263}
264
265/*
266 * Convert BOX2D to BOX3D, using pre-allocated BOX3D as output
267 * Z values are set to NO_Z_VALUE.
268 */
269void
270box2df_to_box3d_p(BOX2DFLOAT4 *box, BOX3D *out)
271{
272        if (box == NULL) return;
273
274        out->xmin = box->xmin;
275        out->ymin = box->ymin;
276
277        out->xmax = box->xmax;
278        out->ymax = box->ymax;
279
280        out->zmin = out->zmax = NO_Z_VALUE;
281}
282
283
284
285/*
286 * Returns a BOX3D that encloses b1 and b2
287 * box3d_union(NULL,A) --> A
288 * box3d_union(A,NULL) --> A
289 * box3d_union(A,B) --> A union B
290 */
291BOX3D *
292box3d_union(BOX3D *b1, BOX3D *b2)
293{
294        BOX3D *result;
295
296        result = lwalloc(sizeof(BOX3D));
297
298        if ( (b1 == NULL) && (b2 == NULL) )
299        {
300                return NULL;
301        }
302
303        if  (b1 == NULL)
304        {
305                /*return b2 */
306                memcpy(result, b2, sizeof(BOX3D));
307                return result;
308        }
309        if (b2 == NULL)
310        {
311                /*return b1 */
312                memcpy(result, b1, sizeof(BOX3D));
313                return result;
314        }
315
316        if (b1->xmin < b2->xmin)
317                result->xmin = b1->xmin;
318        else
319                result->xmin = b2->xmin;
320
321        if (b1->ymin < b2->ymin)
322                result->ymin = b1->ymin;
323        else
324                result->ymin = b2->ymin;
325
326
327        if (b1->xmax > b2->xmax)
328                result->xmax = b1->xmax;
329        else
330                result->xmax = b2->xmax;
331
332        if (b1->ymax > b2->ymax)
333                result->ymax = b1->ymax;
334        else
335                result->ymax = b2->ymax;
336
337        if (b1->zmax > b2->zmax)
338                result->zmax = b1->zmax;
339        else
340                result->zmax = b2->zmax;
341
342        if (b1->zmin > b2->zmin)
343                result->zmin = b1->zmin;
344        else
345                result->zmin = b2->zmin;
346
347        return result;
348}
349
350/* Make given ubox a union of b1 and b2 */
351int
352box3d_union_p(BOX3D *b1, BOX3D *b2, BOX3D *ubox)
353{
354
355        LWDEBUG(2, "box3d_union_p called: (xmin, xmax), (ymin, ymax), (zmin, zmax)");
356        LWDEBUGF(4, "b1: (%.16f, %.16f),(%.16f, %.16f),(%.16f, %.16f)", b1->xmin, b1->xmax, b1->ymin, b1->ymax, b1->zmin, b1->zmax);
357        LWDEBUGF(4, "b2: (%.16f, %.16f),(%.16f, %.16f),(%.16f, %.16f)", b2->xmin, b2->xmax, b2->ymin, b2->ymax, b2->zmin, b2->zmax);
358
359        if ( (b1 == NULL) && (b2 == NULL) )
360        {
361                return 0;
362        }
363
364        if  (b1 == NULL)
365        {
366                memcpy(ubox, b2, sizeof(BOX3D));
367                return 1;
368        }
369        if (b2 == NULL)
370        {
371                memcpy(ubox, b1, sizeof(BOX3D));
372                return 1;
373        }
374
375        if (b1->xmin < b2->xmin)
376                ubox->xmin = b1->xmin;
377        else
378                ubox->xmin = b2->xmin;
379
380        if (b1->ymin < b2->ymin)
381                ubox->ymin = b1->ymin;
382        else
383                ubox->ymin = b2->ymin;
384
385
386        if (b1->xmax > b2->xmax)
387                ubox->xmax = b1->xmax;
388        else
389                ubox->xmax = b2->xmax;
390
391        if (b1->ymax > b2->ymax)
392                ubox->ymax = b1->ymax;
393        else
394                ubox->ymax = b2->ymax;
395
396        if (b1->zmax > b2->zmax)
397                ubox->zmax = b1->zmax;
398        else
399                ubox->zmax = b2->zmax;
400
401        if (b1->zmin < b2->zmin)
402                ubox->zmin = b1->zmin;
403        else
404                ubox->zmin = b2->zmin;
405
406        return 1;
407}
408
409#if 0 /* UNUSED */
410/*
411 * Returns a pointer to internal storage, or NULL
412 * if the serialized form does not have a BBOX.
413 */
414BOX2DFLOAT4 *
415getbox2d_internal(uchar *srl)
416{
417        if (TYPE_HASBBOX(srl[0])) return (BOX2DFLOAT4 *)(srl+1);
418        else return NULL;
419}
420#endif /* UNUSED */
421
422/*
423 * Same as getbox2d, but modifies box instead of returning result on the stack
424 */
425int
426getbox2d_p(uchar *srl, BOX2DFLOAT4 *box)
427{
428        uchar type = srl[0];
429        uchar *loc;
430        BOX3D box3d;
431
432        LWDEBUG(2, "getbox2d_p call");
433
434        loc = srl+1;
435
436        if (lwgeom_hasBBOX(type))
437        {
438                /*woot - this is easy */
439                LWDEBUG(4, "getbox2d_p: has box");
440                memcpy(box, loc, sizeof(BOX2DFLOAT4));
441                return 1;
442        }
443
444        LWDEBUG(4, "getbox2d_p: has no box - computing");
445
446        /* We have to actually compute it! */
447        if ( ! compute_serialized_box3d_p(srl, &box3d) ) return 0;
448
449        LWDEBUGF(4, "getbox2d_p: compute_serialized_box3d returned %p", box3d);
450
451        if ( ! box3d_to_box2df_p(&box3d, box) ) return 0;
452
453        LWDEBUG(4, "getbox2d_p: box3d converted to box2d");
454
455        return 1;
456}
457
458/************************************************************************
459 * POINTARRAY support functions
460 *
461 * TODO: should be moved to ptarray.c probably
462 *
463 ************************************************************************/
464
465/*
466 * Copies a point from the point array into the parameter point
467 * will set point's z=NO_Z_VALUE if pa is 2d
468 * will set point's m=NO_M_VALUE if pa is 3d or 2d
469 *
470 * NOTE: point is a real POINT3D *not* a pointer
471 */
472POINT4D
473getPoint4d(const POINTARRAY *pa, int n)
474{
475        POINT4D result;
476        getPoint4d_p(pa, n, &result);
477        return result;
478}
479
480/*
481 * Copies a point from the point array into the parameter point
482 * will set point's z=NO_Z_VALUE  if pa is 2d
483 * will set point's m=NO_M_VALUE  if pa is 3d or 2d
484 *
485 * NOTE: this will modify the point4d pointed to by 'point'.
486 */
487int
488getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *op)
489{
490        uchar *ptr;
491        int zmflag;
492
493#if PARANOIA_LEVEL > 0
494        if ( ! pa ) lwerror("getPoint4d_p: NULL pointarray");
495
496        if ( (n<0) || (n>=pa->npoints))
497        {
498                lwerror("getPoint4d_p: point offset out of range");
499        }
500#endif
501
502        LWDEBUG(4, "getPoint4d_p called.");
503
504        /* Get a pointer to nth point offset and zmflag */
505        ptr=getPoint_internal(pa, n);
506        zmflag=TYPE_GETZM(pa->dims);
507
508        LWDEBUGF(4, "ptr %p, zmflag %d", ptr, zmflag);
509
510        switch (zmflag)
511        {
512        case 0: /* 2d  */
513                memcpy(op, ptr, sizeof(POINT2D));
514                op->m=NO_M_VALUE;
515                op->z=NO_Z_VALUE;
516                break;
517
518        case 3: /* ZM */
519                memcpy(op, ptr, sizeof(POINT4D));
520                break;
521
522        case 2: /* Z */
523                memcpy(op, ptr, sizeof(POINT3DZ));
524                op->m=NO_M_VALUE;
525                break;
526
527        case 1: /* M */
528                memcpy(op, ptr, sizeof(POINT3DM));
529                op->m=op->z; /* we use Z as temporary storage */
530                op->z=NO_Z_VALUE;
531                break;
532
533        default:
534                lwerror("Unknown ZM flag ??");
535        }
536        return 1;
537
538}
539
540
541
542/*
543 * Copy a point from the point array into the parameter point
544 * will set point's z=NO_Z_VALUE if pa is 2d
545 * NOTE: point is a real POINT3DZ *not* a pointer
546 */
547POINT3DZ
548getPoint3dz(const POINTARRAY *pa, int n)
549{
550        POINT3DZ result;
551        getPoint3dz_p(pa, n, &result);
552        return result;
553}
554
555/*
556 * Copy a point from the point array into the parameter point
557 * will set point's z=NO_Z_VALUE if pa is 2d
558 *
559 * NOTE: point is a real POINT3DZ *not* a pointer
560 */
561POINT3DM
562getPoint3dm(const POINTARRAY *pa, int n)
563{
564        POINT3DM result;
565        getPoint3dm_p(pa, n, &result);
566        return result;
567}
568
569/*
570 * Copy a point from the point array into the parameter point
571 * will set point's z=NO_Z_VALUE if pa is 2d
572 *
573 * NOTE: this will modify the point3dz pointed to by 'point'.
574 */
575int
576getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *op)
577{
578        uchar *ptr;
579
580#if PARANOIA_LEVEL > 0
581        if ( ! pa ) return 0;
582
583        if ( (n<0) || (n>=pa->npoints))
584        {
585                LWDEBUGF(4, "%d out of numpoint range (%d)", n, pa->npoints);
586                return 0; /*error */
587        }
588#endif
589
590        LWDEBUGF(2, "getPoint3dz_p called on array of %d-dimensions / %u pts",
591                 TYPE_NDIMS(pa->dims), pa->npoints);
592
593        /* Get a pointer to nth point offset */
594        ptr=getPoint_internal(pa, n);
595
596        /*
597         * if input POINTARRAY has the Z, it is always
598         * at third position so make a single copy
599         */
600        if ( TYPE_HASZ(pa->dims) )
601        {
602                memcpy(op, ptr, sizeof(POINT3DZ));
603        }
604
605        /*
606         * Otherwise copy the 2d part and initialize
607         * Z to NO_Z_VALUE
608         */
609        else
610        {
611                memcpy(op, ptr, sizeof(POINT2D));
612                op->z=NO_Z_VALUE;
613        }
614
615        return 1;
616
617}
618
619/*
620 * Copy a point from the point array into the parameter point
621 * will set point's m=NO_Z_VALUE if pa has no M
622 *
623 * NOTE: this will modify the point3dm pointed to by 'point'.
624 */
625int
626getPoint3dm_p(const POINTARRAY *pa, int n, POINT3DM *op)
627{
628        uchar *ptr;
629        int zmflag;
630
631#if PARANOIA_LEVEL > 0
632        if ( ! pa ) return 0;
633
634        if ( (n<0) || (n>=pa->npoints))
635        {
636                lwerror("%d out of numpoint range (%d)", n, pa->npoints);
637                return 0; /*error */
638        }
639#endif
640
641        LWDEBUGF(2, "getPoint3dm_p(%d) called on array of %d-dimensions / %u pts",
642                 n, TYPE_NDIMS(pa->dims), pa->npoints);
643
644
645        /* Get a pointer to nth point offset and zmflag */
646        ptr=getPoint_internal(pa, n);
647        zmflag=TYPE_GETZM(pa->dims);
648
649        /*
650         * if input POINTARRAY has the M and NO Z,
651         * we can issue a single memcpy
652         */
653        if ( zmflag == 1 )
654        {
655                memcpy(op, ptr, sizeof(POINT3DM));
656                return 1;
657        }
658
659        /*
660         * Otherwise copy the 2d part and
661         * initialize M to NO_M_VALUE
662         */
663        memcpy(op, ptr, sizeof(POINT2D));
664
665        /*
666         * Then, if input has Z skip it and
667         * copy next double, otherwise initialize
668         * M to NO_M_VALUE
669         */
670        if ( zmflag == 3 )
671        {
672                ptr+=sizeof(POINT3DZ);
673                memcpy(&(op->m), ptr, sizeof(double));
674        }
675        else
676        {
677                op->m=NO_M_VALUE;
678        }
679
680        return 1;
681}
682
683
684/*
685 * Copy a point from the point array into the parameter point
686 * z value (if present) is not returned.
687 *
688 * NOTE: point is a real POINT2D *not* a pointer
689 */
690POINT2D
691getPoint2d(const POINTARRAY *pa, int n)
692{
693        POINT2D result;
694        getPoint2d_p(pa, n, &result);
695        return result;
696}
697
698/*
699 * Copy a point from the point array into the parameter point
700 * z value (if present) is not returned.
701 *
702 * NOTE: this will modify the point2d pointed to by 'point'.
703 */
704int
705getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
706{
707#if PARANOIA_LEVEL > 0
708        if ( ! pa ) return 0;
709
710        if ( (n<0) || (n>=pa->npoints))
711        {
712                lwerror("getPoint2d_p: point offset out of range");
713                return 0; /*error */
714        }
715#endif
716
717        /* this does x,y */
718        memcpy(point, getPoint_internal(pa, n), sizeof(POINT2D));
719        return 1;
720}
721
722/*
723 * set point N to the given value
724 * NOTE that the pointarray can be of any
725 * dimension, the appropriate ordinate values
726 * will be extracted from it
727 *
728 */
729void
730setPoint4d(POINTARRAY *pa, int n, POINT4D *p4d)
731{
732        uchar *ptr=getPoint_internal(pa, n);
733        switch ( TYPE_GETZM(pa->dims) )
734        {
735        case 3:
736                memcpy(ptr, p4d, sizeof(POINT4D));
737                break;
738        case 2:
739                memcpy(ptr, p4d, sizeof(POINT3DZ));
740                break;
741        case 1:
742                memcpy(ptr, p4d, sizeof(POINT2D));
743                ptr+=sizeof(POINT2D);
744                memcpy(ptr, &(p4d->m), sizeof(double));
745                break;
746        case 0:
747                memcpy(ptr, p4d, sizeof(POINT2D));
748                break;
749        }
750}
751
752
753/*
754 * Get a pointer to nth point of a POINTARRAY.
755 * You cannot safely cast this to a real POINT, due to memory alignment
756 * constraints. Use getPoint*_p for that.
757 */
758uchar *
759getPoint_internal(const POINTARRAY *pa, int n)
760{
761        int size;
762
763#if PARANOIA_LEVEL > 0
764        if ( pa == NULL )
765        {
766                lwerror("getPoint got NULL pointarray");
767                return NULL;
768        }
769
770        if ( (n<0) || (n>=pa->npoints))
771        {
772                return NULL; /*error */
773        }
774#endif
775
776        size = pointArray_ptsize(pa);
777
778        return &(pa->serialized_pointlist[size*n]);
779}
780
781
782
783/*
784 * Constructs a POINTARRAY.
785 *
786 * NOTE: points is *not* copied, so be careful about modification
787 * (can be aligned/missaligned).
788 *
789 * NOTE: ndims is descriptive - it describes what type of data 'points'
790 *       points to.  No data conversion is done.
791 */
792POINTARRAY *
793pointArray_construct(uchar *points, char hasz, char hasm, uint32 npoints)
794{
795        POINTARRAY  *pa;
796
797        LWDEBUG(2, "pointArray_construct called.");
798
799        pa = (POINTARRAY*)lwalloc(sizeof(POINTARRAY));
800
801        pa->dims = 0;
802        TYPE_SETZM(pa->dims, hasz?1:0, hasm?1:0);
803        pa->npoints = npoints;
804
805        pa->serialized_pointlist = points;
806
807        LWDEBUGF(4, "pointArray_construct returning %p", pa);
808
809        return pa;
810}
811
812
813/*
814 * Size of point represeneted in the POINTARRAY
815 * 16 for 2d, 24 for 3d, 32 for 4d
816 */
817int
818pointArray_ptsize(const POINTARRAY *pa)
819{
820        LWDEBUGF(2, "pointArray_ptsize: TYPE_NDIMS(pa->dims)=%x",TYPE_NDIMS(pa->dims));
821
822        return sizeof(double)*TYPE_NDIMS(pa->dims);
823}
824
825
826/***************************************************************************
827 * Basic type handling
828 ***************************************************************************/
829
830
831/* Returns true if this type says it has an SRID (S bit set) */
832char
833lwgeom_hasSRID(uchar type)
834{
835        return TYPE_HASSRID(type);
836}
837
838/* Returns either 2,3, or 4 -- 2=2D, 3=3D, 4=4D */
839int
840lwgeom_ndims(uchar type)
841{
842        return TYPE_NDIMS(type);
843}
844
845/* has M ? */
846int lwgeom_hasM(uchar type)
847{
848        return  ( (type & 0x10) >>4);
849}
850
851/* has Z ? */
852int lwgeom_hasZ(uchar type)
853{
854        return  ( (type & 0x20) >>5);
855}
856
857
858/* get base type (ie. POLYGONTYPE) */
859int
860lwgeom_getType(uchar type)
861{
862        LWDEBUGF(2, "lwgeom_getType %d", type);
863
864        return (type & 0x0F);
865}
866
867
868/* Construct a type (hasBOX=false) */
869uchar
870lwgeom_makeType(char hasz, char hasm, char hasSRID, int type)
871{
872        return lwgeom_makeType_full(hasz, hasm, hasSRID, type, 0);
873}
874
875/*
876 * Construct a type
877 * TODO: needs to be expanded to accept explicit MZ type
878 */
879uchar
880lwgeom_makeType_full(char hasz, char hasm, char hasSRID, int type, char hasBBOX)
881{
882        uchar result = (char)type;
883
884        TYPE_SETZM(result, hasz, hasm);
885        TYPE_SETHASSRID(result, hasSRID);
886        TYPE_SETHASBBOX(result, hasBBOX);
887
888        return result;
889}
890
891/* Returns true if there's a bbox in this LWGEOM (B bit set) */
892char
893lwgeom_hasBBOX(uchar type)
894{
895        return TYPE_HASBBOX(type);
896}
897
898/*****************************************************************************
899 * Basic sub-geometry types
900 *****************************************************************************/
901
902/* handle missaligned unsigned int32 data */
903uint32
904lw_get_uint32(const uchar *loc)
905{
906        uint32 result;
907
908        memcpy(&result, loc, sizeof(uint32));
909        return result;
910}
911
912/* handle missaligned signed int32 data */
913int32
914lw_get_int32(const uchar *loc)
915{
916        int32 result;
917
918        memcpy(&result,loc, sizeof(int32));
919        return result;
920}
921
922
923/*************************************************************************
924 *
925 * Multi-geometry support
926 *
927 * Note - for a simple type (ie. point), this will have
928 * sub_geom[0] = serialized_form.
929 *
930 * For multi-geomtries sub_geom[0] will be a few bytes
931 * into the serialized form.
932 *
933 * This function just computes the length of each sub-object and
934 * pre-caches this info.
935 *
936 * For a geometry collection of multi* geometries, you can inspect
937 * the sub-components
938 * as well.
939 */
940LWGEOM_INSPECTED *
941lwgeom_inspect(const uchar *serialized_form)
942{
943        LWGEOM_INSPECTED *result = lwalloc(sizeof(LWGEOM_INSPECTED));
944        uchar typefl = (uchar)serialized_form[0];
945        uchar type;
946        uchar **sub_geoms;
947        const uchar *loc;
948        int     t;
949
950        LWDEBUGF(2, "lwgeom_inspect: serialized@%p", serialized_form);
951
952        if (serialized_form == NULL)
953                return NULL;
954
955        result->serialized_form = serialized_form;
956        result->type = (uchar) serialized_form[0];
957        result->SRID = -1; /* assume */
958
959        type = lwgeom_getType(typefl);
960
961        loc = serialized_form+1;
962
963        if ( lwgeom_hasBBOX(typefl) )
964        {
965                loc += sizeof(BOX2DFLOAT4);
966        }
967
968        if ( lwgeom_hasSRID(typefl) )
969        {
970                result->SRID = lw_get_int32(loc);
971                loc += 4;
972        }
973
974        if ( (type==POINTTYPE) || (type==LINETYPE) || (type==POLYGONTYPE) || (type == CIRCSTRINGTYPE))
975        {
976                /* simple geometry (point/line/polygon/circstring)-- not multi! */
977                result->ngeometries = 1;
978                sub_geoms = (uchar**) lwalloc(sizeof(char*));
979                sub_geoms[0] = (uchar *)serialized_form;
980                result->sub_geoms = (uchar **)sub_geoms;
981                return result;
982        }
983
984        /* its a GeometryCollection or multi* geometry */
985
986        result->ngeometries = lw_get_uint32(loc);
987        loc +=4;
988
989        LWDEBUGF(3, "lwgeom_inspect: geometry is a collection of %d elements",
990                 result->ngeometries);
991
992        if ( ! result->ngeometries ) return result;
993
994        sub_geoms = lwalloc(sizeof(uchar*) * result->ngeometries );
995        result->sub_geoms = sub_geoms;
996        sub_geoms[0] = (uchar *)loc;
997
998        LWDEBUGF(3, "subgeom[0] @ %p (+%d)", sub_geoms[0], sub_geoms[0]-serialized_form);
999
1000        for (t=1; t<result->ngeometries; t++)
1001        {
1002                /* -1 = entire object */
1003                int sub_length = lwgeom_size_subgeom(sub_geoms[t-1], -1);
1004                sub_geoms[t] = sub_geoms[t-1] + sub_length;
1005
1006                LWDEBUGF(3, "subgeom[%d] @ %p (+%d)",
1007                         t, sub_geoms[t], sub_geoms[0]-serialized_form);
1008        }
1009
1010        return result;
1011
1012}
1013
1014
1015/*
1016 * 1st geometry has geom_number = 0
1017 * if the actual sub-geometry isnt a POINT, null is returned (see _gettype()).
1018 * if there arent enough geometries, return null.
1019 * this is fine to call on a point (with geom_num=0),
1020 * multipoint or geometrycollection
1021 */
1022LWPOINT *
1023lwgeom_getpoint(uchar *serialized_form, int geom_number)
1024{
1025        int type = lwgeom_getType((uchar)serialized_form[0]);
1026        uchar *sub_geom;
1027
1028        if ((type == POINTTYPE)  && (geom_number == 0))
1029        {
1030                /* Be nice and do as they want instead of what they say */
1031                return lwpoint_deserialize(serialized_form);
1032        }
1033
1034        if ((type != MULTIPOINTTYPE) && (type != COLLECTIONTYPE) )
1035                return NULL;
1036
1037        sub_geom = lwgeom_getsubgeometry(serialized_form, geom_number);
1038        if (sub_geom == NULL)
1039                return NULL;
1040
1041        type = lwgeom_getType(sub_geom[0]);
1042        if (type != POINTTYPE)
1043                return NULL;
1044
1045        return lwpoint_deserialize(sub_geom);
1046}
1047
1048/*
1049 * 1st geometry has geom_number = 0
1050 * if the actual sub-geometry isnt a POINT, null is returned (see _gettype()).
1051 * if there arent enough geometries, return null.
1052 * this is fine to call on a point (with geom_num=0), multipoint
1053 * or geometrycollection
1054 */
1055LWPOINT *lwgeom_getpoint_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1056{
1057        uchar *sub_geom;
1058        uchar type;
1059
1060        sub_geom = lwgeom_getsubgeometry_inspected(inspected, geom_number);
1061
1062        if (sub_geom == NULL) return NULL;
1063
1064        type = lwgeom_getType(sub_geom[0]);
1065        if (type != POINTTYPE) return NULL;
1066
1067        return lwpoint_deserialize(sub_geom);
1068}
1069
1070
1071/*
1072 * 1st geometry has geom_number = 0
1073 * if the actual geometry isnt a LINE, null is returned (see _gettype()).
1074 * if there arent enough geometries, return null.
1075 * this is fine to call on a line, multiline or geometrycollection
1076 */
1077LWLINE *
1078lwgeom_getline(uchar *serialized_form, int geom_number)
1079{
1080        uchar type = lwgeom_getType( (uchar) serialized_form[0]);
1081        uchar *sub_geom;
1082
1083        if ((type == LINETYPE)  && (geom_number == 0))
1084        {
1085                /* be nice and do as they want instead of what they say */
1086                return lwline_deserialize(serialized_form);
1087        }
1088
1089        if ((type != MULTILINETYPE) && (type != COLLECTIONTYPE) )
1090                return NULL;
1091
1092        sub_geom = lwgeom_getsubgeometry(serialized_form, geom_number);
1093        if (sub_geom == NULL) return NULL;
1094
1095        type = lwgeom_getType((uchar) sub_geom[0]);
1096        if (type != LINETYPE) return NULL;
1097
1098        return lwline_deserialize(sub_geom);
1099}
1100
1101/*
1102 * 1st geometry has geom_number = 0
1103 * if the actual geometry isnt a LINE, null is returned (see _gettype()).
1104 * if there arent enough geometries, return null.
1105 * this is fine to call on a line, multiline or geometrycollection
1106 */
1107LWLINE *
1108lwgeom_getline_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1109{
1110        uchar *sub_geom;
1111        uchar type;
1112
1113        sub_geom = lwgeom_getsubgeometry_inspected(inspected, geom_number);
1114
1115        if (sub_geom == NULL) return NULL;
1116
1117        type = lwgeom_getType((uchar) sub_geom[0]);
1118        if (type != LINETYPE) return NULL;
1119
1120        return lwline_deserialize(sub_geom);
1121}
1122
1123/*
1124 * 1st geometry has geom_number = 0
1125 * if the actual geometry isnt a POLYGON, null is returned (see _gettype()).
1126 * if there arent enough geometries, return null.
1127 * this is fine to call on a polygon, multipolygon or geometrycollection
1128 */
1129LWPOLY *
1130lwgeom_getpoly(uchar *serialized_form, int geom_number)
1131{
1132        uchar type = lwgeom_getType((uchar)serialized_form[0]);
1133        uchar *sub_geom;
1134
1135        if ((type == POLYGONTYPE)  && (geom_number == 0))
1136        {
1137                /* Be nice and do as they want instead of what they say */
1138                return lwpoly_deserialize(serialized_form);
1139        }
1140
1141        if ((type != MULTIPOLYGONTYPE) && (type != COLLECTIONTYPE) )
1142                return NULL;
1143
1144        sub_geom = lwgeom_getsubgeometry(serialized_form, geom_number);
1145        if (sub_geom == NULL) return NULL;
1146
1147        type = lwgeom_getType(sub_geom[0]);
1148        if (type != POLYGONTYPE) return NULL;
1149
1150        return lwpoly_deserialize(sub_geom);
1151}
1152
1153/*
1154 * 1st geometry has geom_number = 0
1155 * if the actual geometry isnt a POLYGON, null is returned (see _gettype()).
1156 * if there arent enough geometries, return null.
1157 * this is fine to call on a polygon, multipolygon or geometrycollection
1158 */
1159LWPOLY *
1160lwgeom_getpoly_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1161{
1162        uchar *sub_geom;
1163        uchar type;
1164
1165        sub_geom = lwgeom_getsubgeometry_inspected(inspected, geom_number);
1166
1167        if (sub_geom == NULL) return NULL;
1168
1169        type = lwgeom_getType(sub_geom[0]);
1170        if (type != POLYGONTYPE) return NULL;
1171
1172        return lwpoly_deserialize(sub_geom);
1173}
1174
1175/*
1176 * 1st geometry has geom_number = 0
1177 * if the actual geometry isnt a CIRCULARSTRING, null is returned (see _gettype()).
1178 * if there arent enough geometries, return null.
1179 * this is fine to call on a circularstring
1180 */
1181LWCIRCSTRING *
1182lwgeom_getcircstring_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1183{
1184        uchar *sub_geom;
1185        uchar type;
1186
1187        sub_geom = lwgeom_getsubgeometry_inspected(inspected, geom_number);
1188
1189        if (sub_geom == NULL) return NULL;
1190
1191        type = lwgeom_getType(sub_geom[0]);
1192        if (type != CIRCSTRINGTYPE) return NULL;
1193
1194        return lwcircstring_deserialize(sub_geom);
1195}
1196
1197/*
1198 * 1st geometry has geom_number = 0
1199 * if there arent enough geometries, return null.
1200 */
1201LWGEOM *lwgeom_getgeom_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1202{
1203        uchar *sub_geom;
1204        uchar type;
1205
1206        sub_geom = lwgeom_getsubgeometry_inspected(inspected, geom_number);
1207
1208        if (sub_geom == NULL) return NULL;
1209
1210        type = lwgeom_getType(sub_geom[0]);
1211
1212        return lwgeom_deserialize(sub_geom);
1213}
1214
1215
1216/*
1217 * This gets the serialized form of a sub-geometry
1218 *
1219 * 1st geometry has geom_number = 0
1220 * if this isnt a multi* geometry, and geom_number ==0 then it returns
1221 * itself.
1222 *
1223 * Returns null on problems.
1224 *
1225 * In the future this is how you would access a muli* portion of a
1226 * geometry collection.
1227 *    GEOMETRYCOLLECTION(MULTIPOINT(0 0, 1 1), LINESTRING(0 0, 1 1))
1228 *   ie. lwgeom_getpoint( lwgeom_getsubgeometry( serialized, 0), 1)
1229 *           --> POINT(1 1)
1230 *
1231 * You can inspect the sub-geometry as well if you wish.
1232 *
1233 */
1234uchar *
1235lwgeom_getsubgeometry(const uchar *serialized_form, int geom_number)
1236{
1237        uchar *result;
1238        /*major cheat!! */
1239        LWGEOM_INSPECTED *inspected = lwgeom_inspect(serialized_form);
1240
1241        result = lwgeom_getsubgeometry_inspected(inspected, geom_number);
1242        lwinspected_release(inspected);
1243        return result;
1244}
1245
1246uchar *
1247lwgeom_getsubgeometry_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1248{
1249        if ((geom_number <0) || (geom_number >= inspected->ngeometries) )
1250        {
1251                lwerror("lwgeom_getsubgeometry_inspected: geom_number out of range");
1252                return NULL;
1253        }
1254
1255        return inspected->sub_geoms[geom_number];
1256}
1257
1258
1259/*
1260 * 1st geometry has geom_number = 0
1261 *  use geom_number = -1 to find the actual type of the serialized form.
1262 *    ie lwgeom_gettype( <'MULTIPOINT(0 0, 1 1)'>, -1)
1263 *                 --> multipoint
1264 *   ie lwgeom_gettype( <'MULTIPOINT(0 0, 1 1)'>, 0)
1265 *                 --> point
1266 * gets the 8bit type of the geometry at location geom_number
1267 */
1268uchar
1269lwgeom_getsubtype(uchar *serialized_form, int geom_number)
1270{
1271        char  result;
1272        /*major cheat!! */
1273        LWGEOM_INSPECTED *inspected = lwgeom_inspect(serialized_form);
1274
1275        result = lwgeom_getsubtype_inspected(inspected, geom_number);
1276        lwinspected_release(inspected);
1277        return result;
1278}
1279
1280uchar
1281lwgeom_getsubtype_inspected(LWGEOM_INSPECTED *inspected, int geom_number)
1282{
1283        if ((geom_number <0) || (geom_number >= inspected->ngeometries) )
1284                return 99;
1285
1286        return inspected->sub_geoms[geom_number][0]; /* 1st byte is type */
1287}
1288
1289
1290/*
1291 * How many sub-geometries are there?
1292 * for point,line,polygon will return 1.
1293 */
1294int
1295lwgeom_getnumgeometries(uchar *serialized_form)
1296{
1297        uchar type = lwgeom_getType((uchar)serialized_form[0]);
1298        uchar *loc;
1299
1300        if ( (type==POINTTYPE) || (type==LINETYPE) || (type==POLYGONTYPE) ||
1301                (type==CIRCSTRINGTYPE) || (type==COMPOUNDTYPE) || (type==CURVEPOLYTYPE) )
1302        {
1303                return 1;
1304        }
1305
1306        loc = serialized_form+1;
1307
1308        if (lwgeom_hasBBOX((uchar) serialized_form[0]))
1309        {
1310                loc += sizeof(BOX2DFLOAT4);
1311        }
1312
1313        if (lwgeom_hasSRID((uchar) serialized_form[0]) )
1314        {
1315                loc += 4;
1316        }
1317        /* its a GeometryCollection or multi* geometry */
1318        return lw_get_uint32(loc);
1319}
1320
1321/*
1322 * How many sub-geometries are there?
1323 * for point,line,polygon will return 1.
1324 */
1325int
1326lwgeom_getnumgeometries_inspected(LWGEOM_INSPECTED *inspected)
1327{
1328        return inspected->ngeometries;
1329}
1330
1331
1332/*
1333 * Set finalType to COLLECTIONTYPE or 0 (0 means choose a best type)
1334 *   (ie. give it 2 points and ask it to be a multipoint)
1335 *  use SRID=-1 for unknown SRID  (will have 8bit type's S = 0)
1336 * all subgeometries must have the same SRID
1337 * if you want to construct an inspected, call this then inspect the result...
1338 */
1339uchar *
1340lwgeom_serialized_construct(int SRID, int finalType, char hasz, char hasm,
1341                            int nsubgeometries, uchar **serialized_subs)
1342{
1343        uint32 *lengths;
1344        int t;
1345        int total_length = 0;
1346        char type = (char)-1;
1347        char this_type = -1;
1348        uchar *result;
1349        uchar *loc;
1350
1351        if (nsubgeometries == 0)
1352                return lwgeom_constructempty(SRID, hasz, hasm);
1353
1354        lengths = lwalloc(sizeof(int32) * nsubgeometries);
1355
1356        for (t=0; t<nsubgeometries; t++)
1357        {
1358                lengths[t] = lwgeom_size_subgeom(serialized_subs[t],-1);
1359                total_length += lengths[t];
1360                this_type = lwgeom_getType((uchar) (serialized_subs[t][0]));
1361                if (type == (char)-1)
1362                {
1363                        type = this_type;
1364                }
1365                else if (type == COLLECTIONTYPE)
1366                {
1367                        /* still a collection type... */
1368                }
1369                else
1370                {
1371                        if ( (this_type == MULTIPOINTTYPE) || (this_type == MULTILINETYPE)  || (this_type == MULTIPOLYGONTYPE) || (this_type == COLLECTIONTYPE) )
1372                        {
1373                                type = COLLECTIONTYPE;
1374                        }
1375                        else
1376                        {
1377                                if ( (this_type == POINTTYPE)  && (type==POINTTYPE) )
1378                                        type=MULTIPOINTTYPE;
1379                                else if ( (this_type == LINETYPE)  && (type==LINETYPE) )
1380                                        type=MULTILINETYPE;
1381                                else if ( (this_type == POLYGONTYPE)  && (type==POLYGONTYPE) )
1382                                        type=MULTIPOLYGONTYPE;
1383                                else if ( (this_type == POLYGONTYPE)  && (type==MULTIPOLYGONTYPE) )
1384                                        ; /* nop */
1385                                else if ( (this_type == LINETYPE)  && (type==MULTILINETYPE) )
1386                                        ; /* nop */
1387                                else if ( (this_type == POINTTYPE)  && (type==MULTIPOINTTYPE) )
1388                                        ; /* nop */
1389                                else
1390                                        type = COLLECTIONTYPE;
1391                        }
1392                }
1393        }
1394
1395        if (type == POINTTYPE)
1396                type = MULTIPOINTTYPE;
1397        if (type == LINETYPE)
1398                type = MULTILINETYPE;
1399        if (type == POINTTYPE)
1400                type = MULTIPOINTTYPE;
1401
1402        if (finalType == COLLECTIONTYPE)
1403                type = COLLECTIONTYPE;
1404
1405        /* now we have a multi* or GEOMETRYCOLLECTION, let's serialize it */
1406
1407        if (SRID != -1)
1408                total_length +=4; /* space for SRID */
1409
1410        total_length +=1 ;   /* main type; */
1411        total_length +=4 ;   /* nsubgeometries */
1412
1413        result = lwalloc(total_length);
1414        result[0] = (uchar) lwgeom_makeType(hasz, hasm, SRID != -1,  type);
1415        if (SRID != -1)
1416        {
1417                memcpy(&result[1],&SRID,4);
1418                loc = result+5;
1419        }
1420        else
1421                loc = result+1;
1422
1423        memcpy(loc,&nsubgeometries,4);
1424        loc +=4;
1425
1426        for (t=0; t<nsubgeometries; t++)
1427        {
1428                memcpy(loc, serialized_subs[t], lengths[t] );
1429                loc += lengths[t] ;
1430        }
1431
1432        lwfree(lengths);
1433        return result;
1434}
1435
1436
1437/*
1438 * Construct the empty geometry (GEOMETRYCOLLECTION(EMPTY)).
1439 *
1440 * Returns serialized form
1441 */
1442uchar *
1443lwgeom_constructempty(int SRID, char hasz, char hasm)
1444{
1445        int size = 0;
1446        uchar *result;
1447        int ngeoms = 0;
1448        uchar *loc;
1449
1450        if (SRID != -1)
1451                size +=4;
1452
1453        size += 5;
1454
1455        result = lwalloc(size);
1456
1457        result[0] = lwgeom_makeType(hasz, hasm, SRID != -1,  COLLECTIONTYPE);
1458        if (SRID != -1)
1459        {
1460                memcpy(&result[1],&SRID,4);
1461                loc = result+5;
1462        }
1463        else
1464                loc = result+1;
1465
1466        memcpy(loc,&ngeoms,4);
1467        return result;
1468}
1469
1470size_t
1471lwgeom_empty_length(int SRID)
1472{
1473        int size = 5;
1474        if ( SRID != 1 ) size += 4;
1475        return size;
1476}
1477
1478/**
1479 * Construct the empty geometry (GEOMETRYCOLLECTION(EMPTY))
1480 * writing it into the provided buffer.
1481 */
1482void
1483lwgeom_constructempty_buf(int SRID, char hasz, char hasm,
1484                          uchar *buf, size_t *retsize)
1485{
1486        int ngeoms = 0;
1487
1488        buf[0] =(uchar) lwgeom_makeType( hasz, hasm, SRID != -1,  COLLECTIONTYPE);
1489        if (SRID != -1)
1490        {
1491                memcpy(&buf[1],&SRID,4);
1492                buf += 5;
1493        }
1494        else
1495                buf += 1;
1496
1497        memcpy(buf, &ngeoms, 4);
1498
1499        if (retsize) *retsize = lwgeom_empty_length(SRID);
1500}
1501
1502/**
1503 * helper function (not for general use)
1504 * find the size a geometry (or a sub-geometry)
1505 * 1st geometry has geom_number = 0
1506 *  use geom_number = -1 to find the actual type of the serialized form.
1507 *    ie lwgeom_gettype( <'MULTIPOINT(0 0, 1 1)'>, -1)
1508 *                 --> size of the multipoint
1509 *   ie lwgeom_gettype( <'MULTIPOINT(0 0, 1 1)'>, 0)
1510 *         --> size of the point
1511 * take a geometry, and find its length
1512 */
1513size_t
1514lwgeom_size(const uchar *serialized_form)
1515{
1516        uchar type = lwgeom_getType((uchar) serialized_form[0]);
1517        int t;
1518        const uchar *loc;
1519        uint32 ngeoms;
1520        int sub_size;
1521        int result = 1; /* type */
1522
1523        LWDEBUG(2, "lwgeom_size called");
1524
1525        if (type == POINTTYPE)
1526        {
1527                LWDEBUG(3, "lwgeom_size: is a point");
1528
1529                return lwgeom_size_point(serialized_form);
1530        }
1531        else if (type == LINETYPE)
1532        {
1533                LWDEBUG(3, "lwgeom_size: is a line");
1534
1535                return lwgeom_size_line(serialized_form);
1536        }
1537        else if (type == CIRCSTRINGTYPE)
1538        {
1539                LWDEBUG(3, "lwgeom_size: is a circularstring");
1540
1541                return lwgeom_size_circstring(serialized_form);
1542        }
1543        else if (type == POLYGONTYPE)
1544        {
1545                LWDEBUG(3, "lwgeom_size: is a polygon");
1546
1547                return lwgeom_size_poly(serialized_form);
1548        }
1549        else if (type == COMPOUNDTYPE)
1550        {
1551                LWDEBUG(3, "lwgeom_size: is a compound curve");
1552        }
1553
1554        if ( type == 0 )
1555        {
1556                lwerror("lwgeom_size called with unknown-typed serialized geometry");
1557                return 0;
1558        }
1559
1560        /*
1561         * Handle all the multi* and geometrycollections the same
1562         *
1563         * NOTE: for a geometry collection of GC of GC of GC we will
1564         *       be recursing...
1565         */
1566
1567        LWDEBUGF(3, "lwgeom_size called on a geoemtry with type %d", type);
1568
1569        loc = serialized_form+1;
1570
1571        if (lwgeom_hasBBOX((uchar) serialized_form[0]))
1572        {
1573                LWDEBUG(3, "lwgeom_size: has bbox");
1574
1575                loc += sizeof(BOX2DFLOAT4);
1576                result +=sizeof(BOX2DFLOAT4);
1577        }
1578
1579        if (lwgeom_hasSRID( (uchar) serialized_form[0]) )
1580        {
1581                LWDEBUG(3, "lwgeom_size: has srid");
1582
1583                result +=4;
1584                loc +=4;
1585        }
1586
1587
1588        ngeoms = lw_get_uint32(loc);
1589        loc +=4;
1590        result += 4; /* numgeoms */
1591
1592        LWDEBUGF(3, "lwgeom_size called on a geoemtry with %d elems (result so far: %d)", ngeoms, result);
1593
1594        for (t=0; t<ngeoms; t++)
1595        {
1596                sub_size = lwgeom_size(loc);
1597
1598                LWDEBUGF(3, " subsize %d", sub_size);
1599
1600                loc += sub_size;
1601                result += sub_size;
1602        }
1603
1604        LWDEBUGF(3, "lwgeom_size returning %d", result);
1605
1606        return result;
1607}
1608
1609size_t
1610lwgeom_size_subgeom(const uchar *serialized_form, int geom_number)
1611{
1612        if (geom_number == -1)
1613        {
1614                return lwgeom_size(serialized_form);
1615        }
1616        return lwgeom_size( lwgeom_getsubgeometry(serialized_form,geom_number));
1617}
1618
1619
1620
1621int
1622lwgeom_size_inspected(const LWGEOM_INSPECTED *inspected, int geom_number)
1623{
1624        return lwgeom_size(inspected->sub_geoms[geom_number]);
1625}
1626
1627int
1628compute_serialized_box3d_p(uchar *srl, BOX3D *out)
1629{
1630        BOX3D *box = compute_serialized_box3d(srl);
1631        if ( ! box ) return 0;
1632        out->xmin = box->xmin;
1633        out->ymin = box->ymin;
1634        out->zmin = box->zmin;
1635        out->xmax = box->xmax;
1636        out->ymax = box->ymax;
1637        out->zmax = box->zmax;
1638        lwfree(box);
1639
1640        return 1;
1641}
1642
1643/**
1644 * Compute bounding box of a serialized LWGEOM, even if it is
1645 * already cached. The computed BOX2DFLOAT4 is stored at
1646 * the given location, the function returns 0 is the geometry
1647 * does not have a bounding box (EMPTY GEOM)
1648 */
1649int
1650compute_serialized_box2d_p(uchar *srl, BOX2DFLOAT4 *out)
1651{
1652        BOX3D *result = compute_serialized_box3d(srl);
1653        if ( ! result ) return 0;
1654        out->xmin = result->xmin;
1655        out->xmax = result->xmax;
1656        out->ymin = result->ymin;
1657        out->ymax = result->ymax;
1658        lwfree(result);
1659
1660        return 1;
1661}
1662
1663/**
1664 * Don't forget to lwfree() result !
1665 */
1666BOX3D *
1667compute_serialized_box3d(uchar *srl)
1668{
1669        int type = lwgeom_getType(srl[0]);
1670        int t;
1671        uchar *loc = srl;
1672        uint32 nelems;
1673        BOX3D *result;
1674        BOX3D b1;
1675        int sub_size;
1676        char nboxes=0;
1677
1678        LWDEBUGF(2, "compute_serialized_box3d called on type %d", type);
1679
1680        loc += 1; /* Move past the 'type' byte. */
1681
1682        if (lwgeom_hasBBOX(srl[0]))
1683        {
1684                loc += sizeof(BOX2DFLOAT4); /* Move past the bbox */
1685        }
1686
1687        if (lwgeom_hasSRID(srl[0]) )
1688        {
1689                loc +=4; /* Move past the SRID */
1690        }
1691
1692        if (type == POINTTYPE)
1693        {
1694                LWPOINT *pt = lwpoint_deserialize(srl);
1695
1696                LWDEBUG(3, "compute_serialized_box3d: lwpoint deserialized");
1697
1698                result = lwpoint_compute_box3d(pt);
1699
1700                LWDEBUG(3, "compute_serialized_box3d: bbox found");
1701
1702                lwpoint_free(pt);
1703                return result;
1704        }
1705
1706        /*
1707        ** For items that have elements (everything except points),
1708        ** nelems == 0 => EMPTY geometry
1709        */
1710        nelems = lw_get_uint32(loc);
1711        if ( nelems == 0 ) return NULL;
1712
1713        if (type == LINETYPE)
1714        {
1715                LWLINE *line = lwline_deserialize(srl);
1716                result = lwline_compute_box3d(line);
1717                lwline_free(line);
1718                return result;
1719
1720        }
1721        else if (type == CIRCSTRINGTYPE)
1722        {
1723                LWCIRCSTRING *curve = lwcircstring_deserialize(srl);
1724                result = lwcircstring_compute_box3d(curve);
1725                lwcircstring_free(curve);
1726                return result;
1727        }
1728        else if (type == POLYGONTYPE)
1729        {
1730                LWPOLY *poly = lwpoly_deserialize(srl);
1731                result = lwpoly_compute_box3d(poly);
1732                lwpoly_free(poly);
1733                return result;
1734        }
1735
1736        if ( ! ( type == MULTIPOINTTYPE || type == MULTILINETYPE ||
1737                 type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ||
1738                 type == COMPOUNDTYPE || type == CURVEPOLYTYPE ||
1739                 type == MULTICURVETYPE || type == MULTISURFACETYPE) )
1740        {
1741                lwnotice("compute_serialized_box3d called on unknown type %d", type);
1742                return NULL;
1743        }
1744
1745        loc += 4;
1746
1747        /* each sub-type */
1748        result = NULL;
1749        for (t=0; t<nelems; t++)
1750        {
1751                if ( compute_serialized_box3d_p(loc, &b1) )
1752                {
1753                        LWDEBUG(3, "Geom %d have box:");
1754#if POSTGIS_DEBUG_LEVEL >= 3
1755                        printBOX3D(&b1);
1756#endif
1757
1758                        if (result)
1759                        {
1760                                nboxes += box3d_union_p(result, &b1, result);
1761                        }
1762                        else
1763                        {
1764                                result = lwalloc(sizeof(BOX3D));
1765                                memcpy(result, &b1, sizeof(BOX3D));
1766                        }
1767                }
1768
1769                sub_size = lwgeom_size(loc);
1770                loc += sub_size;
1771        }
1772
1773        return result;
1774
1775}
1776
1777/****************************************************************
1778 * memory management
1779 *
1780 *  these only delete the memory associated
1781 *  directly with the structure - NOT the stuff pointing into
1782 *  the original de-serialized info
1783 *
1784 ****************************************************************/
1785
1786void
1787lwinspected_release(LWGEOM_INSPECTED *inspected)
1788{
1789        if ( inspected->ngeometries )
1790                lwfree(inspected->sub_geoms);
1791        lwfree(inspected);
1792}
1793
1794
1795
1796
1797
1798
1799/************************************************
1800 * debugging routines
1801 ************************************************/
1802
1803
1804void printBOX3D(BOX3D *box)
1805{
1806        lwnotice("BOX3D: %g %g, %g %g", box->xmin, box->ymin,
1807                 box->xmax, box->ymax);
1808}
1809
1810void printPA(POINTARRAY *pa)
1811{
1812        int t;
1813        POINT4D pt;
1814        char *mflag;
1815
1816
1817        if ( TYPE_HASM(pa->dims) ) mflag = "M";
1818        else mflag = "";
1819
1820        lwnotice("      POINTARRAY%s{", mflag);
1821        lwnotice("                 ndims=%i,   ptsize=%i",
1822                 TYPE_NDIMS(pa->dims), pointArray_ptsize(pa));
1823        lwnotice("                 npoints = %i", pa->npoints);
1824
1825        for (t =0; t<pa->npoints; t++)
1826        {
1827                getPoint4d_p(pa, t, &pt);
1828                if (TYPE_NDIMS(pa->dims) == 2)
1829                {
1830                        lwnotice("                    %i : %lf,%lf",t,pt.x,pt.y);
1831                }
1832                if (TYPE_NDIMS(pa->dims) == 3)
1833                {
1834                        lwnotice("                    %i : %lf,%lf,%lf",t,pt.x,pt.y,pt.z);
1835                }
1836                if (TYPE_NDIMS(pa->dims) == 4)
1837                {
1838                        lwnotice("                    %i : %lf,%lf,%lf,%lf",t,pt.x,pt.y,pt.z,pt.m);
1839                }
1840        }
1841
1842        lwnotice("      }");
1843}
1844
1845void printBYTES(uchar *a, int n)
1846{
1847        int t;
1848        char buff[3];
1849
1850        buff[2] = 0; /* null terminate */
1851
1852        lwnotice(" BYTE ARRAY (n=%i) IN HEX: {", n);
1853        for (t=0; t<n; t++)
1854        {
1855                deparse_hex(a[t], buff);
1856                lwnotice("    %i : %s", t,buff );
1857        }
1858        lwnotice("  }");
1859}
1860
1861
1862void
1863printMULTI(uchar *serialized)
1864{
1865        LWGEOM_INSPECTED *inspected = lwgeom_inspect(serialized);
1866        LWLINE  *line;
1867        LWPOINT *point;
1868        LWPOLY  *poly;
1869        int t;
1870
1871        lwnotice("MULTI* geometry (type = %i), with %i sub-geoms",lwgeom_getType((uchar)serialized[0]), inspected->ngeometries);
1872
1873        for (t=0; t<inspected->ngeometries; t++)
1874        {
1875                lwnotice("      sub-geometry %i:", t);
1876                line = NULL;
1877                point = NULL;
1878                poly = NULL;
1879
1880                line = lwgeom_getline_inspected(inspected,t);
1881                if (line !=NULL)
1882                {
1883                        printLWLINE(line);
1884                }
1885                poly = lwgeom_getpoly_inspected(inspected,t);
1886                if (poly !=NULL)
1887                {
1888                        printLWPOLY(poly);
1889                }
1890                point = lwgeom_getpoint_inspected(inspected,t);
1891                if (point !=NULL)
1892                {
1893                        printPA(point->point);
1894                }
1895        }
1896
1897        lwnotice("end multi*");
1898
1899        lwinspected_release(inspected);
1900}
1901
1902void
1903printType(uchar type)
1904{
1905        lwnotice("type 0x%x ==> hasBBOX=%i, hasSRID=%i, ndims=%i, type=%i",(unsigned int) type, lwgeom_hasBBOX(type), lwgeom_hasSRID(type),lwgeom_ndims(type), lwgeom_getType(type));
1906}
1907
1908/**
1909 * Get the SRID from the LWGEOM.
1910 * None present => -1
1911 */
1912int
1913lwgeom_getsrid(uchar *serialized)
1914{
1915        uchar type = serialized[0];
1916        uchar *loc = serialized+1;
1917
1918        if ( ! lwgeom_hasSRID(type)) return -1;
1919
1920        if (lwgeom_hasBBOX(type))
1921        {
1922                loc += sizeof(BOX2DFLOAT4);
1923        }
1924
1925        return lw_get_int32(loc);
1926}
1927
1928char
1929ptarray_isccw(const POINTARRAY *pa)
1930{
1931        int i;
1932        double area = 0;
1933        POINT2D p1, p2;
1934
1935        for (i=0; i<pa->npoints-1; i++)
1936        {
1937                getPoint2d_p(pa, i, &p1);
1938                getPoint2d_p(pa, i+1, &p2);
1939                area += (p1.y * p2.x) - (p1.x * p2.y);
1940        }
1941        if ( area > 0 ) return 0;
1942        else return 1;
1943}
1944
1945/**
1946 * Returns a BOX2DFLOAT4 that encloses b1 and b2
1947 *
1948 * box2d_union(NULL,A) --> A
1949 * box2d_union(A,NULL) --> A
1950 * box2d_union(A,B) --> A union B
1951 */
1952BOX2DFLOAT4 *
1953box2d_union(BOX2DFLOAT4 *b1, BOX2DFLOAT4 *b2)
1954{
1955        BOX2DFLOAT4 *result;
1956
1957        if ( (b1 == NULL) && (b2 == NULL) )
1958        {
1959                return NULL;
1960        }
1961
1962        result = lwalloc(sizeof(BOX2DFLOAT4));
1963
1964        if  (b1 == NULL)
1965        {
1966                memcpy(result, b2, sizeof(BOX2DFLOAT4));
1967                return result;
1968        }
1969        if (b2 == NULL)
1970        {
1971                memcpy(result, b1, sizeof(BOX2DFLOAT4));
1972                return result;
1973        }
1974
1975        if (b1->xmin < b2->xmin) result->xmin = b1->xmin;
1976        else result->xmin = b2->xmin;
1977
1978        if (b1->ymin < b2->ymin) result->ymin = b1->ymin;
1979        else result->ymin = b2->ymin;
1980
1981        if (b1->xmax > b2->xmax) result->xmax = b1->xmax;
1982        else result->xmax = b2->xmax;
1983
1984        if (b1->ymax > b2->ymax) result->ymax = b1->ymax;
1985        else result->ymax = b2->ymax;
1986
1987        return result;
1988}
1989
1990/**
1991 * ubox may be one of the two args...
1992 * return 1 if done something to ubox, 0 otherwise.
1993 */
1994int
1995box2d_union_p(BOX2DFLOAT4 *b1, BOX2DFLOAT4 *b2, BOX2DFLOAT4 *ubox)
1996{
1997        if ( (b1 == NULL) && (b2 == NULL) )
1998        {
1999                return 0;
2000        }
2001
2002        if  (b1 == NULL)
2003        {
2004                memcpy(ubox, b2, sizeof(BOX2DFLOAT4));
2005                return 1;
2006        }
2007        if (b2 == NULL)
2008        {
2009                memcpy(ubox, b1, sizeof(BOX2DFLOAT4));
2010                return 1;
2011        }
2012
2013        if (b1->xmin < b2->xmin) ubox->xmin = b1->xmin;
2014        else ubox->xmin = b2->xmin;
2015
2016        if (b1->ymin < b2->ymin) ubox->ymin = b1->ymin;
2017        else ubox->ymin = b2->ymin;
2018
2019        if (b1->xmax > b2->xmax) ubox->xmax = b1->xmax;
2020        else ubox->xmax = b2->xmax;
2021
2022        if (b1->ymax > b2->ymax) ubox->ymax = b1->ymax;
2023        else ubox->ymax = b2->ymax;
2024
2025        return 1;
2026}
2027
2028const char *
2029lwgeom_typeflags(uchar type)
2030{
2031        static char flags[5];
2032        int flagno=0;
2033        if ( TYPE_HASZ(type) ) flags[flagno++] = 'Z';
2034        if ( TYPE_HASM(type) ) flags[flagno++] = 'M';
2035        if ( TYPE_HASBBOX(type) ) flags[flagno++] = 'B';
2036        if ( TYPE_HASSRID(type) ) flags[flagno++] = 'S';
2037        flags[flagno] = '\0';
2038
2039        LWDEBUGF(4, "Flags: %s - returning %p", flags, flags);
2040
2041        return flags;
2042}
2043
2044/**
2045 * Given a string with at least 2 chars in it, convert them to
2046 * a byte value.  No error checking done!
2047 */
2048uchar
2049parse_hex(char *str)
2050{
2051        /* do this a little brute force to make it faster */
2052
2053        uchar           result_high = 0;
2054        uchar           result_low = 0;
2055
2056        switch (str[0])
2057        {
2058        case '0' :
2059                result_high = 0;
2060                break;
2061        case '1' :
2062                result_high = 1;
2063                break;
2064        case '2' :
2065                result_high = 2;
2066                break;
2067        case '3' :
2068                result_high = 3;
2069                break;
2070        case '4' :
2071                result_high = 4;
2072                break;
2073        case '5' :
2074                result_high = 5;
2075                break;
2076        case '6' :
2077                result_high = 6;
2078                break;
2079        case '7' :
2080                result_high = 7;
2081                break;
2082        case '8' :
2083                result_high = 8;
2084                break;
2085        case '9' :
2086                result_high = 9;
2087                break;
2088        case 'A' :
2089        case 'a' :
2090                result_high = 10;
2091                break;
2092        case 'B' :
2093        case 'b' :
2094                result_high = 11;
2095                break;
2096        case 'C' :
2097        case 'c' :
2098                result_high = 12;
2099                break;
2100        case 'D' :
2101        case 'd' :
2102                result_high = 13;
2103                break;
2104        case 'E' :
2105        case 'e' :
2106                result_high = 14;
2107                break;
2108        case 'F' :
2109        case 'f' :
2110                result_high = 15;
2111                break;
2112        }
2113        switch (str[1])
2114        {
2115        case '0' :
2116                result_low = 0;
2117                break;
2118        case '1' :
2119                result_low = 1;
2120                break;
2121        case '2' :
2122                result_low = 2;
2123                break;
2124        case '3' :
2125                result_low = 3;
2126                break;
2127        case '4' :
2128                result_low = 4;
2129                break;
2130        case '5' :
2131                result_low = 5;
2132                break;
2133        case '6' :
2134                result_low = 6;
2135                break;
2136        case '7' :
2137                result_low = 7;
2138                break;
2139        case '8' :
2140                result_low = 8;
2141                break;
2142        case '9' :
2143                result_low = 9;
2144                break;
2145        case 'A' :
2146        case 'a' :
2147                result_low = 10;
2148                break;
2149        case 'B' :
2150        case 'b' :
2151                result_low = 11;
2152                break;
2153        case 'C' :
2154        case 'c' :
2155                result_low = 12;
2156                break;
2157        case 'D' :
2158        case 'd' :
2159                result_low = 13;
2160                break;
2161        case 'E' :
2162        case 'e' :
2163                result_low = 14;
2164                break;
2165        case 'F' :
2166        case 'f' :
2167                result_low = 15;
2168                break;
2169        }
2170        return (uchar) ((result_high<<4) + result_low);
2171}
2172
2173
2174/**
2175 * Given one byte, populate result with two byte representing
2176 * the hex number.
2177 *
2178 * Ie. deparse_hex( 255, mystr)
2179 *              -> mystr[0] = 'F' and mystr[1] = 'F'
2180 *
2181 * No error checking done
2182 */
2183void
2184deparse_hex(uchar str, char *result)
2185{
2186        int     input_high;
2187        int  input_low;
2188        static char outchr[]=
2189        {
2190                "0123456789ABCDEF"
2191        };
2192
2193        input_high = (str>>4);
2194        input_low = (str & 0x0F);
2195
2196        result[0] = outchr[input_high];
2197        result[1] = outchr[input_low];
2198
2199}
2200
2201
2202/**
2203 * Find interpolation point I
2204 * between point A and point B
2205 * so that the len(AI) == len(AB)*F
2206 * and I falls on AB segment.
2207 *
2208 * Example:
2209 *
2210 *   F=0.5  :    A----I----B
2211 *   F=1    :    A---------B==I
2212 *   F=0    : A==I---------B
2213 *   F=.2   :    A-I-------B
2214 */
2215void
2216interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F)
2217{
2218#if PARANOIA_LEVEL > 0
2219        double absF=fabs(F);
2220        if ( absF < 0 || absF > 1 )
2221        {
2222                lwerror("interpolate_point4d: invalid F (%g)", F);
2223        }
2224#endif
2225        I->x=A->x+((B->x-A->x)*F);
2226        I->y=A->y+((B->y-A->y)*F);
2227        I->z=A->z+((B->z-A->z)*F);
2228        I->m=A->m+((B->m-A->m)*F);
2229}
Note: See TracBrowser for help on using the repository browser.