ST_Intersects() precision error - returns false instead of true

[cdestigter]

This query returns false, but the box and point intersect:

select st_intersects('POINT(169.69960846592 -46.5061209281002)'::geometry, 'POLYGON((169.699607857174 -46.5061218662,169.699607857174 -46.5061195965597,169.699608806526 -46.5061195965597,169.699608806526 -46.5061218662,169.699607857174 -46.5061218662))'::geometry);
 st_intersects 
---------------
 f
(1 row)

More info:

st_relate | st_intersects | st_disjoint
-----------+---------------+-------------
 0FFFFF212 | f             | f

ST_Relate is correct, ST_Disjoint is correct, ST_Intersects is _wrong_

  I B E
  -----
 I 0 F F
 B F F F
 E 2 1 2

Tested on postgis 1.5.2. postgis_version(): 1.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

Using geos via python it works:

>>> from django.contrib.gis import geos
>>> p = geos.GEOSGeometry('POLYGON((169.699607857174 -46.5061218662,169.699607857174 -46.5061195965597,169.699608806526 -46.5061195965597,169.699608806526 -46.5061218662,169.699607857174 -46.5061218662))')
>>> pt = geos.GEOSGeometry('POINT(169.69960846592 -46.5061209281002)')
>>> pt.intersects(p)
True

Probably a precision error, since this polygon is very small.

​postgis-users thread

[strk]

Not a precision issue, or ST_Relate and ST_Disjoint would also be wrong. Sounds more like a conceptual bug like a shortcut or something. Quick, easy test:

select (a && b) as "&&", st_relate(a, b), st_intersects(a, b), _st_intersects(a,b), st_disjoint(a,b) from
(select
'POINT(169.69960846592 -46.5061209281002)'::geometry
as a,
'POLYGON((169.699607857174 -46.5061218662,169.699607857174
-46.5061195965597,169.699608806526 -46.5061195965597,169.699608806526
-46.5061218662,169.699607857174 -46.5061218662))'::geometry
as b) as foo;

[strk]

point_in_polygon returned -1

[strk]

Distance between the point and the polygon boundary: 3.40605993187637e-07 Being a different path (non-GEOS) the precision issue _might_ indeed have a role here.

[strk]

The FP_CONTAINS_BOTTOM macros always return false, with the consequence of winding number never changing. The computed side (determineSide) is correct.

Nicklas: is that your code ?

[nicklas]

No, It's not.

In distance functions the pt_in_ring function is used for point in polygon tests. I guess that is why ST_Distance returns 0.0. pt_in_ring is a stab-line algorithm implementation that is old. So we also have a winding number algorithm implementation.

Is there a reason to have two different? I have read somewhere that the winding number algorithm is more stable (when bugs are fixed) and the same performance.

In 3D I have rewritten pt_in_ring to project the polygon to 2D and then use stabline algorithm. Maybe it is time top clean up and get winding algorithm working and move to that also in distance?

What is the name on that function?

/Nicklas

[nicklas]

Ok, I see that the name was point in ring.

So we have pt_in_ring and point_in_ring

As far as I understand they are supposed to do the same thing.

here is a good discussion about the two approaches:

​http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm

/Nicklas

[strk]

I understand point_in_ring is the winding one. Do you confirm pt_in_ring is the crossing one instead ?

Surely GEOS is using the crossing algorithm.

There's no unit testing for pt_in_ring, nor for point_in_ring. Also, point_in_ring is still in postgis/ rather than liblwgeom/

[mcayland]

Yes - please can we pick just one algorithm, move it to liblwgeom, fix up the callers and regress it for trunk?

[strk]

I think whoever introduced the new algorithm should take care of that.

[nicklas]

Which one is the new algoritm?

[nicklas]

strk, you are right, pt_in_ring is counting crossings.

As I understand now we have different point in polygon code in 4 places

1) point_in_ring (winding number) 2) pt_in_ring (crossing number) 3) pt_in_ring_3d (winding number) 4) There is also some implementation for geography somewhere (crossing number I think from a discussion some time ago)

The geography implementation have to be separate of course.

All the three others can be more or less merged.

I have not studied winding number enough to understand it. But at the softsurfer link above that is what they recommend.

