-\r
-/**********************************************************************\r
- * $Id: measures.c 5439 2010-03-16 03:13:33Z pramsey $\r
- *\r
- * PostGIS - Spatial Types for PostgreSQL\r
- * http://postgis.refractions.net\r
- * Copyright 2011 Nicklas Avén\r
- *\r
- * This is free software; you can redistribute and/or modify it under\r
- * the terms of the GNU General Public Licence. See the COPYING file.\r
- *\r
- **********************************************************************/\r
-\r
-#include <string.h>\r
-#include <stdlib.h>\r
-\r
-#include "measures3d.h"\r
-#include "lwgeom_log.h"\r
-\r
-\r
-/*------------------------------------------------------------------------------------------------------------\r
-Initializing functions\r
-The functions starting the distance-calculation processses\r
---------------------------------------------------------------------------------------------------------------*/\r
-\r
-/**\r
-Function initializing 3dshortestline and 3dlongestline calculations.\r
-*/\r
-LWGEOM *\r
-lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)\r
-{\r
- double x1,x2,y1,y2, z1, z2;\r
- double initdistance = ( mode == DIST_MIN ? MAXFLOAT : -1.0);\r
- DISTPTS3D thedl;\r
- LWPOINT *lwpoints[2];\r
- LWGEOM *result;\r
-\r
- thedl.mode = mode;\r
- thedl.distance = initdistance;\r
- thedl.tolerance = 0.0;\r
-\r
- LWDEBUG(2, "lw_dist3d_distanceline is called");\r
- if (!lw_dist3d_recursive(lw1, lw2, &thedl))\r
- {\r
- /*should never get here. all cases ought to be error handled earlier*/\r
- lwerror("Some unspecified error.");\r
- result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);\r
- }\r
-\r
- /*if thedl.distance is unchanged there where only empty geometries input*/\r
- if (thedl.distance == initdistance)\r
- {\r
- LWDEBUG(3, "didn't find geometries to measure between, returning null");\r
- result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);\r
- }\r
- else\r
- {\r
- x1=thedl.p1.x;\r
- y1=thedl.p1.y;\r
- z1=thedl.p1.z;\r
- x2=thedl.p2.x;\r
- y2=thedl.p2.y;\r
- z2=thedl.p2.z;\r
-\r
-\r
- lwpoints[0] = lwpoint_make3dz(srid, x1, y1, z1);\r
- lwpoints[1] = lwpoint_make3dz(srid, x2, y2, z2);\r
-\r
- result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);\r
- }\r
-\r
- return result;\r
-}\r
-\r
-/**\r
-Function initializing 3dclosestpoint calculations.\r
-*/\r
-LWGEOM *\r
-lw_dist3d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)\r
-{\r
- double x,y,z;\r
- DISTPTS3D thedl;\r
- double initdistance = MAXFLOAT;\r
- LWGEOM *result;\r
-\r
- thedl.mode = mode;\r
- thedl.distance= initdistance;\r
- thedl.tolerance = 0;\r
-\r
- LWDEBUG(2, "lw_dist3d_distancepoint is called");\r
-\r
- if (!lw_dist3d_recursive(lw1, lw2, &thedl))\r
- {\r
- /*should never get here. all cases ought to be error handled earlier*/\r
- lwerror("Some unspecified error.");\r
- result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);\r
- }\r
-\r
- if (thedl.distance == initdistance)\r
- {\r
- LWDEBUG(3, "didn't find geometries to measure between, returning null");\r
- result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);\r
- }\r
- else\r
- {\r
- x=thedl.p1.x;\r
- y=thedl.p1.y;\r
- z=thedl.p1.z;\r
- result = (LWGEOM *)lwpoint_make3dz(srid, x, y, z);\r
- }\r
-\r
- return result;\r
-}\r
-\r
-\r
-/**\r
-Function initialazing 3d max distance calculation\r
-*/\r
-double\r
-lwgeom_maxdistance3d(LWGEOM *lw1, LWGEOM *lw2)\r
-{\r
- LWDEBUG(2, "lwgeom_maxdistance3d is called");\r
-\r
- return lwgeom_maxdistance3d_tolerance( lw1, lw2, 0.0 );\r
-}\r
-\r
-/**\r
-Function handling 3d max distance calculations and dfyllywithin calculations.\r
-The difference is just the tolerance.\r
-*/\r
-double\r
-lwgeom_maxdistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)\r
-{\r
- /*double thedist;*/\r
- DISTPTS3D thedl;\r
- LWDEBUG(2, "lwgeom_maxdistance3d_tolerance is called");\r
- thedl.mode = DIST_MAX;\r
- thedl.distance= -1;\r
- thedl.tolerance = tolerance;\r
- if (lw_dist3d_recursive(lw1, lw2, &thedl))\r
- {\r
- return thedl.distance;\r
- }\r
- /*should never get here. all cases ought to be error handled earlier*/\r
- lwerror("Some unspecified error.");\r
- return -1;\r
-}\r
-\r
-/**\r
- Function initialazing 3d min distance calculation\r
-*/\r
-double\r
-lwgeom_mindistance3d(LWGEOM *lw1, LWGEOM *lw2)\r
-{\r
- LWDEBUG(2, "lwgeom_mindistance3d is called");\r
- return lwgeom_mindistance3d_tolerance( lw1, lw2, 0.0 );\r
-}\r
-\r
-/**\r
- Function handling 3d min distance calculations and dwithin calculations.\r
- The difference is just the tolerance.\r
-*/\r
-double\r
-lwgeom_mindistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)\r
-{\r
- DISTPTS3D thedl;\r
- LWDEBUG(2, "lwgeom_mindistance3d_tolerance is called");\r
- thedl.mode = DIST_MIN;\r
- thedl.distance= MAXFLOAT;\r
- thedl.tolerance = tolerance;\r
- if (lw_dist3d_recursive(lw1, lw2, &thedl))\r
- {\r
- return thedl.distance;\r
- }\r
- /*should never get here. all cases ought to be error handled earlier*/\r
- lwerror("Some unspecified error.");\r
- return MAXFLOAT;\r
-}\r
-\r
-\r
-/*------------------------------------------------------------------------------------------------------------\r
-End of Initializing functions\r
---------------------------------------------------------------------------------------------------------------*/\r
-\r
-/*------------------------------------------------------------------------------------------------------------\r
-Preprocessing functions\r
-Functions preparing geometries for distance-calculations\r
---------------------------------------------------------------------------------------------------------------*/\r
-\r
-\r
-/**\r
-This is a recursive function delivering every possible combinatin of subgeometries\r
-*/\r
-int lw_dist3d_recursive(const LWGEOM *lwg1,const LWGEOM *lwg2, DISTPTS3D *dl)\r
-{\r
- int i, j;\r
- int n1=1;\r
- int n2=1;\r
- LWGEOM *g1 = NULL;\r
- LWGEOM *g2 = NULL;\r
- LWCOLLECTION *c1 = NULL;\r
- LWCOLLECTION *c2 = NULL;\r
-\r
- LWDEBUGF(2, "lw_dist3d_recursive is called with type1=%d, type2=%d", lwg1->type, lwg2->type);\r
-\r
- if (lwgeom_is_collection(lwg1))\r
- {\r
- LWDEBUG(3, "First geometry is collection");\r
- c1 = lwgeom_as_lwcollection(lwg1);\r
- n1 = c1->ngeoms;\r
- }\r
- if (lwgeom_is_collection(lwg2))\r
- {\r
- LWDEBUG(3, "Second geometry is collection");\r
- c2 = lwgeom_as_lwcollection(lwg2);\r
- n2 = c2->ngeoms;\r
- }\r
-\r
- for ( i = 0; i < n1; i++ )\r
- {\r
-\r
- if (lwgeom_is_collection(lwg1))\r
- {\r
- g1 = c1->geoms[i];\r
- }\r
- else\r
- {\r
- g1 = (LWGEOM*)lwg1;\r
- }\r
-\r
- if (lwgeom_is_empty(g1)) return LW_TRUE;\r
-\r
- if (lwgeom_is_collection(g1))\r
- {\r
- LWDEBUG(3, "Found collection inside first geometry collection, recursing");\r
- if (!lw_dist3d_recursive(g1, lwg2, dl)) return LW_FALSE;\r
- continue;\r
- }\r
- for ( j = 0; j < n2; j++ )\r
- {\r
- if (lwgeom_is_collection(lwg2))\r
- {\r
- g2 = c2->geoms[j];\r
- }\r
- else\r
- {\r
- g2 = (LWGEOM*)lwg2;\r
- }\r
- if (lwgeom_is_collection(g2))\r
- {\r
- LWDEBUG(3, "Found collection inside second geometry collection, recursing");\r
- if (!lw_dist3d_recursive(g1, g2, dl)) return LW_FALSE;\r
- continue;\r
- }\r
-\r
-\r
- /*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/\r
- if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;\r
-\r
-\r
- if (!lw_dist3d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;\r
- if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/\r
- }\r
- }\r
- return LW_TRUE;\r
-}\r
-\r
-\r
-\r
-/**\r
-\r
-This function distributes the brut-force for 3D so far the only type, tasks depending on type\r
-*/\r
-int\r
-lw_dist3d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS3D *dl)\r
-{\r
-\r
- int t1 = lwg1->type;\r
- int t2 = lwg2->type;\r
-\r
- LWDEBUGF(2, "lw_dist3d_distribute_bruteforce is called with typ1=%d, type2=%d", lwg1->type, lwg2->type);\r
-\r
- if ( t1 == POINTTYPE )\r
- {\r
- if ( t2 == POINTTYPE )\r
- {\r
- dl->twisted=1;\r
- return lw_dist3d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);\r
- }\r
- else if ( t2 == LINETYPE )\r
- {\r
- dl->twisted=1;\r
- return lw_dist3d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);\r
- }\r
- else if ( t2 == POLYGONTYPE )\r
- {\r
- dl->twisted=1;\r
- return lw_dist3d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2,dl);\r
- }\r
- else\r
- {\r
- lwerror("Unsupported geometry type: %s", lwtype_name(t2));\r
- return LW_FALSE;\r
- }\r
- }\r
- else if ( t1 == LINETYPE )\r
- {\r
- if ( t2 == POINTTYPE )\r
- {\r
- dl->twisted=(-1);\r
- return lw_dist3d_point_line((LWPOINT *)lwg2,(LWLINE *)lwg1,dl);\r
- }\r
- else if ( t2 == LINETYPE )\r
- {\r
- dl->twisted=1;\r
- return lw_dist3d_line_line((LWLINE *)lwg1,(LWLINE *)lwg2,dl);\r
- }\r
- else if ( t2 == POLYGONTYPE )\r
- {\r
- dl->twisted=1;\r
- return lw_dist3d_line_poly((LWLINE *)lwg1,(LWPOLY *)lwg2,dl);\r
- }\r
- else\r
- {\r
- lwerror("Unsupported geometry type: %s", lwtype_name(t2));\r
- return LW_FALSE;\r
- }\r
- }\r
- else if ( t1 == POLYGONTYPE )\r
- {\r
- if ( t2 == POLYGONTYPE )\r
- {\r
- dl->twisted=1;\r
- return lw_dist3d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2,dl);\r
- }\r
- else if ( t2 == POINTTYPE )\r
- {\r
- dl->twisted=-1;\r
- return lw_dist3d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1,dl);\r
- }\r
- else if ( t2 == LINETYPE )\r
- {\r
- dl->twisted=-1;\r
- return lw_dist3d_line_poly((LWLINE *)lwg2,(LWPOLY *)lwg1,dl);\r
- }\r
- else\r
- {\r
- lwerror("Unsupported geometry type: %s", lwtype_name(t2));\r
- return LW_FALSE;\r
- }\r
- }\r
- else\r
- {\r
- lwerror("Unsupported geometry type: %s", lwtype_name(t1));\r
- return LW_FALSE;\r
- }\r
- /*You shouldn't being able to get here*/\r
- lwerror("unspecified error in function lw_dist3d_distribute_bruteforce");\r
- return LW_FALSE;\r
-}\r
-\r
-\r
-\r
-/*------------------------------------------------------------------------------------------------------------\r
-End of Preprocessing functions\r
---------------------------------------------------------------------------------------------------------------*/\r
-\r
-\r
-/*------------------------------------------------------------------------------------------------------------\r
-Brute force functions\r
-So far the only way to do 3D-calculations\r
---------------------------------------------------------------------------------------------------------------*/\r
-\r
-/**\r
-\r
-point to point calculation\r
-*/\r
-int\r
-lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)\r
-{\r
- LWDEBUG(2, "lw_dist3d_point_point is called");\r
- POINT3DZ p1;\r
- POINT3DZ p2;\r
-\r
- getPoint3dz_p(point1->point, 0, &p1);\r
- getPoint3dz_p(point2->point, 0, &p2);\r
-\r
- return lw_dist3d_pt_pt(&p1, &p2,dl);\r
-}\r
-/**\r
-\r
-point to line calculation\r
-*/\r
-int\r
-lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)\r
-{\r
- LWDEBUG(2, "lw_dist3d_point_line is called");\r
- POINT3DZ p;\r
- POINTARRAY *pa = line->points;\r
-\r
- getPoint3dz_p(point->point, 0, &p);\r
- return lw_dist3d_pt_ptarray(&p, pa, dl);\r
-}\r
-\r
-/**\r
-\r
-Computes point to polygon distance\r
-For mindistance that means:\r
-1)find the plane of the polygon \r
-2)projecting the point to the plane of the polygon \r
-3)finding if that projected point is inside the polygon, if so the distance is measured to that projected point\r
-4) if not in polygon above, check the distance against the boundary of the polygon\r
-for max distance it is always point against boundary\r
-\r
-*/\r
-int\r
-lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)\r
-{\r
- LWDEBUG(2, "lw_dist3d_point_poly is called");\r
- POINT3DZ p, projp;/*projp is "point projected on plane"*/\r
- PLANE3D plane;\r
- getPoint3dz_p(point->point, 0, &p);\r
- \r
- /*If we are lookig for max distance, longestline or dfullywithin*/\r
- if (dl->mode == DIST_MAX)\r
- {\r
- LWDEBUG(3, "looking for maxdistance");\r
- return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl);\r
- }\r
- \r
- /*Find the plane of the polygon, the "holes" have to be on the same plane. so we only care about the boudary*/\r
- if(!define_plane(poly->rings[0], &plane))\r
- return LW_FALSE;\r
- \r
- /*get our point projected on the plane of the polygon*/\r
- project_point_on_plane(&p, &plane, &projp);\r
- \r
- return lw_dist3d_pt_poly(&p, poly,&plane, &projp, dl);\r
-}\r
-\r
-\r
-/**\r
-\r
-line to line calculation\r
-*/\r
-int\r
-lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)\r
-{\r
- LWDEBUG(2, "lw_dist3d_line_line is called");\r
- POINTARRAY *pa1 = line1->points;\r
- POINTARRAY *pa2 = line2->points;\r
-\r
- return lw_dist3d_ptarray_ptarray(pa1, pa2, dl);\r
-}\r
-\r
-/**\r
-\r
-line to polygon calculation\r
-*/\r
-int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)\r
-{\r
- LWDEBUG(2, "lw_dist3d_line_poly is called"); \r
- PLANE3D plane; \r
- \r
- if (dl->mode == DIST_MAX)\r
- {\r
- return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);\r
- }\r
- \r
- if(!define_plane(poly->rings[0], &plane))\r
- return LW_FALSE;\r
- \r
- return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl);\r
-}\r
-\r
-/**\r
-\r
-polygon to polygon calculation\r
-*/\r
-int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)\r
-{ \r
- LWDEBUG(2, "lw_dist3d_poly_poly is called");\r
- PLANE3D plane; \r
- if (dl->mode == DIST_MAX)\r
- {\r
- return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);\r
- }\r
- \r
- if(!define_plane(poly2->rings[0], &plane))\r
- return LW_FALSE;\r
- \r
- /*What we do here is to compare the bondary of one polygon with the other polygon \r
- and then take the second boudary comparing with the first polygon*/\r
- dl->twisted=1;\r
- if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl))\r
- return LW_FALSE;\r
- if(dl->distance==0.0) /*Just check if the answer already is given*/\r
- return LW_TRUE;\r
- \r
- if(!define_plane(poly1->rings[0], &plane))\r
- return LW_FALSE;\r
- 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.*/\r
- return lw_dist3d_ptarray_poly(poly2->rings[0], poly1,&plane, dl);\r
-}\r
-\r
-/**\r
-\r
- * search all the segments of pointarray to see which one is closest to p\r
- * Returns distance between point and pointarray\r
- */\r
-int\r
-lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa,DISTPTS3D *dl)\r
-{\r
- int t;\r
- POINT3DZ start, end;\r
- int twist = dl->twisted;\r
-\r
- LWDEBUG(2, "lw_dist3d_pt_ptarray is called");\r
-\r
- getPoint3dz_p(pa, 0, &start);\r
-\r
- for (t=1; t<pa->npoints; t++)\r
- {\r
- dl->twisted=twist;\r
- getPoint3dz_p(pa, t, &end);\r
- if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;\r
-\r
- if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/\r
- start = end;\r
- }\r
-\r
- return LW_TRUE;\r
-}\r
-\r
-\r
-/**\r
-\r
-If searching for min distance, this one finds the closest point on segment A-B from p.\r
-if searching for max distance it just sends p-A and p-B to pt-pt calculation\r
-*/\r
-int\r
-lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl)\r
-{\r
- POINT3DZ c;\r
- double r;\r
- /*if start==end, then use pt distance */\r
- if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )\r
- {\r
- return lw_dist3d_pt_pt(p,A,dl);\r
- }\r
-\r
-\r
- 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) );\r
-\r
- /*This is for finding the 3Dmaxdistance.\r
- the maxdistance have to be between two vertexes,\r
- compared to mindistance which can be between\r
- tvo vertexes vertex.*/\r
- if (dl->mode == DIST_MAX)\r
- {\r
- if (r>=0.5)\r
- {\r
- return lw_dist3d_pt_pt(p,A,dl);\r
- }\r
- if (r<0.5)\r
- {\r
- return lw_dist3d_pt_pt(p,B,dl);\r
- }\r
- }\r
-\r
- if (r<0) /*If the first vertex A is closest to the point p*/\r
- {\r
- return lw_dist3d_pt_pt(p,A,dl);\r
- }\r
- if (r>1) /*If the second vertex B is closest to the point p*/\r
- {\r
- return lw_dist3d_pt_pt(p,B,dl);\r
- }\r
-\r
- /*else if the point p is closer to some point between a and b\r
- then we find that point and send it to lw_dist3d_pt_pt*/\r
- c.x=A->x + r * (B->x-A->x);\r
- c.y=A->y + r * (B->y-A->y);\r
- c.z=A->z + r * (B->z-A->z);\r
-\r
- return lw_dist3d_pt_pt(p,&c,dl);\r
-}\r
-\r
-\r
-/**\r
-\r
-Compares incomming points and\r
-stores the points closest to each other\r
-or most far away from each other\r
-depending on dl->mode (max or min)\r
-*/\r
-int\r
-lw_dist3d_pt_pt(POINT3DZ *thep1, POINT3DZ *thep2,DISTPTS3D *dl)\r
-{\r
- 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 );\r
- double dx = thep2->x - thep1->x;\r
- double dy = thep2->y - thep1->y;\r
- double dz = thep2->z - thep1->z;\r
- double dist = sqrt ( dx*dx + dy*dy + dz*dz);\r
-\r
- if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/\r
- {\r
- dl->distance = dist;\r
-\r
- 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*/\r
- {\r
- dl->p1 = *thep1;\r
- dl->p2 = *thep2;\r
- }\r
- else\r
- {\r
- dl->p1 = *thep2;\r
- dl->p2 = *thep1;\r
- }\r
- }\r
- return LW_TRUE;\r
-}\r
-\r
-\r
-/**\r
-\r
-Finds all combinationes of segments between two pointarrays\r
-*/\r
-int\r
-lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS3D *dl)\r
-{\r
- LWDEBUGF(2, "lw_dist3d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);\r
- int t,u;\r
- POINT3DZ start, end;\r
- POINT3DZ start2, end2;\r
- int twist = dl->twisted;\r
-\r
-\r
-\r
- 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*/\r
- {\r
- for (t=0; t<l1->npoints; t++) /*for each segment in L1 */\r
- {\r
- getPoint3dz_p(l1, t, &start);\r
- for (u=0; u<l2->npoints; u++) /*for each segment in L2 */\r
- {\r
- getPoint3dz_p(l2, u, &start2);\r
- lw_dist3d_pt_pt(&start,&start2,dl);\r
- LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);\r
- LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",\r
- t, u, dl->distance, dl->tolerance);\r
- }\r
- }\r
- }\r
- else\r
- {\r
- getPoint3dz_p(l1, 0, &start);\r
- for (t=1; t<l1->npoints; t++) /*for each segment in L1 */\r
- {\r
- getPoint3dz_p(l1, t, &end);\r
- getPoint3dz_p(l2, 0, &start2);\r
- for (u=1; u<l2->npoints; u++) /*for each segment in L2 */\r
- {\r
- getPoint3dz_p(l2, u, &end2);\r
- dl->twisted=twist;\r
- lw_dist3d_seg_seg(&start, &end, &start2, &end2,dl);\r
- LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);\r
- LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",\r
- t, u, dl->distance, dl->tolerance);\r
- if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/\r
- start2 = end2;\r
- }\r
- start = end;\r
- }\r
- }\r
- return LW_TRUE;\r
-}\r
-\r
-/**\r
-\r
-Finds the two closest points on two linesegments\r
-*/\r
-int \r
-lw_dist3d_seg_seg(POINT3DZ *s1p1, POINT3DZ *s1p2, POINT3DZ *s2p1, POINT3DZ *s2p2, DISTPTS3D *dl)\r
-{\r
- /*s1p1 and s1p2 are the same point */\r
- if ( ( s1p1->x == s1p2->x) && (s1p1->y == s1p2->y) && (s1p1->z == s1p2->y) )\r
- {\r
- return lw_dist3d_pt_seg(s1p1,s2p1,s2p2,dl);\r
- }\r
- /*s2p1 and s2p2 are the same point */\r
- if ( ( s2p1->x == s2p2->x) && (s2p1->y == s2p2->y) && (s2p1->z == s2p2->y) )\r
- {\r
- dl->twisted= ((dl->twisted) * (-1));\r
- return lw_dist3d_pt_seg(s2p1,s1p1,s1p2,dl);\r
- }\r
- \r
- \r
-/*\r
- Here we use algorithm from softsurfer.com\r
- that can be found here\r
- http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm\r
-*/\r
- \r
- \r
- VECTOR3D v1, v2, vl;\r
- 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*/\r
- POINT3DZ p1, p2;\r
- \r
- if (!get_3dvector_from_points(s1p1, s1p2, &v1))\r
- return LW_FALSE; \r
-\r
- if (!get_3dvector_from_points(s2p1, s2p2, &v2))\r
- return LW_FALSE; \r
-\r
- if (!get_3dvector_from_points(s2p1, s1p1, &vl))\r
- return LW_FALSE; \r
-\r
- double a = DOT(v1,v1);\r
- double b = DOT(v1,v2);\r
- double c = DOT(v2,v2);\r
- double d = DOT(v1,vl);\r
- double e = DOT(v2,vl);\r
- double D = a*c - b*b; \r
-\r
-\r
- if (D <0.000000001) \r
- { /* the lines are almost parallel*/\r
- 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.*/\r
- if(b>c) /* use the largest denominator*/\r
- {\r
- s2k=d/b;\r
- }\r
- else\r
- {\r
- s2k =e/c;\r
- }\r
- }\r
- else \r
- {\r
- s1k = (b*e - c*d) / D;\r
- s2k = (a*e - b*d) / D;\r
- }\r
-\r
- /* 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*/\r
- if(s1k<0.0||s1k>1.0||s2k<0.0||s2k>1.0)\r
- {\r
- if(s1k<0.0) \r
- {\r
-\r
- if (!lw_dist3d_pt_seg(s1p1, s2p1, s2p2, dl))\r
- {\r
- return LW_FALSE;\r
- }\r
- }\r
- if(s1k>1.0)\r
- {\r
-\r
- if (!lw_dist3d_pt_seg(s1p2, s2p1, s2p2, dl))\r
- {\r
- return LW_FALSE;\r
- }\r
- }\r
- if(s2k<0.0)\r
- {\r
- dl->twisted= ((dl->twisted) * (-1));\r
- if (!lw_dist3d_pt_seg(s2p1, s1p1, s1p2, dl))\r
- {\r
- return LW_FALSE;\r
- }\r
- }\r
- if(s2k>1.0)\r
- {\r
- dl->twisted= ((dl->twisted) * (-1));\r
- if (!lw_dist3d_pt_seg(s2p2, s1p1, s1p2, dl))\r
- {\r
- return LW_FALSE;\r
- }\r
- }\r
- }\r
- else\r
- {/*Find the closest point on the edges of both segments*/\r
- p1.x=s1p1->x+s1k*(s1p2->x-s1p1->x);\r
- p1.y=s1p1->y+s1k*(s1p2->y-s1p1->y);\r
- p1.z=s1p1->z+s1k*(s1p2->z-s1p1->z);\r
-\r
- p2.x=s2p1->x+s2k*(s2p2->x-s2p1->x);\r
- p2.y=s2p1->y+s2k*(s2p2->y-s2p1->y);\r
- p2.z=s2p1->z+s2k*(s2p2->z-s2p1->z);\r
-\r
- if (!lw_dist3d_pt_pt(&p1,&p2,dl))/* Send the closest points to point-point calculation*/\r
- {\r
- return LW_FALSE;\r
- }\r
- }\r
- return LW_TRUE;\r
-}\r
-\r
-/**\r
-\r
-Checking if the point projected on the plane of the polygon actually is inside that polygon. \r
-If so the mindistance is between that projected point and our original point.\r
-If not we check from original point to the bounadary.\r
-If the projected point is inside a hole of the polygon we check the distance to the boudary of that hole.\r
-*/\r
-int\r
-lw_dist3d_pt_poly(POINT3DZ *p, LWPOLY *poly, PLANE3D *plane,POINT3DZ *projp, DISTPTS3D *dl)\r
-{ \r
- int i;\r
- \r
- LWDEBUG(2, "lw_dist3d_point_poly called");\r
-\r
- \r
- if(pt_in_ring_3d(projp, poly->rings[0], plane))\r
- {\r
- for (i=1; i<poly->nrings; i++)\r
- {\r
- /* Inside a hole. Distance = pt -> ring */\r
- if ( pt_in_ring_3d(projp, poly->rings[i], plane ))\r
- {\r
- LWDEBUG(3, " inside an hole");\r
- return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);\r
- }\r
- } \r
- \r
- 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*/\r
- }\r
- else\r
- {\r
- 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*/\r
- } \r
- \r
- return LW_TRUE;\r
- \r
-}\r
-\r
-/**\r
-\r
-Computes pointarray to polygon distance\r
-*/\r
-int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3D *dl)\r
-{\r
- \r
-\r
- int i,j,k;\r
- double f, s1, s2;\r
- VECTOR3D projp1_projp2;\r
- POINT3DZ p1, p2,projp1, projp2, intersectionp;\r
- \r
- getPoint3dz_p(pa, 0, &p1);\r
- \r
- s1=project_point_on_plane(&p1, plane, &projp1); /*the sign of s1 tells us on which side of the plane the point is. */\r
- lw_dist3d_pt_poly(&p1, poly, plane,&projp1, dl); \r
- \r
- for (i=1;i<pa->npoints;i++)\r
- { \r
- getPoint3dz_p(pa, i, &p2);\r
- s2=project_point_on_plane(&p2, plane, &projp2); \r
- lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl);\r
- \r
- /*If s1and s2 has different signs that means they are on different sides of the plane of the polygon.\r
- That means that the edge between the points crosses the plane and might intersect with the polygon*/\r
- if((s1*s2)<=0) \r
- {\r
- f=fabs(s1)/(fabs(s1)+fabs(s2)); /*The size of s1 and s2 is the distance from the point to the plane.*/\r
- get_3dvector_from_points(&projp1, &projp2,&projp1_projp2);\r
- \r
- /*get the point where the line segment crosses the plane*/\r
- intersectionp.x=projp1.x+f*projp1_projp2.x;\r
- intersectionp.y=projp1.y+f*projp1_projp2.y;\r
- intersectionp.z=projp1.z+f*projp1_projp2.z;\r
- \r
- int intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/\r
- \r
- if(pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/\r
- {\r
- for (k=1;k<poly->nrings; k++)\r
- {\r
- /* Inside a hole, so no intersection with the polygon*/\r
- if ( pt_in_ring_3d(&intersectionp, poly->rings[k], plane ))\r
- {\r
- intersects=LW_FALSE;\r
- break;\r
- }\r
- } \r
- if(intersects) \r
- {\r
- dl->distance=0.0;\r
- dl->p1.x=intersectionp.x;\r
- dl->p1.y=intersectionp.y;\r
- dl->p1.z=intersectionp.z;\r
- \r
- dl->p2.x=intersectionp.x;\r
- dl->p2.y=intersectionp.y;\r
- dl->p2.z=intersectionp.z;\r
- return LW_TRUE;\r
- \r
- } \r
- } \r
- }\r
- \r
- projp2=projp1;\r
- s2=s1;\r
- p2=p1;\r
- } \r
- \r
- /*check or pointarray against boundary and inner boundaries of the polygon*/\r
- for (j=0;j<poly->nrings;j++)\r
- {\r
- lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);\r
- }\r
- \r
-return LW_TRUE;\r
-} \r
-\r
-\r
-/**\r
-\r
-Here we define the plane of a polygon (boundary pointarray of a polygon)\r
-the plane is stored as a pont in plane (plane.pop) and a normal vector (plane.pv)\r
-*/\r
-int\r
-define_plane(POINTARRAY *pa, PLANE3D *pl)\r
-{\r
- int i,j, numberofvectors, pointsinslice;\r
- POINT3DZ p, p1, p2;\r
-\r
- double sumx=0;\r
- double sumy=0;\r
- double sumz=0;\r
- double vl; /*vector length*/\r
-\r
- VECTOR3D v1, v2, v;\r
- \r
- if((pa->npoints-1)==3) /*Triangle is special case*/\r
- {\r
- pointsinslice=1; \r
- }\r
- else\r
- {\r
- pointsinslice=(int) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/\r
- }\r
- \r
- /*find the avg point*/\r
- for (i=0;i<(pa->npoints-1);i++)\r
- {\r
- getPoint3dz_p(pa, i, &p);\r
- sumx+=p.x;\r
- sumy+=p.y;\r
- sumz+=p.z; \r
- } \r
- pl->pop.x=(sumx/(pa->npoints-1));\r
- pl->pop.y=(sumy/(pa->npoints-1));\r
- pl->pop.z=(sumz/(pa->npoints-1));\r
- \r
- sumx=0;\r
- sumy=0;\r
- sumz=0;\r
- numberofvectors= floor((pa->npoints-1)/pointsinslice); /*the number of vectors we try can be 3, 4 or 5*/\r
- \r
- getPoint3dz_p(pa, 0, &p1);\r
- for (j=pointsinslice;j<pa->npoints;j+=pointsinslice)\r
- {\r
- getPoint3dz_p(pa, j, &p2); \r
- \r
- if (!get_3dvector_from_points(&(pl->pop), &p1, &v1) || !get_3dvector_from_points(&(pl->pop), &p2, &v2))\r
- return LW_FALSE; \r
- /*perpendicular vector is cross product of v1 and v2*/\r
- if (!get_3dcross_product(&v1,&v2, &v))\r
- return LW_FALSE; \r
- vl=VECTORLENGTH(v);\r
- sumx+=(v.x/vl);\r
- sumy+=(v.y/vl);\r
- sumz+=(v.z/vl); \r
- p1=p2;\r
- }\r
- pl->pv.x=(sumx/numberofvectors);\r
- pl->pv.y=(sumy/numberofvectors);\r
- pl->pv.z=(sumz/numberofvectors);\r
- \r
- return 1;\r
-}\r
-\r
-/**\r
-\r
-Finds a point on a plane from where the original point is perpendicular to the plane\r
-*/\r
-double \r
-project_point_on_plane(POINT3DZ *p, PLANE3D *pl, POINT3DZ *p0)\r
-{\r
-/*In our plane definition we have a point on the plane and a normal vektor (pl.pv), perpendicular to the plane\r
-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.\r
-So, we already have a direction from p to find p0, but we don't know the distance.\r
-*/\r
-\r
- VECTOR3D v1;\r
- double f;\r
- \r
- if (!get_3dvector_from_points(&(pl->pop), p, &v1))\r
- return LW_FALSE; \r
- \r
- f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv));\r
- \r
- p0->x=p->x+pl->pv.x*f;\r
- p0->y=p->y+pl->pv.y*f;\r
- p0->z=p->z+pl->pv.z*f; \r
- \r
- return f; \r
-}\r
-\r
-\r
-\r
-\r
-/**\r
- * pt_in_ring_3d(): crossing number test for a point in a polygon\r
- * input: p = a point,\r
- * pa = vertex points of a ring V[n+1] with V[n]=V[0]\r
-* plane=the plane that the vertex points are lying on\r
- * returns: 0 = outside, 1 = inside\r
- *\r
- * Our polygons have first and last point the same,\r
- *\r
-* The difference in 3D variant is that we exclude the dimension that faces the plane least.\r
-* That is the dimension with the highest number in pv\r
- */\r
-int\r
-pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring,PLANE3D *plane)\r
-{\r
- \r
- int cn = 0; /* the crossing number counter */\r
- int i;\r
- POINT3DZ v1, v2;\r
-\r
- POINT3DZ first, last;\r
-\r
- getPoint3dz_p(ring, 0, &first);\r
- getPoint3dz_p(ring, ring->npoints-1, &last);\r
- if ( memcmp(&first, &last, sizeof(POINT3DZ)) )\r
- {\r
- lwerror("pt_in_ring_3d: V[n] != V[0] (%g %g %g!= %g %g %g)",\r
- first.x, first.y, first.z, last.x, last.y, last.z);\r
- return LW_FALSE;\r
- }\r
-\r
- LWDEBUGF(2, "pt_in_ring_3d called with point: %g %g %g", p->x, p->y, p->z);\r
- /* printPA(ring); */\r
-\r
- /* loop through all edges of the polygon */\r
- getPoint3dz_p(ring, 0, &v1);\r
- \r
- \r
- 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*/\r
- {\r
- for (i=0; i<ring->npoints-1; i++)\r
- {\r
- double vt;\r
- getPoint3dz_p(ring, i+1, &v2);\r
-\r
- /* edge from vertex i to vertex i+1 */\r
- if\r
- (\r
- /* an upward crossing */\r
- ((v1.y <= p->y) && (v2.y > p->y))\r
- /* a downward crossing */\r
- || ((v1.y > p->y) && (v2.y <= p->y))\r
- )\r
- {\r
-\r
- vt = (double)(p->y - v1.y) / (v2.y - v1.y);\r
-\r
- /* P.x <intersect */\r
- if (p->x < v1.x + vt * (v2.x - v1.x))\r
- {\r
- /* a valid crossing of y=p.y right of p.x */\r
- ++cn;\r
- }\r
- }\r
- v1 = v2;\r
- }\r
- }\r
- 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*/\r
- {\r
- for (i=0; i<ring->npoints-1; i++)\r
- {\r
- double vt;\r
- getPoint3dz_p(ring, i+1, &v2);\r
-\r
- /* edge from vertex i to vertex i+1 */\r
- if\r
- (\r
- /* an upward crossing */\r
- ((v1.z <= p->z) && (v2.z > p->z))\r
- /* a downward crossing */\r
- || ((v1.z > p->z) && (v2.z <= p->z))\r
- )\r
- {\r
-\r
- vt = (double)(p->z - v1.z) / (v2.z - v1.z);\r
-\r
- /* P.x <intersect */\r
- if (p->x < v1.x + vt * (v2.x - v1.x))\r
- {\r
- /* a valid crossing of y=p.y right of p.x */\r
- ++cn;\r
- }\r
- }\r
- v1 = v2;\r
- }\r
- }\r
- else /*Hopefully we only have the cases where x part of the normal vector is largest left*/\r
- {\r
- for (i=0; i<ring->npoints-1; i++)\r
- {\r
- double vt;\r
- getPoint3dz_p(ring, i+1, &v2);\r
-\r
- /* edge from vertex i to vertex i+1 */\r
- if\r
- (\r
- /* an upward crossing */\r
- ((v1.z <= p->z) && (v2.z > p->z))\r
- /* a downward crossing */\r
- || ((v1.z > p->z) && (v2.z <= p->z))\r
- )\r
- {\r
-\r
- vt = (double)(p->z - v1.z) / (v2.z - v1.z);\r
-\r
- /* P.x <intersect */\r
- if (p->y < v1.y + vt * (v2.y - v1.y))\r
- {\r
- /* a valid crossing of y=p.y right of p.x */\r
- ++cn;\r
- }\r
- }\r
- v1 = v2;\r
- }\r
- }\r
- LWDEBUGF(3, "pt_in_ring_3d returning %d", cn&1);\r
-\r
- return (cn&1); /* 0 if even (out), and 1 if odd (in) */\r
-}\r
-\r
-\r
-\r
-/*------------------------------------------------------------------------------------------------------------\r
-End of Brute force functions\r
---------------------------------------------------------------------------------------------------------------*/\r
-\r
-\r
-\r
+
+/**********************************************************************
+ * $Id: measures.c 5439 2010-03-16 03:13:33Z pramsey $
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2011 Nicklas Avén
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include <string.h>
+#include <stdlib.h>
+
+#include "measures3d.h"
+#include "lwgeom_log.h"
+
+
+/*------------------------------------------------------------------------------------------------------------
+Initializing functions
+The functions starting the distance-calculation processses
+--------------------------------------------------------------------------------------------------------------*/
+
+/**
+Function initializing 3dshortestline and 3dlongestline calculations.
+*/
+LWGEOM *
+lw_dist3d_distanceline(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
+{
+ double x1,x2,y1,y2, z1, z2;
+ double initdistance = ( mode == DIST_MIN ? MAXFLOAT : -1.0);
+ DISTPTS3D thedl;
+ LWPOINT *lwpoints[2];
+ LWGEOM *result;
+
+ thedl.mode = mode;
+ thedl.distance = initdistance;
+ thedl.tolerance = 0.0;
+
+ LWDEBUG(2, "lw_dist3d_distanceline is called");
+ if (!lw_dist3d_recursive(lw1, lw2, &thedl))
+ {
+ /*should never get here. all cases ought to be error handled earlier*/
+ lwerror("Some unspecified error.");
+ result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
+ }
+
+ /*if thedl.distance is unchanged there where only empty geometries input*/
+ if (thedl.distance == initdistance)
+ {
+ LWDEBUG(3, "didn't find geometries to measure between, returning null");
+ result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
+ }
+ else
+ {
+ x1=thedl.p1.x;
+ y1=thedl.p1.y;
+ z1=thedl.p1.z;
+ x2=thedl.p2.x;
+ y2=thedl.p2.y;
+ z2=thedl.p2.z;
+
+
+ lwpoints[0] = lwpoint_make3dz(srid, x1, y1, z1);
+ lwpoints[1] = lwpoint_make3dz(srid, x2, y2, z2);
+
+ result = (LWGEOM *)lwline_from_ptarray(srid, 2, lwpoints);
+ }
+
+ return result;
+}
+
+/**
+Function initializing 3dclosestpoint calculations.
+*/
+LWGEOM *
+lw_dist3d_distancepoint(LWGEOM *lw1, LWGEOM *lw2,int srid,int mode)
+{
+ double x,y,z;
+ DISTPTS3D thedl;
+ double initdistance = MAXFLOAT;
+ LWGEOM *result;
+
+ thedl.mode = mode;
+ thedl.distance= initdistance;
+ thedl.tolerance = 0;
+
+ LWDEBUG(2, "lw_dist3d_distancepoint is called");
+
+ if (!lw_dist3d_recursive(lw1, lw2, &thedl))
+ {
+ /*should never get here. all cases ought to be error handled earlier*/
+ lwerror("Some unspecified error.");
+ result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
+ }
+
+ if (thedl.distance == initdistance)
+ {
+ LWDEBUG(3, "didn't find geometries to measure between, returning null");
+ result = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
+ }
+ else
+ {
+ x=thedl.p1.x;
+ y=thedl.p1.y;
+ z=thedl.p1.z;
+ result = (LWGEOM *)lwpoint_make3dz(srid, x, y, z);
+ }
+
+ return result;
+}
+
+
+/**
+Function initialazing 3d max distance calculation
+*/
+double
+lwgeom_maxdistance3d(LWGEOM *lw1, LWGEOM *lw2)
+{
+ LWDEBUG(2, "lwgeom_maxdistance3d is called");
+
+ return lwgeom_maxdistance3d_tolerance( lw1, lw2, 0.0 );
+}
+
+/**
+Function handling 3d max distance calculations and dfyllywithin calculations.
+The difference is just the tolerance.
+*/
+double
+lwgeom_maxdistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
+{
+ /*double thedist;*/
+ DISTPTS3D thedl;
+ LWDEBUG(2, "lwgeom_maxdistance3d_tolerance is called");
+ thedl.mode = DIST_MAX;
+ thedl.distance= -1;
+ thedl.tolerance = tolerance;
+ if (lw_dist3d_recursive(lw1, lw2, &thedl))
+ {
+ return thedl.distance;
+ }
+ /*should never get here. all cases ought to be error handled earlier*/
+ lwerror("Some unspecified error.");
+ return -1;
+}
+
+/**
+ Function initialazing 3d min distance calculation
+*/
+double
+lwgeom_mindistance3d(LWGEOM *lw1, LWGEOM *lw2)
+{
+ LWDEBUG(2, "lwgeom_mindistance3d is called");
+ return lwgeom_mindistance3d_tolerance( lw1, lw2, 0.0 );
+}
+
+/**
+ Function handling 3d min distance calculations and dwithin calculations.
+ The difference is just the tolerance.
+*/
+double
+lwgeom_mindistance3d_tolerance(LWGEOM *lw1, LWGEOM *lw2, double tolerance)
+{
+ DISTPTS3D thedl;
+ LWDEBUG(2, "lwgeom_mindistance3d_tolerance is called");
+ thedl.mode = DIST_MIN;
+ thedl.distance= MAXFLOAT;
+ thedl.tolerance = tolerance;
+ if (lw_dist3d_recursive(lw1, lw2, &thedl))
+ {
+ return thedl.distance;
+ }
+ /*should never get here. all cases ought to be error handled earlier*/
+ lwerror("Some unspecified error.");
+ return MAXFLOAT;
+}
+
+
+/*------------------------------------------------------------------------------------------------------------
+End of Initializing functions
+--------------------------------------------------------------------------------------------------------------*/
+
+/*------------------------------------------------------------------------------------------------------------
+Preprocessing functions
+Functions preparing geometries for distance-calculations
+--------------------------------------------------------------------------------------------------------------*/
+
+
+/**
+This is a recursive function delivering every possible combinatin of subgeometries
+*/
+int lw_dist3d_recursive(const LWGEOM *lwg1,const LWGEOM *lwg2, DISTPTS3D *dl)
+{
+ int i, j;
+ int n1=1;
+ int n2=1;
+ LWGEOM *g1 = NULL;
+ LWGEOM *g2 = NULL;
+ LWCOLLECTION *c1 = NULL;
+ LWCOLLECTION *c2 = NULL;
+
+ LWDEBUGF(2, "lw_dist3d_recursive is called with type1=%d, type2=%d", lwg1->type, lwg2->type);
+
+ if (lwgeom_is_collection(lwg1))
+ {
+ LWDEBUG(3, "First geometry is collection");
+ c1 = lwgeom_as_lwcollection(lwg1);
+ n1 = c1->ngeoms;
+ }
+ if (lwgeom_is_collection(lwg2))
+ {
+ LWDEBUG(3, "Second geometry is collection");
+ c2 = lwgeom_as_lwcollection(lwg2);
+ n2 = c2->ngeoms;
+ }
+
+ for ( i = 0; i < n1; i++ )
+ {
+
+ if (lwgeom_is_collection(lwg1))
+ {
+ g1 = c1->geoms[i];
+ }
+ else
+ {
+ g1 = (LWGEOM*)lwg1;
+ }
+
+ if (lwgeom_is_empty(g1)) return LW_TRUE;
+
+ if (lwgeom_is_collection(g1))
+ {
+ LWDEBUG(3, "Found collection inside first geometry collection, recursing");
+ if (!lw_dist3d_recursive(g1, lwg2, dl)) return LW_FALSE;
+ continue;
+ }
+ for ( j = 0; j < n2; j++ )
+ {
+ if (lwgeom_is_collection(lwg2))
+ {
+ g2 = c2->geoms[j];
+ }
+ else
+ {
+ g2 = (LWGEOM*)lwg2;
+ }
+ if (lwgeom_is_collection(g2))
+ {
+ LWDEBUG(3, "Found collection inside second geometry collection, recursing");
+ if (!lw_dist3d_recursive(g1, g2, dl)) return LW_FALSE;
+ continue;
+ }
+
+
+ /*If one of geometries is empty, return. True here only means continue searching. False would have stoped the process*/
+ if (lwgeom_is_empty(g1)||lwgeom_is_empty(g2)) return LW_TRUE;
+
+
+ if (!lw_dist3d_distribute_bruteforce(g1, g2, dl)) return LW_FALSE;
+ if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
+ }
+ }
+ return LW_TRUE;
+}
+
+
+
+/**
+
+This function distributes the brut-force for 3D so far the only type, tasks depending on type
+*/
+int
+lw_dist3d_distribute_bruteforce(LWGEOM *lwg1, LWGEOM *lwg2, DISTPTS3D *dl)
+{
+
+ int t1 = lwg1->type;
+ int t2 = lwg2->type;
+
+ LWDEBUGF(2, "lw_dist3d_distribute_bruteforce is called with typ1=%d, type2=%d", lwg1->type, lwg2->type);
+
+ if ( t1 == POINTTYPE )
+ {
+ if ( t2 == POINTTYPE )
+ {
+ dl->twisted=1;
+ return lw_dist3d_point_point((LWPOINT *)lwg1, (LWPOINT *)lwg2, dl);
+ }
+ else if ( t2 == LINETYPE )
+ {
+ dl->twisted=1;
+ return lw_dist3d_point_line((LWPOINT *)lwg1, (LWLINE *)lwg2, dl);
+ }
+ else if ( t2 == POLYGONTYPE )
+ {
+ dl->twisted=1;
+ return lw_dist3d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2,dl);
+ }
+ else
+ {
+ lwerror("Unsupported geometry type: %s", lwtype_name(t2));
+ return LW_FALSE;
+ }
+ }
+ else if ( t1 == LINETYPE )
+ {
+ if ( t2 == POINTTYPE )
+ {
+ dl->twisted=(-1);
+ return lw_dist3d_point_line((LWPOINT *)lwg2,(LWLINE *)lwg1,dl);
+ }
+ else if ( t2 == LINETYPE )
+ {
+ dl->twisted=1;
+ return lw_dist3d_line_line((LWLINE *)lwg1,(LWLINE *)lwg2,dl);
+ }
+ else if ( t2 == POLYGONTYPE )
+ {
+ dl->twisted=1;
+ return lw_dist3d_line_poly((LWLINE *)lwg1,(LWPOLY *)lwg2,dl);
+ }
+ else
+ {
+ lwerror("Unsupported geometry type: %s", lwtype_name(t2));
+ return LW_FALSE;
+ }
+ }
+ else if ( t1 == POLYGONTYPE )
+ {
+ if ( t2 == POLYGONTYPE )
+ {
+ dl->twisted=1;
+ return lw_dist3d_poly_poly((LWPOLY *)lwg1, (LWPOLY *)lwg2,dl);
+ }
+ else if ( t2 == POINTTYPE )
+ {
+ dl->twisted=-1;
+ return lw_dist3d_point_poly((LWPOINT *)lwg2, (LWPOLY *)lwg1,dl);
+ }
+ else if ( t2 == LINETYPE )
+ {
+ dl->twisted=-1;
+ return lw_dist3d_line_poly((LWLINE *)lwg2,(LWPOLY *)lwg1,dl);
+ }
+ else
+ {
+ lwerror("Unsupported geometry type: %s", lwtype_name(t2));
+ return LW_FALSE;
+ }
+ }
+ else
+ {
+ lwerror("Unsupported geometry type: %s", lwtype_name(t1));
+ return LW_FALSE;
+ }
+ /*You shouldn't being able to get here*/
+ lwerror("unspecified error in function lw_dist3d_distribute_bruteforce");
+ return LW_FALSE;
+}
+
+
+
+/*------------------------------------------------------------------------------------------------------------
+End of Preprocessing functions
+--------------------------------------------------------------------------------------------------------------*/
+
+
+/*------------------------------------------------------------------------------------------------------------
+Brute force functions
+So far the only way to do 3D-calculations
+--------------------------------------------------------------------------------------------------------------*/
+
+/**
+
+point to point calculation
+*/
+int
+lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl)
+{
+ LWDEBUG(2, "lw_dist3d_point_point is called");
+ POINT3DZ p1;
+ POINT3DZ p2;
+
+ getPoint3dz_p(point1->point, 0, &p1);
+ getPoint3dz_p(point2->point, 0, &p2);
+
+ return lw_dist3d_pt_pt(&p1, &p2,dl);
+}
+/**
+
+point to line calculation
+*/
+int
+lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl)
+{
+ LWDEBUG(2, "lw_dist3d_point_line is called");
+ POINT3DZ p;
+ POINTARRAY *pa = line->points;
+
+ getPoint3dz_p(point->point, 0, &p);
+ return lw_dist3d_pt_ptarray(&p, pa, dl);
+}
+
+/**
+
+Computes point to polygon distance
+For mindistance that means:
+1)find the plane of the polygon
+2)projecting the point to the plane of the polygon
+3)finding if that projected point is inside the polygon, if so the distance is measured to that projected point
+4) if not in polygon above, check the distance against the boundary of the polygon
+for max distance it is always point against boundary
+
+*/
+int
+lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)
+{
+ LWDEBUG(2, "lw_dist3d_point_poly is called");
+ POINT3DZ p, projp;/*projp is "point projected on plane"*/
+ PLANE3D plane;
+ getPoint3dz_p(point->point, 0, &p);
+
+ /*If we are lookig for max distance, longestline or dfullywithin*/
+ if (dl->mode == DIST_MAX)
+ {
+ LWDEBUG(3, "looking for maxdistance");
+ return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl);
+ }
+
+ /*Find the plane of the polygon, the "holes" have to be on the same plane. so we only care about the boudary*/
+ if(!define_plane(poly->rings[0], &plane))
+ return LW_FALSE;
+
+ /*get our point projected on the plane of the polygon*/
+ project_point_on_plane(&p, &plane, &projp);
+
+ return lw_dist3d_pt_poly(&p, poly,&plane, &projp, dl);
+}
+
+
+/**
+
+line to line calculation
+*/
+int
+lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
+{
+ LWDEBUG(2, "lw_dist3d_line_line is called");
+ POINTARRAY *pa1 = line1->points;
+ POINTARRAY *pa2 = line2->points;
+
+ return lw_dist3d_ptarray_ptarray(pa1, pa2, dl);
+}
+
+/**
+
+line to polygon calculation
+*/
+int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
+{
+ LWDEBUG(2, "lw_dist3d_line_poly is called");
+ PLANE3D plane;
+
+ if (dl->mode == DIST_MAX)
+ {
+ return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);
+ }
+
+ if(!define_plane(poly->rings[0], &plane))
+ return LW_FALSE;
+
+ return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl);
+}
+
+/**
+
+polygon to polygon calculation
+*/
+int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
+{
+ LWDEBUG(2, "lw_dist3d_poly_poly is called");
+ PLANE3D plane;
+ if (dl->mode == DIST_MAX)
+ {
+ return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl);
+ }
+
+ if(!define_plane(poly2->rings[0], &plane))
+ return LW_FALSE;
+
+ /*What we do here is to compare the bondary of one polygon with the other polygon
+ and then take the second boudary comparing with the first polygon*/
+ dl->twisted=1;
+ if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl))
+ return LW_FALSE;
+ if(dl->distance==0.0) /*Just check if the answer already is given*/
+ return LW_TRUE;
+
+ if(!define_plane(poly1->rings[0], &plane))
+ return LW_FALSE;
+ 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.*/
+ return lw_dist3d_ptarray_poly(poly2->rings[0], poly1,&plane, dl);
+}
+
+/**
+
+ * search all the segments of pointarray to see which one is closest to p
+ * Returns distance between point and pointarray
+ */
+int
+lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa,DISTPTS3D *dl)
+{
+ int t;
+ POINT3DZ start, end;
+ int twist = dl->twisted;
+
+ LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
+
+ getPoint3dz_p(pa, 0, &start);
+
+ for (t=1; t<pa->npoints; t++)
+ {
+ dl->twisted=twist;
+ getPoint3dz_p(pa, t, &end);
+ if (!lw_dist3d_pt_seg(p, &start, &end,dl)) return LW_FALSE;
+
+ if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
+ start = end;
+ }
+
+ return LW_TRUE;
+}
+
+
+/**
+
+If searching for min distance, this one finds the closest point on segment A-B from p.
+if searching for max distance it just sends p-A and p-B to pt-pt calculation
+*/
+int
+lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl)
+{
+ POINT3DZ c;
+ double r;
+ /*if start==end, then use pt distance */
+ if ( ( A->x == B->x) && (A->y == B->y) && (A->z == B->z) )
+ {
+ return lw_dist3d_pt_pt(p,A,dl);
+ }
+
+
+ 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) );
+
+ /*This is for finding the 3Dmaxdistance.
+ the maxdistance have to be between two vertexes,
+ compared to mindistance which can be between
+ tvo vertexes vertex.*/
+ if (dl->mode == DIST_MAX)
+ {
+ if (r>=0.5)
+ {
+ return lw_dist3d_pt_pt(p,A,dl);
+ }
+ if (r<0.5)
+ {
+ return lw_dist3d_pt_pt(p,B,dl);
+ }
+ }
+
+ if (r<0) /*If the first vertex A is closest to the point p*/
+ {
+ return lw_dist3d_pt_pt(p,A,dl);
+ }
+ if (r>1) /*If the second vertex B is closest to the point p*/
+ {
+ return lw_dist3d_pt_pt(p,B,dl);
+ }
+
+ /*else if the point p is closer to some point between a and b
+ then we find that point and send it to lw_dist3d_pt_pt*/
+ c.x=A->x + r * (B->x-A->x);
+ c.y=A->y + r * (B->y-A->y);
+ c.z=A->z + r * (B->z-A->z);
+
+ return lw_dist3d_pt_pt(p,&c,dl);
+}
+
+
+/**
+
+Compares incomming points and
+stores the points closest to each other
+or most far away from each other
+depending on dl->mode (max or min)
+*/
+int
+lw_dist3d_pt_pt(POINT3DZ *thep1, POINT3DZ *thep2,DISTPTS3D *dl)
+{
+ 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 );
+ double dx = thep2->x - thep1->x;
+ double dy = thep2->y - thep1->y;
+ double dz = thep2->z - thep1->z;
+ double dist = sqrt ( dx*dx + dy*dy + dz*dz);
+
+ if (((dl->distance - dist)*(dl->mode))>0) /*multiplication with mode to handle mindistance (mode=1) and maxdistance (mode = (-1)*/
+ {
+ dl->distance = dist;
+
+ 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*/
+ {
+ dl->p1 = *thep1;
+ dl->p2 = *thep2;
+ }
+ else
+ {
+ dl->p1 = *thep2;
+ dl->p2 = *thep1;
+ }
+ }
+ return LW_TRUE;
+}
+
+
+/**
+
+Finds all combinationes of segments between two pointarrays
+*/
+int
+lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS3D *dl)
+{
+ LWDEBUGF(2, "lw_dist3d_ptarray_ptarray called (points: %d-%d)",l1->npoints, l2->npoints);
+ int t,u;
+ POINT3DZ start, end;
+ POINT3DZ start2, end2;
+ int twist = dl->twisted;
+
+
+
+ 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*/
+ {
+ for (t=0; t<l1->npoints; t++) /*for each segment in L1 */
+ {
+ getPoint3dz_p(l1, t, &start);
+ for (u=0; u<l2->npoints; u++) /*for each segment in L2 */
+ {
+ getPoint3dz_p(l2, u, &start2);
+ lw_dist3d_pt_pt(&start,&start2,dl);
+ LWDEBUGF(4, "maxdist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
+ LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
+ t, u, dl->distance, dl->tolerance);
+ }
+ }
+ }
+ else
+ {
+ getPoint3dz_p(l1, 0, &start);
+ for (t=1; t<l1->npoints; t++) /*for each segment in L1 */
+ {
+ getPoint3dz_p(l1, t, &end);
+ getPoint3dz_p(l2, 0, &start2);
+ for (u=1; u<l2->npoints; u++) /*for each segment in L2 */
+ {
+ getPoint3dz_p(l2, u, &end2);
+ dl->twisted=twist;
+ lw_dist3d_seg_seg(&start, &end, &start2, &end2,dl);
+ LWDEBUGF(4, "mindist_ptarray_ptarray; seg %i * seg %i, dist = %g\n",t,u,dl->distance);
+ LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",
+ t, u, dl->distance, dl->tolerance);
+ if (dl->distance<=dl->tolerance && dl->mode == DIST_MIN) return LW_TRUE; /*just a check if the answer is already given*/
+ start2 = end2;
+ }
+ start = end;
+ }
+ }
+ return LW_TRUE;
+}
+
+/**
+
+Finds the two closest points on two linesegments
+*/
+int
+lw_dist3d_seg_seg(POINT3DZ *s1p1, POINT3DZ *s1p2, POINT3DZ *s2p1, POINT3DZ *s2p2, DISTPTS3D *dl)
+{
+ /*s1p1 and s1p2 are the same point */
+ if ( ( s1p1->x == s1p2->x) && (s1p1->y == s1p2->y) && (s1p1->z == s1p2->y) )
+ {
+ return lw_dist3d_pt_seg(s1p1,s2p1,s2p2,dl);
+ }
+ /*s2p1 and s2p2 are the same point */
+ if ( ( s2p1->x == s2p2->x) && (s2p1->y == s2p2->y) && (s2p1->z == s2p2->y) )
+ {
+ dl->twisted= ((dl->twisted) * (-1));
+ return lw_dist3d_pt_seg(s2p1,s1p1,s1p2,dl);
+ }
+
+
+/*
+ Here we use algorithm from softsurfer.com
+ that can be found here
+ http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
+*/
+
+
+ VECTOR3D v1, v2, vl;
+ 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*/
+ POINT3DZ p1, p2;
+
+ if (!get_3dvector_from_points(s1p1, s1p2, &v1))
+ return LW_FALSE;
+
+ if (!get_3dvector_from_points(s2p1, s2p2, &v2))
+ return LW_FALSE;
+
+ if (!get_3dvector_from_points(s2p1, s1p1, &vl))
+ return LW_FALSE;
+
+ double a = DOT(v1,v1);
+ double b = DOT(v1,v2);
+ double c = DOT(v2,v2);
+ double d = DOT(v1,vl);
+ double e = DOT(v2,vl);
+ double D = a*c - b*b;
+
+
+ if (D <0.000000001)
+ { /* the lines are almost parallel*/
+ 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.*/
+ if(b>c) /* use the largest denominator*/
+ {
+ s2k=d/b;
+ }
+ else
+ {
+ s2k =e/c;
+ }
+ }
+ else
+ {
+ s1k = (b*e - c*d) / D;
+ s2k = (a*e - b*d) / D;
+ }
+
+ /* 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*/
+ if(s1k<0.0||s1k>1.0||s2k<0.0||s2k>1.0)
+ {
+ if(s1k<0.0)
+ {
+
+ if (!lw_dist3d_pt_seg(s1p1, s2p1, s2p2, dl))
+ {
+ return LW_FALSE;
+ }
+ }
+ if(s1k>1.0)
+ {
+
+ if (!lw_dist3d_pt_seg(s1p2, s2p1, s2p2, dl))
+ {
+ return LW_FALSE;
+ }
+ }
+ if(s2k<0.0)
+ {
+ dl->twisted= ((dl->twisted) * (-1));
+ if (!lw_dist3d_pt_seg(s2p1, s1p1, s1p2, dl))
+ {
+ return LW_FALSE;
+ }
+ }
+ if(s2k>1.0)
+ {
+ dl->twisted= ((dl->twisted) * (-1));
+ if (!lw_dist3d_pt_seg(s2p2, s1p1, s1p2, dl))
+ {
+ return LW_FALSE;
+ }
+ }
+ }
+ else
+ {/*Find the closest point on the edges of both segments*/
+ p1.x=s1p1->x+s1k*(s1p2->x-s1p1->x);
+ p1.y=s1p1->y+s1k*(s1p2->y-s1p1->y);
+ p1.z=s1p1->z+s1k*(s1p2->z-s1p1->z);
+
+ p2.x=s2p1->x+s2k*(s2p2->x-s2p1->x);
+ p2.y=s2p1->y+s2k*(s2p2->y-s2p1->y);
+ p2.z=s2p1->z+s2k*(s2p2->z-s2p1->z);
+
+ if (!lw_dist3d_pt_pt(&p1,&p2,dl))/* Send the closest points to point-point calculation*/
+ {
+ return LW_FALSE;
+ }
+ }
+ return LW_TRUE;
+}
+
+/**
+
+Checking if the point projected on the plane of the polygon actually is inside that polygon.
+If so the mindistance is between that projected point and our original point.
+If not we check from original point to the bounadary.
+If the projected point is inside a hole of the polygon we check the distance to the boudary of that hole.
+*/
+int
+lw_dist3d_pt_poly(POINT3DZ *p, LWPOLY *poly, PLANE3D *plane,POINT3DZ *projp, DISTPTS3D *dl)
+{
+ int i;
+
+ LWDEBUG(2, "lw_dist3d_point_poly called");
+
+
+ if(pt_in_ring_3d(projp, poly->rings[0], plane))
+ {
+ for (i=1; i<poly->nrings; i++)
+ {
+ /* Inside a hole. Distance = pt -> ring */
+ if ( pt_in_ring_3d(projp, poly->rings[i], plane ))
+ {
+ LWDEBUG(3, " inside an hole");
+ return lw_dist3d_pt_ptarray(p, poly->rings[i], dl);
+ }
+ }
+
+ 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*/
+ }
+ else
+ {
+ 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*/
+ }
+
+ return LW_TRUE;
+
+}
+
+/**
+
+Computes pointarray to polygon distance
+*/
+int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3D *dl)
+{
+
+
+ int i,j,k;
+ double f, s1, s2;
+ VECTOR3D projp1_projp2;
+ POINT3DZ p1, p2,projp1, projp2, intersectionp;
+
+ getPoint3dz_p(pa, 0, &p1);
+
+ s1=project_point_on_plane(&p1, plane, &projp1); /*the sign of s1 tells us on which side of the plane the point is. */
+ lw_dist3d_pt_poly(&p1, poly, plane,&projp1, dl);
+
+ for (i=1;i<pa->npoints;i++)
+ {
+ getPoint3dz_p(pa, i, &p2);
+ s2=project_point_on_plane(&p2, plane, &projp2);
+ lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl);
+
+ /*If s1and s2 has different signs that means they are on different sides of the plane of the polygon.
+ That means that the edge between the points crosses the plane and might intersect with the polygon*/
+ if((s1*s2)<=0)
+ {
+ f=fabs(s1)/(fabs(s1)+fabs(s2)); /*The size of s1 and s2 is the distance from the point to the plane.*/
+ get_3dvector_from_points(&projp1, &projp2,&projp1_projp2);
+
+ /*get the point where the line segment crosses the plane*/
+ intersectionp.x=projp1.x+f*projp1_projp2.x;
+ intersectionp.y=projp1.y+f*projp1_projp2.y;
+ intersectionp.z=projp1.z+f*projp1_projp2.z;
+
+ int intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/
+
+ if(pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/
+ {
+ for (k=1;k<poly->nrings; k++)
+ {
+ /* Inside a hole, so no intersection with the polygon*/
+ if ( pt_in_ring_3d(&intersectionp, poly->rings[k], plane ))
+ {
+ intersects=LW_FALSE;
+ break;
+ }
+ }
+ if(intersects)
+ {
+ dl->distance=0.0;
+ dl->p1.x=intersectionp.x;
+ dl->p1.y=intersectionp.y;
+ dl->p1.z=intersectionp.z;
+
+ dl->p2.x=intersectionp.x;
+ dl->p2.y=intersectionp.y;
+ dl->p2.z=intersectionp.z;
+ return LW_TRUE;
+
+ }
+ }
+ }
+
+ projp2=projp1;
+ s2=s1;
+ p2=p1;
+ }
+
+ /*check or pointarray against boundary and inner boundaries of the polygon*/
+ for (j=0;j<poly->nrings;j++)
+ {
+ lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);
+ }
+
+return LW_TRUE;
+}
+
+
+/**
+
+Here we define the plane of a polygon (boundary pointarray of a polygon)
+the plane is stored as a pont in plane (plane.pop) and a normal vector (plane.pv)
+*/
+int
+define_plane(POINTARRAY *pa, PLANE3D *pl)
+{
+ int i,j, numberofvectors, pointsinslice;
+ POINT3DZ p, p1, p2;
+
+ double sumx=0;
+ double sumy=0;
+ double sumz=0;
+ double vl; /*vector length*/
+
+ VECTOR3D v1, v2, v;
+
+ if((pa->npoints-1)==3) /*Triangle is special case*/
+ {
+ pointsinslice=1;
+ }
+ else
+ {
+ pointsinslice=(int) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/
+ }
+
+ /*find the avg point*/
+ for (i=0;i<(pa->npoints-1);i++)
+ {
+ getPoint3dz_p(pa, i, &p);
+ sumx+=p.x;
+ sumy+=p.y;
+ sumz+=p.z;
+ }
+ pl->pop.x=(sumx/(pa->npoints-1));
+ pl->pop.y=(sumy/(pa->npoints-1));
+ pl->pop.z=(sumz/(pa->npoints-1));
+
+ sumx=0;
+ sumy=0;
+ sumz=0;
+ numberofvectors= floor((pa->npoints-1)/pointsinslice); /*the number of vectors we try can be 3, 4 or 5*/
+
+ getPoint3dz_p(pa, 0, &p1);
+ for (j=pointsinslice;j<pa->npoints;j+=pointsinslice)
+ {
+ getPoint3dz_p(pa, j, &p2);
+
+ if (!get_3dvector_from_points(&(pl->pop), &p1, &v1) || !get_3dvector_from_points(&(pl->pop), &p2, &v2))
+ return LW_FALSE;
+ /*perpendicular vector is cross product of v1 and v2*/
+ if (!get_3dcross_product(&v1,&v2, &v))
+ return LW_FALSE;
+ vl=VECTORLENGTH(v);
+ sumx+=(v.x/vl);
+ sumy+=(v.y/vl);
+ sumz+=(v.z/vl);
+ p1=p2;
+ }
+ pl->pv.x=(sumx/numberofvectors);
+ pl->pv.y=(sumy/numberofvectors);
+ pl->pv.z=(sumz/numberofvectors);
+
+ return 1;
+}
+
+/**
+
+Finds a point on a plane from where the original point is perpendicular to the plane
+*/
+double
+project_point_on_plane(POINT3DZ *p, PLANE3D *pl, POINT3DZ *p0)
+{
+/*In our plane definition we have a point on the plane and a normal vektor (pl.pv), perpendicular to the plane
+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.
+So, we already have a direction from p to find p0, but we don't know the distance.
+*/
+
+ VECTOR3D v1;
+ double f;
+
+ if (!get_3dvector_from_points(&(pl->pop), p, &v1))
+ return LW_FALSE;
+
+ f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv));
+
+ p0->x=p->x+pl->pv.x*f;
+ p0->y=p->y+pl->pv.y*f;
+ p0->z=p->z+pl->pv.z*f;
+
+ return f;
+}
+
+
+
+
+/**
+ * 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
+ *
+ * Our polygons have first and last point the same,
+ *
+* The difference in 3D variant is that we exclude the dimension that faces the plane least.
+* That is the dimension with the highest number in pv
+ */
+int
+pt_in_ring_3d(const POINT3DZ *p, const POINTARRAY *ring,PLANE3D *plane)
+{
+
+ int cn = 0; /* the crossing number counter */
+ int i;
+ POINT3DZ v1, v2;
+
+ POINT3DZ first, last;
+
+ getPoint3dz_p(ring, 0, &first);
+ getPoint3dz_p(ring, ring->npoints-1, &last);
+ if ( memcmp(&first, &last, sizeof(POINT3DZ)) )
+ {
+ lwerror("pt_in_ring_3d: V[n] != V[0] (%g %g %g!= %g %g %g)",
+ first.x, first.y, first.z, last.x, last.y, last.z);
+ return LW_FALSE;
+ }
+
+ LWDEBUGF(2, "pt_in_ring_3d called with point: %g %g %g", p->x, p->y, p->z);
+ /* printPA(ring); */
+
+ /* loop through all edges of the polygon */
+ getPoint3dz_p(ring, 0, &v1);
+
+
+ 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*/
+ {
+ for (i=0; i<ring->npoints-1; i++)
+ {
+ double vt;
+ getPoint3dz_p(ring, i+1, &v2);
+
+ /* edge from vertex i to vertex i+1 */
+ if
+ (
+ /* an upward crossing */
+ ((v1.y <= p->y) && (v2.y > p->y))
+ /* a downward crossing */
+ || ((v1.y > p->y) && (v2.y <= p->y))
+ )
+ {
+
+ vt = (double)(p->y - v1.y) / (v2.y - v1.y);
+
+ /* P.x <intersect */
+ if (p->x < v1.x + vt * (v2.x - v1.x))
+ {
+ /* a valid crossing of y=p.y right of p.x */
+ ++cn;
+ }
+ }
+ v1 = v2;
+ }
+ }
+ 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*/
+ {
+ for (i=0; i<ring->npoints-1; i++)
+ {
+ double vt;
+ getPoint3dz_p(ring, i+1, &v2);
+
+ /* edge from vertex i to vertex i+1 */
+ if
+ (
+ /* an upward crossing */
+ ((v1.z <= p->z) && (v2.z > p->z))
+ /* a downward crossing */
+ || ((v1.z > p->z) && (v2.z <= p->z))
+ )
+ {
+
+ vt = (double)(p->z - v1.z) / (v2.z - v1.z);
+
+ /* P.x <intersect */
+ if (p->x < v1.x + vt * (v2.x - v1.x))
+ {
+ /* a valid crossing of y=p.y right of p.x */
+ ++cn;
+ }
+ }
+ v1 = v2;
+ }
+ }
+ else /*Hopefully we only have the cases where x part of the normal vector is largest left*/
+ {
+ for (i=0; i<ring->npoints-1; i++)
+ {
+ double vt;
+ getPoint3dz_p(ring, i+1, &v2);
+
+ /* edge from vertex i to vertex i+1 */
+ if
+ (
+ /* an upward crossing */
+ ((v1.z <= p->z) && (v2.z > p->z))
+ /* a downward crossing */
+ || ((v1.z > p->z) && (v2.z <= p->z))
+ )
+ {
+
+ vt = (double)(p->z - v1.z) / (v2.z - v1.z);
+
+ /* P.x <intersect */
+ if (p->y < v1.y + vt * (v2.y - v1.y))
+ {
+ /* a valid crossing of y=p.y right of p.x */
+ ++cn;
+ }
+ }
+ v1 = v2;
+ }
+ }
+ LWDEBUGF(3, "pt_in_ring_3d returning %d", cn&1);
+
+ return (cn&1); /* 0 if even (out), and 1 if odd (in) */
+}
+
+
+
+/*------------------------------------------------------------------------------------------------------------
+End of Brute force functions
+--------------------------------------------------------------------------------------------------------------*/
+
+
+
projects / postgis / commitdiff