[postgis-tickets] r17053 - Fix undefined behaviour in ST_3DDistance
Raul
raul at rmr.ninja
Thu Nov 22 05:44:47 PST 2018
Author: algunenano
Date: 2018-11-22 05:44:47 -0800 (Thu, 22 Nov 2018)
New Revision: 17053
Modified:
trunk/NEWS
trunk/liblwgeom/cunit/cu_measures.c
trunk/liblwgeom/measures.h
trunk/liblwgeom/measures3d.c
trunk/regress/core/measures.sql
trunk/regress/core/measures_expected
Log:
Fix undefined behaviour in ST_3DDistance
Closes #4246
Closes https://github.com/postgis/postgis/pull/338
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2018-11-22 13:43:39 UTC (rev 17052)
+++ trunk/NEWS 2018-11-22 13:44:47 UTC (rev 17053)
@@ -41,6 +41,7 @@
- #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
Modified: trunk/liblwgeom/cunit/cu_measures.c
===================================================================
--- trunk/liblwgeom/cunit/cu_measures.c 2018-11-22 13:43:39 UTC (rev 17052)
+++ trunk/liblwgeom/cunit/cu_measures.c 2018-11-22 13:44:47 UTC (rev 17053)
@@ -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 @@
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 @@
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_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 @@
{
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);
Modified: trunk/liblwgeom/measures.h
===================================================================
--- trunk/liblwgeom/measures.h 2018-11-22 13:43:39 UTC (rev 17052)
+++ trunk/liblwgeom/measures.h 2018-11-22 13:44:47 UTC (rev 17053)
@@ -34,6 +34,9 @@
*
**********************************************************************/
+#ifndef _MEASURES_H
+#define _MEASURES_H 1
+
#include "liblwgeom_internal.h"
/* for the measure functions*/
@@ -124,4 +127,4 @@
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
Modified: trunk/liblwgeom/measures3d.c
===================================================================
--- trunk/liblwgeom/measures3d.c 2018-11-22 13:43:39 UTC (rev 17052)
+++ trunk/liblwgeom/measures3d.c 2018-11-22 13:44:47 UTC (rev 17053)
@@ -48,7 +48,7 @@
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 @@
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 @@
/*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 @@
}
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,7 +687,8 @@
*/
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)
{
@@ -688,21 +695,37 @@
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 @@
*/
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 @@
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;i<pa->npoints;i++)
{
@@ -1068,10 +1094,14 @@
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 @@
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;
+ const uint32_t POL_BREAKS = 3;
+ uint32_t unique_points = pa->npoints - 1;
+ uint32_t i;
- 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*/
+ /* 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++)
{
- pointsinslice=1;
- }
- else
- {
- pointsinslice=(uint32_t) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/
- }
-
- /*find the avg point*/
- for (i=0;i<(pa->npoints-1);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;j<pa->npoints;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,14 +1239,22 @@
double f;
if (!get_3dvector_from_points(&(pl->pop), p, &v1))
- return LW_FALSE;
+ return 0.0;
- f=-(DOT(pl->pv,v1)/DOT(pl->pv,pl->pv));
+ f = DOT(pl->pv, v1);
+ if (FP_IS_ZERO(f))
+ {
+ /* Point is in the plane */
+ *p0 = *p;
+ return 0;
+ }
- p0->x=p->x+pl->pv.x*f;
- p0->y=p->y+pl->pv.y*f;
- p0->z=p->z+pl->pv.z*f;
+ 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;
+
return f;
}
Modified: trunk/regress/core/measures.sql
===================================================================
--- trunk/regress/core/measures.sql 2018-11-22 13:43:39 UTC (rev 17052)
+++ trunk/regress/core/measures.sql 2018-11-22 13:44:47 UTC (rev 17053)
@@ -221,6 +221,10 @@
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));
Modified: trunk/regress/core/measures_expected
===================================================================
--- trunk/regress/core/measures_expected 2018-11-22 13:43:39 UTC (rev 17052)
+++ trunk/regress/core/measures_expected 2018-11-22 13:44:47 UTC (rev 17053)
@@ -35,6 +35,7 @@
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"
More information about the postgis-tickets
mailing list