Thread: Line intersection point is wrong
While working on removing FP macros of the geometric type operators,
I noticed that line intersection is giving wrong results:
regression=# select '{1,1,1}'::line # '{2,1,0}'::line;
?column?
----------
(1,2)
(1 row)
There are a few other operators relying on this function which are also
returning wrong results. Segfaults should be possible at least on
assert enabled builds, though FP macros makes it hard to get them.
The bug seems to be introduced by 215bc83d on 1997. I couldn't find
any relevant bug report on the archives except Tom Lane's doubtful
comment [1] from 1999. Patch attached to fix the formula.
[1] https://www.postgresql.org/message-id/5537.941498236%40sss.pgh.pa.us
Attachment
Emre Hasegeli <emre@hasegeli.com> writes:
> While working on removing FP macros of the geometric type operators,
> I noticed that line intersection is giving wrong results:
Hm. I don't think I believe the vertical-line cases there either.
They seem to be assuming A = -1 in a vertical line, which would be
true if the line was computed by line_construct_pts, but otherwise
not necessarily.
regression=# select '((1,0),(1,1))'::line;
line
----------
{-1,0,1}
(1 row)
regression=# select '{-1,0,1}'::line # '((0,2),(2,2))'::line;
?column?
----------
(1,2) -- this answer is correct
(1 row)
regression=# select '{-2,0,2}'::line # '((0,2),(2,2))'::line;
?column?
----------
(2,2)
(1 row)
If I haven't lost it completely, {-1,0,1} and {-2,0,2} are
equivalent line values. line_eq certainly thinks so.
Also: your formulation of the general case assumes that
(l1->A * l2->B - l2->A * l1->B) is not zero, which I'm
not entirely convinced of. In principle the line_parallel test
would catch the case, but seeing that that is not exactly how
line_parallel computes its result, roundoff error could bite us
here. I wonder if line_interpt_internal should skip the
line_parallel call and instead do its own tests for zero divide
to detect parallel lines.
regards, tom lane
> Hm. I don't think I believe the vertical-line cases there either.
> They seem to be assuming A = -1 in a vertical line, which would be
> true if the line was computed by line_construct_pts, but otherwise
> not necessarily.
I think we can just remove those cases.
> Also: your formulation of the general case assumes that
> (l1->A * l2->B - l2->A * l1->B) is not zero, which I'm
> not entirely convinced of. In principle the line_parallel test
> would catch the case, but seeing that that is not exactly how
> line_parallel computes its result, roundoff error could bite us
> here. I wonder if line_interpt_internal should skip the
> line_parallel call and instead do its own tests for zero divide
> to detect parallel lines.
If it would do its own test, the result would be inconsistent with ?#
and ?||. It is not hard to get an overflow with the current logic
anyway:
hasegeli=# select '{1,1,-2}'::line # '{1,2,-3}'::line;
?column?
---------------------
(Infinity,Infinity)
(1 row)
Evidently nobody is using these operators. I am going to send a patch
reworking many of those to make them less vulnerable to roundoff
errors for the next release.
Emre Hasegeli <emre@hasegeli.com> writes:
>> Hm. I don't think I believe the vertical-line cases there either.
>> They seem to be assuming A = -1 in a vertical line, which would be
>> true if the line was computed by line_construct_pts, but otherwise
>> not necessarily.
> I think we can just remove those cases.
That would be nice.
>> Also: your formulation of the general case assumes that
>> (l1->A * l2->B - l2->A * l1->B) is not zero, which I'm
>> not entirely convinced of. In principle the line_parallel test
>> would catch the case, but seeing that that is not exactly how
>> line_parallel computes its result, roundoff error could bite us
>> here. I wonder if line_interpt_internal should skip the
>> line_parallel call and instead do its own tests for zero divide
>> to detect parallel lines.
> If it would do its own test, the result would be inconsistent with ?#
> and ?||.
Hmm, well, what if we change line_parallel() so that what it tests
is whether (l1->A * l2->B - l2->A * l1->B) is zero? That seems
simpler and more symmetric than what it does now.
regards, tom lane
... also, a quick review of all the LINE-related code in geo_ops.c
says that it's just *full* of bugs. Aside from the issues already
mentioned:
line_recv fails to reject A=B=0
line_perp can get division by zero
silly coding in line_eq (A and B can't both be zero, no need to look at C)
line_distance ignores l1->A, suspicious
close_pl busted in same way as line_interpt_internal
close_ps ignores tmp->B, suspicious
I suspect that a lot of this code was originally written on the assumption
that vertical and horizontal lines would have B or A exactly -1
respectively. Somewhere along the line that restriction got dropped,
but the code wasn't improved to cope.
regards, tom lane
> Hmm, well, what if we change line_parallel() so that what it tests > is whether (l1->A * l2->B - l2->A * l1->B) is zero? That seems > simpler and more symmetric than what it does now. I was thinking to do so. I would also check for both sides not to overflow, but that can wait. Patch attached. Unfortunately, the change causes too much noise on the regression tests.
Attachment
Emre Hasegeli <emre@hasegeli.com> writes:
>> Hmm, well, what if we change line_parallel() so that what it tests
>> is whether (l1->A * l2->B - l2->A * l1->B) is zero? That seems
>> simpler and more symmetric than what it does now.
> I was thinking to do so. I would also check for both sides not to
> overflow, but that can wait. Patch attached. Unfortunately,
> the change causes too much noise on the regression tests.
I don't think this is right:
+ x = (l1->B * l2->C - l2->B * l1->C) / (l1->A * l2->B - l2->A * l1->B);
+ y = (l1->A * x + l1->C) / l1->B;
You haven't done anything to exclude the possibility that l1->B is zero,
so you could be getting zero-divide in the y calculation.
FWIW, given that you're claiming the calculation is wrong as-is, there
is no particularly good reason to assume that the expected results
embedded in select_views.out are correct :-(. I doubt anyone's ever
tried to verify those independently.
regards, tom lane
> You haven't done anything to exclude the possibility that l1->B is zero, > so you could be getting zero-divide in the y calculation. Ah, yes. We can use l2 when l1->B is zero. I will send another patch on tomorrow addressing other bugs you found.
Emre Hasegeli <emre@hasegeli.com> writes:
>> You haven't done anything to exclude the possibility that l1->B is zero,
>> so you could be getting zero-divide in the y calculation.
> Ah, yes. We can use l2 when l1->B is zero.
After working out the algebra by hand, I think the attached is correct
(and it does pass the regression tests, yay). I also made the
line_parallel and line_perp tests more symmetric and zero-divide-free.
regards, tom lane
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 657bcee..ce796c7 100644
*** a/src/backend/utils/adt/geo_ops.c
--- b/src/backend/utils/adt/geo_ops.c
*************** line_parallel(PG_FUNCTION_ARGS)
*** 1108,1117 ****
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
! if (FPzero(l1->B))
PG_RETURN_BOOL(FPzero(l2->B));
! PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B)));
}
Datum
--- 1108,1119 ----
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
! if (FPzero(l1->A)) /* horizontal? */
! PG_RETURN_BOOL(FPzero(l2->A));
! if (FPzero(l1->B)) /* vertical? */
PG_RETURN_BOOL(FPzero(l2->B));
! PG_RETURN_BOOL(FPeq(l2->A * l1->B, l1->A * l2->B));
}
Datum
*************** line_perp(PG_FUNCTION_ARGS)
*** 1120,1131 ****
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
! if (FPzero(l1->A))
PG_RETURN_BOOL(FPzero(l2->B));
! else if (FPzero(l1->B))
PG_RETURN_BOOL(FPzero(l2->A));
! PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
}
Datum
--- 1122,1133 ----
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
! if (FPzero(l1->A)) /* horizontal? */
PG_RETURN_BOOL(FPzero(l2->B));
! if (FPzero(l1->B)) /* vertical? */
PG_RETURN_BOOL(FPzero(l2->A));
! PG_RETURN_BOOL(FPeq(l2->A * l1->B, -(l1->A * l2->B)));
}
Datum
*************** line_interpt_internal(LINE *l1, LINE *l2
*** 1233,1250 ****
if (FPzero(l1->B)) /* l1 vertical? */
{
! x = l1->C;
! y = (l2->A * x + l2->C);
}
else if (FPzero(l2->B)) /* l2 vertical? */
{
! x = l2->C;
! y = (l1->A * x + l1->C);
}
else
{
! x = (l1->C - l2->C) / (l2->A - l1->A);
! y = (l1->A * x + l1->C);
}
result = point_construct(x, y);
--- 1235,1252 ----
if (FPzero(l1->B)) /* l1 vertical? */
{
! x = -(l1->C / l1->A);
! y = -(l2->A * x + l2->C) / l2->B;
}
else if (FPzero(l2->B)) /* l2 vertical? */
{
! x = -(l2->C / l2->A);
! y = -(l1->A * x + l1->C) / l1->B;
}
else
{
! x = (l1->B * l2->C - l2->B * l1->C) / (l1->A * l2->B - l2->A * l1->B);
! y = -(l1->A * x + l1->C) / l1->B;
}
result = point_construct(x, y);
> ... also, a quick review of all the LINE-related code in geo_ops.c > says that it's just *full* of bugs. Aside from the issues already > mentioned: > > line_recv fails to reject A=B=0 > line_perp can get division by zero > silly coding in line_eq (A and B can't both be zero, no need to look at C) > line_distance ignores l1->A, suspicious > close_pl busted in same way as line_interpt_internal I have tried to fix them. Patch attached. > close_ps ignores tmp->B, suspicious I think we better throw this one away, and implement another one using close_pl().
Attachment
> After working out the algebra by hand, I think the attached is correct
> (and it does pass the regression tests, yay). I also made the
> line_parallel and line_perp tests more symmetric and zero-divide-free.
It looks better to me, though not completely symmetric:
hasegeli=# select '{0.000001,1000,0}'::line ?|| '{0.00009,90000,0}'::line;
?column?
----------
f
(1 row)
hasegeli=# select '{0.00009,90000,0}'::line ?|| '{0.000001,1000,0}'::line;
?column?
----------
t
(1 row)
Emre Hasegeli <emre@hasegeli.com> writes:
>> After working out the algebra by hand, I think the attached is correct
>> (and it does pass the regression tests, yay). I also made the
>> line_parallel and line_perp tests more symmetric and zero-divide-free.
> It looks better to me, though not completely symmetric:
> hasegeli=# select '{0.000001,1000,0}'::line ?|| '{0.00009,90000,0}'::line;
> ?column?
> ----------
> f
> (1 row)
> hasegeli=# select '{0.00009,90000,0}'::line ?|| '{0.000001,1000,0}'::line;
> ?column?
> ----------
> t
> (1 row)
Hmm, that's annoying, although of course the existing code has problems
of the same ilk:
regression=# select '{1000,0.000001,0}'::line ?|| '{90000,0.00009,0}'::line;
?column?
----------
f
(1 row)
regression=# select '{90000,0.00009,0}'::line ?|| '{1000,0.000001,0}'::line;
?column?
----------
t
(1 row)
Basically this comes down to the fuzziness of the FPzero/FPeq tests.
It might be possible to fix this particular problem by changing the FPzero
tests in line_parallel() to exact "== 0" tests; though we would likely
be well advised to still keep the main test as FPeq.
regards, tom lane
> Basically this comes down to the fuzziness of the FPzero/FPeq tests. > > It might be possible to fix this particular problem by changing the FPzero > tests in line_parallel() to exact "== 0" tests; though we would likely > be well advised to still keep the main test as FPeq. That can create whole bunch of different problems, because then the fuzziness wouldn't apply to the values around the zero. I don't think there is any good solution without removing those macros all together. We should just fix what is obviously broken for now.