[geos-devel] Re: [RE]: Intersection of point and linestring. still bad results...

[strk at refractions.net]

ksa, I've verified the problem locally, with postgis 0.9 and head 
branches and GEOS cvs.

The lost of precision is not due to textual conversion, but I'm still
convinced it is due to rounding errors.

Martin, basically the computed intersection point bewteen two
linestrings do not intersect any of the two.

This are a few intersting debugging lines:

 OverlayOp::computeOverlay: li:
 (-271.415,93.57,0)_(-187.14,74.3603,0) (-236.164, 88.6138)_(-238.349,24.5005)
 :  proper

The LineIntersector found proper intersection between those two segments:
POINT(-236.268122207233 85.558597112599)

 OverlayOp::labelIncompleteNode(node (-236.268,85.5585,0) lbl: a:- b:i, 0)
 RobustLineIntersector::computeIntersection((-236.268,85.5585,0), (-236.164,88.6138), (-238.349,24.5005))
  envelope intersects point
  orientationIndex p1p2p==1 p2p1p==-1 <---- CHECK THIS !
  after location set: node (-236.268,85.5585,0) lbl: a:e b:i

The PointLocator is failing to recognize it as 'internal' to geometry 'a':

 LINESTRING(-268.94877508 161.46993124,-268.94877508 161.46993124,-236.16352798 88.61382659,-238.34921112 24.5004545)

Whose last segment is exactly the one on which LineIntersector above
found the proper intersection point.

The orientationIndex computed by CGAlgorithm should be 0 for both p1p2p
and p2p1p for p to be considered ON p1-p2 segment, but it is 1 and -1
respectively.

The RobustDeterminant class uses ==, < and > checks which do not consider
floating point rounding errors. Maybe overriding these operators for
considering rounding errors would give the expected result, but I dunno
how this generalization would impact both speed and precision of GEOS.

Ideas Martin ?

--strk;

On Tue, Nov 30, 2004 at 11:08:50AM +0300, ksa wrote:
> Hi, strk.
> 
> search=# select gid, geom from t_table_46;
>  gid |
> geom
> -----+----------------------------------------------------------------------
> ----------------------------------------------------------------------------
> ----------------------------------------------------------------------------
> ------------------
>    1 |
> SRID=2000;MULTILINESTRING((-346.39495192 -32.22278924,-346.39495192 -32.2227
> 8924,-347.01462094 37.18014143,-261.50029565 23.54742291,-271.41500003
> 93.5700226,-187.14001279 74.36028286,-195.1957101 155.53692498,-118.97642017
> 141.90420646))
>    2 | SRID=2000;MULTILINESTRING((-268.94877508 161.46993124,-268.94877508
> 161.46993124,-236.16352798 88.61382659,-238.34921112 24.5004545))
> (2 rows)
> 
> 
> search=# select setsrid(intersection(t1.geom, t2.geom), 2000) from
> t_table_46 t1, t_table_46 t2 where t1.gid = 1 and t2.gid = 2;
>                       setsrid
> ----------------------------------------------------
>  SRID=2000;POINT(-236.267686750346 85.558502660642)
> (1 row)
> 
> So, I have one point that intersects with both multilinestrings.
> But I still have bad results...
> Pleace, help me...
> 
> 
> search=# select
> intersection('SRID=2000;MULTILINESTRING((-346.39495192 -32.22278924,-346.394
> 95192 -32.22278924,-347.01462094 37.18014143,-261.50029565
> 23.54742291,-271.41500003 93.5700226,-187.14001279 74.36028286,-195.1957101
> 155.53692498,-118.97642017 141.90420646))',
> 'SRID=2000;POINT(-236.267686750346 85.558502660642)');
>            intersection
> -----------------------------------
>  SRID=-1;GEOMETRYCOLLECTION(EMPTY)
> (1 row)
> 
> 
> search=# select
> intersection(transform('SRID=2000;MULTILINESTRING((-346.39495192 -32.2227892
> 4,-346.39495192 -32.22278924,-347.01462094 37.18014143,-261.50029565
> 23.54742291,-271.41500003 93.5700226,-187.14001279 74.36028286,-195.1957101
> 155.53692498,-118.97642017 141.90420646))', 2000),
> transform('SRID=2000;POINT(-236.267686750346 85.558502660642)', 2000));
>            intersection
> -----------------------------------
>  SRID=-1;GEOMETRYCOLLECTION(EMPTY)
> (1 row)
> 
> 
> And why intersection() returns geometry with -1 SRID? All parameters have
> srid=2000.
> 
> search=# select intersection(t1.geom, t2.geom) from t_table_46 t1,
> t_table_46 t2 where t1.gid = 1 and t2.gid = 2;
>                    intersection
> --------------------------------------------------
>  SRID=-1;POINT(-236.267686750346 85.558502660642)
> (1 row)
> 
> ksa.