[postgis-tickets] r17260 - Make lwgeom_mindist3d solid aware. Add TIN support.
Darafei
komzpa at gmail.com
Fri Feb 15 05:05:39 PST 2019
Author: komzpa
Date: 2019-02-15 17:05:39 -0800 (Fri, 15 Feb 2019)
New Revision: 17260
Modified:
trunk/NEWS
trunk/liblwgeom/cunit/cu_algorithm.c
trunk/liblwgeom/cunit/cu_measures.c
trunk/liblwgeom/lwgeom.c
trunk/liblwgeom/lwgeom_api.c
trunk/liblwgeom/lwlinearreferencing.c
trunk/liblwgeom/lwpoly.c
trunk/liblwgeom/measures3d.c
trunk/liblwgeom/measures3d.h
trunk/regress/core/measures.sql
trunk/regress/core/measures_expected
Log:
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
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/NEWS 2019-02-16 01:05:39 UTC (rev 17260)
@@ -71,8 +71,9 @@
- #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
WARNING: If compiling with PostgreSQL+JIT, LLVM >= 6 is required
Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/cunit/cu_algorithm.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -778,6 +778,20 @@
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)
Modified: trunk/liblwgeom/cunit/cu_measures.c
===================================================================
--- trunk/liblwgeom/cunit/cu_measures.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/cunit/cu_measures.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -33,7 +33,8 @@
#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 @@
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 @@
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 @@
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)
Modified: trunk/liblwgeom/lwgeom.c
===================================================================
--- trunk/liblwgeom/lwgeom.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/lwgeom.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -2068,6 +2068,7 @@
return ptarray_startpoint(((LWLINE*)lwgeom)->points, pt);
case POLYGONTYPE:
return lwpoly_startpoint((LWPOLY*)lwgeom, pt);
+ case TINTYPE:
case CURVEPOLYTYPE:
case COMPOUNDTYPE:
case MULTIPOINTTYPE:
@@ -2074,10 +2075,10 @@
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;
}
}
Modified: trunk/liblwgeom/lwgeom_api.c
===================================================================
--- trunk/liblwgeom/lwgeom_api.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/lwgeom_api.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -175,8 +175,6 @@
}
-
-
/*
* 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 @@
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 @@
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 @@
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(" PTARRAY is null pointer!");
+ }
+ else
+ {
+
+ for (t = 0; t < pa->npoints; t++)
{
- lwnotice(" %i : %lf,%lf",t,pt.x,pt.y);
+ 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);
}
- 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(" }");
}
Modified: trunk/liblwgeom/lwlinearreferencing.c
===================================================================
--- trunk/liblwgeom/lwlinearreferencing.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/lwlinearreferencing.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -285,12 +285,6 @@
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 @@
p->m = value;
return;
}
+ lwerror("Cannot set %c ordinate.", ordinate);
+ return;
}
/**
@@ -345,11 +341,10 @@
}
#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 @@
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_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 @@
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 @@
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 @@
static inline LWCOLLECTION *
lwtriangle_clip_to_ordinate_range(const LWTRIANGLE *tri, char ordinate, double from, double to)
{
- LWCOLLECTION *lwgeom_out = NULL;
+ assert(tri);
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);
-
+ 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 @@
case MULTILINETYPE:
case MULTIPOLYGONTYPE:
case COLLECTIONTYPE:
+ case POLYHEDRALSURFACETYPE:
out_col = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)lwin, ordinate, from, to);
break;
default:
Modified: trunk/liblwgeom/lwpoly.c
===================================================================
--- trunk/liblwgeom/lwpoly.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/lwpoly.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -532,19 +532,29 @@
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;
}
Modified: trunk/liblwgeom/measures3d.c
===================================================================
--- trunk/liblwgeom/measures3d.c 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/measures3d.c 2019-02-16 01:05:39 UTC (rev 17260)
@@ -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 @@
{
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 @@
double x, y, z;
DISTPTS3D thedl;
- double initdistance = FLT_MAX;
+ double initdistance = DBL_MAX;
LWGEOM *result;
thedl.mode = mode;
@@ -335,6 +337,103 @@
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 @@
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 @@
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 @@
for (i = 0; i < n1; i++)
{
-
if (lwgeom_is_collection(lwg1))
g1 = c1->geoms[i];
else
@@ -424,8 +529,8 @@
if (lwgeom_is_collection(lwg2))
g2 = c2->geoms[j];
else
+ g2 = (LWGEOM *)lwg2;
- g2 = (LWGEOM *)lwg2;
if (lwgeom_is_collection(g2))
{
LWDEBUG(3, "Found collection inside second geometry collection, recursing");
@@ -455,7 +560,6 @@
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 @@
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 @@
{
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 @@
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,6 +638,11 @@
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));
@@ -530,6 +649,35 @@
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 max distance it is always point against boundary
*/
+
int
lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl)
{
@@ -596,14 +745,11 @@
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 */
@@ -610,16 +756,35 @@
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 @@
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,23 +803,33 @@
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);
+
+ /* if triangle does not define a plane: Return distance line to line */
+ if (!define_plane(tri->points, &plane))
+ return lw_dist3d_ptarray_ptarray(line->points, tri->points, dl);
+
+ return lw_dist3d_ptarray_tri(line->points, tri, &plane, dl);
+}
+
+/** polygon to polygon calculation */
+int
lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl)
{
PLANE3D plane1, plane2;
@@ -663,9 +836,7 @@
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);
- }
planedef1 = define_plane(poly1->rings[0], &plane1);
planedef2 = define_plane(poly2->rings[0], &plane2);
@@ -672,20 +843,17 @@
if (!planedef1 || !planedef2)
{
+ /* Neither polygon define a plane: Return distance line to line */
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 */
+ /* Only poly2 defines a plane: Return distance from line (poly1) to poly2 */
+ else if (!planedef1)
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);
+ else
+ 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
@@ -701,8 +869,89 @@
return lw_dist3d_ptarray_poly(poly2->rings[0], poly1, &plane1, dl);
}
+/** polygon to triangle calculation */
+int
+lw_dist3d_poly_tri(LWPOLY *poly, LWTRIANGLE *tri, DISTPTS3D *dl)
+{
+ PLANE3D plane1, plane2;
+ int planedef1, planedef2;
+
+ if (dl->mode == DIST_MAX)
+ return lw_dist3d_ptarray_ptarray(poly->rings[0], tri->points, dl);
+
+ planedef1 = define_plane(poly->rings[0], &plane1);
+ planedef2 = define_plane(tri->points, &plane2);
+
+ if (!planedef1 || !planedef2)
+ {
+ /* Neither define a plane: Return distance line to line */
+ if (!planedef1 && !planedef2)
+ return lw_dist3d_ptarray_ptarray(poly->rings[0], tri->points, dl);
+
+ /* Only tri defines a plane: Return distance from line (poly) to tri */
+ else if (!planedef1)
+ return lw_dist3d_ptarray_tri(poly->rings[0], tri, &plane2, dl);
+
+ /* Only poly defines a plane: Return distance from line (tri) to poly */
+ else
+ return lw_dist3d_ptarray_poly(tri->points, poly, &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_tri(poly->rings[0], tri, &plane2, dl))
+ return LW_FALSE;
+ if (dl->distance < dl->tolerance) /* Just check if the answer already is given*/
+ return LW_TRUE;
+
+ 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(tri->points, poly, &plane1, dl);
+}
+
+/** triangle to triangle calculation */
+int
+lw_dist3d_tri_tri(LWTRIANGLE *tri1, LWTRIANGLE *tri2, DISTPTS3D *dl)
+{
+ PLANE3D plane1, plane2;
+ int planedef1, planedef2;
+
+ if (dl->mode == DIST_MAX)
+ return lw_dist3d_ptarray_ptarray(tri1->points, tri2->points, dl);
+
+ planedef1 = define_plane(tri1->points, &plane1);
+ planedef2 = define_plane(tri2->points, &plane2);
+
+ if (!planedef1 || !planedef2)
+ {
+ /* Neither define a plane: Return distance line to line */
+ if (!planedef1 && !planedef2)
+ return lw_dist3d_ptarray_ptarray(tri1->points, tri2->points, dl);
+
+ /* Only tri defines a plane: Return distance from line (tri1) to tri2 */
+ else if (!planedef1)
+ return lw_dist3d_ptarray_tri(tri1->points, tri2, &plane2, dl);
+
+ /* Only poly defines a plane: Return distance from line (tri2) to tri1 */
+ else
+ return lw_dist3d_ptarray_tri(tri2->points, tri1, &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_tri(tri1->points, tri2, &plane2, dl))
+ return LW_FALSE;
+ if (dl->distance < dl->tolerance) /* Just check if the answer already is given*/
+ return LW_TRUE;
+
+ 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_tri(tri2->points, tri1, &plane1, dl);
+}
+
/**
-
* search all the segments of pointarray to see which one is closest to p
* Returns distance between point and pointarray
*/
@@ -712,9 +961,9 @@
uint32_t t;
POINT3DZ start, end;
int twist = dl->twisted;
+ if (!pa)
+ return LW_FALSE;
- LWDEBUG(2, "lw_dist3d_pt_ptarray is called");
-
getPoint3dz_p(pa, 0, &start);
for (t = 1; t < pa->npoints; t++)
@@ -1005,28 +1254,37 @@
{
/* Inside a hole. Distance = pt -> ring */
if (pt_in_ring_3d(projp, poly->rings[i], plane))
- {
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*/
+ /* if the projected point is inside the polygon the shortest distance is between that point and the
+ * input point */
+ return lw_dist3d_pt_pt(p, projp, dl);
}
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 boundary instead*/
- }
+ /* if the projected point is outside the polygon we search for the closest distance against the boundary
+ * instead */
+ return lw_dist3d_pt_ptarray(p, poly->rings[0], dl);
return LW_TRUE;
}
-/**
+int
+lw_dist3d_pt_tri(POINT3DZ *p, LWTRIANGLE *tri, PLANE3D *plane, POINT3DZ *projp, DISTPTS3D *dl)
+{
+ if (pt_in_ring_3d(projp, tri->points, plane))
+ /* if the projected point is inside the polygon the shortest distance is between that point and the
+ * input point */
+ return lw_dist3d_pt_pt(p, projp, dl);
+ else
+ /* if the projected point is outside the polygon we search for the closest distance against the boundary
+ * instead */
+ return lw_dist3d_pt_ptarray(p, tri->points, dl);
-Computes pointarray to polygon distance
-*/
+ return LW_TRUE;
+}
+
+/** Computes pointarray to polygon distance */
int
lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, PLANE3D *plane, DISTPTS3D *dl)
{
@@ -1037,8 +1295,8 @@
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. */
+ /* the sign of s1 tells us on which side of the plane the point is. */
+ s1 = project_point_on_plane(&p1, plane, &projp1);
lw_dist3d_pt_poly(&p1, poly, plane, &projp1, dl);
if ((s1 == 0.0) && (dl->distance < dl->tolerance))
return LW_TRUE;
@@ -1052,21 +1310,21 @@
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 and 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.*/
+ /* The size of s1 and s2 is the distance from the point to the plane. */
+ f = fabs(s1) / (fabs(s1) + fabs(s2));
get_3dvector_from_points(&projp1, &projp2, &projp1_projp2);
- /*get the point where the line segment crosses the plane*/
+ /* 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;
- intersects = LW_TRUE; /*We set intersects to true until the opposite is proved*/
+ /* We set intersects to true until the opposite is proved */
+ intersects = LW_TRUE;
if (pt_in_ring_3d(&intersectionp, poly->rings[0], plane)) /*Inside outer ring*/
{
@@ -1099,12 +1357,80 @@
p1 = p2;
}
- /*check or pointarray against boundary and inner boundaries of the polygon*/
+ /* check our 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;
+}
+
+/** Computes pointarray to triangle distance */
+int
+lw_dist3d_ptarray_tri(POINTARRAY *pa, LWTRIANGLE *tri, PLANE3D *plane, DISTPTS3D *dl)
+{
+ uint32_t i;
+ double f, s1, s2;
+ VECTOR3D projp1_projp2;
+ POINT3DZ p1, p2, projp1, projp2, intersectionp;
+
+ getPoint3dz_p(pa, 0, &p1);
+
+ /* the sign of s1 tells us on which side of the plane the point is. */
+ s1 = project_point_on_plane(&p1, plane, &projp1);
+ lw_dist3d_pt_tri(&p1, tri, plane, &projp1, dl);
+ if ((s1 == 0.0) && (dl->distance < dl->tolerance))
+ return LW_TRUE;
+
+ for (i = 1; i < pa->npoints; i++)
{
- lw_dist3d_ptarray_ptarray(pa, poly->rings[j], dl);
+ int intersects;
+ getPoint3dz_p(pa, i, &p2);
+ s2 = project_point_on_plane(&p2, plane, &projp2);
+ lw_dist3d_pt_tri(&p2, tri, plane, &projp2, dl);
+ if ((s2 == 0.0) && (dl->distance < dl->tolerance))
+ return LW_TRUE;
+
+ /* If s1 and s2 has different signs that means they are on different sides of the plane of the triangle.
+ * That means that the edge between the points crosses the plane and might intersect with the triangle
+ */
+ if ((s1 * s2) < 0)
+ {
+ /* The size of s1 and s2 is the distance from the point to the plane. */
+ f = fabs(s1) / (fabs(s1) + fabs(s2));
+ 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;
+
+ /* We set intersects to true until the opposite is proved */
+ intersects = LW_TRUE;
+
+ if (pt_in_ring_3d(&intersectionp, tri->points, plane)) /*Inside outer ring*/
+ {
+ 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;
+ }
+ }
+ }
+
+ projp1 = projp2;
+ s1 = s2;
+ p1 = p2;
}
+ /* check our pointarray against triangle contour */
+ lw_dist3d_ptarray_ptarray(pa, tri->points, dl);
return LW_TRUE;
}
@@ -1121,9 +1447,18 @@
define_plane(POINTARRAY *pa, PLANE3D *pl)
{
const uint32_t POL_BREAKS = 3;
+
+ assert(pa);
+ assert(pa->npoints > 3);
+ if (!pa)
+ return LW_FALSE;
+
uint32_t unique_points = pa->npoints - 1;
uint32_t i;
+ if (pa->npoints < 3)
+ return LW_FALSE;
+
/* Calculate the average point */
pl->pop.x = 0.0;
pl->pop.y = 0.0;
Modified: trunk/liblwgeom/measures3d.h
===================================================================
--- trunk/liblwgeom/measures3d.h 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/liblwgeom/measures3d.h 2019-02-16 01:05:39 UTC (rev 17260)
@@ -19,6 +19,7 @@
**********************************************************************
*
* Copyright 2011 Nicklas Avén
+ * Copyright 2019 Darafei Praliaskouski
*
**********************************************************************/
@@ -76,16 +77,28 @@
int lw_dist3d_pt_ptarray(POINT3DZ *p, POINTARRAY *pa, DISTPTS3D *dl);
int lw_dist3d_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS3D *dl);
int lw_dist3d_point_line(LWPOINT *point, LWLINE *line, DISTPTS3D *dl);
+int lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl);
+int lw_dist3d_point_tri(LWPOINT *point, LWTRIANGLE *tri, DISTPTS3D *dl);
+
int lw_dist3d_line_line(LWLINE *line1, LWLINE *line2, DISTPTS3D *dl);
-int lw_dist3d_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS3D *dl);
int lw_dist3d_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS3D *dl);
+int lw_dist3d_line_tri(LWLINE *line, LWTRIANGLE *tri, DISTPTS3D *dl);
+
int lw_dist3d_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS3D *dl);
-int lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS3D *dl);
-int lw_dist3d_seg_seg(POINT3DZ *A, POINT3DZ *B, POINT3DZ *C, POINT3DZ *D, DISTPTS3D *dl);
+int lw_dist3d_poly_tri(LWPOLY *poly, LWTRIANGLE *tri, DISTPTS3D *dl);
+
+int lw_dist3d_tri_tri(LWTRIANGLE *tri1, LWTRIANGLE *tri2, DISTPTS3D *dl);
+
int lw_dist3d_pt_pt(POINT3DZ *p1, POINT3DZ *p2, DISTPTS3D *dl);
int lw_dist3d_pt_seg(POINT3DZ *p, POINT3DZ *A, POINT3DZ *B, DISTPTS3D *dl);
int lw_dist3d_pt_poly(POINT3DZ *p, LWPOLY *poly, PLANE3D *plane, POINT3DZ *projp, DISTPTS3D *dl);
+int lw_dist3d_pt_tri(POINT3DZ *p, LWTRIANGLE *tri, PLANE3D *plane, POINT3DZ *projp, DISTPTS3D *dl);
+
+int lw_dist3d_seg_seg(POINT3DZ *A, POINT3DZ *B, POINT3DZ *C, POINT3DZ *D, DISTPTS3D *dl);
+
+int lw_dist3d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS3D *dl);
int lw_dist3d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, PLANE3D *plane, DISTPTS3D *dl);
+int lw_dist3d_ptarray_tri(POINTARRAY *pa, LWTRIANGLE *tri, PLANE3D *plane, DISTPTS3D *dl);
double project_point_on_plane(POINT3DZ *p, PLANE3D *pl, POINT3DZ *p0);
int define_plane(POINTARRAY *pa, PLANE3D *pl);
Modified: trunk/regress/core/measures.sql
===================================================================
--- trunk/regress/core/measures.sql 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/regress/core/measures.sql 2019-02-16 01:05:39 UTC (rev 17260)
@@ -274,3 +274,11 @@
select 'length2d_spheroid', ST_Length2DSpheroid('LINESTRING(0 0 0, 0 0 100)'::geometry, 'SPHEROID["GRS_1980",6378137,298.257222101]');
select 'length_spheroid', ST_LengthSpheroid('LINESTRING(0 0 0, 0 0 100)'::geometry, 'SPHEROID["GRS_1980",6378137,298.257222101]');
+
+
+-- Solid intersects solid when contains it
+select '#4278.1', ST_3DIntersects('BOX3D(0 0 0, 1 1 1)'::box3d::geometry, 'BOX3D(-1 -1 -1, 2 2 2)'::box3d::geometry);
+select '#4278.2', ST_3DDistance('BOX3D(0 0 0, 1 1 1)'::box3d::geometry, 'BOX3D(-1 -1 -1, 2 2 2)'::box3d::geometry);
+-- cast to text and back as a way of getting rid of solid flag
+select '#4278.3', ST_3DIntersects('BOX3D(0 0 0, 1 1 1)'::box3d::geometry::text::geometry, 'BOX3D(-1 -1 -1, 2 2 2)'::box3d::geometry::text::geometry);
+select '#4278.4', ST_3DDistance('BOX3D(0 0 0, 1 1 1)'::box3d::geometry::text::geometry, 'BOX3D(-1 -1 -1, 2 2 2)'::box3d::geometry::text::geometry);
Modified: trunk/regress/core/measures_expected
===================================================================
--- trunk/regress/core/measures_expected 2019-02-15 19:01:06 UTC (rev 17259)
+++ trunk/regress/core/measures_expected 2019-02-16 01:05:39 UTC (rev 17260)
@@ -63,3 +63,7 @@
spheroidLength1|85204.52077
length2d_spheroid|100
length_spheroid|100
+#4278.1|t
+#4278.2|0
+#4278.3|f
+#4278.4|1
More information about the postgis-tickets
mailing list