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 for 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.

[pramsey]

This problem is not gathering a lot of effort… my install lacks CGAL and support for ST_MakeSolid(), and if I run it without that, it returns an answer.

pramsey=# 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)))'::geometry
);

 st_3dintersects 
-----------------
 f

If you have a patch, I would be happy to look at it. I'd suggest at a minimum get the offending code into a debugger and see where it's hung up.

[robe]

This wasn't an sfcgal issue, it wa a PostGIS one. ST_3DIntersects (exposed by SFCGAL was removed in PostGIS 3.0)

The answer on my setup with or without sfcgal installed returns false. So I guess the issue has been fixed in PostGIS proper at least to not hang. The only issue that remains is ST_3DIntersects I don't think is Solid aware.

I think we might have a ticket somewhere to reexpose the SFCGAL ST_3DIntersects which is solid aware. If not we should create one.

The ST_3DIntersection function, which IS postgis_sfcgal, seems to return the point fine.

SELECT ST_AsText(ST_3DIntersection('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)))'))

outputs

POINT Z (0 0.5 0.5)

[robe]

Come to think of it, this might have been a GEOS issue. I don't have GEOS 3.8.3 readily handy at the moment. I was testing against GEOS 3.12.

[robe]

Ticketed to reexpose at #5405