line_perp() (?-|) is broken.
I happend to see a strange geometric calcualtion on master/HEAD.
CREATE TABLE t (l1 line, l2 line);
INSERT INTO t
(SELECT line(point(0, 0), point(x, y)), line(point(0,0), point(-y, x))
FROM (SELECT random() x, random() y FROM generate_series(0, 1000)) AS a);
SELECT l1?-|l2 AS is_perp, l1, l2 FROM t;
is_perp | l1 | l2
---------+----------------------------+-----------------------------
f | {0.215980373614968,-1,0} | {-4.63005032940037,-1,0}
f | {1.51653638154567,-1,0} | {-0.659397303070824,-1,0}
f | {0.596861744826871,-1,0} | {-1.67542987746696,-1,0}
...
SELECT l1?-|l2 AS is_perp, l1, l2 FROM t WHERE l1?-|l2;
is_perp | l1 | l2
---------+----+----
(0 rows)
The two lines in a row are always perpendicular, but ?-| doesn't
agree.
line_perp() is working with the following arithmetic.
(l1->A * l2->B) / (l1->B * l2->A) == -1.0
or
(l1->A * l2->B) + (l1->B * l2->A) == 0
If the plus were a minus, it would calculate the cross product of
the two vectors and tell if the two lines are "parallel" or
not... Anyway it doesn't work as expected. At least back to 9.0
has the same expression in the function.
Instead, calculating inner product of the two direction vectors
works as expected.
(l1->A * l2->A) + (l1->B * l2->B) == 0
FPzero checks at the beggining of the function for the purpose of
avoiding div0 is useless with this form of expression.
With the attached patch, we will have the correct result.
is_perp | l1 | l2
---------+----------------------------+-----------------------------
t | {0.215980373614968,-1,0} | {-4.63005032940037,-1,0}
t | {1.51653638154567,-1,0} | {-0.659397303070824,-1,0}
t | {0.596861744826871,-1,0} | {-1.67542987746696,-1,0}
t | {2.14739225724798,-1,0} | {-0.465681105361517,-1,0}
...
SELECT l1?-|l2 AS is_perp, l1, l2 FROM t WHERE NOT l1?-|l2;
is_perp | l1 | l2
---------+----+----
(0 rows)
regards.
Attachments:
fix_line_perp.patchtext/x-patch; charset=us-asciiDownload
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index f57380a..3b12d41 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -1119,12 +1119,7 @@ line_perp(PG_FUNCTION_ARGS)
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));
+ PG_RETURN_BOOL(FPzero(l1->A * l2->A + l1->B * l2->B));
}
Datum
Kyotaro HORIGUCHI <horiguchi.kyotaro@lab.ntt.co.jp> writes:
I happend to see a strange geometric calcualtion on master/HEAD.
...
Instead, calculating inner product of the two direction vectors
works as expected.
(l1->A * l2->A) + (l1->B * l2->B) == 0
This seems to be a strict subset of the changes in Emre Hasgeli's
latest geometric-types patch. I'm disinclined to commit it unless
that patch gets rejected, as it'll just require Emre to rebase again.
However, for either patch, I'm a bit concerned about using FPzero()
on the inner product result. To the extent that the EPSILON bound
has any useful meaning at all, it needs to mean a maximum difference
between two coordinate values. The inner product has units of coordinates
squared, so it seems like EPSILON isn't an appropriate threshold there.
Possibly this objection is pointless, because I'm not at all sure that
the existing code is careful about how it uses FPeq/FPzero; perhaps
we're applying EPSILON to all manner of numbers that don't have units
of the coordinate space. But this won't help make it better.
regards, tom lane
However, for either patch, I'm a bit concerned about using FPzero()
on the inner product result. To the extent that the EPSILON bound
has any useful meaning at all, it needs to mean a maximum difference
between two coordinate values. The inner product has units of coordinates
squared, so it seems like EPSILON isn't an appropriate threshold there.
In this case, I suggest this:
if (FPzero(l1->A))
PG_RETURN_BOOL(FPzero(l2->B));
if (FPzero(l2->A))
PG_RETURN_BOOL(FPzero(l1->B));
if (FPzero(l1->B))
PG_RETURN_BOOL(FPzero(l2->A));
if (FPzero(l2->B))
PG_RETURN_BOOL(FPzero(l1->A));PG_RETURN_BOOL(FPeq((l1->A * l2->A) / (l1->B * l2->B), -1.0));
Possibly this objection is pointless, because I'm not at all sure that
the existing code is careful about how it uses FPeq/FPzero; perhaps
we're applying EPSILON to all manner of numbers that don't have units
of the coordinate space. But this won't help make it better.
The existing code was probably paying attention to this particular
one, but it fails to apply EPSILON meaningfully to many other places.
For example lseg_parallel() and lseg_perp() applies it 2 times to the
same input. First point_sl() compares x coordinates with FPeq(), and
then the returned slopes are compared again with FPeq().
Emre Hasegeli <emre@hasegeli.com> writes:
Possibly this objection is pointless, because I'm not at all sure that
the existing code is careful about how it uses FPeq/FPzero; perhaps
we're applying EPSILON to all manner of numbers that don't have units
of the coordinate space. But this won't help make it better.
The existing code was probably paying attention to this particular
one, but it fails to apply EPSILON meaningfully to many other places.
For example lseg_parallel() and lseg_perp() applies it 2 times to the
same input. First point_sl() compares x coordinates with FPeq(), and
then the returned slopes are compared again with FPeq().
Yeah, comparing a slope to EPSILON sure feels wrong :-(
regards, tom lane
Hello, I'm returning to here sonn.
I regard Emre's work is aiming to refactor (not improve its
functionality of) the code, on the other hand thins one is a
separate bug fix which also should be backpatched.
But, honestrly I'm not sure such small fix would improve the
current situnation of the geometric operators at any rate.
At Sat, 03 Mar 2018 10:43:26 -0500, Tom Lane <tgl@sss.pgh.pa.us> wrote in <31399.1520091806@sss.pgh.pa.us>
Emre Hasegeli <emre@hasegeli.com> writes:
Possibly this objection is pointless, because I'm not at all sure that
the existing code is careful about how it uses FPeq/FPzero; perhaps
we're applying EPSILON to all manner of numbers that don't have units
of the coordinate space. But this won't help make it better.
Agreed.
The existing code was probably paying attention to this particular
one, but it fails to apply EPSILON meaningfully to many other places.
For example lseg_parallel() and lseg_perp() applies it 2 times to the
same input. First point_sl() compares x coordinates with FPeq(), and
then the returned slopes are compared again with FPeq().Yeah, comparing a slope to EPSILON sure feels wrong :-(
It is a quite significant problem to fix and would be
controversial in detail. But, anyway, we are focusing on other
points that are less controversial. Then cook the EPSILON in the
next stage.
regards,
--
Kyotaro Horiguchi
NTT Open Source Software Center