From 03ed7d3458ebe7ef30182aa0fb9618f9df625cb2 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Ra=C3=BAl=20Mar=C3=ADn=20Rodr=C3=ADguez?= Date: Thu, 22 Nov 2018 13:44:47 +0000 Subject: [PATCH] Fix undefined behaviour in ST_3DDistance Closes #4246 Closes https://github.com/postgis/postgis/pull/338 git-svn-id: http://svn.osgeo.org/postgis/trunk@17053 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 1 + liblwgeom/cunit/cu_measures.c | 52 ++++++++- liblwgeom/measures.h | 5 +- liblwgeom/measures3d.c | 187 ++++++++++++++++++++------------- regress/core/measures.sql | 4 + regress/core/measures_expected | 1 + 6 files changed, 172 insertions(+), 78 deletions(-) diff --git a/NEWS b/NEWS index c31ec2ddb..d80306fcb 100644 --- a/NEWS +++ b/NEWS @@ -41,6 +41,7 @@ PostGIS 3.0.0 - #4233, Fix undefined behaviour in gserialized_spgist_picksplit_nd (Raúl Marín) - #4247, Avoid undefined behaviour in next_float functions (Raúl Marín) - #4249, Fix undefined behaviour in raster intersection (Raúl Marín) + - #4246, Fix undefined behaviour in ST_3DDistance (Raúl Marín) PostGIS 2.5.0 2018/09/23 diff --git a/liblwgeom/cunit/cu_measures.c b/liblwgeom/cunit/cu_measures.c index d94b15487..2e93e0012 100644 --- a/liblwgeom/cunit/cu_measures.c +++ b/liblwgeom/cunit/cu_measures.c @@ -19,6 +19,7 @@ #include "liblwgeom_internal.h" #include "cu_tester.h" #include "measures.h" +#include "measures3d.h" #include "lwtree.h" static LWGEOM* lwgeom_from_text(const char *str) @@ -29,9 +30,17 @@ static LWGEOM* lwgeom_from_text(const char *str) return r.geom; } -#define DIST2DTEST(str1, str2, res) do_test_mindistance2d_tolerance(str1, str2, res, __LINE__) +#define DIST2DTEST(str1, str2, res) \ + do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance2d_tolerance) +#define DIST3DTEST(str1, str2, res) \ + do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance3d_tolerance) -static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line) +static void +do_test_mindistance_tolerance(char *in1, + char *in2, + double expected_res, + int line, + double (*distancef)(const LWGEOM *, const LWGEOM *, double)) { LWGEOM *lw1; LWGEOM *lw2; @@ -53,7 +62,7 @@ static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expecte exit(1); } - distance = lwgeom_mindistance2d_tolerance(lw1, lw2, 0.0); + distance = distancef(lw1, lw2, 0.0); lwgeom_free(lw1); lwgeom_free(lw2); @@ -193,6 +202,42 @@ static void test_mindistance2d_tolerance(void) } +static void +test_mindistance3d_tolerance(void) +{ + /* 2D [Z=0] should work just the same */ + DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5); + DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5); + DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0); + DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))", 5.0); + DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0))))", 5.0); + DIST3DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0); + DIST3DTEST("GEOMETRYCOLLECTION(POINT(0 0 0))", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0); + DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0 0)))", + "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))", + 5.0); + DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0 0)))", + "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4 0)))", + 5.0); + DIST3DTEST("LINESTRING(-2 0 0, -0.2 0 0)", "POINT(-2 0 0)", 0); + DIST3DTEST("LINESTRING(-0.2 0 0, -2 0 0)", "POINT(-2 0 0)", 0); + DIST3DTEST("LINESTRING(-1e-8 0 0, -0.2 0 0)", "POINT(-1e-8 0 0)", 0); + DIST3DTEST("LINESTRING(-0.2 0 0, -1e-8 0 0)", "POINT(-1e-8 0 0)", 0); + + /* Tests around intersections */ + DIST3DTEST("LINESTRING(1 0 0 , 2 0 0)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 1.0); + DIST3DTEST("LINESTRING(1 1 1 , 2 1 0)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 0.0); + DIST3DTEST("LINESTRING(1 1 1 , 2 1 1)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 1.0); + /* Line in polygon */ + DIST3DTEST("LINESTRING(1 1 1 , 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 1, 0 0 0))", 0.0); + + /* Line has a point in the same plane as the polygon but isn't the closest*/ + DIST3DTEST("LINESTRING(-10000 -10000 0, 0 0 1)", "POLYGON((0 0 0, 1 0 0, 1 1 0, 0 1 0, 0 0 0))", 1); + + /* This is an invalid polygon since it defines just a line */ + DIST3DTEST("LINESTRING(1 1 1 , 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", 0); +} + static int tree_pt(RECT_NODE *tree, double x, double y) { POINT2D pt; @@ -1404,6 +1449,7 @@ void measures_suite_setup(void) { CU_pSuite suite = CU_add_suite("measures", NULL, NULL); PG_ADD_TEST(suite, test_mindistance2d_tolerance); + PG_ADD_TEST(suite, test_mindistance3d_tolerance); PG_ADD_TEST(suite, test_rect_tree_contains_point); PG_ADD_TEST(suite, test_rect_tree_intersects_tree); PG_ADD_TEST(suite, test_lwgeom_segmentize2d); diff --git a/liblwgeom/measures.h b/liblwgeom/measures.h index 970a7a249..7f723ef06 100644 --- a/liblwgeom/measures.h +++ b/liblwgeom/measures.h @@ -34,6 +34,9 @@ * **********************************************************************/ +#ifndef _MEASURES_H +#define _MEASURES_H 1 + #include "liblwgeom_internal.h" /* for the measure functions*/ @@ -124,4 +127,4 @@ double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3); LWGEOM* lw_dist2d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode); LWGEOM* lw_dist2d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode); - +#endif /* !defined _MEASURES_H */ \ No newline at end of file diff --git a/liblwgeom/measures3d.c b/liblwgeom/measures3d.c index 31e494478..b16e82542 100644 --- a/liblwgeom/measures3d.c +++ b/liblwgeom/measures3d.c @@ -48,7 +48,7 @@ get_3dvector_from_points(POINT3DZ *p1,POINT3DZ *p2, VECTOR3D *v) v->y=p2->y-p1->y; v->z=p2->z-p1->z; - return LW_TRUE; + return (!FP_IS_ZERO(v->x) || !FP_IS_ZERO(v->y) || !FP_IS_ZERO(v->z)); } static inline int @@ -58,7 +58,7 @@ get_3dcross_product(VECTOR3D *v1,VECTOR3D *v2, VECTOR3D *v) v->y=(v1->z*v2->x)-(v1->x*v2->z); v->z=(v1->x*v2->y)-(v1->y*v2->x); - return LW_TRUE; + return (!FP_IS_ZERO(v->x) || !FP_IS_ZERO(v->y) || !FP_IS_ZERO(v->z)); } @@ -632,7 +632,10 @@ lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *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; + { + /* Polygon does not define a plane: Return distance point -> line */ + return lw_dist3d_pt_ptarray(&p, poly->rings[0], dl); + } /*get our point projected on the plane of the polygon*/ project_point_on_plane(&p, &plane, &projp); @@ -670,7 +673,10 @@ int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl) } if(!define_plane(poly->rings[0], &plane)) - return LW_FALSE; + { + /* Polygon does not define a plane: Return distance line to line */ + return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl); + } return lw_dist3d_ptarray_poly(line->points, poly,&plane, dl); } @@ -681,28 +687,45 @@ polygon to polygon calculation */ int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl) { - PLANE3D plane; + PLANE3D plane1, plane2; + int planedef1, planedef2; LWDEBUG(2, "lw_dist3d_poly_poly is called"); 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; + planedef1 = define_plane(poly1->rings[0], &plane1); + planedef2 = define_plane(poly2->rings[0], &plane2); + + if (!planedef1 || !planedef2) + { + if (!planedef1 && !planedef2) + { + /* Neither polygon define a plane: Return distance line to line */ + return lw_dist3d_ptarray_ptarray(poly1->rings[0], poly2->rings[0], dl); + } + + if (!planedef1) + { + /* Only poly2 defines a plane: Return distance from line (poly1) to poly2 */ + return lw_dist3d_ptarray_poly(poly1->rings[0], poly2, &plane2, dl); + } + + /* Only poly1 defines a plane: Return distance from line (poly2) to poly1 */ + return lw_dist3d_ptarray_poly(poly2->rings[0], poly1, &plane1, dl); + } /*What we do here is to compare the boundary of one polygon with the other polygon and then take the second boundary comparing with the first polygon*/ dl->twisted=1; - if(!lw_dist3d_ptarray_poly(poly1->rings[0], poly2,&plane, dl)) + if (!lw_dist3d_ptarray_poly(poly1->rings[0], poly2, &plane2, dl)) return LW_FALSE; - if(dl->distance==0.0) /*Just check if the answer already is given*/ + if (dl->distance < dl->tolerance) /*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 switch the order of geometries we switch "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); + return lw_dist3d_ptarray_poly(poly2->rings[0], poly1, &plane1, dl); } /** @@ -1050,8 +1073,6 @@ Computes pointarray to polygon distance */ int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3D *dl) { - - uint32_t i,j,k; double f, s1, s2; VECTOR3D projp1_projp2; @@ -1059,8 +1080,13 @@ int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3 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. */ + 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); + if ((s1 == 0.0) && (dl->distance < dl->tolerance)) + { + return LW_TRUE; + } for (i=1;inpoints;i++) { @@ -1068,10 +1094,14 @@ int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3 getPoint3dz_p(pa, i, &p2); s2=project_point_on_plane(&p2, plane, &projp2); lw_dist3d_pt_poly(&p2, poly, plane,&projp2, dl); + if ((s2 == 0.0) && (dl->distance < dl->tolerance)) + { + return LW_TRUE; + } /*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) + 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); @@ -1124,72 +1154,73 @@ int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3 return LW_TRUE; } - -/** - -Here we define the plane of a polygon (boundary pointarray of a polygon) -the plane is stored as a point in plane (plane.pop) and a normal vector (plane.pv) -*/ +/* Here we define the plane of a polygon (boundary pointarray of a polygon) + * the plane is stored as a point in plane (plane.pop) and a normal vector (plane.pv) + * + * We are calculating the average point and using it to break the polygon into + * POL_BREAKS (3) parts. Then we calculate the normal of those 3 vectors and + * use its average to define the normal of the plane.This is done to minimize + * the error introduced by FP precision + * We use [3] breaks to contemplate the special case of 3d triangles + */ int define_plane(POINTARRAY *pa, PLANE3D *pl) { - uint32_t 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=(uint32_t) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/ - } + const uint32_t POL_BREAKS = 3; + uint32_t unique_points = pa->npoints - 1; + uint32_t i; - /*find the avg point*/ - for (i=0;i<(pa->npoints-1);i++) + /* Calculate the average point */ + pl->pop.x = 0.0; + pl->pop.y = 0.0; + pl->pop.z = 0.0; + for (i = 0; i < unique_points; i++) { + POINT3DZ p; getPoint3dz_p(pa, i, &p); - sumx+=p.x; - sumy+=p.y; - sumz+=p.z; + pl->pop.x += p.x; + pl->pop.y += p.y; + pl->pop.z += 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*/ + pl->pop.x /= unique_points; + pl->pop.y /= unique_points; + pl->pop.z /= unique_points; - getPoint3dz_p(pa, 0, &p1); - for (j=pointsinslice;jnpoints;j+=pointsinslice) + pl->pv.x = 0.0; + pl->pv.y = 0.0; + pl->pv.z = 0.0; + for (i = 0; i < POL_BREAKS; i++) { - getPoint3dz_p(pa, j, &p2); + POINT3DZ point1, point2; + VECTOR3D v1, v2, vp; + uint32_t n1, n2; + n1 = i * unique_points / POL_BREAKS; + n2 = n1 + unique_points / POL_BREAKS; /* May overflow back to the first / last point */ - 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; + if (n1 == n2) + continue; + + getPoint3dz_p(pa, n1, &point1); + if (!get_3dvector_from_points(&pl->pop, &point1, &v1)) + continue; + + getPoint3dz_p(pa, n2, &point2); + if (!get_3dvector_from_points(&pl->pop, &point2, &v2)) + continue; + + if (get_3dcross_product(&v1, &v2, &vp)) + { + /* If the cross product isn't zero these 3 points aren't collinear + * We divide by its lengthsq to normalize the additions */ + double vl = vp.x * vp.x + vp.y * vp.y + vp.z * vp.z; + pl->pv.x += vp.x / vl; + pl->pv.y += vp.y / vl; + pl->pv.z += vp.z / vl; + } } - pl->pv.x=(sumx/numberofvectors); - pl->pv.y=(sumy/numberofvectors); - pl->pv.z=(sumz/numberofvectors); - return 1; + return (!FP_IS_ZERO(pl->pv.x) || !FP_IS_ZERO(pl->pv.y) || !FP_IS_ZERO(pl->pv.z)); } /** @@ -1208,13 +1239,21 @@ So, we already have a direction from p to find p0, but we don't know the distanc double f; if (!get_3dvector_from_points(&(pl->pop), p, &v1)) - return LW_FALSE; + return 0.0; + + f = DOT(pl->pv, v1); + if (FP_IS_ZERO(f)) + { + /* Point is in the plane */ + *p0 = *p; + return 0; + } - f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv)); + f = -f / 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; + 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; } diff --git a/regress/core/measures.sql b/regress/core/measures.sql index a093cb421..93964e3d6 100644 --- a/regress/core/measures.sql +++ b/regress/core/measures.sql @@ -221,6 +221,10 @@ SELECT '3dDistancetest6', ST_3DDistance(a,b) FROM ( SELECT 'LINESTRING(1 1 1 , 2 2 2)'::geometry as a, 'POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))'::geometry as b) as foo; +SELECT '3dDistancetest7', + ST_3DDistance(a,b) FROM ( + SELECT 'LINESTRING(1 1 1 , 2 2 2)'::geometry as a, 'POLYGON((0 0 0, 2 2 2, 3 3 1, 0 0 0))'::geometry as b) as foo; + -- 3D mixed dimmentionality #2034 --closestpoint with 2d as first point and 3d as second select st_astext(st_3dclosestpoint('linestring(0 0,1 1,2 0)'::geometry, 'linestring(0 2 3, 3 2 3)'::geometry)); diff --git a/regress/core/measures_expected b/regress/core/measures_expected index 65bb8b237..e11c9734d 100644 --- a/regress/core/measures_expected +++ b/regress/core/measures_expected @@ -35,6 +35,7 @@ distancepoly6|0|32.5269119345812|LINESTRING(2 2,2 2)|LINESTRING(2 2,2 2)|LINESTR 3dDistancetest4|2|10.0498756211209|t|f|LINESTRING(1 1 3,1 1 1)|POINT(1 1 3)|LINESTRING(5 7 8,1 1 1) 3dDistancetest5|2|10|t|f|LINESTRING(5 0 5,5 2 5)|POINT(5 0 5)|LINESTRING(11 0 5,5 0 13) 3dDistancetest6|0 +3dDistancetest7|0 NOTICE: One or both of the geometries is missing z-value. The unknown z-value will be regarded as "any value" POINT Z (1 1 3) NOTICE: One or both of the geometries is missing z-value. The unknown z-value will be regarded as "any value" -- 2.40.0