ST_Intersection give results that are inconsistent with ST_Distance

[nbvfgh]

SELECT ST_AsText(ST_Intersection(geom.b, geom.a))AS intersection 
    FROM (SELECT ST_GeomFromText('POLYGON ((2 0, 0 2, 0 0, 2 0))') As a,
    ST_GeomFromText('LINESTRING(0 3, 0.4 1.6, 2 3)') As b) As geom;

-- result : {LINESTRING EMPTY}

And when calculating the intersection of POINT (0.4 1.6) and LINESTRING (2 0, 0 2), we obtain an empty set, but POINT (0.4 1.6) is common sense on LINESTRING (2 0, 0 2), but the intersection result is an empty set.

SELECT ST_AsText(ST_Intersection(geom.b, geom.a))AS intersection 
    FROM (SELECT ST_GeomFromText('LINESTRING(2 0, 0 2)') As a,
    ST_GeomFromText('POINT(0.4 1.6)') As b) As geom;

-- expected: {POINT(0.4 1.6)}
-- actual  : {POINT EMPTY}

Then I made the following attempts:

SELECT ST_AsText(ST_Intersection(geom.b, geom.a))AS intersection 
    FROM (SELECT ST_GeomFromText('POLYGON ((2 0, 0 2, 2 2, 2 0))') As a,
    ST_GeomFromText('POINT(0.4 1.6)') As b) As geom;

-- result : {POINT(0.4 1.6)}

SELECT ST_AsText(ST_Intersection(geom.b, geom.a))AS intersection 
    FROM (SELECT ST_GeomFromText('POLYGON ((2 0, 0 2, 0 0, 2 0))') As a,
    ST_GeomFromText('POINT(0.4 1.6)') As b) As geom;

-- result : {POINT EMPTY}

But when calculating their distances separately using ST_Distance, the result obtained is 0.

SELECT ST_Distance(geom.a, geom.c), ST_Distance(geom.b, geom.c)
    FROM (SELECT ST_GeomFromText('POLYGON ((2 0, 0 2, 0 0, 2 0))') As a,
    ST_GeomFromText('LINESTRING(0 3, 0.4 1.6, 2 3)') As b,
    ST_GeomFromText('POINT(0.4 1.6)') As c) As geom;

-- result: {0, 0}

Version Info：GEOS="3.13.0-CAPI-1.19.0"

[nbvfgh]

This seems to be caused by the inability to accurately represent geometric objects

[dr-jts]

This seems to be caused by the inability to accurately represent geometric objects

Correct. Note that 0.4 and 1.6 can't be represented as exact binary fractions (due to the presence of a factor of 5). Unfortunately, the extended-precision arithmetic used for the point-line orientation predicates is a little too exact in this case.

One possible fix is to compute the point-line orientation using regular FP math first, and if the result is exact, use that result. Only when the point is not on the line is the extended-precision math used (to be more precise about which side of the line the point lies on).

In the future the topological predicates may support a distance tolerance, which would solve this problem as well.

[dr-jts]

Here's a repro in GEOS:

geosop -a "POLYGON ((2 0, 0 2, 0 0, 2 0))" -b "POINT(0.4 1.6)" intersection

POINT EMPTY

[nbvfgh]

This seems to be caused by the inability to accurately represent geometric objects

Correct. Note that 0.4 and 1.6 can't be represented as exact binary fractions (due to the presence of a factor of 5). Unfortunately, the extended-precision arithmetic used for the point-line orientation predicates is a little too exact in this case.

One possible fix is to compute the point-line orientation using regular FP math first, and if the result is exact, use that result. Only when the point is not on the line is the extended-precision math used (to be more precise about which side of the line the point lies on).

In the future the topological predicates may support a distance tolerance, which would solve this problem as well.

Yes, using a distance tolerance can effectively avoid such problems caused by discretization, but the tolerance needs to be controlled within an appropriate range.
Thank you for your detailed explaination！