Opened 2 years ago

Closed 11 months ago

Last modified 11 months ago

#5144 closed defect (wontfix)

3d intersects in postgis hangs on simple polyhedron

Reported by: lbartoletti Owned by: pramsey
Priority: medium Milestone: PostGIS 3.4.0
Component: postgis Version: 3.2.x
Keywords: Cc:

Description

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

Change History (15)

comment:1 by lbartoletti, 2 years ago

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

comment:2 by lbartoletti, 2 years ago

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

comment:3 by lbartoletti, 2 years ago

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

comment:4 by lbartoletti, 2 years ago

Tested with geos and sfcgal backend.

Last edited 2 years ago by lbartoletti (previous) (diff)

comment:5 by lbartoletti, 2 years ago

Component: sfcgalpostgis
Owner: changed from robe to pramsey

comment:6 by robe, 2 years ago

Milestone: PostGIS SFCGALPostGIS 3.2.2

comment:7 by robe, 2 years ago

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) );

comment:8 by robe, 22 months ago

Milestone: PostGIS 3.2.2PostGIS 3.2.3

comment:9 by pramsey, 21 months ago

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 :)

comment:10 by robe, 21 months ago

Milestone: PostGIS 3.2.3PostGIS 3.4.0

comment:11 by robe, 21 months ago

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.

comment:12 by pramsey, 11 months ago

Resolution: wontfix
Status: newclosed

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.

comment:13 by robe, 11 months ago

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)

comment:14 by robe, 11 months ago

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.

comment:15 by robe, 11 months ago

Ticketed to reexpose at #5405

Note: See TracTickets for help on using tickets.