Commit 3abc9a25 authored by Darafei Praliaskouski's avatar Darafei Praliaskouski

Make lwgeom_mindist3d solid aware. Add TIN support.

Thanks Hugo Mercier for suggesting test cases.

Closes #4278
Closes https://github.com/postgis/postgis/pull/373


git-svn-id: http://svn.osgeo.org/postgis/trunk@17260 b70326c6-7e19-0410-871a-916f4a2858ee
parent 646cbc01
Pipeline #47770447 passed with stage
in 29 minutes and 24 seconds
......@@ -71,7 +71,8 @@ PostGIS 3.0.0
- #4313, #4307, PostgreSQL 12 compatibility (Laurenz Albe, Raúl Marín)
- #4299, #4304, ST_GeneratePoints is now VOLATILE. IMMUTABLE version with
seed parameter added. (Mike Taves)
- #4278, ST_3DDistance and ST_3DIntersects now support Solid TIN and Solid
POLYHEDRALSURFACE
PostGIS 2.5.0
2018/09/23
......
......@@ -778,6 +778,20 @@ test_lwtriangle_clip(void)
lwfree(ewkt);
lwcollection_free(c);
lwgeom_free(g);
g = lwgeom_from_wkt(
"TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(g, 'Z', 0.0, DBL_MAX, 0);
ewkt = lwgeom_to_ewkt((LWGEOM *)c);
// printf("c = %s\n", ewkt);
ASSERT_STRING_EQUAL(
ewkt,
"TIN(((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 0 0,-1 -1 2)),((-1 -1 2,-1 0 0,-1 -1 0,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,0 2 0,-1 2 2)),((-1 2 2,0 2 0,-1 2 0,-1 2 2)),((-1 0 0,-1 2 2,-1 2 0,-1 0 0)),((-1 -1 2,-1 -1 0,1 -1 0,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 1 0,2 2 2)),((2 2 2,2 1 0,2 2 0,2 2 2)),((0 2 0,2 2 2,2 2 0,0 2 0)),((-1 -1 2,1 -1 0,2 -1 0,-1 -1 2)),((-1 -1 2,2 -1 0,2 -1 2,-1 -1 2)),((2 1 0,2 -1 2,2 -1 0,2 1 0)))");
lwfree(ewkt);
lwcollection_free(c);
lwgeom_free(g);
}
static void test_lwmline_clip(void)
......
......@@ -33,7 +33,8 @@ static LWGEOM* lwgeom_from_text(const char *str)
#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)
do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance3d_tolerance);\
do_test_mindistance_tolerance(str2, str1, res, __LINE__, lwgeom_mindistance3d_tolerance)
static void
do_test_mindistance_tolerance(char *in1,
......@@ -62,6 +63,9 @@ do_test_mindistance_tolerance(char *in1,
exit(1);
}
FLAGS_SET_SOLID(lw1->flags, 1);
FLAGS_SET_SOLID(lw2->flags, 1);
distance = distancef(lw1, lw2, 0.0);
lwgeom_free(lw1);
lwgeom_free(lw2);
......@@ -228,6 +232,15 @@ test_mindistance3d_tolerance(void)
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);
/* Same but triangles */
DIST3DTEST("LINESTRING(1 0 0 , 2 0 0)", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 1.0);
DIST3DTEST("LINESTRING(1 1 1 , 2 1 0)", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 0.0);
DIST3DTEST("LINESTRING(1 1 1 , 2 1 1)", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 1.0);
/* Triangle to triangle*/
DIST3DTEST("TRIANGLE((-1 1 0, -2 1 0, -1 2 0, -1 1 0))", "TRIANGLE((1 1 0, 2 1 0, 1 2 0, 1 1 0))", 2.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);
......@@ -235,7 +248,58 @@ test_mindistance3d_tolerance(void)
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);
DIST3DTEST("LINESTRING(1 1 1, 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", 0);
DIST3DTEST("TRIANGLE((1 1 1, 2 2 2, 3 3 3, 1 1 1))", "POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", 0);
DIST3DTEST("POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", "TRIANGLE((1 1 1, 2 2 2, 3 3 3, 1 1 1))", 0);
DIST3DTEST("TRIANGLE((0 0 0, 2 2 2, 3 3 3, 0 0 0))", "LINESTRING(1 1 1, 2 2 2)", 0);
/* A box in a box: two solids, one inside another */
DIST3DTEST(
"POLYHEDRALSURFACE Z (((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)))",
"POLYHEDRALSURFACE Z (((-1 -1 -1,-1 2 -1,2 2 -1,2 -1 -1,-1 -1 -1)),((-1 -1 2,2 -1 2,2 2 2,-1 2 2,-1 -1 2)),((-1 -1 -1,-1 -1 2,-1 2 2,-1 2 -1,-1 -1 -1)),((2 -1 -1,2 2 -1,2 2 2,2 -1 2,2 -1 -1)),((-1 -1 -1,2 -1 -1,2 -1 2,-1 -1 2,-1 -1 -1)),((-1 2 -1,-1 2 2,2 2 2,2 2 -1,-1 2 -1)))",
0);
/* A box in a box with a hat: two solids, one inside another, Z ray up is hitting hat */
DIST3DTEST(
"POLYHEDRALSURFACE Z (((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)))",
"POLYHEDRALSURFACE Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
0);
/* Same but as TIN */
DIST3DTEST(
"TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
"POLYHEDRALSURFACE Z (((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 1,1 0 1,1 1 1,0 1 1,0 0 1)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((1 0 0,1 1 0,1 1 1,1 0 1,1 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 1 1,1 1 1,1 1 0,0 1 0)))",
0);
/* Same but both are TIN */
DIST3DTEST(
"TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
"TIN Z (((0 0 0,0 1 0,1 1 0,0 0 0)),((1 0 0,0 0 0,1 1 0,1 0 0)),((0 1 1,1 0 1,1 1 1,0 1 1)),((0 1 1,0 0 1,1 0 1,0 1 1)),((0 0 0,0 0 1,0 1 1,0 0 0)),((0 1 0,0 0 0,0 1 1,0 1 0)),((1 0 1,1 1 0,1 1 1,1 0 1)),((1 0 1,1 0 0,1 1 0,1 0 1)),((0 0 1,1 0 0,1 0 1,0 0 1)),((0 0 1,0 0 0,1 0 0,0 0 1)),((0 1 0,0 1 1,1 1 1,0 1 0)),((1 1 0,0 1 0,1 1 1,1 1 0)))",
0);
/* Point inside TIN */
DIST3DTEST(
"TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
"POINT(0 0 0)",
0);
/* A point hits vertical Z edge */
DIST3DTEST(
"POLYHEDRALSURFACE Z (((0 -1 1,-1 -1 1,-1 -1 -1,0 -1 -1,1 -1 -1,0 -1 2,0 -1 1)),((0 1 1,0 1 2,1 1 -1,0 1 -1,-1 1 -1,-1 1 1,0 1 1)),((0 -1 1,0 1 1,-1 1 1,-1 -1 1,0 -1 1)),((-1 -1 1,-1 1 1,-1 1 -1,-1 -1 -1,-1 -1 1)),((-1 -1 -1,-1 1 -1,0 1 -1,0 -1 -1,-1 -1 -1)),((0 -1 -1,0 1 -1,1 1 -1,1 -1 -1,0 -1 -1)),((1 -1 -1,1 1 -1,0 1 2,0 -1 2,1 -1 -1)),((0 -1 2,0 1 2,0 1 1,0 -1 1,0 -1 2)))",
"POINT(0 0 0)",
0);
/* A point in the middle of a hole of extruded polygon */
DIST3DTEST(
"POLYHEDRALSURFACE Z (((-3 -3 0,-3 3 0,3 3 0,3 -3 0,-3 -3 0),(-1 -1 0,1 -1 0,1 1 0,-1 1 0,-1 -1 0)),((-3 -3 3,3 -3 3,3 3 3,-3 3 3,-3 -3 3),(-1 -1 3,-1 1 3,1 1 3,1 -1 3,-1 -1 3)),((-3 -3 0,-3 -3 3,-3 3 3,-3 3 0,-3 -3 0)),((-3 3 0,-3 3 3,3 3 3,3 3 0,-3 3 0)),((3 3 0,3 3 3,3 -3 3,3 -3 0,3 3 0)),((3 -3 0,3 -3 3,-3 -3 3,-3 -3 0,3 -3 0)),((-1 -1 0,-1 -1 3,1 -1 3,1 -1 0,-1 -1 0)),((1 -1 0,1 -1 3,1 1 3,1 1 0,1 -1 0)),((1 1 0,1 1 3,-1 1 3,-1 1 0,1 1 0)),((-1 1 0,-1 1 3,-1 -1 3,-1 -1 0,-1 1 0)))",
"POINT(0 0 1)",
1);
/* A point at the face of a hole of extruded polygon */
DIST3DTEST(
"POLYHEDRALSURFACE Z (((-3 -3 0,-3 3 0,3 3 0,3 -3 0,-3 -3 0),(-1 -1 0,1 -1 0,1 1 0,-1 1 0,-1 -1 0)),((-3 -3 3,3 -3 3,3 3 3,-3 3 3,-3 -3 3),(-1 -1 3,-1 1 3,1 1 3,1 -1 3,-1 -1 3)),((-3 -3 0,-3 -3 3,-3 3 3,-3 3 0,-3 -3 0)),((-3 3 0,-3 3 3,3 3 3,3 3 0,-3 3 0)),((3 3 0,3 3 3,3 -3 3,3 -3 0,3 3 0)),((3 -3 0,3 -3 3,-3 -3 3,-3 -3 0,3 -3 0)),((-1 -1 0,-1 -1 3,1 -1 3,1 -1 0,-1 -1 0)),((1 -1 0,1 -1 3,1 1 3,1 1 0,1 -1 0)),((1 1 0,1 1 3,-1 1 3,-1 1 0,1 1 0)),((-1 1 0,-1 1 3,-1 -1 3,-1 -1 0,-1 1 0)))",
"POINT(1 1 2)",
0);
}
static int tree_pt(RECT_NODE *tree, double x, double y)
......
......@@ -2068,16 +2068,17 @@ lwgeom_startpoint(const LWGEOM *lwgeom, POINT4D *pt)
return ptarray_startpoint(((LWLINE*)lwgeom)->points, pt);
case POLYGONTYPE:
return lwpoly_startpoint((LWPOLY*)lwgeom, pt);
case TINTYPE:
case CURVEPOLYTYPE:
case COMPOUNDTYPE:
case MULTIPOINTTYPE:
case MULTILINETYPE:
case MULTIPOLYGONTYPE:
case COLLECTIONTYPE:
case POLYHEDRALSURFACETYPE:
return lwcollection_startpoint((LWCOLLECTION*)lwgeom, pt);
default:
lwerror("int: unsupported geometry type: %s",
lwtype_name(lwgeom->type));
lwerror("lwgeom_startpoint: unsupported geometry type: %s", lwtype_name(lwgeom->type));
return LW_FAILURE;
}
}
......
......@@ -175,8 +175,6 @@ getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *op)
}
/*
* Copy a point from the point array into the parameter point
* will set point's z=NO_Z_VALUE if pa is 2d
......@@ -221,6 +219,7 @@ getPoint3dz_p(const POINTARRAY *pa, uint32_t n, POINT3DZ *op)
return 0;
}
assert(n < pa->npoints);
if ( n>=pa->npoints )
{
lwnotice("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
......@@ -276,7 +275,7 @@ getPoint3dm_p(const POINTARRAY *pa, uint32_t n, POINT3DM *op)
if ( n>=pa->npoints )
{
lwnotice("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
lwerror("%s [%d] called with n=%d and npoints=%d", __FILE__, __LINE__, n, pa->npoints);
return 0;
}
......@@ -497,23 +496,24 @@ void printPA(POINTARRAY *pa)
FLAGS_NDIMS(pa->flags), ptarray_point_size(pa));
lwnotice(" npoints = %i", pa->npoints);
for (t =0; t<pa->npoints; t++)
if (!pa)
{
getPoint4d_p(pa, t, &pt);
if (FLAGS_NDIMS(pa->flags) == 2)
{
lwnotice(" %i : %lf,%lf",t,pt.x,pt.y);
}
if (FLAGS_NDIMS(pa->flags) == 3)
{
lwnotice(" %i : %lf,%lf,%lf",t,pt.x,pt.y,pt.z);
}
if (FLAGS_NDIMS(pa->flags) == 4)
lwnotice(" PTARRAY is null pointer!");
}
else
{
for (t = 0; t < pa->npoints; t++)
{
lwnotice(" %i : %lf,%lf,%lf,%lf",t,pt.x,pt.y,pt.z,pt.m);
getPoint4d_p(pa, t, &pt);
if (FLAGS_NDIMS(pa->flags) == 2)
lwnotice(" %i : %lf,%lf", t, pt.x, pt.y);
if (FLAGS_NDIMS(pa->flags) == 3)
lwnotice(" %i : %lf,%lf,%lf", t, pt.x, pt.y, pt.z);
if (FLAGS_NDIMS(pa->flags) == 4)
lwnotice(" %i : %lf,%lf,%lf,%lf", t, pt.x, pt.y, pt.z, pt.m);
}
}
lwnotice(" }");
}
......
......@@ -285,12 +285,6 @@ lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
return;
}
if (!(ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M'))
{
lwerror("Cannot set %c ordinate.", ordinate);
return;
}
switch ( ordinate )
{
case 'X':
......@@ -306,6 +300,8 @@ lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
p->m = value;
return;
}
lwerror("Cannot set %c ordinate.", ordinate);
return;
}
/**
......@@ -345,11 +341,10 @@ point_interpolate(const POINT4D *p1,
}
#endif
proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
proportion = (interpolation_value - p1_value) / (p2_value - p1_value);
for (i = 0; i < 4; i++)
{
double newordinate = 0.0;
if (dims[i] == 'Z' && !hasz)
continue;
if (dims[i] == 'M' && !hasm)
......@@ -358,6 +353,7 @@ point_interpolate(const POINT4D *p1,
lwpoint_set_ordinate(p, dims[i], interpolation_value);
else
{
double newordinate = 0.0;
p1_value = lwpoint_get_ordinate(p1, dims[i]);
p2_value = lwpoint_get_ordinate(p2, dims[i]);
newordinate = p1_value + proportion * (p2_value - p1_value);
......@@ -478,7 +474,7 @@ ptarray_clamp_to_ordinate_range(const POINTARRAY *ipa, char ordinate, double fro
{
ptarray_append_point(opa, &p2, LW_FALSE);
}
else if (p1out == p2out && p1out != 0) /* both invisible */
else if (p1out == p2out && p1out != 0) /* both invisible on the same side */
{
/* skip */
}
......@@ -728,16 +724,12 @@ lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, do
static inline LWCOLLECTION *
lwpoly_clip_to_ordinate_range(const LWPOLY *poly, char ordinate, double from, double to)
{
LWCOLLECTION *lwgeom_out = NULL;
uint32_t i, nrings;
assert(poly);
char hasz = FLAGS_GET_Z(poly->flags), hasm = FLAGS_GET_M(poly->flags);
LWPOLY *poly_res = lwpoly_construct_empty(poly->srid, hasz, hasm);
LWCOLLECTION *lwgeom_out = lwcollection_construct_empty(MULTIPOLYGONTYPE, poly->srid, hasz, hasm);
assert(poly);
lwgeom_out = lwcollection_construct_empty(MULTIPOLYGONTYPE, poly->srid, hasz, hasm);
nrings = poly->nrings;
for (i = 0; i < nrings; i++)
for (uint32_t i = 0; i < poly->nrings; i++)
{
/* Ret number of points */
POINTARRAY *pa = ptarray_clamp_to_ordinate_range(poly->rings[i], ordinate, from, to, LW_TRUE);
......@@ -751,7 +743,10 @@ lwpoly_clip_to_ordinate_range(const LWPOLY *poly, char ordinate, double from, do
break;
}
}
lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, (LWGEOM *)poly_res);
if (poly_res->nrings > 0)
lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, (LWGEOM *)poly_res);
else
lwpoly_free(poly_res);
return lwgeom_out;
}
......@@ -762,12 +757,9 @@ lwpoly_clip_to_ordinate_range(const LWPOLY *poly, char ordinate, double from, do
static inline LWCOLLECTION *
lwtriangle_clip_to_ordinate_range(const LWTRIANGLE *tri, char ordinate, double from, double to)
{
LWCOLLECTION *lwgeom_out = NULL;
char hasz = FLAGS_GET_Z(tri->flags), hasm = FLAGS_GET_M(tri->flags);
assert(tri);
lwgeom_out = lwcollection_construct_empty(TINTYPE, tri->srid, hasz, hasm);
char hasz = FLAGS_GET_Z(tri->flags), hasm = FLAGS_GET_M(tri->flags);
LWCOLLECTION *lwgeom_out = lwcollection_construct_empty(TINTYPE, tri->srid, hasz, hasm);
POINTARRAY *pa = ptarray_clamp_to_ordinate_range(tri->points, ordinate, from, to, LW_TRUE);
if (pa->npoints >= 4)
......@@ -870,6 +862,7 @@ lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, do
case MULTILINETYPE:
case MULTIPOLYGONTYPE:
case COLLECTIONTYPE:
case POLYHEDRALSURFACETYPE:
out_col = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)lwin, ordinate, from, to);
break;
default:
......
......@@ -532,19 +532,29 @@ int
lwpoly_contains_point(const LWPOLY *poly, const POINT2D *pt)
{
uint32_t i;
int t;
if ( lwpoly_is_empty(poly) )
return LW_FALSE;
return LW_OUTSIDE;
if ( ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE )
return LW_FALSE;
t = ptarray_contains_point(poly->rings[0], pt);
for ( i = 1; i < poly->nrings; i++ )
if (t == LW_INSIDE)
{
if ( ptarray_contains_point(poly->rings[i], pt) == LW_INSIDE )
return LW_FALSE;
for (i = 1; i < poly->nrings; i++)
{
t = ptarray_contains_point(poly->rings[i], pt);
if (t == LW_INSIDE)
return LW_OUTSIDE;
if (t == LW_BOUNDARY)
{
return LW_BOUNDARY;
}
}
return LW_INSIDE;
}
return LW_TRUE;
else
return t;
}
......@@ -19,6 +19,7 @@
**********************************************************************
*
* Copyright 2011 Nicklas Avén
* Copyright 2019 Darafei Praliaskouski
*
**********************************************************************/
......@@ -26,6 +27,7 @@
#include <stdlib.h>
#include "measures3d.h"
#include "lwrandom.h"
#include "lwgeom_log.h"
static inline int
......@@ -96,7 +98,7 @@ lw_dist3d_distanceline(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode)
{
LWDEBUG(2, "lw_dist3d_distanceline is called");
double x1, x2, y1, y2, z1, z2, x, y;
double initdistance = (mode == DIST_MIN ? FLT_MAX : -1.0);
double initdistance = (mode == DIST_MIN ? DBL_MAX : -1.0);
DISTPTS3D thedl;
LWPOINT *lwpoints[2];
LWGEOM *result;
......@@ -201,7 +203,7 @@ lw_dist3d_distancepoint(const LWGEOM *lw1, const LWGEOM *lw2, int srid, int mode
double x, y, z;
DISTPTS3D thedl;
double initdistance = FLT_MAX;
double initdistance = DBL_MAX;
LWGEOM *result;
thedl.mode = mode;
......@@ -335,6 +337,103 @@ lwgeom_mindistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
return lwgeom_mindistance3d_tolerance(lw1, lw2, 0.0);
}
static inline int
gbox_contains_3d(const GBOX *g1, const GBOX *g2)
{
return (g2->xmin >= g1->xmin) && (g2->xmax <= g1->xmax) && (g2->ymin >= g1->ymin) && (g2->ymax <= g1->ymax) &&
(g2->zmin >= g1->zmin) && (g2->zmax <= g1->zmax);
}
static inline int
lwgeom_solid_contains_lwgeom(const LWGEOM *solid, const LWGEOM *g)
{
const GBOX *b1;
const GBOX *b2;
/* If solid isn't solid it can't contain anything */
if (!FLAGS_GET_SOLID(solid->flags))
return LW_FALSE;
b1 = lwgeom_get_bbox(solid);
b2 = lwgeom_get_bbox(g);
/* If box won't contain box, shape won't contain shape */
if (!gbox_contains_3d(b1, b2))
return LW_FALSE;
else /* Raycast upwards in Z */
{
POINT4D pt;
/* We'll skew copies if we're not lucky */
LWGEOM *solid_copy = lwgeom_clone_deep(solid);
LWGEOM *g_copy = lwgeom_clone_deep(g);
while (LW_TRUE)
{
uint8_t is_boundary = LW_FALSE;
uint8_t is_inside = LW_FALSE;
/* take first point */
if (!lwgeom_startpoint(g_copy, &pt))
return LW_FALSE;
/* get part of solid that is upwards */
LWCOLLECTION *c = lwgeom_clip_to_ordinate_range(solid_copy, 'Z', pt.z, DBL_MAX, 0);
for (uint32_t i = 0; i < c->ngeoms; i++)
{
if (c->geoms[i]->type == POLYGONTYPE)
{
/* 3D raycast along Z is 2D point in polygon */
int t = lwpoly_contains_point((LWPOLY *)c->geoms[i], (POINT2D *)&pt);
if (t == LW_INSIDE)
is_inside = !is_inside;
else if (t == LW_BOUNDARY)
{
is_boundary = LW_TRUE;
break;
}
}
else if (c->geoms[i]->type == TRIANGLETYPE)
{
/* 3D raycast along Z is 2D point in polygon */
LWTRIANGLE *tri = (LWTRIANGLE *)c->geoms[i];
int t = ptarray_contains_point(tri->points, (POINT2D *)&pt);
if (t == LW_INSIDE)
is_inside = !is_inside;
else if (t == LW_BOUNDARY)
{
is_boundary = LW_TRUE;
break;
}
}
}
lwcollection_free(c);
lwgeom_free(solid_copy);
lwgeom_free(g_copy);
if (is_boundary)
{
/* randomly skew a bit and restart*/
double cx = lwrandom_uniform() - 0.5;
double cy = lwrandom_uniform() - 0.5;
AFFINE aff = {1, 0, cx, 0, 1, cy, 0, 0, 1, 0, 0, 0};
solid_copy = lwgeom_clone_deep(solid);
lwgeom_affine(solid_copy, &aff);
g_copy = lwgeom_clone_deep(g);
lwgeom_affine(g_copy, &aff);
continue;
}
return is_inside;
}
}
}
/**
Function handling 3d min distance calculations and dwithin calculations.
The difference is just the tolerance.
......@@ -342,6 +441,7 @@ lwgeom_mindistance3d(const LWGEOM *lw1, const LWGEOM *lw2)
double
lwgeom_mindistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
{
assert(tolerance >= 0);
if (!lwgeom_has_z(lw1) || !lwgeom_has_z(lw2))
{
lwnotice(
......@@ -350,17 +450,23 @@ lwgeom_mindistance3d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tole
return lwgeom_mindistance2d_tolerance(lw1, lw2, tolerance);
}
DISTPTS3D thedl;
LWDEBUG(2, "lwgeom_mindistance3d_tolerance is called");
thedl.mode = DIST_MIN;
thedl.distance = FLT_MAX;
thedl.distance = DBL_MAX;
thedl.tolerance = tolerance;
if (lw_dist3d_recursive(lw1, lw2, &thedl))
{
if (thedl.distance <= tolerance)
return thedl.distance;
if (lwgeom_solid_contains_lwgeom(lw1, lw2) || lwgeom_solid_contains_lwgeom(lw2, lw1))
return 0;
return thedl.distance;
}
/* should never get here. all cases ought to be error handled earlier */
lwerror("Some unspecified error.");
return FLT_MAX;
return DBL_MAX;
}
/*------------------------------------------------------------------------------------------------------------
......@@ -403,7 +509,6 @@ lw_dist3d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
for (i = 0; i < n1; i++)
{
if (lwgeom_is_collection(lwg1))
g1 = c1->geoms[i];
else
......@@ -424,8 +529,8 @@ lw_dist3d_recursive(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
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");
......@@ -455,7 +560,6 @@ This function distributes the brute-force for 3D so far the only type, tasks dep
int
lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3D *dl)
{
int t1 = lwg1->type;
int t2 = lwg2->type;
......@@ -478,6 +582,11 @@ lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3
dl->twisted = 1;
return lw_dist3d_point_poly((LWPOINT *)lwg1, (LWPOLY *)lwg2, dl);
}
else if (t2 == TRIANGLETYPE)
{
dl->twisted = 1;
return lw_dist3d_point_tri((LWPOINT *)lwg1, (LWTRIANGLE *)lwg2, dl);
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
......@@ -488,7 +597,7 @@ lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3
{
if (t2 == POINTTYPE)
{
dl->twisted = (-1);
dl->twisted = -1;
return lw_dist3d_point_line((LWPOINT *)lwg2, (LWLINE *)lwg1, dl);
}
else if (t2 == LINETYPE)
......@@ -501,6 +610,11 @@ lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3
dl->twisted = 1;
return lw_dist3d_line_poly((LWLINE *)lwg1, (LWPOLY *)lwg2, dl);
}
else if (t2 == TRIANGLETYPE)
{
dl->twisted = 1;
return lw_dist3d_line_tri((LWLINE *)lwg1, (LWTRIANGLE *)lwg2, dl);
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
......@@ -524,12 +638,46 @@ lw_dist3d_distribute_bruteforce(const LWGEOM *lwg1, const LWGEOM *lwg2, DISTPTS3
dl->twisted = -1;
return lw_dist3d_line_poly((LWLINE *)lwg2, (LWPOLY *)lwg1, dl);
}
else if (t2 == TRIANGLETYPE)
{
dl->twisted = 1;
return lw_dist3d_poly_tri((LWPOLY *)lwg1, (LWTRIANGLE *)lwg2, dl);
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
return LW_FALSE;
}
}
else if (t1 == TRIANGLETYPE)
{
if (t2 == POLYGONTYPE)
{
dl->twisted = -1;
return lw_dist3d_poly_tri((LWPOLY *)lwg2, (LWTRIANGLE *)lwg1, dl);
}
else if (t2 == POINTTYPE)
{
dl->twisted = -1;
return lw_dist3d_point_tri((LWPOINT *)lwg2, (LWTRIANGLE *)lwg1, dl);
}
else if (t2 == LINETYPE)
{
dl->twisted = -1;
return lw_dist3d_line_tri((LWLINE *)lwg2, (LWTRIANGLE *)lwg1, dl);
}
else if (t2 == TRIANGLETYPE)
{
dl->twisted = 1;
return lw_dist3d_tri_tri((LWTRIANGLE *)lwg1, (LWTRIANGLE *)lwg2, dl);
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t2));
return LW_FALSE;
}
}
else
{
lwerror("Unsupported geometry type: %s", lwtype_name(t1));
......@@ -588,6 +736,7 @@ For mindistance that means:
for max distance it is always point against boundary
*/
int
lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)
{
......@@ -596,30 +745,46 @@ lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)
LWDEBUG(2, "lw_dist3d_point_poly is called");
getPoint3dz_p(point->point, 0, &p);
/*If we are looking for max distance, longestline or dfullywithin*/
/* If we are looking 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*/
/* 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))
{
/* 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*/
/* 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);
}
/**
/* point to triangle calculation */
int
lw_dist3d_point_tri(LWPOINT *point, LWTRIANGLE *tri, DISTPTS3D *dl)
{
POINT3DZ p, projp; /*projp is "point projected on plane"*/
PLANE3D plane;
getPoint3dz_p(point->point, 0, &p);
line to line calculation
*/
/* If we are looking for max distance, longestline or dfullywithin */
if (dl->mode == DIST_MAX)
return lw_dist3d_pt_ptarray(&p, tri->points, dl);
/* If triangle does not define a plane, treat it as a line */
if (!define_plane(tri->points, &plane))
return lw_dist3d_pt_ptarray(&p, tri->points, dl);
/* Get our point projected on the plane of triangle */
project_point_on_plane(&p, &plane, &projp);
return lw_dist3d_pt_tri(&p, tri, &plane, &projp, dl);
}
/** line to line calculation */
int
lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
{
......@@ -630,9 +795,7 @@ lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl)
return lw_dist3d_ptarray_ptarray(pa1, pa2, dl);
}
/**
line to polygon calculation
*/
/** line to polygon calculation */
int
lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
{
......@@ -640,22 +803,32 @@ lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl)
LWDEBUG(2, "lw_dist3d_line_poly is called");
if (dl->mode == DIST_MAX)
{
return lw_dist3d_ptarray_ptarray(line->points, poly->rings[0], dl);
}
/* if polygon does not define a plane: Return distance line to line */
if (!define_plane(poly->rings[0], &plane))
{
/* 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);
}
/**
polygon to polygon calculation
*/
/** line to triangle calculation */
int
lw_dist3d_line_tri(LWLINE *line, LWTRIANGLE *tri, DISTPTS3D *dl)
{
PLANE3D plane;
if (dl->mode == DIST_MAX)
return lw_dist3d_ptarray_ptarray(line->points, tri->points, dl);