[geos-devel] [postgis-devel] SNapPointToLine : GEOS : where is the code for st_intersects between line and point?
Rémi Cura
remi.cura at gmail.com
Thu Nov 7 05:36:26 PST 2013
Hey,
after some testing it appears that simply evaluating the = with precision
in mind won't suffice.
The core of the problem is that the determinant algorithm geos uses seems
to be designed to robustly say left or right of a line, but not robustly on
the line.
So an idea may be to offset the line to the left and the right with the
precision limit distance, and use the robust determinant to test if the
point is right of the left and left of the right (between).
Can I have some answers on this please? Is it doable, ?
Cheers,
Rémi-C
2013/11/6 Rémi Cura <remi.cura at gmail.com>
> Hey Nicklas,
>
> I agree tolerance is even more complicated than precision. This is a
> completely different matter as when defining tolerance you go into fuzzy
> logic.
>
> Here I'm just talking about all computing being "true" for the double
> precision world.
>
> Finding a point on a line with double precision is not a problem, you can
> just do computing with better precision and round.
> The problem is this point would be on the line (for a double world), but
> not on the line for st_intersects.
>
> I found experimentally that using st_intersect, you have point spacing
> following the prime decomposition of the slope , which is very bad news in
> most of case.
>
> The issue seems to be in robust_determinant computing, at a point, for a
> point M being on line P1 P2,
> you must have
> (x2-x1)*(ym-y1) = (y2-y1)*(xm-x1)
> by translating you can make x1=y1=0 with x2t and y2t being the x2 y2
> translated.
> x2t*ym = y2t * xm
> Now the "equal" should be evaluated with precision in mind :
> Ina typical case of a small slope : (slope = 0.001), we could take
> precision(x2t=1000m)= 10^(-15+4), precision(y2t=1m)= 10^(-15+1);
>
> In this case the "=" should be evaluated with a 10^(-11) precision (is it
> sufficient? too much??)
>
> Cheers,
> Rémi-C
>
>
>
> 2013/11/6 Nicklas Avén <nicklas.aven at jordogskog.no>
>
>> Hallo
>>
>> Sorry, I was wrong. It was only the point in polygon test that is done
>> in the PostGIS code, see file postgis/lwgeom_geos.c row 2602
>>
>> I thought I remembered that also point vs linestring was handled in the
>> postgis code.
>>
>> I think it was last time we discussed this problem that I created the
>> tolerance discussion wiki-page. But it is not very easy to handle
>> correct.
>>
>> Maybe the point should be projected to the line and see if the projected
>> point is equal to the unprojected point and then declare it
>> intersecting. But I doubt that is very robust when comparing projected
>> points from different machines on different platforms. We have seen
>> differences in the last decimal before.
>>
>>
>> Your idea, if I understand it right to find a point that actually is on
>> the line and is representable with double precission coorinates can
>> result in points far away from where it is supposed to be.
>>
>> I think that in most cases with not horizontal or vertical lines you
>> will just get the vertex points and nothing else.
>>
>> Representable values on not horizontal or vertical lines is actually
>> quite rare.
>>
>> Best Regards
>>
>> Nicklas Avén
>>
>>
>>
>> On Wed, 2013-11-06 at 11:15 +0100, Rémi Cura wrote:
>> > Hey Nicklas,
>> > thanks for your answer.
>> >
>> > I'm simply trying to solve a problem in current ST_Intersects behavior
>> > regarding a point and a line. IN most of the case, a point on a
>> > line(to double limit precision) is not on the line for st_intersects !
>> >
>> >
>> > it works as intended in very very few cases (ie, all the points that
>> > intersects a line is a very sparse set and is not a line at all !).
>> >
>> >
>> > Here is a simple illustration : http://hpics.li/c3485b5 represents a
>> > line with max zoom : it is all the point that follow the line equation
>> > and are rounded to the(16 or 17 digits)
>> > http://hpics.li/f30aeb7 now the big points represents the point that
>> > are on the line for ST_Intersect, in a favorable case (bad case could
>> > mean no points at all).
>> >
>> >
>> >
>> > This is without speaking of custom precision model, ie working with
>> > the postgis 'infinite' precision (double precision, 15 to 17 digits)
>> >
>> >
>> > I tried several workaround (like several random walk along the line to
>> > find a point that st_intersects like,
>> > and also reverse engineer to find which points are seen on the line by
>> > st_intersects (related to prime factorization of slope it appears) ).
>> >
>> >
>> > The issue is there is a flaw in ST_Intersects and in bad cases there
>> > is no workaround.
>> >
>> >
>> > Could you elaborate about 'native' functions please?
>> >
>> >
>> > Cheers,
>> >
>> > Rémi-C
>> >
>> >
>> >
>> >
>> > 2013/11/6 Nicklas Avén <nicklas.aven at jordogskog.no>
>> > I think postgis has native functions for st_intersects when at
>> > least one of the geometries is a point. And there is two paths
>> > in the code. One that prepares the "not point" geometry if it
>> > is to be used multiple times and one that doesn't.
>> >
>> >
>> > I do not get the idea do far. Is it a tolerance you are
>> > implementing?
>> >
>> >
>> >
>> >
>> >
>> >
>> > Best regards
>> > Nicklad
>> >
>> >
>> > Skickat från min Samsung Mobil
>> >
>> >
>> >
>> > -------- Originalmeddelande --------
>> > Från: Rémi Cura <remi.cura at gmail.com>
>> > Datum: 06-11-2013 8:59 (GMT+01:00)
>> > Till: PostGIS Development Discussion
>> > <postgis-devel at lists.osgeo.org>
>> > Rubrik: [postgis-devel] SNapPointToLine : GEOS : where is the
>> > code for st_intersects between line and point?
>> >
>> >
>> > Hey,
>> > I've been working on a ST_SnapPointToLine for the past week,
>> > by reverse engineering when points are considered to be on
>> > line by st_intersect.
>> >
>> >
>> > The (classic) problem is purely a numerical issue, and I
>> > suspect a problem in design of st_intersects.
>> >
>> >
>> >
>> >
>> > So far my solution doesn't work for all line, can someone
>> > please point me in the right direction in GEOS so I can see if
>> > something can be improved regarding precision.
>> >
>> >
>> > I know where the API is, but I'm looking for the exact part
>> > that decide if a point and a line intersects.
>> >
>> >
>> > Cheers,
>> >
>> >
>> > Rémi-C
>> >
>> > _______________________________________________
>> > postgis-devel mailing list
>> > postgis-devel at lists.osgeo.org
>> > http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel
>> >
>> >
>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geos-devel/attachments/20131107/313bbdad/attachment.html>
More information about the geos-devel
mailing list