My reason for rewriting the iterations in 3d-version even if it is the same thing happening as 2d-version is that it gets quite messy if I project the polygon to, say x-z plane instead of x-y since x and y is hard coded in the functions. If we can write the function so it just takes coord1 and coord2 instead of x and y it can be used more widely.

How is that done in a robust way.

In other words, is this possible to do in some efficient robust way?

coord1=x; pt.coord1;

instead of pt.x;

a macro can do that?

Will it be to ugly (unefficient) or is it worth it to share more code?

/Nicklas

[strk]

Let's try not to make things even more confusing than they are already. First of all I'd move _all_ implementations in a single file, under liblwgeom.

lwalgorithm.{c,h} would seem a sensible candidate for this.

I've seen pt_in_ring in liblwgeom/measures.c and point_in_ring in postgis/lwgeom_functions_analitic.c ... no idea where the others are.

I was assuming point_in_ring was the new one (but could be wrong)

[strk]

I've found this interesting comment:

commit ac4fc4ee809bfcc5247f67a776d5deb7ac5fd768
Author: mleslie <mleslie@b70326c6-7e19-0410-871a-916f4a2858ee>
Date:   Thu Aug 2 19:57:38 2007 +0000

    Removed a call to the deprecated point_in_ring function.
    
    git-svn-id: http://svn.osgeo.org/postgis/trunk@2677 b70326c6-7e19-0410-871a-916f4a2858ee

I'm hoping it referred to the same point_in_ring, as that would mean we should be focusing on the one in liblwgeom already...

[strk]

Forget previous comment:

-                if(point_in_ring_deprecated(root[i], &pt) == 1)
+                if(point_in_ring(root[i], &pt) == 1)

[strk]

pt_in_ring_2d is surely pretty old, anyway:

commit db0ccc046165153cc1f2db84775d2f8f483aa541
Author: strk <strk@b70326c6-7e19-0410-871a-916f4a2858ee>
Date:   Mon Nov 29 16:37:12 2004 +0000

    Fixed a bug in pt_in_ring_2d.
    
    
    git-svn-id: http://svn.osgeo.org/postgis/trunk@1123 b70326c6-7e19-0410-871a-916f4a2858ee

[strk]

I've tried switching point_in_polygon to use the crossing pt_in_ring_2d and got failures with within and contains:

--- regress_ogc_expected        2011-01-04 09:54:58.000000000 +0100
+++ /tmp/pgis_reg_12408/test_39_out     2011-02-24 14:04:39.000000000 +0100
@@ -10,9 +10,9 @@
 crosses|f
 crosses|t
 within100|t
-within101|f
+within101|t
 within102|f
-within103|f
+within103|t
 within104|f
 within105|t
 within106|t
@@ -45,9 +45,9 @@
 intersects155|t
 intersects156|f
 contains100|t
-contains101|f
+contains101|t
 contains102|f
-contains103|f
+contains103|t
 contains104|f
 contains105|t
 contains106|t

Those tests might need a double-check for correctness, but anyway both 'within' and 'contains' make a distinction between _inside_ and _on_the_boundary_, which pt_in_ring_2d doesnt' do.

[strk]

contains101, contains103, within101 and within103 are indeed all point-on-boundary cases.

[nicklas]

Ok, does this mean that ST_Within, ST_Contains, ST_Intersects and so on uses native PostGIS functions for simplier cases where a point is one of the geometries?

As you say that must be the reason that the new point in polygon is implemented to make that distinction.

So, the wisest probably is to sit down and get a grip of the winding algorithm and review it. As discussed in the link above there is significant performance differences depending on how the counting winding is implemented.

It would be nice to have the fastest and most stable way of doing it of course.

I think I can look at it in a few days, if it is not clear to someone else what should be done, already now.

And I suggest we leave 3d out of this. Maybe it should be rewritten to the same algorithm but I think it has to be separate code.

/Nicklas

[strk]

Yes, within and contains (and intersects) all have a short-circuit using the postgis algorithm when one operand is a polygon (or multipolygon) and the other is a point.

Given the fact the same operations would use PreparedGeometries I'm wondering if that shortcut has any real advantage... After all the same operation against prepared geometries should be faster (as the prepared polygon would be indexed).

[strk]

about the boundary vs. interior: if the pt_in_ring_2d is really currently unused, we might also fix that to have the -1,0,+1 returns. It'd be trivial to support the boundary case.

