PostGIS  2.3.8dev-r@@SVN_REVISION@@
measures3d.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright 2011 Nicklas Avén
22  *
23  **********************************************************************/
24 
25 
26 /**********************************************************************
27  *
28  * PostGIS - Spatial Types for PostgreSQL
29  * http://postgis.net
30  * Copyright 2011 Nicklas Avén
31  *
32  * This is free software; you can redistribute and/or modify it under
33  * the terms of the GNU General Public Licence. See the COPYING file.
34  *
35  **********************************************************************/
36 
37 #include <string.h>
38 #include <stdlib.h>
39 
40 #include "measures3d.h"
41 #include "lwgeom_log.h"
42 
43 
44 static inline int
46 {
47  v->x=p2->x-p1->x;
48  v->y=p2->y-p1->y;
49  v->z=p2->z-p1->z;
50 
51  return LW_TRUE;
52 }
53 
54 static inline int
56 {
57  v->x=(v1->y*v2->z)-(v1->z*v2->y);
58  v->y=(v1->z*v2->x)-(v1->x*v2->z);
59  v->z=(v1->x*v2->y)-(v1->y*v2->x);
60 
61  return LW_TRUE;
62 }
63 
64 
70 static
71 LWGEOM* create_v_line(const LWGEOM *lwgeom,double x, double y, int srid)
72 {
73 
74  LWPOINT *lwpoints[2];
75  GBOX gbox;
76  int rv = lwgeom_calculate_gbox(lwgeom, &gbox);
77 
78  if ( rv == LW_FAILURE )
79  return NULL;
80 
81  lwpoints[0] = lwpoint_make3dz(srid, x, y, gbox.zmin);
82  lwpoints[1] = lwpoint_make3dz(srid, x, y, gbox.zmax);
83 
84  return (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
85 }
86 
87 LWGEOM *
88 lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
89 {
90  return lw_dist3d_distanceline(lw1, lw2, lw1->srid, DIST_MIN);
91 }
92 
93 LWGEOM *
95 {
96  return lw_dist3d_distanceline(lw1, lw2, lw1->srid, DIST_MAX);
97 }
98 
99 LWGEOM *
100 lwgeom_closest_point_3d(const LWGEOM *lw1, const LWGEOM *lw2)
101 {
102  return lw_dist3d_distancepoint(lw1, lw2, lw1->srid, DIST_MIN);
103 }
104 
105 
109 LWGEOM *
110 lw_dist3d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
111 {
112  LWDEBUG(2, "lw_dist3d_distanceline is called");
113  double x1,x2,y1,y2, z1, z2, x, y;
114  double initdistance = ( mode == DIST_MIN ? FLT_MAX : -1.0);
115  DISTPTS3D thedl;
116  LWPOINT *lwpoints[2];
117  LWGEOM *result;
118 
119  thedl.mode = mode;
120  thedl.distance = initdistance;
121  thedl.tolerance = 0.0;
122 
123  /*Check if we really have 3D geoemtries*/
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 geoemtries*/
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 stoped 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  /*You shouldn't being able to get here*/
563  lwerror("unspecified error in function lw_dist3d_distribute_bruteforce");
564  return LW_FALSE;
565 }
566 
567 
568 
569 /*------------------------------------------------------------------------------------------------------------
570 End of Preprocessing functions
571 --------------------------------------------------------------------------------------------------------------*/
572 
573 
574 /*------------------------------------------------------------------------------------------------------------
575 Brute force functions
576 So far the only way to do 3D-calculations
577 --------------------------------------------------------------------------------------------------------------*/
578 
583 int
585 {
586  POINT3DZ p1;
587  POINT3DZ p2;
588  LWDEBUG(2, "lw_dist3d_point_point is called");
589 
590  getPoint3dz_p(point1->point, 0, &p1);
591  getPoint3dz_p(point2->point, 0, &p2);
592 
593  return lw_dist3d_pt_pt(&p1, &p2,dl);
594 }
599 int
601 {
602  POINT3DZ p;
603  POINTARRAY *pa = line->points;
604  LWDEBUG(2, "lw_dist3d_point_line is called");
605 
606  getPoint3dz_p(point->point, 0, &p);
607  return lw_dist3d_pt_ptarray(&p, pa, dl);
608 }
609 
621 int
623 {
624  POINT3DZ p, projp;/*projp is "point projected on plane"*/
625  PLANE3D plane;
626  LWDEBUG(2, "lw_dist3d_point_poly is called");
627  getPoint3dz_p(point->point, 0, &p);
628 
629  /*If we are lookig for max distance, longestline or dfullywithin*/
630  if (dl->mode == DIST_MAX)
631  {
632  LWDEBUG(3, "looking for maxdistance");
633  return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl);
634  }
635 
636  /*Find the plane of the polygon, the "holes" have to be on the same plane. so we only care about the boudary*/
637  if(!define_plane(poly->rings[0], &plane))
638  return LW_FALSE;
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  return LW_FALSE;
677 
678  return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl);
679 }
680 
685 int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
686 {
687  PLANE3D plane;
688  LWDEBUG(2, "lw_dist3d_poly_poly is called");
689  if (dl->mode == DIST_MAX)
690  {
691  return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
692  }
693 
694  if(!define_plane(poly2->rings[0], &plane))
695  return LW_FALSE;
696 
697  /*What we do here is to compare the bondary of one polygon with the other polygon
698  and then take the second boudary comparing with the first polygon*/
699  dl->twisted=1;
700  if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl))
701  return LW_FALSE;
702  if(dl->distance==0.0) /*Just check if the answer already is given*/
703  return LW_TRUE;
704 
705  if(!define_plane(poly1->rings[0], &plane))
706  return LW_FALSE;
707  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.*/
708  return lw_dist3d_ptarray_poly(poly2->rings[0], poly1,&plane, dl);
709 }
710 
716 int
718 {
719  int t;
720  POINT3DZ start, end;
721  int twist = dl->twisted;
722 
723  LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
724 
725  getPoint3dz_p(pa, 0, &start);
726 
727  for (t=1; t<pa->npoints; t++)
728  {
729  dl->twisted=twist;
730  getPoint3dz_p(pa, t, &end);
731  if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
732 
733  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
734  start = end;
735  }
736 
737  return LW_TRUE;
738 }
739 
740 
746 int
748 {
749  POINT3DZ c;
750  double r;
751  /*if start==end, then use pt distance */
752  if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )
753  {
754  return lw_dist3d_pt_pt(p,A,dl);
755  }
756 
757 
758  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) );
759 
760  /*This is for finding the 3Dmaxdistance.
761  the maxdistance have to be between two vertexes,
762  compared to mindistance which can be between
763  tvo vertexes vertex.*/
764  if (dl->mode == DIST_MAX)
765  {
766  if (r>=0.5)
767  {
768  return lw_dist3d_pt_pt(p,A,dl);
769  }
770  if (r<0.5)
771  {
772  return lw_dist3d_pt_pt(p,B,dl);
773  }
774  }
775 
776  if (r<0) /*If the first vertex A is closest to the point p*/
777  {
778  return lw_dist3d_pt_pt(p,A,dl);
779  }
780  if (r>1) /*If the second vertex B is closest to the point p*/
781  {
782  return lw_dist3d_pt_pt(p,B,dl);
783  }
784 
785  /*else if the point p is closer to some point between a and b
786  then we find that point and send it to lw_dist3d_pt_pt*/
787  c.x=A->x + r * (B->x-A->x);
788  c.y=A->y + r * (B->y-A->y);
789  c.z=A->z + r * (B->z-A->z);
790 
791  return lw_dist3d_pt_pt(p,&c,dl);
792 }
793 
794 double
796 {
797  double dx = p2->x - p1->x;
798  double dy = p2->y - p1->y;
799  double dz = p2->z - p1->z;
800  return sqrt ( dx*dx + dy*dy + dz*dz);
801 }
802 
803 
811 int
813 {
814  double dx = thep2->x - thep1->x;
815  double dy = thep2->y - thep1->y;
816  double dz = thep2->z - thep1->z;
817  double dist = sqrt ( dx*dx + dy*dy + dz*dz);
818  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 );
819 
820  if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
821  {
822  dl->distance = dist;
823 
824  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*/
825  {
826  dl->p1 = *thep1;
827  dl->p2 = *thep2;
828  }
829  else
830  {
831  dl->p1 = *thep2;
832  dl->p2 = *thep1;
833  }
834  }
835  return LW_TRUE;
836 }
837 
838 
843 int
845 {
846  int t,u;
847  POINT3DZ start, end;
848  POINT3DZ start2, end2;
849  int twist = dl->twisted;
850  LWDEBUGF(2, "lw_dist3d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
851 
852 
853 
854  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*/
855  {
856  for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
857  {
858  getPoint3dz_p(l1, t, &start);
859  for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
860  {
861  getPoint3dz_p(l2, u, &start2);
862  lw_dist3d_pt_pt(&start,&start2,dl);
863  LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
864  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
865  t, u, dl->distance, dl->tolerance);
866  }
867  }
868  }
869  else
870  {
871  getPoint3dz_p(l1, 0, &start);
872  for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
873  {
874  getPoint3dz_p(l1, t, &end);
875  getPoint3dz_p(l2, 0, &start2);
876  for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
877  {
878  getPoint3dz_p(l2, u, &end2);
879  dl->twisted=twist;
880  lw_dist3d_seg_seg(&start, &end, &start2, &end2,dl);
881  LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
882  LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
883  t, u, dl->distance, dl->tolerance);
884  if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
885  start2 = end2;
886  }
887  start = end;
888  }
889  }
890  return LW_TRUE;
891 }
892 
897 int
899 {
900  VECTOR3D v1, v2, vl;
901  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*/
902  POINT3DZ p1, p2;
903  double a, b, c, d, e, D;
904 
905  /*s1p1 and s1p2 are the same point */
906  if ( ( s1p1->x == s1p2->x) && (s1p1->y == s1p2->y) && (s1p1->z == s1p2->z) )
907  {
908  return lw_dist3d_pt_seg(s1p1,s2p1,s2p2,dl);
909  }
910  /*s2p1 and s2p2 are the same point */
911  if ( ( s2p1->x == s2p2->x) && (s2p1->y == s2p2->y) && (s2p1->z == s2p2->z) )
912  {
913  dl->twisted= ((dl->twisted) * (-1));
914  return lw_dist3d_pt_seg(s2p1,s1p1,s1p2,dl);
915  }
916 
917 /*
918  Here we use algorithm from softsurfer.com
919  that can be found here
920  http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
921 */
922 
923  if (!get_3dvector_from_points(s1p1, s1p2, &v1))
924  return LW_FALSE;
925 
926  if (!get_3dvector_from_points(s2p1, s2p2, &v2))
927  return LW_FALSE;
928 
929  if (!get_3dvector_from_points(s2p1, s1p1, &vl))
930  return LW_FALSE;
931 
932  a = DOT(v1,v1);
933  b = DOT(v1,v2);
934  c = DOT(v2,v2);
935  d = DOT(v1,vl);
936  e = DOT(v2,vl);
937  D = a*c - b*b;
938 
939 
940  if (D <0.000000001)
941  { /* the lines are almost parallel*/
942  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.*/
943  if(b>c) /* use the largest denominator*/
944  {
945  s2k=d/b;
946  }
947  else
948  {
949  s2k =e/c;
950  }
951  }
952  else
953  {
954  s1k = (b*e - c*d) / D;
955  s2k = (a*e - b*d) / D;
956  }
957 
958  /* 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*/
959  if(s1k<0.0||s1k>1.0||s2k<0.0||s2k>1.0)
960  {
961  if(s1k<0.0)
962  {
963 
964  if (!lw_dist3d_pt_seg(s1p1, s2p1, s2p2, dl))
965  {
966  return LW_FALSE;
967  }
968  }
969  if(s1k>1.0)
970  {
971 
972  if (!lw_dist3d_pt_seg(s1p2, s2p1, s2p2, dl))
973  {
974  return LW_FALSE;
975  }
976  }
977  if(s2k<0.0)
978  {
979  dl->twisted= ((dl->twisted) * (-1));
980  if (!lw_dist3d_pt_seg(s2p1, s1p1, s1p2, dl))
981  {
982  return LW_FALSE;
983  }
984  }
985  if(s2k>1.0)
986  {
987  dl->twisted= ((dl->twisted) * (-1));
988  if (!lw_dist3d_pt_seg(s2p2, s1p1, s1p2, dl))
989  {
990  return LW_FALSE;
991  }
992  }
993  }
994  else
995  {/*Find the closest point on the edges of both segments*/
996  p1.x=s1p1->x+s1k*(s1p2->x-s1p1->x);
997  p1.y=s1p1->y+s1k*(s1p2->y-s1p1->y);
998  p1.z=s1p1->z+s1k*(s1p2->z-s1p1->z);
999 
1000  p2.x=s2p1->x+s2k*(s2p2->x-s2p1->x);
1001  p2.y=s2p1->y+s2k*(s2p2->y-s2p1->y);
1002  p2.z=s2p1->z+s2k*(s2p2->z-s2p1->z);
1003 
1004  if (!lw_dist3d_pt_pt(&p1,&p2,dl))/* Send the closest points to point-point calculation*/
1005  {
1006  return LW_FALSE;
1007  }
1008  }
1009  return LW_TRUE;
1010 }
1011 
1019 int
1021 {
1022  int i;
1023 
1024  LWDEBUG(2, "lw_dist3d_point_poly called");
1025 
1026 
1027  if(pt_in_ring_3d(projp, poly->rings[0], plane))
1028  {
1029  for (i=1; i<poly->nrings; i++)
1030  {
1031  /* Inside a hole. Distance = pt -> ring */
1032  if ( pt_in_ring_3d(projp, poly->rings[i], plane ))
1033  {
1034  LWDEBUG(3, " inside an hole");
1035  return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);
1036  }
1037  }
1038 
1039  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*/
1040  }
1041  else
1042  {
1043  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*/
1044  }
1045 
1046  return LW_TRUE;
1047 
1048 }
1049 
1055 {
1056 
1057 
1058  int i,j,k;
1059  double f, s1, s2;
1060  VECTOR3D projp1_projp2;
1061  POINT3DZ p1, p2,projp1, projp2, intersectionp;
1062 
1063  getPoint3dz_p(pa, 0, &p1);
1064 
1065  s1=project_point_on_plane(&p1, plane, &projp1); /*the sign of s1 tells us on which side of the plane the point is. */
1066  lw_dist3d_pt_poly(&p1, poly, plane,&projp1, dl);
1067 
1068  for (i=1;i<pa->npoints;i++)
1069  {
1070  int intersects;
1071  getPoint3dz_p(pa, i, &p2);
1072  s2=project_point_on_plane(&p2, plane, &projp2);
1073  lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl);
1074 
1075  /*If s1and s2 has different signs that means they are on different sides of the plane of the polygon.
1076  That means that the edge between the points crosses the plane and might intersect with the polygon*/
1077  if((s1*s2)<=0)
1078  {
1079  f=fabs(s1)/(fabs(s1)+fabs(s2)); /*The size of s1 and s2 is the distance from the point to the plane.*/
1080  get_3dvector_from_points(&projp1, &projp2,&projp1_projp2);
1081 
1082  /*get the point where the line segment crosses the plane*/
1083  intersectionp.x=projp1.x+f*projp1_projp2.x;
1084  intersectionp.y=projp1.y+f*projp1_projp2.y;
1085  intersectionp.z=projp1.z+f*projp1_projp2.z;
1086 
1087  intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/
1088 
1089  if(pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/
1090  {
1091  for (k=1;k<poly->nrings; k++)
1092  {
1093  /* Inside a hole, so no intersection with the polygon*/
1094  if ( pt_in_ring_3d(&intersectionp, poly->rings[k], plane ))
1095  {
1096  intersects=LW_FALSE;
1097  break;
1098  }
1099  }
1100  if(intersects)
1101  {
1102  dl->distance=0.0;
1103  dl->p1.x=intersectionp.x;
1104  dl->p1.y=intersectionp.y;
1105  dl->p1.z=intersectionp.z;
1106 
1107  dl->p2.x=intersectionp.x;
1108  dl->p2.y=intersectionp.y;
1109  dl->p2.z=intersectionp.z;
1110  return LW_TRUE;
1111 
1112  }
1113  }
1114  }
1115 
1116  projp1=projp2;
1117  s1=s2;
1118  p1=p2;
1119  }
1120 
1121  /*check or pointarray against boundary and inner boundaries of the polygon*/
1122  for (j=0;j<poly->nrings;j++)
1123  {
1124  lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);
1125  }
1126 
1127 return LW_TRUE;
1128 }
1129 
1130 
1136 int
1138 {
1139  int i,j, numberofvectors, pointsinslice;
1140  POINT3DZ p, p1, p2;
1141 
1142  double sumx=0;
1143  double sumy=0;
1144  double sumz=0;
1145  double vl; /*vector length*/
1146 
1147  VECTOR3D v1, v2, v;
1148 
1149  if((pa->npoints-1)==3) /*Triangle is special case*/
1150  {
1151  pointsinslice=1;
1152  }
1153  else
1154  {
1155  pointsinslice=(int) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/
1156  }
1157 
1158  /*find the avg point*/
1159  for (i=0;i<(pa->npoints-1);i++)
1160  {
1161  getPoint3dz_p(pa, i, &p);
1162  sumx+=p.x;
1163  sumy+=p.y;
1164  sumz+=p.z;
1165  }
1166  pl->pop.x=(sumx/(pa->npoints-1));
1167  pl->pop.y=(sumy/(pa->npoints-1));
1168  pl->pop.z=(sumz/(pa->npoints-1));
1169 
1170  sumx=0;
1171  sumy=0;
1172  sumz=0;
1173  numberofvectors= floor((pa->npoints-1)/pointsinslice); /*the number of vectors we try can be 3, 4 or 5*/
1174 
1175  getPoint3dz_p(pa, 0, &p1);
1176  for (j=pointsinslice;j<pa->npoints;j+=pointsinslice)
1177  {
1178  getPoint3dz_p(pa, j, &p2);
1179 
1180  if (!get_3dvector_from_points(&(pl->pop), &p1, &v1) || !get_3dvector_from_points(&(pl->pop), &p2, &v2))
1181  return LW_FALSE;
1182  /*perpendicular vector is cross product of v1 and v2*/
1183  if (!get_3dcross_product(&v1,&v2, &v))
1184  return LW_FALSE;
1185  vl=VECTORLENGTH(v);
1186  sumx+=(v.x/vl);
1187  sumy+=(v.y/vl);
1188  sumz+=(v.z/vl);
1189  p1=p2;
1190  }
1191  pl->pv.x=(sumx/numberofvectors);
1192  pl->pv.y=(sumy/numberofvectors);
1193  pl->pv.z=(sumz/numberofvectors);
1194 
1195  return 1;
1196 }
1197 
1202 double
1204 {
1205 /*In our plane definition we have a point on the plane and a normal vektor (pl.pv), perpendicular to the plane
1206 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.
1207 So, we already have a direction from p to find p0, but we don't know the distance.
1208 */
1209 
1210  VECTOR3D v1;
1211  double f;
1212 
1213  if (!get_3dvector_from_points(&(pl->pop), p, &v1))
1214  return LW_FALSE;
1215 
1216  f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv));
1217 
1218  p0->x=p->x+pl->pv.x*f;
1219  p0->y=p->y+pl->pv.y*f;
1220  p0->z=p->z+pl->pv.z*f;
1221 
1222  return f;
1223 }
1224 
1225 
1226 
1227 
1240 int
1241 pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring,PLANE3D *plane)
1242 {
1243 
1244  int cn = 0; /* the crossing number counter */
1245  int i;
1246  POINT3DZ v1, v2;
1247 
1248  POINT3DZ first, last;
1249 
1250  getPoint3dz_p(ring, 0, &first);
1251  getPoint3dz_p(ring, ring->npoints-1, &last);
1252  if ( memcmp(&first, &last, sizeof(POINT3DZ)) )
1253  {
1254  lwerror("pt_in_ring_3d: V[n] != V[0] (%g %g %g!= %g %g %g)",
1255  first.x, first.y, first.z, last.x, last.y, last.z);
1256  return LW_FALSE;
1257  }
1258 
1259  LWDEBUGF(2, "pt_in_ring_3d called with point: %g %g %g", p->x, p->y, p->z);
1260  /* printPA(ring); */
1261 
1262  /* loop through all edges of the polygon */
1263  getPoint3dz_p(ring, 0, &v1);
1264 
1265 
1266  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*/
1267  {
1268  for (i=0; i<ring->npoints-1; i++)
1269  {
1270  double vt;
1271  getPoint3dz_p(ring, i+1, &v2);
1272 
1273  /* edge from vertex i to vertex i+1 */
1274  if
1275  (
1276  /* an upward crossing */
1277  ((v1.y <= p->y) && (v2.y > p->y))
1278  /* a downward crossing */
1279  || ((v1.y > p->y) && (v2.y <= p->y))
1280  )
1281  {
1282 
1283  vt = (double)(p->y - v1.y) / (v2.y - v1.y);
1284 
1285  /* P.x <intersect */
1286  if (p->x < v1.x + vt * (v2.x - v1.x))
1287  {
1288  /* a valid crossing of y=p.y right of p.x */
1289  ++cn;
1290  }
1291  }
1292  v1 = v2;
1293  }
1294  }
1295  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*/
1296  {
1297  for (i=0; i<ring->npoints-1; i++)
1298  {
1299  double vt;
1300  getPoint3dz_p(ring, i+1, &v2);
1301 
1302  /* edge from vertex i to vertex i+1 */
1303  if
1304  (
1305  /* an upward crossing */
1306  ((v1.z <= p->z) && (v2.z > p->z))
1307  /* a downward crossing */
1308  || ((v1.z > p->z) && (v2.z <= p->z))
1309  )
1310  {
1311 
1312  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1313 
1314  /* P.x <intersect */
1315  if (p->x < v1.x + vt * (v2.x - v1.x))
1316  {
1317  /* a valid crossing of y=p.y right of p.x */
1318  ++cn;
1319  }
1320  }
1321  v1 = v2;
1322  }
1323  }
1324  else /*Hopefully we only have the cases where x part of the normal vector is largest left*/
1325  {
1326  for (i=0; i<ring->npoints-1; i++)
1327  {
1328  double vt;
1329  getPoint3dz_p(ring, i+1, &v2);
1330 
1331  /* edge from vertex i to vertex i+1 */
1332  if
1333  (
1334  /* an upward crossing */
1335  ((v1.z <= p->z) && (v2.z > p->z))
1336  /* a downward crossing */
1337  || ((v1.z > p->z) && (v2.z <= p->z))
1338  )
1339  {
1340 
1341  vt = (double)(p->z - v1.z) / (v2.z - v1.z);
1342 
1343  /* P.x <intersect */
1344  if (p->y < v1.y + vt * (v2.y - v1.y))
1345  {
1346  /* a valid crossing of y=p.y right of p.x */
1347  ++cn;
1348  }
1349  }
1350  v1 = v2;
1351  }
1352  }
1353  LWDEBUGF(3, "pt_in_ring_3d returning %d", cn&1);
1354 
1355  return (cn&1); /* 0 if even (out), and 1 if odd (in) */
1356 }
1357 
1358 
1359 
1360 /*------------------------------------------------------------------------------------------------------------
1361 End of Brute force functions
1362 --------------------------------------------------------------------------------------------------------------*/
1363 
1364 
1365 
#define LINETYPE
Definition: liblwgeom.h:85
double z
Definition: liblwgeom.h:333
double y
Definition: liblwgeom.h:333
double distance
Definition: measures3d.h:41
LWGEOM * lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing closestpoint calculations.
Definition: measures.c:131
int lwgeom_is_collection(const LWGEOM *lwgeom)
Determine whether a LWGEOM can contain sub-geometries or not.
Definition: lwgeom.c:1004
int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
line to polygon calculation
Definition: measures3d.c:665
int lw_dist2d_comp(const LWGEOM *lw1, const LWGEOM *lw2, DISTPTS *dl)
This function just deserializes geometries Bboxes is not checked here since it is the subgeometries b...
Definition: measures.c:245
double x
Definition: liblwgeom.h:333
char * r
Definition: cu_in_wkt.c:24
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:89
void lwfree(void *mem)
Definition: lwutil.c:242
double z
Definition: measures3d.h:51
double y
Definition: liblwgeom.h:339
double lwgeom_maxdistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling max distance calculations and dfyllywithin calculations.
Definition: measures.c:181
POINT3DZ p2
Definition: measures3d.h:43
int npoints
Definition: liblwgeom.h:370
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:213
#define POLYGONTYPE
Definition: liblwgeom.h:86
int mode
Definition: measures.h:51
LWLINE * lwline_from_ptarray(int srid, uint32_t npoints, LWPOINT **points)
Definition: lwline.c:243
int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
polygon to polygon calculation
Definition: measures3d.c:685
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:795
POINT2D * p1
Definition: lwtree.h:33
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
double x
Definition: liblwgeom.h:339
VECTOR3D pv
Definition: measures3d.h:58
POINTARRAY * point
Definition: liblwgeom.h:410
Structure used in distance-calculations.
Definition: measures3d.h:39
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
int32_t srid
Definition: liblwgeom.h:398
#define VECTORLENGTH(v)
Definition: measures3d.h:32
int lw_dist3d_seg_seg(POINT3DZ *s1p1, POINT3DZ *s1p2, POINT3DZ *s2p1, POINT3DZ *s2p2, DISTPTS3D *dl)
Finds the two closest points on two linesegments.
Definition: measures3d.c:898
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:1203
POINT3DZ p1
Definition: measures3d.h:42
double y
Definition: measures3d.h:51
POINT2D * p2
Definition: lwtree.h:34
#define LW_FAILURE
Definition: liblwgeom.h:78
LWGEOM * lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing shortestline and longestline calculations.
Definition: measures.c:84
POINT2D p1
Definition: measures.h:49
double z
Definition: liblwgeom.h:339
#define DOT(u, v)
Definition: measures3d.h:31
double tolerance
Definition: measures.h:53
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
Calculate bounding box of a geometry, automatically taking into account whether it is cartesian or ge...
Definition: lwgeom.c:665
LWGEOM * lwgeom_closest_point_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:100
double x
Definition: liblwgeom.h:327
int lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
This function distributes the brute-force for 3D so far the only type, tasks depending on type...
Definition: measures3d.c:480
double zmax
Definition: liblwgeom.h:296
const char * lwtype_name(uint8_t type)
Return the type name string associated with a type number (e.g.
Definition: lwutil.c:216
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:1020
#define DIST_MIN
Definition: measures.h:41
int getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *point)
Definition: lwgeom_api.c:332
LWGEOM * lw_dist3d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing 3dshortestline and 3dlongestline calculations.
Definition: measures3d.c:110
#define LW_FALSE
Definition: liblwgeom.h:76
int lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)
Computes point to polygon distance For mindistance that means: 1)find the plane of the polygon 2)proj...
Definition: measures3d.c:622
Datum intersects(PG_FUNCTION_ARGS)
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
int lw_dist3d_pt_pt(POINT3DZ *thep1, POINT3DZ *thep2, DISTPTS3D *dl)
Compares incomming points and stores the points closest to each other or most far away from each othe...
Definition: measures3d.c:812
LWGEOM ** geoms
Definition: liblwgeom.h:508
LWGEOM * lwgeom_closest_line_3d(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures3d.c:88
LWGEOM * lw_dist3d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
Function initializing 3dclosestpoint calculations.
Definition: measures3d.c:214
int twisted
Definition: measures3d.h:45
POINTARRAY ** rings
Definition: liblwgeom.h:456
LWPOINT * lwpoint_make3dz(int srid, double x, double y, double z)
Definition: lwpoint.c:155
static int get_3dvector_from_points(POINT3DZ *p1, POINT3DZ *p2, VECTOR3D *v)
Definition: measures3d.c:45
POINT2D p2
Definition: measures.h:50
int nrings
Definition: liblwgeom.h:454
int mode
Definition: measures3d.h:44
double y
Definition: liblwgeom.h:327
LWGEOM * lwgeom_furthest_line_3d(LWGEOM *lw1, LWGEOM *lw2)
Definition: measures3d.c:94
double lwgeom_mindistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling 3d min distance calculations and dwithin calculations.
Definition: measures3d.c:363
POINT3DZ pop
Definition: measures3d.h:57
int lw_dist3d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
This is a recursive function delivering every possible combination of subgeometries.
Definition: measures3d.c:399
double lwgeom_maxdistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling 3d max distance calculations and dfullywithin calculations.
Definition: measures3d.c:326
int lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)
point to line calculation
Definition: measures3d.c:600
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:1137
double tolerance
Definition: measures3d.h:46
int pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring, PLANE3D *plane)
pt_in_ring_3d(): crossing number test for a point in a polygon input: p = a point, pa = vertex points of a ring V[n+1] with V[n]=V[0] plane=the plane that the vertex points are lying on returns: 0 = outside, 1 = inside
Definition: measures3d.c:1241
double distance
Definition: measures.h:48
#define DIST_MAX
Definition: measures.h:40
LWCOLLECTION * lwgeom_as_lwcollection(const LWGEOM *lwgeom)
Definition: lwgeom.c:156
double zmin
Definition: liblwgeom.h:295
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
int lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)
point to point calculation
Definition: measures3d.c:584
double lwgeom_maxdistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 3d max distance calculation.
Definition: measures3d.c:314
int lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl)
If searching for min distance, this one finds the closest point on segment A-B from p...
Definition: measures3d.c:747
uint8_t type
Definition: liblwgeom.h:395
double x
Definition: measures3d.h:51
static LWGEOM * create_v_line(const LWGEOM *lwgeom, double x, double y, int srid)
This function is used to create a vertical line used for cases where one if the geometries lacks z-va...
Definition: measures3d.c:71
double lwgeom_mindistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 3d min distance calculation.
Definition: measures3d.c:352
int lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS3D *dl)
Finds all combinationes of segments between two pointarrays.
Definition: measures3d.c:844
int lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
line to line calculation
Definition: measures3d.c:652
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:1310
Structure used in distance-calculations.
Definition: measures.h:46
LWCOLLECTION * lwcollection_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwcollection.c:94
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
int lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa, DISTPTS3D *dl)
search all the segments of pointarray to see which one is closest to p Returns distance between point...
Definition: measures3d.c:717
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, PLANE3D *plane, DISTPTS3D *dl)
Computes pointarray to polygon distance.
Definition: measures3d.c:1054
#define COLLECTIONTYPE
Definition: liblwgeom.h:90
POINTARRAY * points
Definition: liblwgeom.h:421
static int get_3dcross_product(VECTOR3D *v1, VECTOR3D *v2, VECTOR3D *v)
Definition: measures3d.c:55