PostGIS  2.5.1dev-r@@SVN_REVISION@@
measures3d.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright 2011 Nicklas Avén
22  *
23  **********************************************************************/
24 
25 
26 /**********************************************************************
27  *
28  * PostGIS - Spatial Types for PostgreSQL
29  * http://postgis.net
30  * Copyright 2011 Nicklas Avén
31  *
32  * This is free software; you can redistribute and/or modify it under
33  * the terms of the GNU General Public Licence. See the COPYING file.
34  *
35  **********************************************************************/
36 
37 #include <string.h>
38 #include <stdlib.h>
39 
40 #include "measures3d.h"
41 #include "lwgeom_log.h"
42 
43 
44 static inline int
46 {
47  v->x=p2->x-p1->x;
48  v->y=p2->y-p1->y;
49  v->z=p2->z-p1->z;
50 
51  return LW_TRUE;
52 }
53 
54 static inline int
56 {
57  v->x=(v1->y*v2->z)-(v1->z*v2->y);
58  v->y=(v1->z*v2->x)-(v1->x*v2->z);
59  v->z=(v1->x*v2->y)-(v1->y*v2->x);
60 
61  return LW_TRUE;
62 }
63 
64 
70 static
71 LWGEOM* create_v_line(const LWGEOM *lwgeom,double x, double y, int srid)
72 {
73 
74  LWPOINT *lwpoints[2];
75  GBOX gbox;
76  int rv = lwgeom_calculate_gbox(lwgeom, &gbox);
77 
78  if ( rv == LW_FAILURE )
79  return NULL;
80 
81  lwpoints[0] = lwpoint_make3dz(srid, x, y, gbox.zmin);
82  lwpoints[1] = lwpoint_make3dz(srid, x, y, gbox.zmax);
83 
84  return (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
85 }
86 
87 LWGEOM *
88 lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
89 {
90  return lw_dist3d_distanceline(lw1, lw2, lw1->srid, DIST_MIN);
91 }
92 
93 LWGEOM *
95 {
96  return lw_dist3d_distanceline(lw1, lw2, lw1->srid, DIST_MAX);
97 }
98 
99 LWGEOM *
100 lwgeom_closest_point_3d(const LWGEOM *lw1, const LWGEOM *lw2)
101 {
102  return lw_dist3d_distancepoint(lw1, lw2, lw1->srid, DIST_MIN);
103 }
104 
105 
109 LWGEOM *
110 lw_dist3d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
111 {
112  LWDEBUG(2, "lw_dist3d_distanceline is called");
113  double x1,x2,y1,y2, z1, z2, x, y;
114  double initdistance = ( mode == DIST_MIN ? FLT_MAX : -1.0);
115  DISTPTS3D thedl;
116  LWPOINT *lwpoints[2];
117  LWGEOM *result;
118 
119  thedl.mode = mode;
120  thedl.distance = initdistance;
121  thedl.tolerance = 0.0;
122 
123  /*Check if we really have 3D geometries*/
124  /*If not, send it to 2D-calculations which will give the same result*/
125  /*as an infinite z-value at one or two of the geometries*/
126  if(!lwgeom_has_z(lw1) || !lwgeom_has_z(lw2))
127  {
128 
129  lwnotice("One or both of the geometries is missing z-value. The unknown z-value will be regarded as \"any value\"");
130 
131  if(!lwgeom_has_z(lw1) && !lwgeom_has_z(lw2))
132  return lw_dist2d_distanceline(lw1, lw2, srid, mode);
133 
134  DISTPTS thedl2d;
135  thedl2d.mode = mode;
136  thedl2d.distance = initdistance;
137  thedl2d.tolerance = 0.0;
138  if (!lw_dist2d_comp( lw1,lw2,&thedl2d))
139  {
140  /*should never get here. all cases ought to be error handled earlier*/
141  lwerror("Some unspecified error.");
142  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
143  }
144  LWGEOM *vertical_line;
145  if(!lwgeom_has_z(lw1))
146  {
147  x=thedl2d.p1.x;
148  y=thedl2d.p1.y;
149 
150  vertical_line = create_v_line(lw2,x,y,srid);
151  if (!lw_dist3d_recursive(vertical_line, lw2, &thedl))
152  {
153  /*should never get here. all cases ought to be error handled earlier*/
154  lwfree(vertical_line);
155  lwerror("Some unspecified error.");
156  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
157  }
158  lwfree(vertical_line);
159  }
160  if(!lwgeom_has_z(lw2))
161  {
162  x=thedl2d.p2.x;
163  y=thedl2d.p2.y;
164 
165  vertical_line = create_v_line(lw1,x,y,srid);
166  if (!lw_dist3d_recursive(lw1, vertical_line, &thedl))
167  {
168  /*should never get here. all cases ought to be error handled earlier*/
169  lwfree(vertical_line);
170  lwerror("Some unspecified error.");
171  return (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
172  }
173  lwfree(vertical_line);
174  }
175 
176  }
177  else
178  {
179  if (!lw_dist3d_recursive(lw1, lw2, &thedl))
180  {
181  /*should never get here. all cases ought to be error handled earlier*/
182  lwerror("Some unspecified error.");
183  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
184  }
185  }
186  /*if thedl.distance is unchanged there where only empty geometries input*/
187  if (thedl.distance == initdistance)
188  {
189  LWDEBUG(3, "didn't find geometries to measure between, returning null");
190  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
191  }
192  else
193  {
194  x1=thedl.p1.x;
195  y1=thedl.p1.y;
196  z1=thedl.p1.z;
197  x2=thedl.p2.x;
198  y2=thedl.p2.y;
199  z2=thedl.p2.z;
200 
201  lwpoints[0] = lwpoint_make3dz(srid, x1, y1, z1);
202  lwpoints[1] = lwpoint_make3dz(srid, x2, y2, z2);
203 
204  result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
205  }
206 
207  return result;
208 }
209 
213 LWGEOM *
214 lw_dist3d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
215 {
216 
217  double x,y,z;
218  DISTPTS3D thedl;
219  double initdistance = FLT_MAX;
220  LWGEOM *result;
221 
222  thedl.mode = mode;
223  thedl.distance= initdistance;
224  thedl.tolerance = 0;
225 
226  LWDEBUG(2, "lw_dist3d_distancepoint is called");
227 
228  /*Check if we really have 3D geometries*/
229  /*If not, send it to 2D-calculations which will give the same result*/
230  /*as an infinite z-value at one or two of the geometries*/
231  if(!lwgeom_has_z(lw1) || !lwgeom_has_z(lw2))
232  {
233  lwnotice("One or both of the geometries is missing z-value. The unknown z-value will be regarded as \"any value\"");
234 
235  if(!lwgeom_has_z(lw1) && !lwgeom_has_z(lw2))
236  return lw_dist2d_distancepoint(lw1, lw2, srid, mode);
237 
238 
239  DISTPTS thedl2d;
240  thedl2d.mode = mode;
241  thedl2d.distance = initdistance;
242  thedl2d.tolerance = 0.0;
243  if (!lw_dist2d_comp( lw1,lw2,&thedl2d))
244  {
245  /*should never get here. all cases ought to be error handled earlier*/
246  lwerror("Some unspecified error.");
247  return (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
248  }
249 
250  LWGEOM *vertical_line;
251  if(!lwgeom_has_z(lw1))
252  {
253  x=thedl2d.p1.x;
254  y=thedl2d.p1.y;
255 
256  vertical_line = create_v_line(lw2,x,y,srid);
257  if (!lw_dist3d_recursive(vertical_line, lw2, &thedl))
258  {
259  /*should never get here. all cases ought to be error handled earlier*/
260  lwfree(vertical_line);
261  lwerror("Some unspecified error.");
262  return (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
263  }
264  lwfree(vertical_line);
265  }
266 
267  if(!lwgeom_has_z(lw2))
268  {
269  x=thedl2d.p2.x;
270  y=thedl2d.p2.y;
271 
272  vertical_line = create_v_line(lw1,x,y,srid);
273  if (!lw_dist3d_recursive(lw1, vertical_line, &thedl))
274  {
275  /*should never get here. all cases ought to be error handled earlier*/
276  lwfree(vertical_line);
277  lwerror("Some unspecified error.");
278  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
279  }
280  lwfree(vertical_line);
281  }
282 
283  }
284  else
285  {
286  if (!lw_dist3d_recursive(lw1, lw2, &thedl))
287  {
288  /*should never get here. all cases ought to be error handled earlier*/
289  lwerror("Some unspecified error.");
290  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
291  }
292  }
293  if (thedl.distance == initdistance)
294  {
295  LWDEBUG(3, "didn't find geometries to measure between, returning null");
296  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
297  }
298  else
299  {
300  x=thedl.p1.x;
301  y=thedl.p1.y;
302  z=thedl.p1.z;
303  result = (LWGEOM *)lwpoint_make3dz(srid, x, y, z);
304  }
305 
306  return result;
307 }
308 
309 
313 double
314 lwgeom_maxdistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
315 {
316  LWDEBUG(2, "lwgeom_maxdistance3d is called");
317 
318  return lwgeom_maxdistance3d_tolerance( lw1, lw2, 0.0 );
319 }
320 
325 double
326 lwgeom_maxdistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
327 {
328  if(!lwgeom_has_z(lw1) || !lwgeom_has_z(lw2))
329  {
330  lwnotice("One or both of the geometries is missing z-value. The unknown z-value will be regarded as \"any value\"");
331  return lwgeom_maxdistance2d_tolerance(lw1, lw2, tolerance);
332  }
333  /*double thedist;*/
334  DISTPTS3D thedl;
335  LWDEBUG(2, "lwgeom_maxdistance3d_tolerance is called");
336  thedl.mode = DIST_MAX;
337  thedl.distance= -1;
338  thedl.tolerance = tolerance;
339  if (lw_dist3d_recursive(lw1, lw2, &thedl))
340  {
341  return thedl.distance;
342  }
343  /*should never get here. all cases ought to be error handled earlier*/
344  lwerror("Some unspecified error.");
345  return -1;
346 }
347 
351 double
352 lwgeom_mindistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
353 {
354  LWDEBUG(2, "lwgeom_mindistance3d is called");
355  return lwgeom_mindistance3d_tolerance( lw1, lw2, 0.0 );
356 }
357 
362 double
363 lwgeom_mindistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
364 {
365  if(!lwgeom_has_z(lw1) || !lwgeom_has_z(lw2))
366  {
367  lwnotice("One or both of the geometries is missing z-value. The unknown z-value will be regarded as \"any value\"");
368 
369  return lwgeom_mindistance2d_tolerance(lw1, lw2, tolerance);
370  }
371  DISTPTS3D thedl;
372  LWDEBUG(2, "lwgeom_mindistance3d_tolerance is called");
373  thedl.mode = DIST_MIN;
374  thedl.distance= FLT_MAX;
375  thedl.tolerance = tolerance;
376  if (lw_dist3d_recursive(lw1, lw2, &thedl))
377  {
378  return thedl.distance;
379  }
380  /*should never get here. all cases ought to be error handled earlier*/
381  lwerror("Some unspecified error.");
382  return FLT_MAX;
383 }
384 
385 
386 /*------------------------------------------------------------------------------------------------------------
387 End of Initializing functions
388 --------------------------------------------------------------------------------------------------------------*/
389 
390 /*------------------------------------------------------------------------------------------------------------
391 Preprocessing functions
392 Functions preparing geometries for distance-calculations
393 --------------------------------------------------------------------------------------------------------------*/
394 
395 
399 int lw_dist3d_recursive(const LWGEOM *lwg1,const LWGEOM *lwg2, DISTPTS3D *dl)
400 {
401  int i, j;
402  int n1=1;
403  int n2=1;
404  LWGEOM *g1 = NULL;
405  LWGEOM *g2 = NULL;
406  LWCOLLECTION *c1 = NULL;
407  LWCOLLECTION *c2 = NULL;
408 
409  LWDEBUGF(2, "lw_dist3d_recursive is called with type1=%d, type2=%d", lwg1->type, lwg2->type);
410 
411  if (lwgeom_is_collection(lwg1))
412  {
413  LWDEBUG(3, "First geometry is collection");
414  c1 = lwgeom_as_lwcollection(lwg1);
415  n1 = c1->ngeoms;
416  }
417  if (lwgeom_is_collection(lwg2))
418  {
419  LWDEBUG(3, "Second geometry is collection");
420  c2 = lwgeom_as_lwcollection(lwg2);
421  n2 = c2->ngeoms;
422  }
423 
424  for ( i = 0; i < n1; i++ )
425  {
426 
427  if (lwgeom_is_collection(lwg1))
428  {
429  g1 = c1->geoms[i];
430  }
431  else
432  {
433  g1 = (LWGEOM*)lwg1;
434  }
435 
436  if (lwgeom_is_empty(g1)) return LW_TRUE;
437 
438  if (lwgeom_is_collection(g1))
439  {
440  LWDEBUG(3, "Found collection inside first geometry collection, recursing");
441  if (!lw_dist3d_recursive(g1, lwg2, dl)) return LW_FALSE;
442  continue;
443  }
444  for ( j = 0; j < n2; j++ )
445  {
446  if (lwgeom_is_collection(lwg2))
447  {
448  g2 = c2->geoms[j];
449  }
450  else
451  {
452  g2 = (LWGEOM*)lwg2;
453  }
454  if (lwgeom_is_collection(g2))
455  {
456  LWDEBUG(3, "Found collection inside second geometry collection, recursing");
457  if (!lw_dist3d_recursive(g1, g2, dl)) return LW_FALSE;
458  continue;
459  }
460 
461 
462  /*If one of geometries is empty, return. True here only means continue searching. False would have stopped the process*/
463  if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
464 
465 
466  if (!lw_dist3d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
467  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
468  }
469  }
470  return LW_TRUE;
471 }
472 
473 
474 
479 int
481 {
482 
483  int t1 = lwg1->type;
484  int t2 = lwg2->type;
485 
486  LWDEBUGF(2, "lw_dist3d_distribute_bruteforce is called with typ1=%d, type2=%d", lwg1->type, lwg2->type);
487 
488  if ( t1 == POINTTYPE )
489  {
490  if ( t2 == POINTTYPE )
491  {
492  dl->twisted=1;
493  return lw_dist3d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
494  }
495  else if ( t2 == LINETYPE )
496  {
497  dl->twisted=1;
498  return lw_dist3d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
499  }
500  else if ( t2 == POLYGONTYPE )
501  {
502  dl->twisted=1;
503  return lw_dist3d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2,dl);
504  }
505  else
506  {
507  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
508  return LW_FALSE;
509  }
510  }
511  else if ( t1 == LINETYPE )
512  {
513  if ( t2 == POINTTYPE )
514  {
515  dl->twisted=(-1);
516  return lw_dist3d_point_line((LWPOINT *)lwg2,(LWLINE *)lwg1,dl);
517  }
518  else if ( t2 == LINETYPE )
519  {
520  dl->twisted=1;
521  return lw_dist3d_line_line((LWLINE *)lwg1,(LWLINE *)lwg2,dl);
522  }
523  else if ( t2 == POLYGONTYPE )
524  {
525  dl->twisted=1;
526  return lw_dist3d_line_poly((LWLINE *)lwg1,(LWPOLY *)lwg2,dl);
527  }
528  else
529  {
530  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
531  return LW_FALSE;
532  }
533  }
534  else if ( t1 == POLYGONTYPE )
535  {
536  if ( t2 == POLYGONTYPE )
537  {
538  dl->twisted=1;
539  return lw_dist3d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2,dl);
540  }
541  else if ( t2 == POINTTYPE )
542  {
543  dl->twisted=-1;
544  return lw_dist3d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1,dl);
545  }
546  else if ( t2 == LINETYPE )
547  {
548  dl->twisted=-1;
549  return lw_dist3d_line_poly((LWLINE *)lwg2,(LWPOLY *)lwg1,dl);
550  }
551  else
552  {
553  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
554  return LW_FALSE;
555  }
556  }
557  else
558  {
559  lwerror("Unsupported geometry type: %s", lwtype_name(t1));
560  return LW_FALSE;
561  }
562 }
563 
564 
565 
566 /*------------------------------------------------------------------------------------------------------------
567 End of Preprocessing functions
568 --------------------------------------------------------------------------------------------------------------*/
569 
570 
571 /*------------------------------------------------------------------------------------------------------------
572 Brute force functions
573 So far the only way to do 3D-calculations
574 --------------------------------------------------------------------------------------------------------------*/
575 
580 int
582 {
583  POINT3DZ p1;
584  POINT3DZ p2;
585  LWDEBUG(2, "lw_dist3d_point_point is called");
586 
587  getPoint3dz_p(point1->point, 0, &p1);
588  getPoint3dz_p(point2->point, 0, &p2);
589 
590  return lw_dist3d_pt_pt(&p1, &p2,dl);
591 }
596 int
598 {
599  POINT3DZ p;
600  POINTARRAY *pa = line->points;
601  LWDEBUG(2, "lw_dist3d_point_line is called");
602 
603  getPoint3dz_p(point->point, 0, &p);
604  return lw_dist3d_pt_ptarray(&p, pa, dl);
605 }
606 
618 int
620 {
621  POINT3DZ p, projp;/*projp is "point projected on plane"*/
622  PLANE3D plane;
623  LWDEBUG(2, "lw_dist3d_point_poly is called");
624  getPoint3dz_p(point->point, 0, &p);
625 
626  /*If we are looking for max distance, longestline or dfullywithin*/
627  if (dl->mode == DIST_MAX)
628  {
629  LWDEBUG(3, "looking for maxdistance");
630  return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl);
631  }
632 
633  /*Find the plane of the polygon, the "holes" have to be on the same plane. so we only care about the boudary*/
634  if(!define_plane(poly->rings[0], &plane))
635  return LW_FALSE;
636 
637  /*get our point projected on the plane of the polygon*/
638  project_point_on_plane(&p, &plane, &projp);
639 
640  return lw_dist3d_pt_poly(&p, poly,&plane, &projp, dl);
641 }
642 
643 
648 int
650 {
651  POINTARRAY *pa1 = line1->points;
652  POINTARRAY *pa2 = line2->points;
653  LWDEBUG(2, "lw_dist3d_line_line is called");
654 
655  return lw_dist3d_ptarray_ptarray(pa1, pa2, dl);
656 }
657 
663 {
664  PLANE3D plane;
665  LWDEBUG(2, "lw_dist3d_line_poly is called");
666 
667  if (dl->mode == DIST_MAX)
668  {
669  return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);
670  }
671 
672  if(!define_plane(poly->rings[0], &plane))
673  return LW_FALSE;
674 
675  return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl);
676 }
677 
682 int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
683 {
684  PLANE3D plane;
685  LWDEBUG(2, "lw_dist3d_poly_poly is called");
686  if (dl->mode == DIST_MAX)
687  {
688  return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
689  }
690 
691  if(!define_plane(poly2->rings[0], &plane))
692  return LW_FALSE;
693 
694  /*What we do here is to compare the boundary of one polygon with the other polygon
695  and then take the second boundary comparing with the first polygon*/
696  dl->twisted=1;
697  if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl))
698  return LW_FALSE;
699  if(dl->distance==0.0) /*Just check if the answer already is given*/
700  return LW_TRUE;
701 
702  if(!define_plane(poly1->rings[0], &plane))
703  return LW_FALSE;
704  dl->twisted=-1; /*because we switch the order of geometries we switch "twisted" to -1 which will give the right order of points in shortest line.*/
705  return lw_dist3d_ptarray_poly(poly2->rings[0], poly1,&plane, dl);
706 }
707 
713 int
715 {
716  uint32_t t;
717  POINT3DZ start, end;
718  int twist = dl->twisted;
719 
720  LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
721 
722  getPoint3dz_p(pa, 0, &start);
723 
724  for (t=1; t<pa->npoints; t++)
725  {
726  dl->twisted=twist;
727  getPoint3dz_p(pa, t, &end);
728  if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
729 
730  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
731  start = end;
732  }
733 
734  return LW_TRUE;
735 }
736 
737 
743 int
745 {
746  POINT3DZ c;
747  double r;
748  /*if start==end, then use pt distance */
749  if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )
750  {
751  return lw_dist3d_pt_pt(p,A,dl);
752  }
753 
754 
755  r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) + ( p->z-A->z) * (B->z-A->z) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y)+(B->z-A->z)*(B->z-A->z) );
756 
757  /*This is for finding the 3Dmaxdistance.
758  the maxdistance have to be between two vertexes,
759  compared to mindistance which can be between
760  two vertexes vertex.*/
761  if (dl->mode == DIST_MAX)
762  {
763  if (r>=0.5)
764  {
765  return lw_dist3d_pt_pt(p,A,dl);
766  }
767  if (r<0.5)
768  {
769  return lw_dist3d_pt_pt(p,B,dl);
770  }
771  }
772 
773  if (r<0) /*If the first vertex A is closest to the point p*/
774  {
775  return lw_dist3d_pt_pt(p,A,dl);
776  }
777  if (r>1) /*If the second vertex B is closest to the point p*/
778  {
779  return lw_dist3d_pt_pt(p,B,dl);
780  }
781 
782  /*else if the point p is closer to some point between a and b
783  then we find that point and send it to lw_dist3d_pt_pt*/
784  c.x=A->x + r * (B->x-A->x);
785  c.y=A->y + r * (B->y-A->y);
786  c.z=A->z + r * (B->z-A->z);
787 
788  return lw_dist3d_pt_pt(p,&c,dl);
789 }
790 
791 double
792 distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
793 {
794  double dx = p2->x - p1->x;
795  double dy = p2->y - p1->y;
796  double dz = p2->z - p1->z;
797  return sqrt ( dx*dx + dy*dy + dz*dz);
798 }
799 
800 
808 int
810 {
811  double dx = thep2->x - thep1->x;
812  double dy = thep2->y - thep1->y;
813  double dz = thep2->z - thep1->z;
814  double dist = sqrt ( dx*dx + dy*dy + dz*dz);
815  LWDEBUGF(2, "lw_dist3d_pt_pt called (with points: p1.x=%f, p1.y=%f,p1.z=%f,p2.x=%f, p2.y=%f,p2.z=%f)",thep1->x,thep1->y,thep1->z,thep2->x,thep2->y,thep2->z );
816 
817  if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
818  {
819  dl->distance = dist;
820 
821  if (dl->twisted>0) /*To get the points in right order. twisted is updated between 1 and (-1) every time the order is changed earlier in the chain*/
822  {
823  dl->p1 = *thep1;
824  dl->p2 = *thep2;
825  }
826  else
827  {
828  dl->p1 = *thep2;
829  dl->p2 = *thep1;
830  }
831  }
832  return LW_TRUE;
833 }
834 
835 
840 int
842 {
843  uint32_t t,u;
844  POINT3DZ start, end;
845  POINT3DZ start2, end2;
846  int twist = dl->twisted;
847  LWDEBUGF(2, "lw_dist3d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
848 
849 
850 
851  if (dl->mode == DIST_MAX)/*If we are searching for maxdistance we go straight to point-point calculation since the maxdistance have to be between two vertexes*/
852  {
853  for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
854  {
855  getPoint3dz_p(l1, t, &start);
856  for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
857  {
858  getPoint3dz_p(l2, u, &start2);
859  lw_dist3d_pt_pt(&start,&start2,dl);
860  LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
861  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
862  t, u, dl->distance, dl->tolerance);
863  }
864  }
865  }
866  else
867  {
868  getPoint3dz_p(l1, 0, &start);
869  for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
870  {
871  getPoint3dz_p(l1, t, &end);
872  getPoint3dz_p(l2, 0, &start2);
873  for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
874  {
875  getPoint3dz_p(l2, u, &end2);
876  dl->twisted=twist;
877  lw_dist3d_seg_seg(&start, &end, &start2, &end2,dl);
878  LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
879  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
880  t, u, dl->distance, dl->tolerance);
881  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
882  start2 = end2;
883  }
884  start = end;
885  }
886  }
887  return LW_TRUE;
888 }
889 
894 int
896 {
897  VECTOR3D v1, v2, vl;
898  double s1k, s2k; /*two variables representing where on Line 1 (s1k) and where on Line 2 (s2k) a connecting line between the two lines is perpendicular to both lines*/
899  POINT3DZ p1, p2;
900  double a, b, c, d, e, D;
901 
902  /*s1p1 and s1p2 are the same point */
903  if ( ( s1p1->x == s1p2->x) && (s1p1->y == s1p2->y) && (s1p1->z == s1p2->z) )
904  {
905  return lw_dist3d_pt_seg(s1p1,s2p1,s2p2,dl);
906  }
907  /*s2p1 and s2p2 are the same point */
908  if ( ( s2p1->x == s2p2->x) && (s2p1->y == s2p2->y) && (s2p1->z == s2p2->z) )
909  {
910  dl->twisted= ((dl->twisted) * (-1));
911  return lw_dist3d_pt_seg(s2p1,s1p1,s1p2,dl);
912  }
913 
914 /*
915  Here we use algorithm from softsurfer.com
916  that can be found here
917  http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
918 */
919 
920  if (!get_3dvector_from_points(s1p1, s1p2, &v1))
921  return LW_FALSE;
922 
923  if (!get_3dvector_from_points(s2p1, s2p2, &v2))
924  return LW_FALSE;
925 
926  if (!get_3dvector_from_points(s2p1, s1p1, &vl))
927  return LW_FALSE;
928 
929  a = DOT(v1,v1);
930  b = DOT(v1,v2);
931  c = DOT(v2,v2);
932  d = DOT(v1,vl);
933  e = DOT(v2,vl);
934  D = a*c - b*b;
935 
936 
937  if (D <0.000000001)
938  { /* the lines are almost parallel*/
939  s1k = 0.0; /*If the lines are parallel we try by using the startpoint of first segment. If that gives a projected point on the second line outside segment 2 it wil be found that s2k is >1 or <0.*/
940  if(b>c) /* use the largest denominator*/
941  {
942  s2k=d/b;
943  }
944  else
945  {
946  s2k =e/c;
947  }
948  }
949  else
950  {
951  s1k = (b*e - c*d) / D;
952  s2k = (a*e - b*d) / D;
953  }
954 
955  /* Now we check if the projected closest point on the infinite lines is outside our segments. If so the combinations with start and end points will be tested*/
956  if(s1k<0.0||s1k>1.0||s2k<0.0||s2k>1.0)
957  {
958  if(s1k<0.0)
959  {
960 
961  if (!lw_dist3d_pt_seg(s1p1, s2p1, s2p2, dl))
962  {
963  return LW_FALSE;
964  }
965  }
966  if(s1k>1.0)
967  {
968 
969  if (!lw_dist3d_pt_seg(s1p2, s2p1, s2p2, dl))
970  {
971  return LW_FALSE;
972  }
973  }
974  if(s2k<0.0)
975  {
976  dl->twisted= ((dl->twisted) * (-1));
977  if (!lw_dist3d_pt_seg(s2p1, s1p1, s1p2, dl))
978  {
979  return LW_FALSE;
980  }
981  }
982  if(s2k>1.0)
983  {
984  dl->twisted= ((dl->twisted) * (-1));
985  if (!lw_dist3d_pt_seg(s2p2, s1p1, s1p2, dl))
986  {
987  return LW_FALSE;
988  }
989  }
990  }
991  else
992  {/*Find the closest point on the edges of both segments*/
993  p1.x=s1p1->x+s1k*(s1p2->x-s1p1->x);
994  p1.y=s1p1->y+s1k*(s1p2->y-s1p1->y);
995  p1.z=s1p1->z+s1k*(s1p2->z-s1p1->z);
996 
997  p2.x=s2p1->x+s2k*(s2p2->x-s2p1->x);
998  p2.y=s2p1->y+s2k*(s2p2->y-s2p1->y);
999  p2.z=s2p1->z+s2k*(s2p2->z-s2p1->z);
1000 
1001  if (!lw_dist3d_pt_pt(&p1,&p2,dl))/* Send the closest points to point-point calculation*/
1002  {
1003  return LW_FALSE;
1004  }
1005  }
1006  return LW_TRUE;
1007 }
1008 
1016 int
1018 {
1019  uint32_t i;
1020 
1021  LWDEBUG(2, "lw_dist3d_point_poly called");
1022 
1023 
1024  if(pt_in_ring_3d(projp, poly->rings[0], plane))
1025  {
1026  for (i=1; i<poly->nrings; i++)
1027  {
1028  /* Inside a hole. Distance = pt -> ring */
1029  if ( pt_in_ring_3d(projp, poly->rings[i], plane ))
1030  {
1031  LWDEBUG(3, " inside an hole");
1032  return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);
1033  }
1034  }
1035 
1036  return lw_dist3d_pt_pt(p,projp,dl);/* If the projected point is inside the polygon the shortest distance is between that point and the inputed point*/
1037  }
1038  else
1039  {
1040  return lw_dist3d_pt_ptarray(p, poly->rings[0], dl); /*If the projected point is outside the polygon we search for the closest distance against the boundary instead*/
1041  }
1042 
1043  return LW_TRUE;
1044 
1045 }
1046 
1052 {
1053 
1054 
1055  uint32_t i,j,k;
1056  double f, s1, s2;
1057  VECTOR3D projp1_projp2;
1058  POINT3DZ p1, p2,projp1, projp2, intersectionp;
1059 
1060  getPoint3dz_p(pa, 0, &p1);
1061 
1062  s1=project_point_on_plane(&p1, plane, &projp1); /*the sign of s1 tells us on which side of the plane the point is. */
1063  lw_dist3d_pt_poly(&p1, poly, plane,&projp1, dl);
1064 
1065  for (i=1;i<pa->npoints;i++)
1066  {
1067  int intersects;
1068  getPoint3dz_p(pa, i, &p2);
1069  s2=project_point_on_plane(&p2, plane, &projp2);
1070  lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl);
1071 
1072  /*If s1and s2 has different signs that means they are on different sides of the plane of the polygon.
1073  That means that the edge between the points crosses the plane and might intersect with the polygon*/
1074  if((s1*s2)<=0)
1075  {
1076  f=fabs(s1)/(fabs(s1)+fabs(s2)); /*The size of s1 and s2 is the distance from the point to the plane.*/
1077  get_3dvector_from_points(&projp1, &projp2,&projp1_projp2);
1078 
1079  /*get the point where the line segment crosses the plane*/
1080  intersectionp.x=projp1.x+f*projp1_projp2.x;
1081  intersectionp.y=projp1.y+f*projp1_projp2.y;
1082  intersectionp.z=projp1.z+f*projp1_projp2.z;
1083 
1084  intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/
1085 
1086  if(pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/
1087  {
1088  for (k=1;k<poly->nrings; k++)
1089  {
1090  /* Inside a hole, so no intersection with the polygon*/
1091  if ( pt_in_ring_3d(&intersectionp, poly->rings[k], plane ))
1092  {
1093  intersects=LW_FALSE;
1094  break;
1095  }
1096  }
1097  if(intersects)
1098  {
1099  dl->distance=0.0;
1100  dl->p1.x=intersectionp.x;
1101  dl->p1.y=intersectionp.y;
1102  dl->p1.z=intersectionp.z;
1103 
1104  dl->p2.x=intersectionp.x;
1105  dl->p2.y=intersectionp.y;
1106  dl->p2.z=intersectionp.z;
1107  return LW_TRUE;
1108 
1109  }
1110  }
1111  }
1112 
1113  projp1=projp2;
1114  s1=s2;
1115  p1=p2;
1116  }
1117 
1118  /*check or pointarray against boundary and inner boundaries of the polygon*/
1119  for (j=0;j<poly->nrings;j++)
1120  {
1121  lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);
1122  }
1123 
1124 return LW_TRUE;
1125 }
1126 
1127 
1133 int
1135 {
1136  uint32_t i,j, numberofvectors, pointsinslice;
1137  POINT3DZ p, p1, p2;
1138 
1139  double sumx=0;
1140  double sumy=0;
1141  double sumz=0;
1142  double vl; /*vector length*/
1143 
1144  VECTOR3D v1, v2, v;
1145 
1146  if((pa->npoints-1)==3) /*Triangle is special case*/
1147  {
1148  pointsinslice=1;
1149  }
1150  else
1151  {
1152  pointsinslice=(uint32_t) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/
1153  }
1154 
1155  /*find the avg point*/
1156  for (i=0;i<(pa->npoints-1);i++)
1157  {
1158  getPoint3dz_p(pa, i, &p);
1159  sumx+=p.x;
1160  sumy+=p.y;
1161  sumz+=p.z;
1162  }
1163  pl->pop.x=(sumx/(pa->npoints-1));
1164  pl->pop.y=(sumy/(pa->npoints-1));
1165  pl->pop.z=(sumz/(pa->npoints-1));
1166 
1167  sumx=0;
1168  sumy=0;
1169  sumz=0;
1170  numberofvectors= floor((pa->npoints-1)/pointsinslice); /*the number of vectors we try can be 3, 4 or 5*/
1171 
1172  getPoint3dz_p(pa, 0, &p1);
1173  for (j=pointsinslice;j<pa->npoints;j+=pointsinslice)
1174  {
1175  getPoint3dz_p(pa, j, &p2);
1176 
1177  if (!get_3dvector_from_points(&(pl->pop), &p1, &v1) || !get_3dvector_from_points(&(pl->pop), &p2, &v2))
1178  return LW_FALSE;
1179  /*perpendicular vector is cross product of v1 and v2*/
1180  if (!get_3dcross_product(&v1,&v2, &v))
1181  return LW_FALSE;
1182  vl=VECTORLENGTH(v);
1183  sumx+=(v.x/vl);
1184  sumy+=(v.y/vl);
1185  sumz+=(v.z/vl);
1186  p1=p2;
1187  }
1188  pl->pv.x=(sumx/numberofvectors);
1189  pl->pv.y=(sumy/numberofvectors);
1190  pl->pv.z=(sumz/numberofvectors);
1191 
1192  return 1;
1193 }
1194 
1199 double
1201 {
1202 /*In our plane definition we have a point on the plane and a normal vector (pl.pv), perpendicular to the plane
1203 this vector will be parallel to the line between our inputted point above the plane and the point we are searching for on the plane.
1204 So, we already have a direction from p to find p0, but we don't know the distance.
1205 */
1206 
1207  VECTOR3D v1;
1208  double f;
1209 
1210  if (!get_3dvector_from_points(&(pl->pop), p, &v1))
1211  return LW_FALSE;
1212 
1213  f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv));
1214 
1215  p0->x=p->x+pl->pv.x*f;
1216  p0->y=p->y+pl->pv.y*f;
1217  p0->z=p->z+pl->pv.z*f;
1218 
1219  return f;
1220 }
1221 
1222 
1223 
1224 
1237 int
1238 pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring,PLANE3D *plane)
1239 {
1240 
1241  uint32_t cn = 0; /* the crossing number counter */
1242  uint32_t i;
1243  POINT3DZ v1, v2;
1244 
1245  POINT3DZ first, last;
1246 
1247  getPoint3dz_p(ring, 0, &first);
1248  getPoint3dz_p(ring, ring->npoints-1, &last);
1249  if ( memcmp(&first, &last, sizeof(POINT3DZ)) )
1250  {
1251  lwerror("pt_in_ring_3d: V[n] != V[0] (%g %g %g!= %g %g %g)",
1252  first.x, first.y, first.z, last.x, last.y, last.z);
1253  return LW_FALSE;
1254  }
1255 
1256  LWDEBUGF(2, "pt_in_ring_3d called with point: %g %g %g", p->x, p->y, p->z);
1257  /* printPA(ring); */
1258 
1259  /* loop through all edges of the polygon */
1260  getPoint3dz_p(ring, 0, &v1);
1261 
1262 
1263  if(fabs(plane->pv.z)>=fabs(plane->pv.x)&&fabs(plane->pv.z)>=fabs(plane->pv.y)) /*If the z vector of the normal vector to the plane is larger than x and y vector we project the ring to the xy-plane*/
1264  {
1265  for (i=0; i<ring->npoints-1; i++)
1266  {
1267  double vt;
1268  getPoint3dz_p(ring, i+1, &v2);
1269 
1270  /* edge from vertex i to vertex i+1 */
1271  if
1272  (
1273  /* an upward crossing */
1274  ((v1.y <= p->y) && (v2.y > p->y))
1275  /* a downward crossing */
1276  || ((v1.y > p->y) && (v2.y <= p->y))
1277  )
1278  {
1279 
1280  vt = (double)(p->y - v1.y) / (v2.y - v1.y);
1281 
1282  /* P.x <intersect */
1283  if (p->x < v1.x + vt * (v2.x - v1.x))
1284  {
1285  /* a valid crossing of y=p.y right of p.x */
1286  ++cn;
1287  }
1288  }
1289  v1 = v2;
1290  }
1291  }
1292  else if(fabs(plane->pv.y)>=fabs(plane->pv.x)&&fabs(plane->pv.y)>=fabs(plane->pv.z)) /*If the y vector of the normal vector to the plane is larger than x and z vector we project the ring to the xz-plane*/
1293  {
1294  for (i=0; i<ring->npoints-1; i++)
1295  {
1296  double vt;
1297  getPoint3dz_p(ring, i+1, &v2);
1298 
1299  /* edge from vertex i to vertex i+1 */
1300  if
1301  (
1302  /* an upward crossing */
1303  ((v1.z <= p->z) && (v2.z > p->z))
1304  /* a downward crossing */
1305  || ((v1.z > p->z) && (v2.z <= p->z))
1306  )
1307  {
1308 
1309  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1310 
1311  /* P.x <intersect */
1312  if (p->x < v1.x + vt * (v2.x - v1.x))
1313  {
1314  /* a valid crossing of y=p.y right of p.x */
1315  ++cn;
1316  }
1317  }
1318  v1 = v2;
1319  }
1320  }
1321  else /*Hopefully we only have the cases where x part of the normal vector is largest left*/
1322  {
1323  for (i=0; i<ring->npoints-1; i++)
1324  {
1325  double vt;
1326  getPoint3dz_p(ring, i+1, &v2);
1327 
1328  /* edge from vertex i to vertex i+1 */
1329  if
1330  (
1331  /* an upward crossing */
1332  ((v1.z <= p->z) && (v2.z > p->z))
1333  /* a downward crossing */
1334  || ((v1.z > p->z) && (v2.z <= p->z))
1335  )
1336  {
1337 
1338  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1339 
1340  /* P.x <intersect */
1341  if (p->y < v1.y + vt * (v2.y - v1.y))
1342  {
1343  /* a valid crossing of y=p.y right of p.x */
1344  ++cn;
1345  }
1346  }
1347  v1 = v2;
1348  }
1349  }
1350  LWDEBUGF(3, "pt_in_ring_3d returning %d", cn&1);
1351 
1352  return (cn&1); /* 0 if even (out), and 1 if odd (in) */
1353 }
1354 
1355 
1356 
1357 /*------------------------------------------------------------------------------------------------------------
1358 End of Brute force functions
1359 --------------------------------------------------------------------------------------------------------------*/
1360 
1361 
1362 
#define LINETYPE
Definition: liblwgeom.h:85
double z
Definition: liblwgeom.h:336
double y
Definition: liblwgeom.h:336
double distance
Definition: measures3d.h:41
LWGEOM * lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing closestpoint calculations.
Definition: measures.c:131
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1085
int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
line to polygon calculation
Definition: measures3d.c:662
int lw_dist2d_comp(const LWGEOM *lw1, const LWGEOM *lw2, DISTPTS *dl)
This function just deserializes geometries Bboxes is not checked here since it is the subgeometries b...
Definition: measures.c:245
double x
Definition: liblwgeom.h:336
char * r
Definition: cu_in_wkt.c:24
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
double z
Definition: measures3d.h:51
double y
Definition: liblwgeom.h:342
double lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling max distance calculations and dfullywithin calculations.
Definition: measures.c:181
POINT3DZ p2
Definition: measures3d.h:43
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:213
#define POLYGONTYPE
Definition: liblwgeom.h:86
int mode
Definition: measures.h:51
LWLINE * lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points)
Definition: lwline.c:237
int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
polygon to polygon calculation
Definition: measures3d.c:682
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:792
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
double x
Definition: liblwgeom.h:342
VECTOR3D pv
Definition: measures3d.h:58
uint32_t ngeoms
Definition: liblwgeom.h:509
POINTARRAY * point
Definition: liblwgeom.h:413
Structure used in distance-calculations.
Definition: measures3d.h:39
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:930
int32_t srid
Definition: liblwgeom.h:401
#define VECTORLENGTH(v)
Definition: measures3d.h:32
int lw_dist3d_seg_seg(POINT3DZ *s1p1, POINT3DZ *s1p2, POINT3DZ *s2p1, POINT3DZ *s2p2, DISTPTS3D *dl)
Finds the two closest points on two linesegments.
Definition: measures3d.c:895
double project_point_on_plane(POINT3DZ *p, PLANE3D *pl, POINT3DZ *p0)
Finds a point on a plane from where the original point is perpendicular to the plane.
Definition: measures3d.c:1200
POINT3DZ p1
Definition: measures3d.h:42
double y
Definition: measures3d.h:51
uint32_t nrings
Definition: liblwgeom.h:457
#define LW_FAILURE
Definition: liblwgeom.h:78
LWGEOM * lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing shortestline and longestline calculations.
Definition: measures.c:84
POINT2D p1
Definition: measures.h:49
double z
Definition: liblwgeom.h:342
#define DOT(u, v)
Definition: measures3d.h:31
unsigned int uint32_t
Definition: uthash.h:78
double tolerance
Definition: measures.h:53
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:746
LWGEOM * lwgeom_closest_point_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:100
double x
Definition: liblwgeom.h:330
int lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
This function distributes the brute-force for 3D so far the only type, tasks depending on type...
Definition: measures3d.c:480
double zmax
Definition: liblwgeom.h:299
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
int lw_dist3d_pt_poly(POINT3DZ *p, LWPOLY *poly, PLANE3D *plane, POINT3DZ *projp, DISTPTS3D *dl)
Checking if the point projected on the plane of the polygon actually is inside that polygon...
Definition: measures3d.c:1017
#define DIST_MIN
Definition: measures.h:41
LWGEOM * lw_dist3d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing 3dshortestline and 3dlongestline calculations.
Definition: measures3d.c:110
#define LW_FALSE
Definition: liblwgeom.h:76
int lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)
Computes point to polygon distance For mindistance that means: 1)find the plane of the polygon 2)proj...
Definition: measures3d.c:619
Datum intersects(PG_FUNCTION_ARGS)
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
int lw_dist3d_pt_pt(POINT3DZ *thep1, POINT3DZ *thep2, DISTPTS3D *dl)
Compares incoming points and stores the points closest to each other or most far away from each other...
Definition: measures3d.c:809
LWGEOM ** geoms
Definition: liblwgeom.h:511
LWGEOM * lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:88
LWGEOM * lw_dist3d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing 3dclosestpoint calculations.
Definition: measures3d.c:214
int twisted
Definition: measures3d.h:45
POINTARRAY ** rings
Definition: liblwgeom.h:459
LWPOINT * lwpoint_make3dz(int srid, double x, double y, double z)
Definition: lwpoint.c:173
static int get_3dvector_from_points(POINT3DZ *p1, POINT3DZ *p2, VECTOR3D *v)
Definition: measures3d.c:45
POINT2D p2
Definition: measures.h:50
int mode
Definition: measures3d.h:44
double y
Definition: liblwgeom.h:330
LWGEOM * lwgeom_furthest_line_3d(LWGEOM *lw1, LWGEOM *lw2)
Definition: measures3d.c:94
double lwgeom_mindistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling 3d min distance calculations and dwithin calculations.
Definition: measures3d.c:363
POINT3DZ pop
Definition: measures3d.h:57
int lw_dist3d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
This is a recursive function delivering every possible combination of subgeometries.
Definition: measures3d.c:399
double lwgeom_maxdistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling 3d max distance calculations and dfullywithin calculations.
Definition: measures3d.c:326
int lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)
point to line calculation
Definition: measures3d.c:597
int define_plane(POINTARRAY *pa, PLANE3D *pl)
Here we define the plane of a polygon (boundary pointarray of a polygon) the plane is stored as a poi...
Definition: measures3d.c:1134
int getPoint3dz_p(const POINTARRAY *pa, uint32_t n, POINT3DZ *point)
Definition: lwgeom_api.c:205
double tolerance
Definition: measures3d.h:46
int pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring, PLANE3D *plane)
pt_in_ring_3d(): crossing number test for a point in a polygon input: p = a point, pa = vertex points of a ring V[n+1] with V[n]=V[0] plane=the plane that the vertex points are lying on returns: 0 = outside, 1 = inside
Definition: measures3d.c:1238
double distance
Definition: measures.h:48
#define DIST_MAX
Definition: measures.h:40
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:224
double zmin
Definition: liblwgeom.h:298
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
int lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)
point to point calculation
Definition: measures3d.c:581
double lwgeom_maxdistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 3d max distance calculation.
Definition: measures3d.c:314
int lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl)
If searching for min distance, this one finds the closest point on segment A-B from p...
Definition: measures3d.c:744
uint8_t type
Definition: liblwgeom.h:398
double x
Definition: measures3d.h:51
static LWGEOM * create_v_line(const LWGEOM *lwgeom, double x, double y, int srid)
This function is used to create a vertical line used for cases where one if the geometries lacks z-va...
Definition: measures3d.c:71
double lwgeom_mindistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 3d min distance calculation.
Definition: measures3d.c:352
int lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS3D *dl)
Finds all combinations of segments between two pointarrays.
Definition: measures3d.c:841
int lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
line to line calculation
Definition: measures3d.c:649
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1393
Structure used in distance-calculations.
Definition: measures.h:46
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
int lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa, DISTPTS3D *dl)
search all the segments of pointarray to see which one is closest to p Returns distance between point...
Definition: measures3d.c:714
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, PLANE3D *plane, DISTPTS3D *dl)
Computes pointarray to polygon distance.
Definition: measures3d.c:1051
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
POINTARRAY * points
Definition: liblwgeom.h:424
static int get_3dcross_product(VECTOR3D *v1, VECTOR3D *v2, VECTOR3D *v)
Definition: measures3d.c:55
uint32_t npoints
Definition: liblwgeom.h:373