[postgis-users] Infinite loop in st_intersects - because of incorrect data out of st_transform?

[Mark Cave-Ayland]

On 28/02/11 15:18, Magnus Hagander wrote:

> Hi!
>
> Running the following query locks up postgis completely (in
> geos::algorithm::RobustDeterminant):
>
> SELECT st_intersects(somegeometry,
> '0103000020E61000000100000005000000737979F3DDCC2CC0F92154F9E7534540000000000000F07F000000000000F07F8F806E993F7E55C0304B29FFEA8554400634E8D1DD424540B5FEE6A37FCD4540737979F3DDCC2CC0F92154F9E7534540'::geometry)
>
> I believe this is because there are infinite values in that geometry:
>
> # select ST_AsText('0103000020E61000000100000005000000737979F3DDCC2CC0F92154F9E7534540000000000000F07F000000000000F07F8F806E993F7E55C0304B29FFEA8554400634E8D1DD424540B5FEE6A37FCD4540737979F3DDCC2CC0F92154F9E7534540'::geometry);
>
>    st_astext
> ---------------------------------------------------------------------------------------------------------------------------------------------------------------
>   POLYGON((-14.4001308522972 42.6555167828373,inf inf,-85.9726317957995
> 82.0924680617579,42.5223944076352 43.6054577711015,-14.4001308522972
> 42.6555167828373))
> (1 row)
>
>
> ISTM that this should either be rejected as an invalid geometry, or at
> least not hang....

Hi Magnus,

Hmmmm - I can definitely reproduce this on trunk with GEOS 3.2 series. 
The backtrace from inside GEOS looks like this:

(gdb) bt
#0  0x00007fb17c672525 in 
geos::algorithm::RobustDeterminant::signOfDet2x2 
(x1=-nan(0x8000000000000), y1=-nan(0x8000000000000), 
x2=-nan(0x8000000000000), y2=-nan(0x8000000000000))
     at RobustDeterminant.cpp:175
#1  0x00007fb17c6670a0 in geos::algorithm::CGAlgorithms::isCCW 
(ring=<value optimized out>) at CGAlgorithms.cpp:163
#2  0x00007fb17c6a0dae in geos::geomgraph::GeometryGraph::addPolygonRing 
(this=0x1a27970, lr=0x1a26560, cwLeft=2, cwRight=0) at GeometryGraph.cpp:262
#3  0x00007fb17c6a0f8e in geos::geomgraph::GeometryGraph::addPolygon 
(this=0x1a27970, p=0x1a265e0) at GeometryGraph.cpp:288
#4  0x00007fb17c6a16b5 in geos::geomgraph::GeometryGraph::add 
(this=0x1a27970, g=0x1a265e0) at GeometryGraph.cpp:183
#5  0x00007fb17c6a1a54 in GeometryGraph (this=0x1a27970, newArgIndex=1, 
newParentGeom=0x1a265e0, bnr=...) at GeometryGraph.cpp:522
#6  0x00007fb17c6d74bb in GeometryGraphOperation (this=<value optimized 
out>, g0=0x1a26480, g1=0x1a265e0) at GeometryGraphOperation.cpp:59
#7  0x00007fb17c701e16 in RelateOp (this=0xfff80000, g0=0xfff80000, 
g1=0x400) at RelateOp.cpp:56
#8  0x00007fb17c701e7f in geos::operation::relate::RelateOp::relate 
(a=<value optimized out>, b=<value optimized out>) at RelateOp.cpp:42
#9  0x00007fb17c67e549 in geos::geom::Geometry::intersects 
(this=0x1a26480, g=0x1a265e0) at Geometry.cpp:371
#10 0x00007fb17cefa047 in GEOSIntersects_r (extHandle=0x1a1d480, 
g1=0x7ff, g2=0x400) at geos_ts_c.cpp:284


Obviously the robust determinant is not as robust as it's name would 
suggest ;)  AFAICT taking a quick look at the GEOS code, I don't think 
it is able to detect the NaNs within x1/y1/x2/y2 as an exit condition 
for the main loop, but someone more familiar with the GEOS code should 
verify this.

> Now, this geometry actually comes back as the result of st_transform,
> which suggests to me that there is *also* a bug in st_transform, which
> should never output an invalid geometry. It is the result of:
>
> st_transform('0103000020CD0B0000010000000500000099E6673CA2FC2DC1AD7BF5ED45CC534199E6673CA2FC2DC1D7BDFAF6D28960415A06E670D7E94B41D7BDFAF6D28960415A06E670D7E94B41AD7BF5ED45CC534199E6673CA2FC2DC1AD7BF5ED45CC5341',
> 4326)
>
> And that geometry is valid before it goes into st_transform, from what
> I can tell:
> # select st_astext('0103000020CD0B0000010000000500000099E6673CA2FC2DC1AD7BF5ED45CC534199E6673CA2FC2DC1D7BDFAF6D28960415A06E670D7E94B41D7BDFAF6D28960415A06E670D7E94B41AD7BF5ED45CC534199E6673CA2FC2DC1AD7BF5ED45CC5341'::geometry);
>
>          st_astext
>
> ----------------------------------------------------------------------------------------------------------------------------------
> ------------------------------------------
>   POLYGON((-982609.1179802 5189911.7181081,-982609.1179802
> 8670871.7181081,3658670.8820198 8670871.7181081,3658670.8820198
> 5189911.
> 7181081,-982609.1179802 5189911.7181081))
> (1 row)
>
>
>
> To me this looks like bugs "below" postgis - it's certainly below the
> PostgreSQL level that I know well :-) Not sure if it can be worked
> around/fixed at the postgis level or if it needs to be done in the
> lower level libraries, but either way it shouldn't behave this way :-)

Yeah. The second question here is how did the (Inf Inf) point get into 
PostGIS anyway, which I think is a separate bug related to this.

Looking here 
(http://trac.osgeo.org/postgis/browser/trunk/postgis/lwgeom_transform.c) 
at transform_point(), if PROJ.4 returns an error then *pj_errno_ref 
should be set to non-zero but even though the projection result is a 
point at infinity, the error code is still returning success.

I think the key thing to knowing whether this part is a PROJ.4 bug or a 
PostGIS bug would be to clarify whether PROJ.4's *pj_errno_ref should be 
non-zero if a projected point appears at infinity - if yes, it's a 
PROJ.4 bug and if no, it's a PostGIS bug because we should detect this 
condition and fail the transformation with a suitable error. Frank - 
what's your take on this?


ATB,

Mark.

-- 
Mark Cave-Ayland - Senior Technical Architect
PostgreSQL - PostGIS
Sirius Corporation plc - control through freedom
http://www.siriusit.co.uk
t: +44 870 608 0063

Sirius Labs: http://www.siriusit.co.uk/labs