libgeos/libgeos — pull request #30 by Remi-C (2013-11-22 16:17:27+00:00)
Title: This upgrade the detection of an intersection between a line and a point

... by adding a tolerance over digits from 13 to 17

The problem is the test to know is point is left or right of a line is precise, but the test to know if a point is on a line doesn't work well.e found to intersect the line.
(an example of such a line : "LINESTRING(1447.451535 534.348973,1447.89016 543.902917)" )
In previous version, some line had no point at all (ie : no points were on the line, except the points on extremities.).
Most of the time ( p>(1-1/10^6) ), computing a point on a line using ST_ClosestPoint or ST_LineInterpolatePoint would give a point which was not on the line according to the test.

So instead, we create 2 points offsetted by an epsilon (corresponding to the 13 digits of the coordinates), and then check if the left offsetted point is left of the line, and the right offseted point is right of the line, or if any is on the line.

NOTE : now we have a probability of error < 1/10^6 for coordinates up to digits (i didn't tested more). This means there is no warranty that a line will have points at all , even if odds are very low that a line would have no points. Snatching vertex to grid would certainly solve the problem if such a line was encountered.

The postgis functions that have been improved are :
ST_Union ;ST_UnaryUnion ;ST_Intersection; ST_Relate, ST_Substring, ST_BuildArea ...

![image](https://f.cloud.github.com/assets/4725794/1601775/577f5062-538e-11e3-9c68-a1113a91ed7f.png)

The plpgsql code to test the changes regarding the new behavior of st_intersects is :

---building a test table with random lines
SELECT setseed(0.192169354297221); --settign seed so every execution gives same result
DROP TABLE IF EXISTS temp_test_precision_geos_serie_max;
CREATE TABLE temp_test_precision_geos_serie_max AS
WITH point_number AS (
SELECT power(10,prec.p-1)::bigint As pprec, 1000000 AS pnum, trunc(random() ::numeric,prec.p) AS rand, prec.p AS prec
FROM (SELECT 8 aS p) AS prec
),
serie AS (
SELECT s1 AS s
FROM point_number,generate_series(1,pnum) AS s1
),
points AS(
SELECT ST_MakePoint(
trunc((s*random()/pnum )::numeric,prec)*pnum
,trunc((s*random()/pnum )::numeric,prec)*pnum )AS p1
,ST_MakePoint(
trunc((s*random()/pnum )::numeric,prec)*pnum
,trunc((s*random()/pnum )::numeric,prec)*pnum ) AS p2
FROM point_number,serie
),
line AS (
SELECT *, St_MakeLine(p1,p2) AS line
FROM points
)
SELECT row_number() over() AS uid
,*
, ST_LineInterpolatePoint(
line, rand) AS ipoint
FROM point_number,line
ORDER BY p1 ASC , p2 ASC;

--testing for each line
SELECT setseed(0.192169354297221);
DROP TABLE IF EXISTS temp_test_precision_postgis;
CREATE TABLE temp_test_precision_postgis AS
with test AS (
SELECT uid, p1, p2, line, ipoint
FROM temp_test_precision_geos_serie_max
--LIMIT 1
),
closest AS (
SELECT * , ST_ClosestPoint(line,ST_Translate(ipoint, ST_X(ipoint)*round(random()::numeric,10),ST_Y(ipoint)*round(random()::numeric,10))) AS cpoint
FROM test
),
verif AS (
SELECT St_Intersects(line,ipoint) AS intersect_line_interpolated_point
, St_Intersects(line,cpoint) AS intersect_line_closest_point
, ST_Relate(line, cpoint) AS relate_line_closest_point
, line, ipoint
FROM closest
)
SELECT row_number() OVEr () AS id_qgis,*
FROM verif
WHERE intersect_line_interpolated_point=false OR intersect_line_closest_point=FALSE ;