[nicklas]

Unused?

No, it is used in 2d distance calculations.

That is no big problem, but all cases (just a very few) have to be rewritten from

if (pt_in_ring_2d(..., ...))

to

if (pt_in_ring_2d(..., ...) >= 0)

or something like that.

[strk]

Maybe using 2 for the boundary would help reducing the rewrite needs

[nicklas]

Yes, and also changing "outside polygon" from -1 to 0, but I think that will be more rewriting in within, intersects and som on than doing the changes in distance. But it looks like the algoritm in counting windings actually returns 0 for outside polygon and then it is chenged:

01168         if (wn == 0)
01169                 return -1;
01170         return 1;

But, don't you think from what yo said earlier that it maybe is better to just send all of it to GEOS? At least until maybe we get some sort of indexing of vertex-points in PostGIS like Paul has sketched. If that things gets implemented I guess at east ST_Intersects maybe can be a very fast native PostGIS function. Maybe also contains and Within, I don't know.

Is there some link describing the prepared geometries? How, are they cached?

/Nicklas

[strk]

It'd take some profiling to really tell speed differences. But surely it's much easier to disable the short-cut than to do whatever we proposed so far. If not else it'd give us the opportunity to put something in the regression testsuite for this bug, so it's easier to optimize later w/out breaking it again.

[strk]

btw, is this confirmed to happen againt 1.5 ? I only tried againt trunk.

[nicklas]

I see in doxygen that FP_CONTAINS_BOTTOM also is referenced by point_in_ring_rtree(). When is that used?

In doxygen I can not follow it further than point_in_polygon_rtree, which doesn't seems to be referenced by any function. Or at least no c-function in PostGIS I guess it means?

/Nicklas

[strk]

It's in point_in_ring, around line 1143

FP_CONTAINS_BOTTOM(A, X, B) supposedly checks that A <= X < B, within tolerance of 1e-12

[strk]

NOTICE: FP_CONTAINS_BOTTOM(-46.506119596559699, -46.506120928100202, -46.506121866199997):0

The above seems to be wrong answer, doesn't it ? I see no fabs() involved in FP_CONTAINS_BOTTOM implementation either

[strk]

uhm, the answer to previous comment is not wrong indeed.

[strk]

This is interesting:

NOTICE:  Neither left rising nor rigth falling
NOTICE:  Neither left rising nor rigth falling
NOTICE:  < 1e-12...
NOTICE:  < 1e-12...

There's no check at all performed for the "within tolerance" cases

[strk]

It seems I've found the culprit. Testing in progress.

[strk]

So, this seems to be fixed with r6857. It should be back-ported to 1.5 (no idea about 1.4). I didn't update the NEWS file as if it'll get into next 1.5 it would not be a NEWS for 2.0.0...

Could anyone else take care of the backport and NEWS file please ?

@robe: history keyword added :)

[strk]

actally, I'll backport. I happen tohave that tree around

[nicklas]

Nice catch :-)

When we are here maybe we should also do the enhancement proposed in the link above and split the different tests and do them in the right order?

Now we check for side and upward/downward even if the segment is strict above or below the point which is meaningless.

The would also be more readable if we split it.

But that should only go for trunk of course.

/Nicklas

[nicklas]

and move it to liblwgeom and also use it for distance.

I'll make a new ticket for that.

/Nicklas

[strk]

+1 for moving to liblwgeom +1 for putting under unit testing

wouldn't change anything until things are regress-tested.

Oh, while porting to 1.5 I've seen the _same_ bug is in point_in_ring_tree. Should I blindly fix that as well or leave alone ? Did you mention it is unused ?

[nicklas]

Well, I don't know if it really is unused. It just has no links in dOxygen. I don't know if that means it is totally unused.

I think you can blindly do the same thing there. It really seems to have the same problem.

Can it be that this function is called when creating the spatial index. I mean pointed at from sql and referenced by postgresql? Then it wouldn't show in doxygen.

[strk]

I've committed r6858 in 1.5 branch, including the point_in_ring_rtree. I'll fix _rtree in trunk and call me out of this.

[strk]

r6859 fixes point_in_ring_rtree in trunk as well.

[cdestigter]

Wow - that was uber quick! Thanks guys :)