[postgis-tickets] r17050 - Fix undefined behaviour in ST_3DDistance
Raul
raul at rmr.ninja
Thu Nov 22 05:40:39 PST 2018
Author: algunenano
Date: 2018-11-22 05:40:38 -0800 (Thu, 22 Nov 2018)
New Revision: 17050
Modified:
branches/2.3/NEWS
branches/2.3/liblwgeom/cunit/cu_measures.c
branches/2.3/liblwgeom/measures.h
branches/2.3/liblwgeom/measures3d.c
branches/2.3/regress/measures.sql
branches/2.3/regress/measures_expected
Log:
Fix undefined behaviour in ST_3DDistance
References #4246
Modified: branches/2.3/NEWS
===================================================================
--- branches/2.3/NEWS 2018-11-22 07:33:52 UTC (rev 17049)
+++ branches/2.3/NEWS 2018-11-22 13:40:38 UTC (rev 17050)
@@ -18,6 +18,7 @@
- #4326, Allocate enough memory in gidx_to_string (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.3.7
2018/04/06
Modified: branches/2.3/liblwgeom/cunit/cu_measures.c
===================================================================
--- branches/2.3/liblwgeom/cunit/cu_measures.c 2018-11-22 07:33:52 UTC (rev 17049)
+++ branches/2.3/liblwgeom/cunit/cu_measures.c 2018-11-22 13:40:38 UTC (rev 17050)
@@ -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 void test_rect_tree_contains_point(void)
{
LWPOLY *poly;
@@ -1279,6 +1324,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: branches/2.3/liblwgeom/measures.h
===================================================================
--- branches/2.3/liblwgeom/measures.h 2018-11-22 07:33:52 UTC (rev 17049)
+++ branches/2.3/liblwgeom/measures.h 2018-11-22 13:40:38 UTC (rev 17050)
@@ -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: branches/2.3/liblwgeom/measures3d.c
===================================================================
--- branches/2.3/liblwgeom/measures3d.c 2018-11-22 07:33:52 UTC (rev 17049)
+++ branches/2.3/liblwgeom/measures3d.c 2018-11-22 13:40:38 UTC (rev 17050)
@@ -47,8 +47,8 @@
v->x=p2->x-p1->x;
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));
}
@@ -635,7 +635,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);
@@ -673,7 +676,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);
}
@@ -683,29 +689,46 @@
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;
-
- /*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*/
+
+ 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 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);
+
+ 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, &plane1, dl);
}
/**
@@ -1053,18 +1076,21 @@
*/
int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly,PLANE3D *plane, DISTPTS3D *dl)
{
-
-
- int i,j,k;
+ uint32_t 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);
+ if ((s1 == 0.0) && (dl->distance < dl->tolerance))
+ {
+ return LW_TRUE;
+ }
- 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++)
{
int intersects;
@@ -1071,6 +1097,10 @@
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*/
@@ -1128,71 +1158,73 @@
}
-/**
-
-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)
-*/
+/* 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)
{
- int 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;
+ POINT3DZ p;
+ getPoint3dz_p(pa, i, &p);
+ pl->pop.x += p.x;
+ pl->pop.y += p.y;
+ pl->pop.z += p.z;
}
- else
+
+ pl->pop.x /= unique_points;
+ pl->pop.y /= unique_points;
+ pl->pop.z /= unique_points;
+
+ pl->pv.x = 0.0;
+ pl->pv.y = 0.0;
+ pl->pv.z = 0.0;
+ for (i = 0; i < POL_BREAKS; i++)
{
- pointsinslice=(int) floor((pa->npoints-1)/4); /*divide the pointarray into 4 slices*/
+ 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 (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;
+ }
}
-
- /*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;
+
+ return (!FP_IS_ZERO(pl->pv.x) || !FP_IS_ZERO(pl->pv.y) || !FP_IS_ZERO(pl->pv.z));
}
/**
@@ -1202,24 +1234,32 @@
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.
+/*In our plane definition we have a point on the plane and a normal vector (pl.pv), perpendicular to the plane
+this vector will be parallel 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;
+ return 0.0;
+
+ f = DOT(pl->pv, v1);
+ if (FP_IS_ZERO(f))
+ {
+ /* Point is in the plane */
+ *p0 = *p;
+ return 0;
+ }
+
+ 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: branches/2.3/regress/measures.sql
===================================================================
--- branches/2.3/regress/measures.sql 2018-11-22 07:33:52 UTC (rev 17049)
+++ branches/2.3/regress/measures.sql 2018-11-22 13:40:38 UTC (rev 17050)
@@ -239,6 +239,10 @@
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: branches/2.3/regress/measures_expected
===================================================================
--- branches/2.3/regress/measures_expected 2018-11-22 07:33:52 UTC (rev 17049)
+++ branches/2.3/regress/measures_expected 2018-11-22 13:40:38 UTC (rev 17050)
@@ -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