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