The use of tolerance in function point_in_ring is not robust.

[nicklas]

I have taken a closer look at function point_in_ring because of my task in #846.

There is two independent problems. the row numbers references the file lwgeom_functions_analytic.c

1) on line 1133 we find this

        /* a point on the boundary of a ring is not contained. */
                if (fabs(side) < 1e-12)
                {

the problem here is that the calculated values "side" is not only related to the distance between the point and the segment but also to the length of the segment. That results in that the first example here returns false and the second true even if the point is on the exact same distance from the polygon in both cases.

select st_intersects('POINT(0.5 0.5000000000001)'::geometry, 'POLYGON((0 0, 10 10, 1 0, 0 0))'::geometry);


select st_intersects('POINT(0.5 0.500000000001)'::geometry, 'POLYGON((0 0, 1 1, 1 0, 0 0))'::geometry);

2) the second problem is the construct to first check on which side of the segment the point is with a tolerance and then deciding if the point is within the endpoints without a tolerance. The second test will make the tolerance worthless if the segment is vertical or horizontal. That is actually the reason why the bug #845 showed.

there is also a third problem with this tolerance but that doesn't cause anything wrong I think after twisting it around many times. But the use of tolerance in the macros in the begining of lilwgeom.h that looks like this

#define FP_TOLERANCE 1e-12
#define FP_LT(A, B) (((A) + FP_TOLERANCE) < (B))
#define FP_LTEQ(A, B) (((A) - FP_TOLERANCE) <= (B))
#define FP_CONTAINS_BOTTOM(A, X, B) (FP_LTEQ(A, X) && FP_LT(X, B))

makes no difference in the result in any case. It causes only a small shift in the whole calculation. If it would have made any difference I am quite sure it would also have been one or many cases where the tolerance would have caused wrong answer. But now I think it is only confusing.

The last issue I would like to just remove the tolerance part, but I suspect it is there for some reason. Does anyone remember? the regression tests passes on my computer if I remove it.

The two first problems is worse. I guess we need some sort of tolerance to be sure to catch the boundary-cases. But this is not the right way, I think I have found. Both because it gives an unpredictable tolerance as shown in 1) and second because there is no tolerance at all in vertical and horizontal segments.

I attach an image showing the point and polygon from #845 to illustrate. The distance betwen the point and the left boundary is about 3e-7 m.

/Nicklas

[strk]

Ouch, this one hurts, we want it fixed in 2.0 I guess, and also in 1.5 if affected. Is this a GEOS/JTS bug or all internal ? Why are we using internal (if we do) for _ST_Intersects ?

[strk]

Ok, I checked, it is PostGIS implementation, not GEOS's.

GEOS correctly returns FALSE for both cases. The point is slightly above the polygon.

Any reason not to disable the short-circuit ?

[nicklas]

By removing, do you mean removing all the _rtree stuff as chodgson suggests in #846? He states that the point in poly test for _st_intersects is necessary to save for non cached geometries. If so, the problem will not go away, if I remember right. I have not looked at the code now.

But I guess geos could take care of non cached cases too, then we maybe can just get rid of the whole problem. And closing #846 as a bonus :-)

But is it only _ST_Intersects that is using those functions? I should take a look in doxygen, but not now, sorry.

[strk]

I mean removing the short-circuit in ST_Intersects:

         * short-circuit 2: if the geoms are a point and a polygon,
         * call the point_outside_polygon function.

Which would solve the bug exposed by the st_intersects calls in the original description, thus giving us at least a working testcase.

I don't see any other use for point_in_ring from the ST_Intersects path.

There is a failing testcase under raster when I try, but it isn't easy for me to tell if it is a false positive (ie: a bogus testcase). Maybe raster people may take a look ?

--- rt_spatial_relationship_expected    2011-11-25 17:13:23.000000000 +0100
+++ /tmp/pgis_reg_7435/test_68_out      2012-01-19 09:18:29.000000000 +0100
@@ -18,7 +18,7 @@
 test 1.1|2|33|f
 test 1.1|2|34|t
 test 1.1|2|35|t
-test 1.1|2|36|t
+test 1.1|2|36|f
 test 1.1|2|41|f
 test 1.1|2|42|t
 test 1.1|2|43|t
@@ -58,7 +58,7 @@
 test 1.3|2|33|f
 test 1.3|2|34|t
 test 1.3|2|35|t
-test 1.3|2|36|t
+test 1.3|2|36|f
 test 1.3|2|41|f
 test 1.3|2|42|t
 test 1.3|2|43|t
@@ -153,6 +153,7 @@
 test 2.3|2|33|GEOMETRYCOLLECTION EMPTY|
 test 2.3|2|34|POINT(-75.5518328537098 49.2814585505576)|4
 test 2.3|2|35|POINT(-75.5523150760919 49.2818643977074)|4
+test 2.3|2|36|GEOMETRYCOLLECTION EMPTY|
 test 2.3|2|41|GEOMETRYCOLLECTION EMPTY|
 test 2.3|2|42|LINESTRING(-75.5533120424461 49.2823793618213,-75.5531972042327 49.2824942000347)|1
 test 2.3|2|43|LINESTRING(-75.5529884291587 49.2811479840607,-75.5522356128797 49.2815620330142)|3

[strk]

Note: the case in the original description does _not_ hit the cache, thus calls point_in_polygon which gets to point_in_ring.

[strk]

Alright, I removed the tolerance with r8873. Nothing fails and we now have the original case under regression testing (partial,as I couldn't handle to trigger the cached version from test yet).

Not sure we want to backport this to 1.5

[strk]

r8874 improves the testcase to also hit the caches.

[strk]

Backported in 1.5 as r8875. Nothing more to do here, moving to #846