PostGIS  2.1.10dev-r@@SVN_REVISION@@
measures3d.c
Go to the documentation of this file.
1 
2 /**********************************************************************
3  * $Id: measures.c 5439 2010-03-16 03:13:33Z pramsey $
4  *
5  * PostGIS - Spatial Types for PostgreSQL
6  * http://postgis.net
7  * Copyright 2011 Nicklas Avén
8  *
9  * This is free software; you can redistribute and/or modify it under
10  * the terms of the GNU General Public Licence. See the COPYING file.
11  *
12  **********************************************************************/
13 
14 #include <string.h>
15 #include <stdlib.h>
16 
17 #include "measures3d.h"
18 #include "lwgeom_log.h"
19 
20 
21 /*------------------------------------------------------------------------------------------------------------
22 Initializing functions
23 The functions starting the distance-calculation processses
24 --------------------------------------------------------------------------------------------------------------*/
25 
29 LWGEOM *
30 lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
31 {
32  double x1,x2,y1,y2, z1, z2;
33  double initdistance = ( mode == DIST_MIN ? MAXFLOAT : -1.0);
34  DISTPTS3D thedl;
35  LWPOINT *lwpoints[2];
36  LWGEOM *result;
37 
38  thedl.mode = mode;
39  thedl.distance = initdistance;
40  thedl.tolerance = 0.0;
41 
42  LWDEBUG(2, "lw_dist3d_distanceline is called");
43  if (!lw_dist3d_recursive(lw1, lw2, &thedl))
44  {
45  /*should never get here. all cases ought to be error handled earlier*/
46  lwerror("Some unspecified error.");
47  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
48  }
49 
50  /*if thedl.distance is unchanged there where only empty geometries input*/
51  if (thedl.distance == initdistance)
52  {
53  LWDEBUG(3, "didn't find geometries to measure between, returning null");
54  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
55  }
56  else
57  {
58  x1=thedl.p1.x;
59  y1=thedl.p1.y;
60  z1=thedl.p1.z;
61  x2=thedl.p2.x;
62  y2=thedl.p2.y;
63  z2=thedl.p2.z;
64 
65 
66  lwpoints[0] = lwpoint_make3dz(srid, x1, y1, z1);
67  lwpoints[1] = lwpoint_make3dz(srid, x2, y2, z2);
68 
69  result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
70  }
71 
72  return result;
73 }
74 
78 LWGEOM *
79 lw_dist3d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
80 {
81  double x,y,z;
82  DISTPTS3D thedl;
83  double initdistance = MAXFLOAT;
84  LWGEOM *result;
85 
86  thedl.mode = mode;
87  thedl.distance= initdistance;
88  thedl.tolerance = 0;
89 
90  LWDEBUG(2, "lw_dist3d_distancepoint is called");
91 
92  if (!lw_dist3d_recursive(lw1, lw2, &thedl))
93  {
94  /*should never get here. all cases ought to be error handled earlier*/
95  lwerror("Some unspecified error.");
96  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
97  }
98 
99  if (thedl.distance == initdistance)
100  {
101  LWDEBUG(3, "didn't find geometries to measure between, returning null");
102  result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
103  }
104  else
105  {
106  x=thedl.p1.x;
107  y=thedl.p1.y;
108  z=thedl.p1.z;
109  result = (LWGEOM *)lwpoint_make3dz(srid, x, y, z);
110  }
111 
112  return result;
113 }
114 
115 
119 double
121 {
122  LWDEBUG(2, "lwgeom_maxdistance3d is called");
123 
124  return lwgeom_maxdistance3d_tolerance( lw1, lw2, 0.0 );
125 }
126 
131 double
132 lwgeom_maxdistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
133 {
134  /*double thedist;*/
135  DISTPTS3D thedl;
136  LWDEBUG(2, "lwgeom_maxdistance3d_tolerance is called");
137  thedl.mode = DIST_MAX;
138  thedl.distance= -1;
139  thedl.tolerance = tolerance;
140  if (lw_dist3d_recursive(lw1, lw2, &thedl))
141  {
142  return thedl.distance;
143  }
144  /*should never get here. all cases ought to be error handled earlier*/
145  lwerror("Some unspecified error.");
146  return -1;
147 }
148 
152 double
154 {
155  LWDEBUG(2, "lwgeom_mindistance3d is called");
156  return lwgeom_mindistance3d_tolerance( lw1, lw2, 0.0 );
157 }
158 
163 double
164 lwgeom_mindistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
165 {
166  DISTPTS3D thedl;
167  LWDEBUG(2, "lwgeom_mindistance3d_tolerance is called");
168  thedl.mode = DIST_MIN;
169  thedl.distance= MAXFLOAT;
170  thedl.tolerance = tolerance;
171  if (lw_dist3d_recursive(lw1, lw2, &thedl))
172  {
173  return thedl.distance;
174  }
175  /*should never get here. all cases ought to be error handled earlier*/
176  lwerror("Some unspecified error.");
177  return MAXFLOAT;
178 }
179 
180 
181 /*------------------------------------------------------------------------------------------------------------
182 End of Initializing functions
183 --------------------------------------------------------------------------------------------------------------*/
184 
185 /*------------------------------------------------------------------------------------------------------------
186 Preprocessing functions
187 Functions preparing geometries for distance-calculations
188 --------------------------------------------------------------------------------------------------------------*/
189 
190 
194 int lw_dist3d_recursive(const LWGEOM *lwg1,const LWGEOM *lwg2, DISTPTS3D *dl)
195 {
196  int i, j;
197  int n1=1;
198  int n2=1;
199  LWGEOM *g1 = NULL;
200  LWGEOM *g2 = NULL;
201  LWCOLLECTION *c1 = NULL;
202  LWCOLLECTION *c2 = NULL;
203 
204  LWDEBUGF(2, "lw_dist3d_recursive is called with type1=%d, type2=%d", lwg1->type, lwg2->type);
205 
206  if (lwgeom_is_collection(lwg1))
207  {
208  LWDEBUG(3, "First geometry is collection");
209  c1 = lwgeom_as_lwcollection(lwg1);
210  n1 = c1->ngeoms;
211  }
212  if (lwgeom_is_collection(lwg2))
213  {
214  LWDEBUG(3, "Second geometry is collection");
215  c2 = lwgeom_as_lwcollection(lwg2);
216  n2 = c2->ngeoms;
217  }
218 
219  for ( i = 0; i < n1; i++ )
220  {
221 
222  if (lwgeom_is_collection(lwg1))
223  {
224  g1 = c1->geoms[i];
225  }
226  else
227  {
228  g1 = (LWGEOM*)lwg1;
229  }
230 
231  if (lwgeom_is_empty(g1)) return LW_TRUE;
232 
233  if (lwgeom_is_collection(g1))
234  {
235  LWDEBUG(3, "Found collection inside first geometry collection, recursing");
236  if (!lw_dist3d_recursive(g1, lwg2, dl)) return LW_FALSE;
237  continue;
238  }
239  for ( j = 0; j < n2; j++ )
240  {
241  if (lwgeom_is_collection(lwg2))
242  {
243  g2 = c2->geoms[j];
244  }
245  else
246  {
247  g2 = (LWGEOM*)lwg2;
248  }
249  if (lwgeom_is_collection(g2))
250  {
251  LWDEBUG(3, "Found collection inside second geometry collection, recursing");
252  if (!lw_dist3d_recursive(g1, g2, dl)) return LW_FALSE;
253  continue;
254  }
255 
256 
257  /*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/
258  if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
259 
260 
261  if (!lw_dist3d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
262  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
263  }
264  }
265  return LW_TRUE;
266 }
267 
268 
269 
274 int
276 {
277 
278  int t1 = lwg1->type;
279  int t2 = lwg2->type;
280 
281  LWDEBUGF(2, "lw_dist3d_distribute_bruteforce is called with typ1=%d, type2=%d", lwg1->type, lwg2->type);
282 
283  if ( t1 == POINTTYPE )
284  {
285  if ( t2 == POINTTYPE )
286  {
287  dl->twisted=1;
288  return lw_dist3d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
289  }
290  else if ( t2 == LINETYPE )
291  {
292  dl->twisted=1;
293  return lw_dist3d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
294  }
295  else if ( t2 == POLYGONTYPE )
296  {
297  dl->twisted=1;
298  return lw_dist3d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2,dl);
299  }
300  else
301  {
302  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
303  return LW_FALSE;
304  }
305  }
306  else if ( t1 == LINETYPE )
307  {
308  if ( t2 == POINTTYPE )
309  {
310  dl->twisted=(-1);
311  return lw_dist3d_point_line((LWPOINT *)lwg2,(LWLINE *)lwg1,dl);
312  }
313  else if ( t2 == LINETYPE )
314  {
315  dl->twisted=1;
316  return lw_dist3d_line_line((LWLINE *)lwg1,(LWLINE *)lwg2,dl);
317  }
318  else if ( t2 == POLYGONTYPE )
319  {
320  dl->twisted=1;
321  return lw_dist3d_line_poly((LWLINE *)lwg1,(LWPOLY *)lwg2,dl);
322  }
323  else
324  {
325  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
326  return LW_FALSE;
327  }
328  }
329  else if ( t1 == POLYGONTYPE )
330  {
331  if ( t2 == POLYGONTYPE )
332  {
333  dl->twisted=1;
334  return lw_dist3d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2,dl);
335  }
336  else if ( t2 == POINTTYPE )
337  {
338  dl->twisted=-1;
339  return lw_dist3d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1,dl);
340  }
341  else if ( t2 == LINETYPE )
342  {
343  dl->twisted=-1;
344  return lw_dist3d_line_poly((LWLINE *)lwg2,(LWPOLY *)lwg1,dl);
345  }
346  else
347  {
348  lwerror("Unsupported geometry type: %s", lwtype_name(t2));
349  return LW_FALSE;
350  }
351  }
352  else
353  {
354  lwerror("Unsupported geometry type: %s", lwtype_name(t1));
355  return LW_FALSE;
356  }
357  /*You shouldn't being able to get here*/
358  lwerror("unspecified error in function lw_dist3d_distribute_bruteforce");
359  return LW_FALSE;
360 }
361 
362 
363 
364 /*------------------------------------------------------------------------------------------------------------
365 End of Preprocessing functions
366 --------------------------------------------------------------------------------------------------------------*/
367 
368 
369 /*------------------------------------------------------------------------------------------------------------
370 Brute force functions
371 So far the only way to do 3D-calculations
372 --------------------------------------------------------------------------------------------------------------*/
373 
378 int
380 {
381  POINT3DZ p1;
382  POINT3DZ p2;
383  LWDEBUG(2, "lw_dist3d_point_point is called");
384 
385  getPoint3dz_p(point1->point, 0, &p1);
386  getPoint3dz_p(point2->point, 0, &p2);
387 
388  return lw_dist3d_pt_pt(&p1, &p2,dl);
389 }
394 int
396 {
397  POINT3DZ p;
398  POINTARRAY *pa = line->points;
399  LWDEBUG(2, "lw_dist3d_point_line is called");
400 
401  getPoint3dz_p(point->point, 0, &p);
402  return lw_dist3d_pt_ptarray(&p, pa, dl);
403 }
404 
416 int
418 {
419  POINT3DZ p, projp;/*projp is "point projected on plane"*/
420  PLANE3D plane;
421  LWDEBUG(2, "lw_dist3d_point_poly is called");
422  getPoint3dz_p(point->point, 0, &p);
423 
424  /*If we are lookig for max distance, longestline or dfullywithin*/
425  if (dl->mode == DIST_MAX)
426  {
427  LWDEBUG(3, "looking for maxdistance");
428  return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl);
429  }
430 
431  /*Find the plane of the polygon, the "holes" have to be on the same plane. so we only care about the boudary*/
432  if(!define_plane(poly->rings[0], &plane))
433  return LW_FALSE;
434 
435  /*get our point projected on the plane of the polygon*/
436  project_point_on_plane(&p, &plane, &projp);
437 
438  return lw_dist3d_pt_poly(&p, poly,&plane, &projp, dl);
439 }
440 
441 
446 int
448 {
449  POINTARRAY *pa1 = line1->points;
450  POINTARRAY *pa2 = line2->points;
451  LWDEBUG(2, "lw_dist3d_line_line is called");
452 
453  return lw_dist3d_ptarray_ptarray(pa1, pa2, dl);
454 }
455 
461 {
462  PLANE3D plane;
463  LWDEBUG(2, "lw_dist3d_line_poly is called");
464 
465  if (dl->mode == DIST_MAX)
466  {
467  return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);
468  }
469 
470  if(!define_plane(poly->rings[0], &plane))
471  return LW_FALSE;
472 
473  return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl);
474 }
475 
480 int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
481 {
482  PLANE3D plane;
483  LWDEBUG(2, "lw_dist3d_poly_poly is called");
484  if (dl->mode == DIST_MAX)
485  {
486  return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
487  }
488 
489  if(!define_plane(poly2->rings[0], &plane))
490  return LW_FALSE;
491 
492  /*What we do here is to compare the bondary of one polygon with the other polygon
493  and then take the second boudary comparing with the first polygon*/
494  dl->twisted=1;
495  if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl))
496  return LW_FALSE;
497  if(dl->distance==0.0) /*Just check if the answer already is given*/
498  return LW_TRUE;
499 
500  if(!define_plane(poly1->rings[0], &plane))
501  return LW_FALSE;
502  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.*/
503  return lw_dist3d_ptarray_poly(poly2->rings[0], poly1,&plane, dl);
504 }
505 
511 int
513 {
514  int t;
515  POINT3DZ start, end;
516  int twist = dl->twisted;
517 
518  LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
519 
520  getPoint3dz_p(pa, 0, &start);
521 
522  for (t=1; t<pa->npoints; t++)
523  {
524  dl->twisted=twist;
525  getPoint3dz_p(pa, t, &end);
526  if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
527 
528  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
529  start = end;
530  }
531 
532  return LW_TRUE;
533 }
534 
535 
541 int
543 {
544  POINT3DZ c;
545  double r;
546  /*if start==end, then use pt distance */
547  if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )
548  {
549  return lw_dist3d_pt_pt(p,A,dl);
550  }
551 
552 
553  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) );
554 
555  /*This is for finding the 3Dmaxdistance.
556  the maxdistance have to be between two vertexes,
557  compared to mindistance which can be between
558  tvo vertexes vertex.*/
559  if (dl->mode == DIST_MAX)
560  {
561  if (r>=0.5)
562  {
563  return lw_dist3d_pt_pt(p,A,dl);
564  }
565  if (r<0.5)
566  {
567  return lw_dist3d_pt_pt(p,B,dl);
568  }
569  }
570 
571  if (r<0) /*If the first vertex A is closest to the point p*/
572  {
573  return lw_dist3d_pt_pt(p,A,dl);
574  }
575  if (r>1) /*If the second vertex B is closest to the point p*/
576  {
577  return lw_dist3d_pt_pt(p,B,dl);
578  }
579 
580  /*else if the point p is closer to some point between a and b
581  then we find that point and send it to lw_dist3d_pt_pt*/
582  c.x=A->x + r * (B->x-A->x);
583  c.y=A->y + r * (B->y-A->y);
584  c.z=A->z + r * (B->z-A->z);
585 
586  return lw_dist3d_pt_pt(p,&c,dl);
587 }
588 
589 
597 int
599 {
600  double dx = thep2->x - thep1->x;
601  double dy = thep2->y - thep1->y;
602  double dz = thep2->z - thep1->z;
603  double dist = sqrt ( dx*dx + dy*dy + dz*dz);
604  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 );
605 
606  if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
607  {
608  dl->distance = dist;
609 
610  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*/
611  {
612  dl->p1 = *thep1;
613  dl->p2 = *thep2;
614  }
615  else
616  {
617  dl->p1 = *thep2;
618  dl->p2 = *thep1;
619  }
620  }
621  return LW_TRUE;
622 }
623 
624 
629 int
631 {
632  int t,u;
633  POINT3DZ start, end;
634  POINT3DZ start2, end2;
635  int twist = dl->twisted;
636  LWDEBUGF(2, "lw_dist3d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
637 
638 
639 
640  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*/
641  {
642  for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
643  {
644  getPoint3dz_p(l1, t, &start);
645  for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
646  {
647  getPoint3dz_p(l2, u, &start2);
648  lw_dist3d_pt_pt(&start,&start2,dl);
649  LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
650  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
651  t, u, dl->distance, dl->tolerance);
652  }
653  }
654  }
655  else
656  {
657  getPoint3dz_p(l1, 0, &start);
658  for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
659  {
660  getPoint3dz_p(l1, t, &end);
661  getPoint3dz_p(l2, 0, &start2);
662  for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
663  {
664  getPoint3dz_p(l2, u, &end2);
665  dl->twisted=twist;
666  lw_dist3d_seg_seg(&start, &end, &start2, &end2,dl);
667  LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
668  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
669  t, u, dl->distance, dl->tolerance);
670  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
671  start2 = end2;
672  }
673  start = end;
674  }
675  }
676  return LW_TRUE;
677 }
678 
683 int
685 {
686  VECTOR3D v1, v2, vl;
687  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*/
688  POINT3DZ p1, p2;
689  double a, b, c, d, e, D;
690 
691  /*s1p1 and s1p2 are the same point */
692  if ( ( s1p1->x == s1p2->x) && (s1p1->y == s1p2->y) && (s1p1->z == s1p2->z) )
693  {
694  return lw_dist3d_pt_seg(s1p1,s2p1,s2p2,dl);
695  }
696  /*s2p1 and s2p2 are the same point */
697  if ( ( s2p1->x == s2p2->x) && (s2p1->y == s2p2->y) && (s2p1->z == s2p2->z) )
698  {
699  dl->twisted= ((dl->twisted) * (-1));
700  return lw_dist3d_pt_seg(s2p1,s1p1,s1p2,dl);
701  }
702 
703 /*
704  Here we use algorithm from softsurfer.com
705  that can be found here
706  http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
707 */
708 
709  if (!get_3dvector_from_points(s1p1, s1p2, &v1))
710  return LW_FALSE;
711 
712  if (!get_3dvector_from_points(s2p1, s2p2, &v2))
713  return LW_FALSE;
714 
715  if (!get_3dvector_from_points(s2p1, s1p1, &vl))
716  return LW_FALSE;
717 
718  a = DOT(v1,v1);
719  b = DOT(v1,v2);
720  c = DOT(v2,v2);
721  d = DOT(v1,vl);
722  e = DOT(v2,vl);
723  D = a*c - b*b;
724 
725 
726  if (D <0.000000001)
727  { /* the lines are almost parallel*/
728  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.*/
729  if(b>c) /* use the largest denominator*/
730  {
731  s2k=d/b;
732  }
733  else
734  {
735  s2k =e/c;
736  }
737  }
738  else
739  {
740  s1k = (b*e - c*d) / D;
741  s2k = (a*e - b*d) / D;
742  }
743 
744  /* 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*/
745  if(s1k<0.0||s1k>1.0||s2k<0.0||s2k>1.0)
746  {
747  if(s1k<0.0)
748  {
749 
750  if (!lw_dist3d_pt_seg(s1p1, s2p1, s2p2, dl))
751  {
752  return LW_FALSE;
753  }
754  }
755  if(s1k>1.0)
756  {
757 
758  if (!lw_dist3d_pt_seg(s1p2, s2p1, s2p2, dl))
759  {
760  return LW_FALSE;
761  }
762  }
763  if(s2k<0.0)
764  {
765  dl->twisted= ((dl->twisted) * (-1));
766  if (!lw_dist3d_pt_seg(s2p1, s1p1, s1p2, dl))
767  {
768  return LW_FALSE;
769  }
770  }
771  if(s2k>1.0)
772  {
773  dl->twisted= ((dl->twisted) * (-1));
774  if (!lw_dist3d_pt_seg(s2p2, s1p1, s1p2, dl))
775  {
776  return LW_FALSE;
777  }
778  }
779  }
780  else
781  {/*Find the closest point on the edges of both segments*/
782  p1.x=s1p1->x+s1k*(s1p2->x-s1p1->x);
783  p1.y=s1p1->y+s1k*(s1p2->y-s1p1->y);
784  p1.z=s1p1->z+s1k*(s1p2->z-s1p1->z);
785 
786  p2.x=s2p1->x+s2k*(s2p2->x-s2p1->x);
787  p2.y=s2p1->y+s2k*(s2p2->y-s2p1->y);
788  p2.z=s2p1->z+s2k*(s2p2->z-s2p1->z);
789 
790  if (!lw_dist3d_pt_pt(&p1,&p2,dl))/* Send the closest points to point-point calculation*/
791  {
792  return LW_FALSE;
793  }
794  }
795  return LW_TRUE;
796 }
797 
805 int
807 {
808  int i;
809 
810  LWDEBUG(2, "lw_dist3d_point_poly called");
811 
812 
813  if(pt_in_ring_3d(projp, poly->rings[0], plane))
814  {
815  for (i=1; i<poly->nrings; i++)
816  {
817  /* Inside a hole. Distance = pt -> ring */
818  if ( pt_in_ring_3d(projp, poly->rings[i], plane ))
819  {
820  LWDEBUG(3, " inside an hole");
821  return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);
822  }
823  }
824 
825  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*/
826  }
827  else
828  {
829  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*/
830  }
831 
832  return LW_TRUE;
833 
834 }
835 
841 {
842 
843 
844  int i,j,k;
845  double f, s1, s2;
846  VECTOR3D projp1_projp2;
847  POINT3DZ p1, p2,projp1, projp2, intersectionp;
848 
849  getPoint3dz_p(pa, 0, &p1);
850 
851  s1=project_point_on_plane(&p1, plane, &projp1); /*the sign of s1 tells us on which side of the plane the point is. */
852  lw_dist3d_pt_poly(&p1, poly, plane,&projp1, dl);
853 
854  for (i=1;i<pa->npoints;i++)
855  {
856  int intersects;
857  getPoint3dz_p(pa, i, &p2);
858  s2=project_point_on_plane(&p2, plane, &projp2);
859  lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl);
860 
861  /*If s1and s2 has different signs that means they are on different sides of the plane of the polygon.
862  That means that the edge between the points crosses the plane and might intersect with the polygon*/
863  if((s1*s2)<=0)
864  {
865  f=fabs(s1)/(fabs(s1)+fabs(s2)); /*The size of s1 and s2 is the distance from the point to the plane.*/
866  get_3dvector_from_points(&projp1, &projp2,&projp1_projp2);
867 
868  /*get the point where the line segment crosses the plane*/
869  intersectionp.x=projp1.x+f*projp1_projp2.x;
870  intersectionp.y=projp1.y+f*projp1_projp2.y;
871  intersectionp.z=projp1.z+f*projp1_projp2.z;
872 
873  intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/
874 
875  if(pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/
876  {
877  for (k=1;k<poly->nrings; k++)
878  {
879  /* Inside a hole, so no intersection with the polygon*/
880  if ( pt_in_ring_3d(&intersectionp, poly->rings[k], plane ))
881  {
882  intersects=LW_FALSE;
883  break;
884  }
885  }
886  if(intersects)
887  {
888  dl->distance=0.0;
889  dl->p1.x=intersectionp.x;
890  dl->p1.y=intersectionp.y;
891  dl->p1.z=intersectionp.z;
892 
893  dl->p2.x=intersectionp.x;
894  dl->p2.y=intersectionp.y;
895  dl->p2.z=intersectionp.z;
896  return LW_TRUE;
897 
898  }
899  }
900  }
901 
902  projp1=projp2;
903  s1=s2;
904  p1=p2;
905  }
906 
907  /*check or pointarray against boundary and inner boundaries of the polygon*/
908  for (j=0;j<poly->nrings;j++)
909  {
910  lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);
911  }
912 
913 return LW_TRUE;
914 }
915 
916 
922 int
924 {
925  int i,j, numberofvectors, pointsinslice;
926  POINT3DZ p, p1, p2;
927 
928  double sumx=0;
929  double sumy=0;
930  double sumz=0;
931  double vl; /*vector length*/
932 
933  VECTOR3D v1, v2, v;
934 
935  if((pa->npoints-1)==3) /*Triangle is special case*/
936  {
937  pointsinslice=1;
938  }
939  else
940  {
941  pointsinslice=(int) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/
942  }
943 
944  /*find the avg point*/
945  for (i=0;i<(pa->npoints-1);i++)
946  {
947  getPoint3dz_p(pa, i, &p);
948  sumx+=p.x;
949  sumy+=p.y;
950  sumz+=p.z;
951  }
952  pl->pop.x=(sumx/(pa->npoints-1));
953  pl->pop.y=(sumy/(pa->npoints-1));
954  pl->pop.z=(sumz/(pa->npoints-1));
955 
956  sumx=0;
957  sumy=0;
958  sumz=0;
959  numberofvectors= floor((pa->npoints-1)/pointsinslice); /*the number of vectors we try can be 3, 4 or 5*/
960 
961  getPoint3dz_p(pa, 0, &p1);
962  for (j=pointsinslice;j<pa->npoints;j+=pointsinslice)
963  {
964  getPoint3dz_p(pa, j, &p2);
965 
966  if (!get_3dvector_from_points(&(pl->pop), &p1, &v1) || !get_3dvector_from_points(&(pl->pop), &p2, &v2))
967  return LW_FALSE;
968  /*perpendicular vector is cross product of v1 and v2*/
969  if (!get_3dcross_product(&v1,&v2, &v))
970  return LW_FALSE;
971  vl=VECTORLENGTH(v);
972  sumx+=(v.x/vl);
973  sumy+=(v.y/vl);
974  sumz+=(v.z/vl);
975  p1=p2;
976  }
977  pl->pv.x=(sumx/numberofvectors);
978  pl->pv.y=(sumy/numberofvectors);
979  pl->pv.z=(sumz/numberofvectors);
980 
981  return 1;
982 }
983 
988 double
990 {
991 /*In our plane definition we have a point on the plane and a normal vektor (pl.pv), perpendicular to the plane
992 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.
993 So, we already have a direction from p to find p0, but we don't know the distance.
994 */
995 
996  VECTOR3D v1;
997  double f;
998 
999  if (!get_3dvector_from_points(&(pl->pop), p, &v1))
1000  return LW_FALSE;
1001 
1002  f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv));
1003 
1004  p0->x=p->x+pl->pv.x*f;
1005  p0->y=p->y+pl->pv.y*f;
1006  p0->z=p->z+pl->pv.z*f;
1007 
1008  return f;
1009 }
1010 
1011 
1012 
1013 
1026 int
1027 pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring,PLANE3D *plane)
1028 {
1029 
1030  int cn = 0; /* the crossing number counter */
1031  int i;
1032  POINT3DZ v1, v2;
1033 
1034  POINT3DZ first, last;
1035 
1036  getPoint3dz_p(ring, 0, &first);
1037  getPoint3dz_p(ring, ring->npoints-1, &last);
1038  if ( memcmp(&first, &last, sizeof(POINT3DZ)) )
1039  {
1040  lwerror("pt_in_ring_3d: V[n] != V[0] (%g %g %g!= %g %g %g)",
1041  first.x, first.y, first.z, last.x, last.y, last.z);
1042  return LW_FALSE;
1043  }
1044 
1045  LWDEBUGF(2, "pt_in_ring_3d called with point: %g %g %g", p->x, p->y, p->z);
1046  /* printPA(ring); */
1047 
1048  /* loop through all edges of the polygon */
1049  getPoint3dz_p(ring, 0, &v1);
1050 
1051 
1052  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*/
1053  {
1054  for (i=0; i<ring->npoints-1; i++)
1055  {
1056  double vt;
1057  getPoint3dz_p(ring, i+1, &v2);
1058 
1059  /* edge from vertex i to vertex i+1 */
1060  if
1061  (
1062  /* an upward crossing */
1063  ((v1.y <= p->y) && (v2.y > p->y))
1064  /* a downward crossing */
1065  || ((v1.y > p->y) && (v2.y <= p->y))
1066  )
1067  {
1068 
1069  vt = (double)(p->y - v1.y) / (v2.y - v1.y);
1070 
1071  /* P.x <intersect */
1072  if (p->x < v1.x + vt * (v2.x - v1.x))
1073  {
1074  /* a valid crossing of y=p.y right of p.x */
1075  ++cn;
1076  }
1077  }
1078  v1 = v2;
1079  }
1080  }
1081  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*/
1082  {
1083  for (i=0; i<ring->npoints-1; i++)
1084  {
1085  double vt;
1086  getPoint3dz_p(ring, i+1, &v2);
1087 
1088  /* edge from vertex i to vertex i+1 */
1089  if
1090  (
1091  /* an upward crossing */
1092  ((v1.z <= p->z) && (v2.z > p->z))
1093  /* a downward crossing */
1094  || ((v1.z > p->z) && (v2.z <= p->z))
1095  )
1096  {
1097 
1098  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1099 
1100  /* P.x <intersect */
1101  if (p->x < v1.x + vt * (v2.x - v1.x))
1102  {
1103  /* a valid crossing of y=p.y right of p.x */
1104  ++cn;
1105  }
1106  }
1107  v1 = v2;
1108  }
1109  }
1110  else /*Hopefully we only have the cases where x part of the normal vector is largest left*/
1111  {
1112  for (i=0; i<ring->npoints-1; i++)
1113  {
1114  double vt;
1115  getPoint3dz_p(ring, i+1, &v2);
1116 
1117  /* edge from vertex i to vertex i+1 */
1118  if
1119  (
1120  /* an upward crossing */
1121  ((v1.z <= p->z) && (v2.z > p->z))
1122  /* a downward crossing */
1123  || ((v1.z > p->z) && (v2.z <= p->z))
1124  )
1125  {
1126 
1127  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1128 
1129  /* P.x <intersect */
1130  if (p->y < v1.y + vt * (v2.y - v1.y))
1131  {
1132  /* a valid crossing of y=p.y right of p.x */
1133  ++cn;
1134  }
1135  }
1136  v1 = v2;
1137  }
1138  }
1139  LWDEBUGF(3, "pt_in_ring_3d returning %d", cn&1);
1140 
1141  return (cn&1); /* 0 if even (out), and 1 if odd (in) */
1142 }
1143 
1144 
1145 
1146 /*------------------------------------------------------------------------------------------------------------
1147 End of Brute force functions
1148 --------------------------------------------------------------------------------------------------------------*/
1149 
1150 
1151 
#define LINETYPE
Definition: liblwgeom.h:61
double z
Definition: liblwgeom.h:290
double y
Definition: liblwgeom.h:290
double distance
Definition: measures3d.h:26
double lwgeom_mindistance3d(LWGEOM *lw1, LWGEOM *lw2)
Function initialazing 3d min distance calculation.
Definition: measures3d.c:153
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:947
int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
line to polygon calculation
Definition: measures3d.c:460
double x
Definition: liblwgeom.h:290
char * r
Definition: cu_in_wkt.c:25
double z
Definition: measures3d.h:36
POINT3DZ p2
Definition: measures3d.h:28
int npoints
Definition: liblwgeom.h:327
#define DIST_MAX
#define POLYGONTYPE
Definition: liblwgeom.h:62
#define DIST_MIN
LWLINE * lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points)
Definition: lwline.c:213
int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
polygon to polygon calculation
Definition: measures3d.c:480
POINT2D * p1
Definition: lwtree.h:12
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
LWGEOM * lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2, int srid, int mode)
Function initializing 3dshortestline and 3dlongestline calculations.
Definition: measures3d.c:30
VECTOR3D pv
Definition: measures3d.h:43
double lwgeom_mindistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
Function handling 3d min distance calculations and dwithin calculations.
Definition: measures3d.c:164
POINTARRAY * point
Definition: liblwgeom.h:367
Structure used in distance-calculations.
Definition: measures3d.h:24
char ** result
Definition: liblwgeom.h:218
#define VECTORLENGTH(v)
Definition: measures3d.h:17
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:684
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:989
POINT3DZ p1
Definition: measures3d.h:27
double y
Definition: measures3d.h:36
POINT2D * p2
Definition: lwtree.h:13
#define DOT(u, v)
Definition: measures3d.h:16
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:164
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:806
double lwgeom_maxdistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
Function handling 3d max distance calculations and dfyllywithin calculations.
Definition: measures3d.c:132
int getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *point)
Definition: lwgeom_api.c:305
#define LW_FALSE
Definition: liblwgeom.h:52
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:417
Datum intersects(PG_FUNCTION_ARGS)
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:51
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:598
LWGEOM ** geoms
Definition: liblwgeom.h:465
int twisted
Definition: measures3d.h:30
int lw_dist3d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS3D *dl)
This function distributes the brut-force for 3D so far the only type, tasks depending on type...
Definition: measures3d.c:275
POINTARRAY ** rings
Definition: liblwgeom.h:413
LWPOINT * lwpoint_make3dz(int srid, double x, double y, double z)
Definition: lwpoint.c:140
int nrings
Definition: liblwgeom.h:411
int mode
Definition: measures3d.h:29
tuple x
Definition: pixval.py:53
POINT3DZ pop
Definition: measures3d.h:42
int lw_dist3d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
This is a recursive function delivering every possible combinatin of subgeometries.
Definition: measures3d.c:194
double lwgeom_maxdistance3d(LWGEOM *lw1, LWGEOM *lw2)
Function initialazing 3d max distance calculation.
Definition: measures3d.c:120
LWGEOM * lw_dist3d_distancepoint(LWGEOM *lw1, LWGEOM *lw2, int srid, int mode)
Function initializing 3dclosestpoint calculations.
Definition: measures3d.c:79
int lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)
point to line calculation
Definition: measures3d.c:395
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:923
double tolerance
Definition: measures3d.h:31
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:1027
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:143
int get_3dvector_from_points(POINT3DZ *p1, POINT3DZ *p2, VECTOR3D *v)
Definition: measures3d.h:85
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
int lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)
point to point calculation
Definition: measures3d.c:379
#define MAXFLOAT
Largest float value.
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:542
double x
Definition: measures3d.h:36
uint8_t type
Definition: liblwgeom.h:352
int lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS3D *dl)
Finds all combinationes of segments between two pointarrays.
Definition: measures3d.c:630
int lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
line to line calculation
Definition: measures3d.c:447
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:1229
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:81
int get_3dcross_product(VECTOR3D *v1, VECTOR3D *v2, VECTOR3D *v)
Definition: measures3d.h:95
#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:512
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:840
#define COLLECTIONTYPE
Definition: liblwgeom.h:66
POINTARRAY * points
Definition: liblwgeom.h:378