ST_3DIntersects says 'no' despite geometries intersecting

[tilt]

This ticket is a follow up on a question on stackexchgange: ​https://gis.stackexchange.com/questions/293619/st-3dintersects-says-no-despite-geometries-intersecting?noredirect=1#comment471093_293619

Considering the following query I expect both geomb and geomc to 3D Intersect with geoma. According to postgis, only geomc intersects. Unfortunately I can't come up with a more simplified polygon as an example, this situation seems largely to depend on the configuration of the polygon.

WITH data AS (
SELECT 
ST_SetSrid(ST_GeometryFromText('
POLYGON Z ((
    122409.416 489410.194 6.64975853832649,
    122404.759 489412.266 3.76021877729423,
    122407.18 489417.708 3.82360012934401,
    122405.377 489418.511 2.70471051793069,
    122406.754 489421.606 2.74081838637221,
    122408.218 489420.954 3.64933000361367,
    122409.263 489423.303 3.67668436361495,
    122407.787 489423.959 2.76102260948305,
    122409.313 489427.389 2.80101513611236,
    122411.091 489426.598 3.90419939305872,
    122411.172 489426.781 3.9061156669998,
    122409.45 489427.547 2.83769648670792,
    122410.905 489430.819 2.87547799988339,
    122412.594 489430.067 3.92357220158512,
    122413.511 489432.131 3.94697470225165,
    122411.723 489432.926 2.83768486408258,
    122413.229 489436.312 2.87694266543155,
    122414.931 489435.555 3.93292998257874,
    122415.345 489436.484 3.94412204281057,
    122420.157 489434.343 6.9298432176791,
    122409.416 489410.194 6.64975853832649
))
'),28992) as geoma,


ST_Extrude(ST_SetSrid(ST_GeometryFromText('
LINESTRING Z (
    122412.982 489429.097 0,
    122414.559 489432.685 0
)
'),28992),0,0,10) as geomb,
ST_Extrude(ST_SetSrid(ST_GeometryFromText('
LINESTRING Z (
    122414.559 489432.685 0,
    122412.756 489433.477 0
)
'),28992),0,0,10) as geomc
)


SELECT 
ST_3DIntersects(geoma,geomb) ab,
ST_3DIntersects(geoma,geomc) ac
FROM data;

Results in: f;t

I guess I should look at sfcgal for the culprit.

Please note: the given geoma above is non-planar because of rounding when exporting to text. In the database it shows ST_IsPlanar(geoma) as true and gives the same results. On top of that, in the database running ST_3DIntersection on geoma vs. geomb and geoma vs. geomc both give a good result. This is not replicable here because of non-planarity.

The situation is perhaps better described in the attached picture, where geoma is blue, geomb is red and geomc is green. The geometries in the image have been Tesselated for display purpose only.

[tilt]

Picture showing how the planes are intersecting in 3D space

[tilt]

For completeness:

SELECT postgis_full_version();
"POSTGIS="2.5.0 r16836" [EXTENSION] PGSQL="100" GEOS="3.7.0-CAPI-1.11.0 673b9939" SFCGAL="1.3.5" PROJ="Rel. 5.2.0, September 15th, 2018" GDAL="GDAL 2.4.0dev-4b763dd896-dirty, released 2018/11/19" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" RASTER"

[komzpa]

ST_3DIntersects in GEOS backend is a plug for ST_3DDWithin(geom, geom, 0):

​https://github.com/postgis/postgis/blob/dd281b1e57b590c0339271cedbf5feaf84049cc4/postgis/lwgeom_backend_api.c#L153

We are going to remove it during #4258 and (probably) rename to ST_3DSurfaceIntersects. May it happen that surfaces of ac intersect, and surfaces of ab don't? If yes, then it works as expected codewise, and we'll untangle this in 3.0 by renaming function.

[tilt]

No. They both have intersecting surfaces (see attached pic). Also, I didn't manage to isolate an easier example (out of hundreds of intersections) than this one making me think it is due to some exception that is not supposed to happen.

[nicklas]

You are running this function with geos backend I goes (which means running it native PostGIS in this case).

If you try with SFCGal backend it gives a hint about the problem:

NOTICE:  During intersects_3d(A,B) :
NOTICE:    with A: POLYGON((4205955507539673/34359738368 4204001555118031/8589934592 7486962518827739/1125899906842624,8411590988476187/68719476736 8408038706925011/17179869184 4233629971063459/1125899906842624,2102939339582341/17179869184 4204066099886555/8589934
NOTICE:     and B: POLYHEDRALSURFACE(((8412156068733387/68719476736 8408327861303247/17179869184 0/1,8412264439348199/68719476736 8408389502673879/17179869184 0/1,8412264439348199/68719476736 8408389502673879/17179869184 10/1,8412156068733387/68719476736 840832786


ERROR:  Polygon is invalid : points don't lie in the same plane : POLYGON((4205955507539673/34359738368 4204001555118031/8589934592 7486962518827739/1125899906842624,8411590988476187/68719476736 8408038706925011/17179869184 4233629971063459/1125899906842624,21029

********** Error **********

ERROR: Polygon is invalid : points don't lie in the same plane : POLYGON((4205955507539673/34359738368 4204001555118031/8589934592 7486962518827739/1125899906842624,8411590988476187/68719476736 8408038706925011/17179869184 4233629971063459/1125899906842624,21029
SQL state: XX000

There is no real check in the postgis code if the polygon is coplanar. But it tries to calculate an "average plane". From the illustration you attached I think that calculation should manage this situation, but it is a little bit difficult so see what the polygons really looks like (and I have no good viewer by hand). But if for instance some part of the polygon is wrapped on the backside of the erst of the polygon this "average calculation" will give a very strange answer and the calculation will be totally wrong.

[tilt]

The non-planarity is due to rounding effects from the ST_AsText (see the 'note' in the ticket). The original polygons are all planar.

You might like to use ​https://github.com/geondev/PostGIS3DExplorer for visualizing.

[nicklas]

Ok, I will take a look

PostGIS3DExplorer seems unusable for Debian though :-)

[nicklas]

Sorry for the delay. It took me some time to get the time to understand all this again. It was many years ago.

I have found the problem. It is when defining the plane. The plane of a polygon is stored aas a point on the polygon and a point perpendicular to that point from the plane. To make the calculations robust it calculates this perpendicular point from some samples. But what happens in some cases is that the average point calculation gets zero because some points gets on the upper side of the plane and others on the down side. This can happen when the line is going both concave and convex as your did. So the simple fix is to add fabs to the plane definition.

But, just some weeks ago Algunenano changed the 3D calculation code quite a lot, in #4546 and then this problem is very much reduced by not calculating an averge point this way, and instead just add to the point. This means that the perpendular point no more defines a normal vector perpendular to the plane, but that doesn't matter. In theory it could still be possible to create a polygon that gets an undefined plane for the same reason as before, but the chance is minimal.

So I leave this without any changes.

[nicklas]

Fixed with #4246