PostGIS  2.5.7dev-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 (!FP_IS_ZERO(v->x) || !FP_IS_ZERO(v->y) || !FP_IS_ZERO(v->z));
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 (!FP_IS_ZERO(v->x) || !FP_IS_ZERO(v->y) || !FP_IS_ZERO(v->z));
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  {
636  /* Polygon does not define a plane: Return distance point -> line */
637  return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl);
638  }
639 
640  /*get our point projected on the plane of the polygon*/
641  project_point_on_plane(&p, &plane, &projp);
642 
643  return lw_dist3d_pt_poly(&p, poly,&plane, &projp, dl);
644 }
645 
646 
651 int
653 {
654  POINTARRAY *pa1 = line1->points;
655  POINTARRAY *pa2 = line2->points;
656  LWDEBUG(2, "lw_dist3d_line_line is called");
657 
658  return lw_dist3d_ptarray_ptarray(pa1, pa2, dl);
659 }
660 
666 {
667  PLANE3D plane;
668  LWDEBUG(2, "lw_dist3d_line_poly is called");
669 
670  if (dl->mode == DIST_MAX)
671  {
672  return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);
673  }
674 
675  if(!define_plane(poly->rings[0], &plane))
676  {
677  /* Polygon does not define a plane: Return distance line to line */
678  return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);
679  }
680 
681  return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl);
682 }
683 
688 int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
689 {
690  PLANE3D plane1, plane2;
691  int planedef1, planedef2;
692  LWDEBUG(2, "lw_dist3d_poly_poly is called");
693  if (dl->mode == DIST_MAX)
694  {
695  return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
696  }
697 
698  planedef1 = define_plane(poly1->rings[0], &plane1);
699  planedef2 = define_plane(poly2->rings[0], &plane2);
700 
701  if (!planedef1 || !planedef2)
702  {
703  if (!planedef1 && !planedef2)
704  {
705  /* Neither polygon define a plane: Return distance line to line */
706  return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
707  }
708 
709  if (!planedef1)
710  {
711  /* Only poly2 defines a plane: Return distance from line (poly1) to poly2 */
712  return lw_dist3d_ptarray_poly(poly1->rings[0], poly2, &plane2, dl);
713  }
714 
715  /* Only poly1 defines a plane: Return distance from line (poly2) to poly1 */
716  return lw_dist3d_ptarray_poly(poly2->rings[0], poly1, &plane1, dl);
717  }
718 
719  /*What we do here is to compare the boundary of one polygon with the other polygon
720  and then take the second boundary comparing with the first polygon*/
721  dl->twisted=1;
722  if (!lw_dist3d_ptarray_poly(poly1->rings[0], poly2, &plane2, dl))
723  return LW_FALSE;
724  if (dl->distance < dl->tolerance) /*Just check if the answer already is given*/
725  return LW_TRUE;
726 
727  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.*/
728  return lw_dist3d_ptarray_poly(poly2->rings[0], poly1, &plane1, dl);
729 }
730 
736 int
738 {
739  uint32_t t;
740  POINT3DZ start, end;
741  int twist = dl->twisted;
742 
743  LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
744 
745  getPoint3dz_p(pa, 0, &start);
746 
747  for (t=1; t<pa->npoints; t++)
748  {
749  dl->twisted=twist;
750  getPoint3dz_p(pa, t, &end);
751  if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
752 
753  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
754  start = end;
755  }
756 
757  return LW_TRUE;
758 }
759 
760 
766 int
768 {
769  POINT3DZ c;
770  double r;
771  /*if start==end, then use pt distance */
772  if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )
773  {
774  return lw_dist3d_pt_pt(p,A,dl);
775  }
776 
777 
778  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) );
779 
780  /*This is for finding the 3Dmaxdistance.
781  the maxdistance have to be between two vertexes,
782  compared to mindistance which can be between
783  two vertexes vertex.*/
784  if (dl->mode == DIST_MAX)
785  {
786  if (r>=0.5)
787  {
788  return lw_dist3d_pt_pt(p,A,dl);
789  }
790  if (r<0.5)
791  {
792  return lw_dist3d_pt_pt(p,B,dl);
793  }
794  }
795 
796  if (r<0) /*If the first vertex A is closest to the point p*/
797  {
798  return lw_dist3d_pt_pt(p,A,dl);
799  }
800  if (r>1) /*If the second vertex B is closest to the point p*/
801  {
802  return lw_dist3d_pt_pt(p,B,dl);
803  }
804 
805  /*else if the point p is closer to some point between a and b
806  then we find that point and send it to lw_dist3d_pt_pt*/
807  c.x=A->x + r * (B->x-A->x);
808  c.y=A->y + r * (B->y-A->y);
809  c.z=A->z + r * (B->z-A->z);
810 
811  return lw_dist3d_pt_pt(p,&c,dl);
812 }
813 
814 double
815 distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
816 {
817  double dx = p2->x - p1->x;
818  double dy = p2->y - p1->y;
819  double dz = p2->z - p1->z;
820  return sqrt ( dx*dx + dy*dy + dz*dz);
821 }
822 
823 
831 int
833 {
834  double dx = thep2->x - thep1->x;
835  double dy = thep2->y - thep1->y;
836  double dz = thep2->z - thep1->z;
837  double dist = sqrt ( dx*dx + dy*dy + dz*dz);
838  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 );
839 
840  if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
841  {
842  dl->distance = dist;
843 
844  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*/
845  {
846  dl->p1 = *thep1;
847  dl->p2 = *thep2;
848  }
849  else
850  {
851  dl->p1 = *thep2;
852  dl->p2 = *thep1;
853  }
854  }
855  return LW_TRUE;
856 }
857 
858 
863 int
865 {
866  uint32_t t,u;
867  POINT3DZ start, end;
868  POINT3DZ start2, end2;
869  int twist = dl->twisted;
870  LWDEBUGF(2, "lw_dist3d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
871 
872 
873 
874  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*/
875  {
876  for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
877  {
878  getPoint3dz_p(l1, t, &start);
879  for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
880  {
881  getPoint3dz_p(l2, u, &start2);
882  lw_dist3d_pt_pt(&start,&start2,dl);
883  LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
884  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
885  t, u, dl->distance, dl->tolerance);
886  }
887  }
888  }
889  else
890  {
891  getPoint3dz_p(l1, 0, &start);
892  for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
893  {
894  getPoint3dz_p(l1, t, &end);
895  getPoint3dz_p(l2, 0, &start2);
896  for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
897  {
898  getPoint3dz_p(l2, u, &end2);
899  dl->twisted=twist;
900  lw_dist3d_seg_seg(&start, &end, &start2, &end2,dl);
901  LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
902  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
903  t, u, dl->distance, dl->tolerance);
904  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
905  start2 = end2;
906  }
907  start = end;
908  }
909  }
910  return LW_TRUE;
911 }
912 
917 int
919 {
920  VECTOR3D v1, v2, vl;
921  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*/
922  POINT3DZ p1, p2;
923  double a, b, c, d, e, D;
924 
925  /*s1p1 and s1p2 are the same point */
926  if ( ( s1p1->x == s1p2->x) && (s1p1->y == s1p2->y) && (s1p1->z == s1p2->z) )
927  {
928  return lw_dist3d_pt_seg(s1p1,s2p1,s2p2,dl);
929  }
930  /*s2p1 and s2p2 are the same point */
931  if ( ( s2p1->x == s2p2->x) && (s2p1->y == s2p2->y) && (s2p1->z == s2p2->z) )
932  {
933  dl->twisted= ((dl->twisted) * (-1));
934  return lw_dist3d_pt_seg(s2p1,s1p1,s1p2,dl);
935  }
936 
937 /*
938  Here we use algorithm from softsurfer.com
939  that can be found here
940  http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
941 */
942 
943  if (!get_3dvector_from_points(s1p1, s1p2, &v1))
944  return LW_FALSE;
945 
946  if (!get_3dvector_from_points(s2p1, s2p2, &v2))
947  return LW_FALSE;
948 
949  if (!get_3dvector_from_points(s2p1, s1p1, &vl))
950  return LW_FALSE;
951 
952  a = DOT(v1,v1);
953  b = DOT(v1,v2);
954  c = DOT(v2,v2);
955  d = DOT(v1,vl);
956  e = DOT(v2,vl);
957  D = a*c - b*b;
958 
959 
960  if (D <0.000000001)
961  { /* the lines are almost parallel*/
962  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.*/
963  if(b>c) /* use the largest denominator*/
964  {
965  s2k=d/b;
966  }
967  else
968  {
969  s2k =e/c;
970  }
971  }
972  else
973  {
974  s1k = (b*e - c*d) / D;
975  s2k = (a*e - b*d) / D;
976  }
977 
978  /* 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*/
979  if(s1k<0.0||s1k>1.0||s2k<0.0||s2k>1.0)
980  {
981  if(s1k<0.0)
982  {
983 
984  if (!lw_dist3d_pt_seg(s1p1, s2p1, s2p2, dl))
985  {
986  return LW_FALSE;
987  }
988  }
989  if(s1k>1.0)
990  {
991 
992  if (!lw_dist3d_pt_seg(s1p2, s2p1, s2p2, dl))
993  {
994  return LW_FALSE;
995  }
996  }
997  if(s2k<0.0)
998  {
999  dl->twisted= ((dl->twisted) * (-1));
1000  if (!lw_dist3d_pt_seg(s2p1, s1p1, s1p2, dl))
1001  {
1002  return LW_FALSE;
1003  }
1004  }
1005  if(s2k>1.0)
1006  {
1007  dl->twisted= ((dl->twisted) * (-1));
1008  if (!lw_dist3d_pt_seg(s2p2, s1p1, s1p2, dl))
1009  {
1010  return LW_FALSE;
1011  }
1012  }
1013  }
1014  else
1015  {/*Find the closest point on the edges of both segments*/
1016  p1.x=s1p1->x+s1k*(s1p2->x-s1p1->x);
1017  p1.y=s1p1->y+s1k*(s1p2->y-s1p1->y);
1018  p1.z=s1p1->z+s1k*(s1p2->z-s1p1->z);
1019 
1020  p2.x=s2p1->x+s2k*(s2p2->x-s2p1->x);
1021  p2.y=s2p1->y+s2k*(s2p2->y-s2p1->y);
1022  p2.z=s2p1->z+s2k*(s2p2->z-s2p1->z);
1023 
1024  if (!lw_dist3d_pt_pt(&p1,&p2,dl))/* Send the closest points to point-point calculation*/
1025  {
1026  return LW_FALSE;
1027  }
1028  }
1029  return LW_TRUE;
1030 }
1031 
1039 int
1041 {
1042  uint32_t i;
1043 
1044  LWDEBUG(2, "lw_dist3d_point_poly called");
1045 
1046 
1047  if(pt_in_ring_3d(projp, poly->rings[0], plane))
1048  {
1049  for (i=1; i<poly->nrings; i++)
1050  {
1051  /* Inside a hole. Distance = pt -> ring */
1052  if ( pt_in_ring_3d(projp, poly->rings[i], plane ))
1053  {
1054  LWDEBUG(3, " inside an hole");
1055  return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);
1056  }
1057  }
1058 
1059  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*/
1060  }
1061  else
1062  {
1063  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*/
1064  }
1065 
1066  return LW_TRUE;
1067 
1068 }
1069 
1075 {
1076  uint32_t i,j,k;
1077  double f, s1, s2;
1078  VECTOR3D projp1_projp2;
1079  POINT3DZ p1, p2,projp1, projp2, intersectionp;
1080 
1081  getPoint3dz_p(pa, 0, &p1);
1082 
1084  &p1, plane, &projp1); /*the sign of s1 tells us on which side of the plane the point is. */
1085  lw_dist3d_pt_poly(&p1, poly, plane,&projp1, dl);
1086  if ((s1 == 0.0) && (dl->distance < dl->tolerance))
1087  {
1088  return LW_TRUE;
1089  }
1090 
1091  for (i=1;i<pa->npoints;i++)
1092  {
1093  int intersects;
1094  getPoint3dz_p(pa, i, &p2);
1095  s2=project_point_on_plane(&p2, plane, &projp2);
1096  lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl);
1097  if ((s2 == 0.0) && (dl->distance < dl->tolerance))
1098  {
1099  return LW_TRUE;
1100  }
1101 
1102  /*If s1and s2 has different signs that means they are on different sides of the plane of the polygon.
1103  That means that the edge between the points crosses the plane and might intersect with the polygon*/
1104  if ((s1 * s2) < 0)
1105  {
1106  f=fabs(s1)/(fabs(s1)+fabs(s2)); /*The size of s1 and s2 is the distance from the point to the plane.*/
1107  get_3dvector_from_points(&projp1, &projp2,&projp1_projp2);
1108 
1109  /*get the point where the line segment crosses the plane*/
1110  intersectionp.x=projp1.x+f*projp1_projp2.x;
1111  intersectionp.y=projp1.y+f*projp1_projp2.y;
1112  intersectionp.z=projp1.z+f*projp1_projp2.z;
1113 
1114  intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/
1115 
1116  if(pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/
1117  {
1118  for (k=1;k<poly->nrings; k++)
1119  {
1120  /* Inside a hole, so no intersection with the polygon*/
1121  if ( pt_in_ring_3d(&intersectionp, poly->rings[k], plane ))
1122  {
1124  break;
1125  }
1126  }
1127  if(intersects)
1128  {
1129  dl->distance=0.0;
1130  dl->p1.x=intersectionp.x;
1131  dl->p1.y=intersectionp.y;
1132  dl->p1.z=intersectionp.z;
1133 
1134  dl->p2.x=intersectionp.x;
1135  dl->p2.y=intersectionp.y;
1136  dl->p2.z=intersectionp.z;
1137  return LW_TRUE;
1138 
1139  }
1140  }
1141  }
1142 
1143  projp1=projp2;
1144  s1=s2;
1145  p1=p2;
1146  }
1147 
1148  /*check or pointarray against boundary and inner boundaries of the polygon*/
1149  for (j=0;j<poly->nrings;j++)
1150  {
1151  lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);
1152  }
1153 
1154 return LW_TRUE;
1155 }
1156 
1157 /* Here we define the plane of a polygon (boundary pointarray of a polygon)
1158  * the plane is stored as a point in plane (plane.pop) and a normal vector (plane.pv)
1159  *
1160  * We are calculating the average point and using it to break the polygon into
1161  * POL_BREAKS (3) parts. Then we calculate the normal of those 3 vectors and
1162  * use its average to define the normal of the plane.This is done to minimize
1163  * the error introduced by FP precision
1164  * We use [3] breaks to contemplate the special case of 3d triangles
1165  */
1166 int
1168 {
1169  const uint32_t POL_BREAKS = 3;
1170  uint32_t unique_points = pa->npoints - 1;
1171  uint32_t i;
1172 
1173  /* Calculate the average point */
1174  pl->pop.x = 0.0;
1175  pl->pop.y = 0.0;
1176  pl->pop.z = 0.0;
1177  for (i = 0; i < unique_points; i++)
1178  {
1179  POINT3DZ p;
1180  getPoint3dz_p(pa, i, &p);
1181  pl->pop.x += p.x;
1182  pl->pop.y += p.y;
1183  pl->pop.z += p.z;
1184  }
1185 
1186  pl->pop.x /= unique_points;
1187  pl->pop.y /= unique_points;
1188  pl->pop.z /= unique_points;
1189 
1190  pl->pv.x = 0.0;
1191  pl->pv.y = 0.0;
1192  pl->pv.z = 0.0;
1193  for (i = 0; i < POL_BREAKS; i++)
1194  {
1195  POINT3DZ point1, point2;
1196  VECTOR3D v1, v2, vp;
1197  uint32_t n1, n2;
1198  n1 = i * unique_points / POL_BREAKS;
1199  n2 = n1 + unique_points / POL_BREAKS; /* May overflow back to the first / last point */
1200 
1201  if (n1 == n2)
1202  continue;
1203 
1204  getPoint3dz_p(pa, n1, &point1);
1205  if (!get_3dvector_from_points(&pl->pop, &point1, &v1))
1206  continue;
1207 
1208  getPoint3dz_p(pa, n2, &point2);
1209  if (!get_3dvector_from_points(&pl->pop, &point2, &v2))
1210  continue;
1211 
1212  if (get_3dcross_product(&v1, &v2, &vp))
1213  {
1214  /* If the cross product isn't zero these 3 points aren't collinear
1215  * We divide by its lengthsq to normalize the additions */
1216  double vl = vp.x * vp.x + vp.y * vp.y + vp.z * vp.z;
1217  pl->pv.x += vp.x / vl;
1218  pl->pv.y += vp.y / vl;
1219  pl->pv.z += vp.z / vl;
1220  }
1221  }
1222 
1223  return (!FP_IS_ZERO(pl->pv.x) || !FP_IS_ZERO(pl->pv.y) || !FP_IS_ZERO(pl->pv.z));
1224 }
1225 
1230 double
1232 {
1233 /*In our plane definition we have a point on the plane and a normal vector (pl.pv), perpendicular to the plane
1234 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.
1235 So, we already have a direction from p to find p0, but we don't know the distance.
1236 */
1237 
1238  VECTOR3D v1;
1239  double f;
1240 
1241  if (!get_3dvector_from_points(&(pl->pop), p, &v1))
1242  return 0.0;
1243 
1244  f = DOT(pl->pv, v1);
1245  if (FP_IS_ZERO(f))
1246  {
1247  /* Point is in the plane */
1248  *p0 = *p;
1249  return 0;
1250  }
1251 
1252  f = -f / DOT(pl->pv, pl->pv);
1253 
1254  p0->x = p->x + pl->pv.x * f;
1255  p0->y = p->y + pl->pv.y * f;
1256  p0->z = p->z + pl->pv.z * f;
1257 
1258  return f;
1259 }
1260 
1261 
1262 
1263 
1276 int
1277 pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring,PLANE3D *plane)
1278 {
1279 
1280  uint32_t cn = 0; /* the crossing number counter */
1281  uint32_t i;
1282  POINT3DZ v1, v2;
1283 
1284  POINT3DZ first, last;
1285 
1286  getPoint3dz_p(ring, 0, &first);
1287  getPoint3dz_p(ring, ring->npoints-1, &last);
1288  if ( memcmp(&first, &last, sizeof(POINT3DZ)) )
1289  {
1290  lwerror("pt_in_ring_3d: V[n] != V[0] (%g %g %g!= %g %g %g)",
1291  first.x, first.y, first.z, last.x, last.y, last.z);
1292  return LW_FALSE;
1293  }
1294 
1295  LWDEBUGF(2, "pt_in_ring_3d called with point: %g %g %g", p->x, p->y, p->z);
1296  /* printPA(ring); */
1297 
1298  /* loop through all edges of the polygon */
1299  getPoint3dz_p(ring, 0, &v1);
1300 
1301 
1302  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*/
1303  {
1304  for (i=0; i<ring->npoints-1; i++)
1305  {
1306  double vt;
1307  getPoint3dz_p(ring, i+1, &v2);
1308 
1309  /* edge from vertex i to vertex i+1 */
1310  if
1311  (
1312  /* an upward crossing */
1313  ((v1.y <= p->y) && (v2.y > p->y))
1314  /* a downward crossing */
1315  || ((v1.y > p->y) && (v2.y <= p->y))
1316  )
1317  {
1318 
1319  vt = (double)(p->y - v1.y) / (v2.y - v1.y);
1320 
1321  /* P.x <intersect */
1322  if (p->x < v1.x + vt * (v2.x - v1.x))
1323  {
1324  /* a valid crossing of y=p.y right of p.x */
1325  ++cn;
1326  }
1327  }
1328  v1 = v2;
1329  }
1330  }
1331  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*/
1332  {
1333  for (i=0; i<ring->npoints-1; i++)
1334  {
1335  double vt;
1336  getPoint3dz_p(ring, i+1, &v2);
1337 
1338  /* edge from vertex i to vertex i+1 */
1339  if
1340  (
1341  /* an upward crossing */
1342  ((v1.z <= p->z) && (v2.z > p->z))
1343  /* a downward crossing */
1344  || ((v1.z > p->z) && (v2.z <= p->z))
1345  )
1346  {
1347 
1348  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1349 
1350  /* P.x <intersect */
1351  if (p->x < v1.x + vt * (v2.x - v1.x))
1352  {
1353  /* a valid crossing of y=p.y right of p.x */
1354  ++cn;
1355  }
1356  }
1357  v1 = v2;
1358  }
1359  }
1360  else /*Hopefully we only have the cases where x part of the normal vector is largest left*/
1361  {
1362  for (i=0; i<ring->npoints-1; i++)
1363  {
1364  double vt;
1365  getPoint3dz_p(ring, i+1, &v2);
1366 
1367  /* edge from vertex i to vertex i+1 */
1368  if
1369  (
1370  /* an upward crossing */
1371  ((v1.z <= p->z) && (v2.z > p->z))
1372  /* a downward crossing */
1373  || ((v1.z > p->z) && (v2.z <= p->z))
1374  )
1375  {
1376 
1377  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1378 
1379  /* P.x <intersect */
1380  if (p->y < v1.y + vt * (v2.y - v1.y))
1381  {
1382  /* a valid crossing of y=p.y right of p.x */
1383  ++cn;
1384  }
1385  }
1386  v1 = v2;
1387  }
1388  }
1389  LWDEBUGF(3, "pt_in_ring_3d returning %d", cn&1);
1390 
1391  return (cn&1); /* 0 if even (out), and 1 if odd (in) */
1392 }
1393 
1394 
1395 
1396 /*------------------------------------------------------------------------------------------------------------
1397 End of Brute force functions
1398 --------------------------------------------------------------------------------------------------------------*/
1399 
1400 
1401 
char * r
Definition: cu_in_wkt.c:24
#define LW_FALSE
Definition: liblwgeom.h:77
#define COLLECTIONTYPE
Definition: liblwgeom.h:91
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
#define LW_FAILURE
Definition: liblwgeom.h:79
#define LINETYPE
Definition: liblwgeom.h:86
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:930
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:85
int getPoint3dz_p(const POINTARRAY *pa, uint32_t n, POINT3DZ *point)
Definition: lwgeom_api.c:215
void lwfree(void *mem)
Definition: lwutil.c:244
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:213
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1085
#define POLYGONTYPE
Definition: liblwgeom.h:87
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:218
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
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
LWLINE * lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points)
Definition: lwline.c:237
LWPOINT * lwpoint_make3dz(int srid, double x, double y, double z)
Definition: lwpoint.c:173
double lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling max distance calculations and dfullywithin calculations.
Definition: measures.c:181
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:224
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
#define FP_IS_ZERO(A)
Datum intersects(PG_FUNCTION_ARGS)
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
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,...
Definition: measures3d.c:1277
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:737
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_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_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)
point to point calculation
Definition: measures3d.c:581
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:918
LWGEOM * lwgeom_furthest_line_3d(LWGEOM *lw1, LWGEOM *lw2)
Definition: measures3d.c:94
LWGEOM * lwgeom_closest_point_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:100
int define_plane(POINTARRAY *pa, PLANE3D *pl)
Definition: measures3d.c:1167
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
static int get_3dvector_from_points(POINT3DZ *p1, POINT3DZ *p2, VECTOR3D *v)
Definition: measures3d.c:45
static int get_3dcross_product(VECTOR3D *v1, VECTOR3D *v2, VECTOR3D *v)
Definition: measures3d.c:55
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:1040
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:815
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:832
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:767
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
int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
polygon to polygon calculation
Definition: measures3d.c:688
int lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
line to line calculation
Definition: measures3d.c:652
int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, PLANE3D *plane, DISTPTS3D *dl)
Computes pointarray to polygon distance.
Definition: measures3d.c:1074
LWGEOM * lw_dist3d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing 3dshortestline and 3dlongestline calculations.
Definition: measures3d.c:110
double lwgeom_maxdistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 3d max distance calculation.
Definition: measures3d.c:314
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
int lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)
point to line calculation
Definition: measures3d.c:597
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:1231
int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
line to polygon calculation
Definition: measures3d.c:665
LWGEOM * lw_dist3d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing 3dclosestpoint calculations.
Definition: measures3d.c:214
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
LWGEOM * lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:88
int lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS3D *dl)
Finds all combinations of segments between two pointarrays.
Definition: measures3d.c:864
double lwgeom_mindistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 3d min distance calculation.
Definition: measures3d.c:352
#define DOT(u, v)
Definition: measures3d.h:31
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
LWGEOM * lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing shortestline and longestline calculations.
Definition: measures.c:84
LWGEOM * lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing closestpoint calculations.
Definition: measures.c:131
#define DIST_MIN
Definition: measures.h:44
#define DIST_MAX
Definition: measures.h:43
POINT3DZ p2
Definition: measures3d.h:43
int twisted
Definition: measures3d.h:45
POINT3DZ p1
Definition: measures3d.h:42
double distance
Definition: measures3d.h:41
int mode
Definition: measures3d.h:44
double tolerance
Definition: measures3d.h:46
Structure used in distance-calculations.
Definition: measures3d.h:40
POINT2D p1
Definition: measures.h:52
POINT2D p2
Definition: measures.h:53
double tolerance
Definition: measures.h:56
int mode
Definition: measures.h:54
double distance
Definition: measures.h:51
Structure used in distance-calculations.
Definition: measures.h:50
double zmax
Definition: liblwgeom.h:300
double zmin
Definition: liblwgeom.h:299
uint32_t ngeoms
Definition: liblwgeom.h:510
LWGEOM ** geoms
Definition: liblwgeom.h:512
uint8_t type
Definition: liblwgeom.h:399
int32_t srid
Definition: liblwgeom.h:402
POINTARRAY * points
Definition: liblwgeom.h:425
POINTARRAY * point
Definition: liblwgeom.h:414
POINTARRAY ** rings
Definition: liblwgeom.h:460
uint32_t nrings
Definition: liblwgeom.h:458
POINT3DZ pop
Definition: measures3d.h:57
VECTOR3D pv
Definition: measures3d.h:58
double y
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:331
double z
Definition: liblwgeom.h:337
double x
Definition: liblwgeom.h:337
double y
Definition: liblwgeom.h:337
double z
Definition: liblwgeom.h:343
double x
Definition: liblwgeom.h:343
double y
Definition: liblwgeom.h:343
uint32_t npoints
Definition: liblwgeom.h:374
double z
Definition: measures3d.h:51
double x
Definition: measures3d.h:51
double y
Definition: measures3d.h:51
unsigned int uint32_t
Definition: uthash.h:78