[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
Wed Nov 6 06:11:10 PST 2013
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/20131106/0d908002/attachment-0001.html>
More information about the geos-devel
mailing list