Potential bug in the ST_3DIntersects function.

[cuteDen-ECNU]

Potential bug in the ST_3DIntersects function.

29/9/2023

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 =\begin{pmatrix} -0.9248934130614574 & 0.2534855979991293 & 0.2834029394388001 \ -0.2259466894528656 & -0.9658860720787614 & 0.12653927963108094 \ 0.3058108369577812 & 0.05300139027692074 & 0.9506158975253333 \end{pmatrix}$

$M*M^T = E$

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;

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

[joyemang33]

How about using this equation $MM^{T} = I$ and telling them $M$ is an orthogonal matrix? I think it may be more clear. A visualized figure for Line a1 and Line a2 will also be good I think.

[cuteDen-ECNU]

bugtracker link: https://trac.osgeo.org/postgis/ticket/5557