Potential bug in the ST_3DIntersects function.

[Wenjing]

Potential bug in the ST_3DIntersects function.

​visualization report

Geometry Description

Line a('LINESTRING(5 4 0, 5 0 0)') intersects Line b ('LINESTRING(6 4 0, 4 2 0)') at Point c(’POINT(5 3 0)’).

We affine Line a and Line b using the same orthogonal matrix M and get Line a1 and Line b1 respectively.

M = [[-0.9248934130614574, 0.2534855979991293, 0.2834029394388001], 
[-0.2259466894528656, -0.9658860720787614, 0.12653927963108094 ],  
[0.3058108369577812, 0.05300139027692074, 0.9506158975253333]]

Expected Behavior

Line a1 and Line b1 should intersect at Point c1 (affine Point c using M) because we only transform the coordinates.

Point c1 is actually the intersection point of Line a1 and Line b1.

SELECT ST_Distance(Foo.c1, Foo.a1) As c1_to_a1, ST_Distance(Foo.c1, Foo.b1) As c1_to_b1 
    FROM (
    SELECT ST_GeomFromText('LINESTRING(5 4 0, 5 0 0)', 0) As a
    ,ST_Affine(ST_GeomFromText('LINESTRING(5 4 0, 5 0 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As a1
    ,ST_GeomFromText('LINESTRING(6 4 0, 4 2 0)', 0) As b
    ,ST_Affine(ST_GeomFromText('LINESTRING(6 4 0, 4 2 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As b1
    ,ST_GeomFromText('POINT(5 3 0)', 0) As c
    ,ST_Affine(ST_GeomFromText('POINT(5 3 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As c1
    ) As Foo;

c1_to_a1 |       c1_to_b1        
----------+-----------------------
        0 | 4.440892098500626e-16
(1 row)

Actual Behavior

PostGIS thinks that Line a and Line b intersect, and doesn’t think that Line a1 and Line b1 intersect.

SELECT ST_3DIntersects(Foo.a, Foo.b) As a_Intersects_b
, ST_3DIntersects(Foo.a1, Foo.b1) As a1_Intersects_b1
    FROM (
    SELECT ST_GeomFromText('LINESTRING(5 4 0, 5 0 0)', 0) As a
    ,ST_Affine(ST_GeomFromText('LINESTRING(5 4 0, 5 0 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As a1
    ,ST_GeomFromText('LINESTRING(6 4 0, 4 2 0)', 0) As b
    ,ST_Affine(ST_GeomFromText('LINESTRING(6 4 0, 4 2 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As b1
    ) As Foo;

-- SQL result:

a_intersects_b | a1_intersects_b1 
----------------+------------------
 t              | f
(1 row)

[mdavis]

The issue is due to the inherent inaccuracy of finite-precision floating point arithmetic. There's a lot of arithmetic being carried out, so the likelihood of 2 transformed intersecting lines NOT intersecting in 3D is very high.

You can't rely on topological relationships being preserved through floating-point operations. This is the essence of the well-known issue of geometrical robustness failures. The solution is usually either to use exact arithmetic (often impractical) or use tolerance values (not a panacea, but very often works).

[robe]

In regard to using tolerance, try using ST_3DWithin instead. You can treat it like 3DIntersects with tolerance.

SELECT ST_3DDWithin(Foo.a, Foo.b,0) As a_Intersects_b
, ST_3DDWithin(Foo.a1, Foo.b1,0.000001) As a1_Intersects_b1
    FROM (
    SELECT ST_GeomFromText('LINESTRING(5 4 0, 5 0 0)', 0) As a
    ,ST_Affine(ST_GeomFromText('LINESTRING(5 4 0, 5 0 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As a1
    ,ST_GeomFromText('LINESTRING(6 4 0, 4 2 0)', 0) As b
    ,ST_Affine(ST_GeomFromText('LINESTRING(6 4 0, 4 2 0)', 0), -0.9248934130614574, 0.2534855979991293, 0.2834029394388001, -0.2259466894528656, -0.9658860720787614, 0.12653927963108094, 0.3058108369577812, 0.05300139027692074, 0.9506158975253333, 0, 0, 0) As b1
    ) As Foo;

But realistically I don't think there is anything we can do to fix the ST_3DIntersects, we have similar issues in ST_Intersects. e.g. the ST_ClosestPoint on a 2d line almost never intersects the 2d line.

[robe]

Added reference to ST_3DDWithin on ST_3DDIntersects In f62ef66/git:

[strk]

Commit [f62ef66b699e9e85e5a6798f1118e571c6b4116c/git] was really for #5557, just had a wrong reference to this ticket. I've added a note to amend that: ​https://git.osgeo.org/gitea/postgis/postgis/commit/f62ef66b699e9e85e5a6798f1118e571c6b4116c