PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwgeom_api.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * Copyright 2001-2006 Refractions Research Inc.
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 
14 #include "liblwgeom_internal.h"
15 #include "lwgeom_log.h"
16 
17 #include <float.h>
18 #include <stdio.h>
19 #include <errno.h>
20 #include <assert.h>
21 
22 /*
23  * Lower this to reduce integrity checks
24  */
25 #define PARANOIA_LEVEL 1
26 
27 
28 /**********************************************************************
29  * BOX routines
30  *
31  * returns the float thats very close to the input, but <=
32  * handles the funny differences in float4 and float8 reps.
33  **********************************************************************/
34 
35 typedef union
36 {
37  float value;
38  uint32_t word;
40 
41 #define GET_FLOAT_WORD(i,d) \
42  do { \
43  ieee_float_shape_type gf_u; \
44  gf_u.value = (d); \
45  (i) = gf_u.word; \
46  } while (0)
47 
48 
49 #define SET_FLOAT_WORD(d,i) \
50  do { \
51  ieee_float_shape_type sf_u; \
52  sf_u.word = (i); \
53  (d) = sf_u.value; \
54  } while (0)
55 
56 
57 /*
58  * Returns the next smaller or next larger float
59  * from x (in direction of y).
60  */
61 static float
62 nextafterf_custom(float x, float y)
63 {
64  int hx,hy,ix,iy;
65 
66  GET_FLOAT_WORD(hx,x);
67  GET_FLOAT_WORD(hy,y);
68  ix = hx&0x7fffffff; /* |x| */
69  iy = hy&0x7fffffff; /* |y| */
70 
71  if ((ix>0x7f800000) || /* x is nan */
72  (iy>0x7f800000)) /* y is nan */
73  return x+y;
74  if (x==y) return y; /* x=y, return y */
75  if (ix==0)
76  {
77  /* x == 0 */
78  SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
79  y = x*x;
80  if (y==x) return y;
81  else return x; /* raise underflow flag */
82  }
83  if (hx>=0)
84  {
85  /* x > 0 */
86  if (hx>hy)
87  {
88  /* x > y, x -= ulp */
89  hx -= 1;
90  }
91  else
92  {
93  /* x < y, x += ulp */
94  hx += 1;
95  }
96  }
97  else
98  {
99  /* x < 0 */
100  if (hy>=0||hx>hy)
101  {
102  /* x < y, x -= ulp */
103  hx -= 1;
104  }
105  else
106  {
107  /* x > y, x += ulp */
108  hx += 1;
109  }
110  }
111  hy = hx&0x7f800000;
112  if (hy>=0x7f800000) return x+x; /* overflow */
113  if (hy<0x00800000)
114  {
115  /* underflow */
116  y = x*x;
117  if (y!=x)
118  {
119  /* raise underflow flag */
120  SET_FLOAT_WORD(y,hx);
121  return y;
122  }
123  }
124  SET_FLOAT_WORD(x,hx);
125  return x;
126 }
127 
128 
129 float next_float_down(double d)
130 {
131  float result = d;
132 
133  if ( ((double) result) <=d)
134  return result;
135 
136  return nextafterf_custom(result, result - 1000000);
137 
138 }
139 
140 /*
141  * Returns the float thats very close to the input, but >=.
142  * handles the funny differences in float4 and float8 reps.
143  */
144 float
145 next_float_up(double d)
146 {
147  float result = d;
148 
149  if ( ((double) result) >=d)
150  return result;
151 
152  return nextafterf_custom(result, result + 1000000);
153 }
154 
155 
156 /*
157  * Returns the double thats very close to the input, but <.
158  * handles the funny differences in float4 and float8 reps.
159  */
160 double
162 {
163  double result = d;
164 
165  if ( result < d)
166  return result;
167 
168  return nextafterf_custom(result, result - 1000000);
169 }
170 
171 /*
172  * Returns the double thats very close to the input, but >
173  * handles the funny differences in float4 and float8 reps.
174  */
175 double
177 {
178  double result = d;
179 
180  if ( result > d)
181  return result;
182 
183  return nextafterf_custom(result, result + 1000000);
184 }
185 
186 
187 /************************************************************************
188  * POINTARRAY support functions
189  *
190  * TODO: should be moved to ptarray.c probably
191  *
192  ************************************************************************/
193 
194 /*
195  * Copies a point from the point array into the parameter point
196  * will set point's z=NO_Z_VALUE if pa is 2d
197  * will set point's m=NO_M_VALUE if pa is 3d or 2d
198  *
199  * NOTE: point is a real POINT3D *not* a pointer
200  */
201 POINT4D
202 getPoint4d(const POINTARRAY *pa, int n)
203 {
204  POINT4D result;
205  getPoint4d_p(pa, n, &result);
206  return result;
207 }
208 
209 /*
210  * Copies a point from the point array into the parameter point
211  * will set point's z=NO_Z_VALUE if pa is 2d
212  * will set point's m=NO_M_VALUE if pa is 3d or 2d
213  *
214  * NOTE: this will modify the point4d pointed to by 'point'.
215  */
216 int
217 getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *op)
218 {
219  uint8_t *ptr;
220  int zmflag;
221 
222 #if PARANOIA_LEVEL > 0
223  if ( ! pa ) lwerror("getPoint4d_p: NULL pointarray");
224 
225  if ( (n<0) || (n>=pa->npoints))
226  {
227  lwerror("getPoint4d_p: point offset out of range");
228  }
229 #endif
230 
231  LWDEBUG(4, "getPoint4d_p called.");
232 
233  /* Get a pointer to nth point offset and zmflag */
234  ptr=getPoint_internal(pa, n);
235  zmflag=FLAGS_GET_ZM(pa->flags);
236 
237  LWDEBUGF(4, "ptr %p, zmflag %d", ptr, zmflag);
238 
239  switch (zmflag)
240  {
241  case 0: /* 2d */
242  memcpy(op, ptr, sizeof(POINT2D));
243  op->m=NO_M_VALUE;
244  op->z=NO_Z_VALUE;
245  break;
246 
247  case 3: /* ZM */
248  memcpy(op, ptr, sizeof(POINT4D));
249  break;
250 
251  case 2: /* Z */
252  memcpy(op, ptr, sizeof(POINT3DZ));
253  op->m=NO_M_VALUE;
254  break;
255 
256  case 1: /* M */
257  memcpy(op, ptr, sizeof(POINT3DM));
258  op->m=op->z; /* we use Z as temporary storage */
259  op->z=NO_Z_VALUE;
260  break;
261 
262  default:
263  lwerror("Unknown ZM flag ??");
264  }
265  return 1;
266 
267 }
268 
269 
270 
271 /*
272  * Copy a point from the point array into the parameter point
273  * will set point's z=NO_Z_VALUE if pa is 2d
274  * NOTE: point is a real POINT3DZ *not* a pointer
275  */
276 POINT3DZ
277 getPoint3dz(const POINTARRAY *pa, int n)
278 {
280  getPoint3dz_p(pa, n, &result);
281  return result;
282 }
283 
284 /*
285  * Copy a point from the point array into the parameter point
286  * will set point's z=NO_Z_VALUE if pa is 2d
287  *
288  * NOTE: point is a real POINT3DZ *not* a pointer
289  */
290 POINT3DM
291 getPoint3dm(const POINTARRAY *pa, int n)
292 {
294  getPoint3dm_p(pa, n, &result);
295  return result;
296 }
297 
298 /*
299  * Copy a point from the point array into the parameter point
300  * will set point's z=NO_Z_VALUE if pa is 2d
301  *
302  * NOTE: this will modify the point3dz pointed to by 'point'.
303  */
304 int
305 getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *op)
306 {
307  uint8_t *ptr;
308 
309 #if PARANOIA_LEVEL > 0
310  if ( ! pa ) return 0;
311 
312  if ( (n<0) || (n>=pa->npoints))
313  {
314  LWDEBUGF(4, "%d out of numpoint range (%d)", n, pa->npoints);
315  return 0; /*error */
316  }
317 #endif
318 
319  LWDEBUGF(2, "getPoint3dz_p called on array of %d-dimensions / %u pts",
320  FLAGS_NDIMS(pa->flags), pa->npoints);
321 
322  /* Get a pointer to nth point offset */
323  ptr=getPoint_internal(pa, n);
324 
325  /*
326  * if input POINTARRAY has the Z, it is always
327  * at third position so make a single copy
328  */
329  if ( FLAGS_GET_Z(pa->flags) )
330  {
331  memcpy(op, ptr, sizeof(POINT3DZ));
332  }
333 
334  /*
335  * Otherwise copy the 2d part and initialize
336  * Z to NO_Z_VALUE
337  */
338  else
339  {
340  memcpy(op, ptr, sizeof(POINT2D));
341  op->z=NO_Z_VALUE;
342  }
343 
344  return 1;
345 
346 }
347 
348 /*
349  * Copy a point from the point array into the parameter point
350  * will set point's m=NO_Z_VALUE if pa has no M
351  *
352  * NOTE: this will modify the point3dm pointed to by 'point'.
353  */
354 int
355 getPoint3dm_p(const POINTARRAY *pa, int n, POINT3DM *op)
356 {
357  uint8_t *ptr;
358  int zmflag;
359 
360 #if PARANOIA_LEVEL > 0
361  if ( ! pa ) return 0;
362 
363  if ( (n<0) || (n>=pa->npoints))
364  {
365  lwerror("%d out of numpoint range (%d)", n, pa->npoints);
366  return 0; /*error */
367  }
368 #endif
369 
370  LWDEBUGF(2, "getPoint3dm_p(%d) called on array of %d-dimensions / %u pts",
371  n, FLAGS_NDIMS(pa->flags), pa->npoints);
372 
373 
374  /* Get a pointer to nth point offset and zmflag */
375  ptr=getPoint_internal(pa, n);
376  zmflag=FLAGS_GET_ZM(pa->flags);
377 
378  /*
379  * if input POINTARRAY has the M and NO Z,
380  * we can issue a single memcpy
381  */
382  if ( zmflag == 1 )
383  {
384  memcpy(op, ptr, sizeof(POINT3DM));
385  return 1;
386  }
387 
388  /*
389  * Otherwise copy the 2d part and
390  * initialize M to NO_M_VALUE
391  */
392  memcpy(op, ptr, sizeof(POINT2D));
393 
394  /*
395  * Then, if input has Z skip it and
396  * copy next double, otherwise initialize
397  * M to NO_M_VALUE
398  */
399  if ( zmflag == 3 )
400  {
401  ptr+=sizeof(POINT3DZ);
402  memcpy(&(op->m), ptr, sizeof(double));
403  }
404  else
405  {
406  op->m=NO_M_VALUE;
407  }
408 
409  return 1;
410 }
411 
412 
413 /*
414  * Copy a point from the point array into the parameter point
415  * z value (if present) is not returned.
416  *
417  * NOTE: point is a real POINT2D *not* a pointer
418  */
419 POINT2D
420 getPoint2d(const POINTARRAY *pa, int n)
421 {
422  const POINT2D *result;
423  result = getPoint2d_cp(pa, n);
424  return *result;
425 }
426 
427 /*
428  * Copy a point from the point array into the parameter point
429  * z value (if present) is not returned.
430  *
431  * NOTE: this will modify the point2d pointed to by 'point'.
432  */
433 int
434 getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
435 {
436 #if PARANOIA_LEVEL > 0
437  if ( ! pa ) return 0;
438 
439  if ( (n<0) || (n>=pa->npoints))
440  {
441  lwerror("getPoint2d_p: point offset out of range");
442  return 0; /*error */
443  }
444 #endif
445 
446  /* this does x,y */
447  memcpy(point, getPoint_internal(pa, n), sizeof(POINT2D));
448  return 1;
449 }
450 
457 const POINT2D*
458 getPoint2d_cp(const POINTARRAY *pa, int n)
459 {
460  if ( ! pa ) return 0;
461 
462  if ( (n<0) || (n>=pa->npoints))
463  {
464  lwerror("getPoint2D_const_p: point offset out of range");
465  return 0; /*error */
466  }
467 
468  return (const POINT2D*)getPoint_internal(pa, n);
469 }
470 
471 const POINT3DZ*
472 getPoint3dz_cp(const POINTARRAY *pa, int n)
473 {
474  if ( ! pa ) return 0;
475 
476  if ( ! FLAGS_GET_Z(pa->flags) )
477  {
478  lwerror("getPoint3dz_cp: no Z coordinates in point array");
479  return 0; /*error */
480  }
481 
482  if ( (n<0) || (n>=pa->npoints))
483  {
484  lwerror("getPoint3dz_cp: point offset out of range");
485  return 0; /*error */
486  }
487 
488  return (const POINT3DZ*)getPoint_internal(pa, n);
489 }
490 
491 
492 
493 /*
494  * set point N to the given value
495  * NOTE that the pointarray can be of any
496  * dimension, the appropriate ordinate values
497  * will be extracted from it
498  *
499  */
500 void
501 ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
502 {
503  uint8_t *ptr;
504  assert(n >= 0 && n < pa->npoints);
505  ptr=getPoint_internal(pa, n);
506  switch ( FLAGS_GET_ZM(pa->flags) )
507  {
508  case 3:
509  memcpy(ptr, p4d, sizeof(POINT4D));
510  break;
511  case 2:
512  memcpy(ptr, p4d, sizeof(POINT3DZ));
513  break;
514  case 1:
515  memcpy(ptr, p4d, sizeof(POINT2D));
516  ptr+=sizeof(POINT2D);
517  memcpy(ptr, &(p4d->m), sizeof(double));
518  break;
519  case 0:
520  memcpy(ptr, p4d, sizeof(POINT2D));
521  break;
522  }
523 }
524 
525 
526 
527 
528 /*****************************************************************************
529  * Basic sub-geometry types
530  *****************************************************************************/
531 
532 /* handle missaligned uint32_t32 data */
533 uint32_t
534 lw_get_uint32_t(const uint8_t *loc)
535 {
536  uint32_t result;
537 
538  memcpy(&result, loc, sizeof(uint32_t));
539  return result;
540 }
541 
542 /* handle missaligned signed int32_t data */
543 int32_t
544 lw_get_int32_t(const uint8_t *loc)
545 {
546  int32_t result;
547 
548  memcpy(&result,loc, sizeof(int32_t));
549  return result;
550 }
551 
552 
553 /************************************************
554  * debugging routines
555  ************************************************/
556 
557 void printBOX3D(BOX3D *box)
558 {
559  lwnotice("BOX3D: %g %g, %g %g", box->xmin, box->ymin,
560  box->xmax, box->ymax);
561 }
562 
564 {
565  int t;
566  POINT4D pt;
567  char *mflag;
568 
569 
570  if ( FLAGS_GET_M(pa->flags) ) mflag = "M";
571  else mflag = "";
572 
573  lwnotice(" POINTARRAY%s{", mflag);
574  lwnotice(" ndims=%i, ptsize=%i",
576  lwnotice(" npoints = %i", pa->npoints);
577 
578  for (t =0; t<pa->npoints; t++)
579  {
580  getPoint4d_p(pa, t, &pt);
581  if (FLAGS_NDIMS(pa->flags) == 2)
582  {
583  lwnotice(" %i : %lf,%lf",t,pt.x,pt.y);
584  }
585  if (FLAGS_NDIMS(pa->flags) == 3)
586  {
587  lwnotice(" %i : %lf,%lf,%lf",t,pt.x,pt.y,pt.z);
588  }
589  if (FLAGS_NDIMS(pa->flags) == 4)
590  {
591  lwnotice(" %i : %lf,%lf,%lf,%lf",t,pt.x,pt.y,pt.z,pt.m);
592  }
593  }
594 
595  lwnotice(" }");
596 }
597 
598 
603 uint8_t
604 parse_hex(char *str)
605 {
606  /* do this a little brute force to make it faster */
607 
608  uint8_t result_high = 0;
609  uint8_t result_low = 0;
610 
611  switch (str[0])
612  {
613  case '0' :
614  result_high = 0;
615  break;
616  case '1' :
617  result_high = 1;
618  break;
619  case '2' :
620  result_high = 2;
621  break;
622  case '3' :
623  result_high = 3;
624  break;
625  case '4' :
626  result_high = 4;
627  break;
628  case '5' :
629  result_high = 5;
630  break;
631  case '6' :
632  result_high = 6;
633  break;
634  case '7' :
635  result_high = 7;
636  break;
637  case '8' :
638  result_high = 8;
639  break;
640  case '9' :
641  result_high = 9;
642  break;
643  case 'A' :
644  case 'a' :
645  result_high = 10;
646  break;
647  case 'B' :
648  case 'b' :
649  result_high = 11;
650  break;
651  case 'C' :
652  case 'c' :
653  result_high = 12;
654  break;
655  case 'D' :
656  case 'd' :
657  result_high = 13;
658  break;
659  case 'E' :
660  case 'e' :
661  result_high = 14;
662  break;
663  case 'F' :
664  case 'f' :
665  result_high = 15;
666  break;
667  }
668  switch (str[1])
669  {
670  case '0' :
671  result_low = 0;
672  break;
673  case '1' :
674  result_low = 1;
675  break;
676  case '2' :
677  result_low = 2;
678  break;
679  case '3' :
680  result_low = 3;
681  break;
682  case '4' :
683  result_low = 4;
684  break;
685  case '5' :
686  result_low = 5;
687  break;
688  case '6' :
689  result_low = 6;
690  break;
691  case '7' :
692  result_low = 7;
693  break;
694  case '8' :
695  result_low = 8;
696  break;
697  case '9' :
698  result_low = 9;
699  break;
700  case 'A' :
701  case 'a' :
702  result_low = 10;
703  break;
704  case 'B' :
705  case 'b' :
706  result_low = 11;
707  break;
708  case 'C' :
709  case 'c' :
710  result_low = 12;
711  break;
712  case 'D' :
713  case 'd' :
714  result_low = 13;
715  break;
716  case 'E' :
717  case 'e' :
718  result_low = 14;
719  break;
720  case 'F' :
721  case 'f' :
722  result_low = 15;
723  break;
724  }
725  return (uint8_t) ((result_high<<4) + result_low);
726 }
727 
728 
738 void
739 deparse_hex(uint8_t str, char *result)
740 {
741  int input_high;
742  int input_low;
743  static char outchr[]=
744  {
745  "0123456789ABCDEF"
746  };
747 
748  input_high = (str>>4);
749  input_low = (str & 0x0F);
750 
751  result[0] = outchr[input_high];
752  result[1] = outchr[input_low];
753 
754 }
755 
756 
770 void
772 {
773 #if PARANOIA_LEVEL > 0
774  double absF=fabs(F);
775  if ( absF < 0 || absF > 1 )
776  {
777  lwerror("interpolate_point4d: invalid F (%g)", F);
778  }
779 #endif
780  I->x=A->x+((B->x-A->x)*F);
781  I->y=A->y+((B->y-A->y)*F);
782  I->z=A->z+((B->z-A->z)*F);
783  I->m=A->m+((B->m-A->m)*F);
784 }
785 
786 
787 
double x
Definition: liblwgeom.h:308
double z
Definition: liblwgeom.h:290
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:501
uint32_t lw_get_uint32_t(const uint8_t *loc)
Definition: lwgeom_api.c:534
double m
Definition: liblwgeom.h:308
int npoints
Definition: liblwgeom.h:327
double next_double_up(float d)
Definition: lwgeom_api.c:176
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
int32_t lw_get_int32_t(const uint8_t *loc)
Definition: lwgeom_api.c:544
float next_float_down(double d)
Definition: lwgeom_api.c:129
#define FLAGS_GET_ZM(flags)
Definition: liblwgeom.h:119
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
double ymin
Definition: liblwgeom.h:233
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a pointer into the POINTARRAY serialized_ptlist, suitable for reading from.
Definition: lwgeom_api.c:458
char ** result
Definition: liblwgeom.h:218
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:54
uint8_t parse_hex(char *str)
Given a string with at least 2 chars in it, convert them to a byte value.
Definition: lwgeom_api.c:604
double xmin
Definition: liblwgeom.h:233
double m
Definition: liblwgeom.h:302
uint8_t flags
Definition: liblwgeom.h:325
void printBOX3D(BOX3D *box)
Definition: lwgeom_api.c:557
uint8_t * getPoint_internal(const POINTARRAY *pa, int n)
Definition: ptarray.c:1645
int ptarray_point_size(const POINTARRAY *pa)
Definition: ptarray.c:41
static float nextafterf_custom(float x, float y)
Definition: lwgeom_api.c:62
#define FLAGS_GET_Z(flags)
Macros for manipulating the 'flags' byte.
Definition: liblwgeom.h:106
double z
Definition: liblwgeom.h:308
tuple x
Definition: pixval.py:53
POINT3DM getPoint3dm(const POINTARRAY *pa, int n)
Definition: lwgeom_api.c:291
#define NO_Z_VALUE
#define SET_FLOAT_WORD(d, i)
Definition: lwgeom_api.c:49
int getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *op)
Definition: lwgeom_api.c:305
double xmax
Definition: liblwgeom.h:234
const POINT3DZ * getPoint3dz_cp(const POINTARRAY *pa, int n)
Returns a POINT3DZ pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:472
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *op)
Definition: lwgeom_api.c:217
POINT4D getPoint4d(const POINTARRAY *pa, int n)
Definition: lwgeom_api.c:202
#define FLAGS_GET_M(flags)
Definition: liblwgeom.h:107
int getPoint3dm_p(const POINTARRAY *pa, int n, POINT3DM *op)
Definition: lwgeom_api.c:355
#define NO_M_VALUE
void interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F)
Find interpolation point I between point A and point B so that the len(AI) == len(AB)*F and I falls o...
Definition: lwgeom_api.c:771
void deparse_hex(uint8_t str, char *result)
Given one byte, populate result with two byte representing the hex number.
Definition: lwgeom_api.c:739
float next_float_up(double d)
Definition: lwgeom_api.c:145
void printPA(POINTARRAY *pa)
Definition: lwgeom_api.c:563
double next_double_down(float d)
Definition: lwgeom_api.c:161
double y
Definition: liblwgeom.h:308
double ymax
Definition: liblwgeom.h:234
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FLAGS_NDIMS(flags)
Definition: liblwgeom.h:118
POINT2D getPoint2d(const POINTARRAY *pa, int n)
Definition: lwgeom_api.c:420
tuple y
Definition: pixval.py:54
POINT3DZ getPoint3dz(const POINTARRAY *pa, int n)
Definition: lwgeom_api.c:277
#define GET_FLOAT_WORD(i, d)
Definition: lwgeom_api.c:41