3d intersects in postgis hangs on simple polyhedron

[lbartoletti]

This postgis query hangs without returning,

SELECT ST_3DIntersects('POINT(0 0.5 0.5)', ST_MakeSolid('POLYHEDRALSURFACE Z (((2 0.5 0.5, 0 0 1, 0 0 0, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 1, 0 0 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 0, 0 1 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 0 0, 0 1 0, 2 0.5 0.5)), ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)))'))

The polyhedron is a pyramid with it's base being the y-z unit square and apex in (2, 0.5, 0.5). The problematic point (found by happenstance) seems to be precicely (0, 0.5, 0.5), located on the base of the rotated pyramid.

Tested both on

postgis 3.2 with SFCGAL 1.3.8 (official postgis docker image)

as well as this environment:

POSTGIS="3.3.0dev 3.2.0-783-g042894569" \[EXTENSION\] PGSQL="110" GEOS="3.8.3-CAPI-1.13.4" SFCGAL="SFCGAL 1.4.1, CGAL 5.3, BOOST 1.78.0" PROJ="7.1.1" GDAL="GDAL 3.2.3, released 2021/04/27" LIBXML="2.9.9" LIBJSON="0.12" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" RASTER , PostgreSQL 11.2 on x86_64-pc-mingw64, compiled by gcc.exe (x86_64-posix-seh-rev0, Built by MinGW-W64 project) 8.1.0, 64-bit

[lbartoletti]

Nota:

SELECT ST_3DIntersects('POINT(0 0.5 0.5)', 'POLYHEDRALSURFACE Z (((2 0.5 0.5, 0 0 1, 0 0 0, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 1, 0 0 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 0, 0 1 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 0 0, 0 1 0, 2 0.5 0.5)), ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)))')

does not hang, but returns false

[lbartoletti]

BTW, it's ok on pure SFCGAL side (c++ and python):

Python example:

import pysfcgal.sfcgal as sfcgal
geom = sfcgal.read_wkt('POLYHEDRALSURFACE Z (((2 0.5 0.5, 0 0 1, 0 0 0, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 1, 0 0 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 0, 0 1 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 0 0, 0 1 0, 2 0.5 0.5)), ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)))')
geom_p = sfcgal.read_wkt('POINT( 0 0.5 0.5 )')
geom_p.intersects_3d(geom)
# True
geom.intersects_3d(geom_p)
# True

[lbartoletti]

So:

• Why ST_3DIntersects('POINT(0 0.5 0.5)', 'POLYHEDRALSURFACE Z (... returns false? I guess it's a precision problem (fix to float)
• Why ST_3DIntersects hangs with ST_MakeSolid

[lbartoletti]

Tested with geos and sfcgal backend.

[robe]

This appears to be still an issue with GEOS 3.11. I think ST_3DIntersects is completely native PostGIS / GEOS now.

— this returns false

SELECT ST_3DIntersects('POINT(0 0.5 0.5)'::geometry, 'POLYHEDRALSURFACE Z (((2 0.5 0.5, 0 0 1, 0 0 0, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 1, 0 0 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 0, 0 1 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 0 0, 0 1 0, 2 0.5 0.5)), ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)))'::geometry)

— this goes on forever

SELECT ST_3DIntersects('POINT(0 0.5 0.5)'::geometry, ST_MakeSolid('POLYHEDRALSURFACE Z (((2 0.5 0.5, 0 0 1, 0 0 0, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 1, 0 0 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 1 0, 0 1 1, 2 0.5 0.5)), ((2 0.5 0.5, 0 0 0, 0 1 0, 2 0.5 0.5)), ((0 0 0, 0 0 1, 0 1 1, 0 1 0, 0 0 0)))'::geometry) );

[pramsey]

GEOS does not do 3DIntersects, I'm pretty sure it's just a wrapper on 3DDistance == 0 which is 100% native to PostGIS. Not that I really want to try and fix it

[robe]

As Loic mentioned, I think at least for a basic fix, bring back the old SFCGAL ST_3DIntersects, but named CG_3DIntersects or whatever sweet prefix we want to come up with for SFCGAL functions where the two backends could overlap in purpose.