Re: Query on distance between 2D line segments

<hoffmann@xxxxxxxxxxxx> wrote in message
news:1121192252.735978.91450@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
> Show your data set where my method fails.
<and>
> Maybe I'm wrong here and there, typos for example, but
> if my concepts are wrong then you should prove this,
> or don't attack me.

Impatience and paranoia. How sad. I post when I have the
time away from my regular job. Sometimes that time is not
immediately available. Also, you made a claim, I disagreed
with it, you said RTFM, and I said post details to support
your claim. Attack? I do not think so. If you do not like it
when folks question what you post, consider defending your
post in a manner different than RTFM. Or consider not
posting at all here.

Well, I RTFM'd. Posting your implementations in Pascal, and as
a PDF, make it tedious to convert to something you can compile
in a modern language. Regardless, I implemented a C++ version
which duplicates what you have. Some testing convinced me that
I ported the code correctly.

To show one problem, consider an example in 2D, so in my discussion
think of the 3D case with z = 0. Consider two line segments of equal
length L and whose midpoints are at the origin. Let the unit-length
direction of the first segment be U1 and the unit-length direction of the
second segment be U2. Let A be a small angle and set U2 to the vector
U1 rotated counterclockwise by A radians. Define
P1 = -0.5*L*U1, Q1 = 0.5*L*U1, P2 = -0.5*L*U2, Q2 = 0.5*L*U2
I chose L = 10000, A = 1e-08, and generated U1.x and U1.y in [-1,1]
using a random number generator and set U1.z to zero, and then normalized
U1 to be unit length. U2 was computed from it by rotation.

Your distance calculator produces D0 = 4 and D1 = D2 = 0.5 and the
values SA = SB = 1/8. The distance is 3.75e-05, which appears to be
close to the true value of zero. However, your closest points are
X1 = (3719.891010, -474.247692, 0)
X2 = (3719.891015, -474.247655, 0)
Each of these points is *far away from the true value (0,0,0)*.

The problem is that D0 = 4 occurs because of catastrophic cancellation
in the expression a11*a22 - a12*a12. It really should be D0 = 1 so
that SA = SB = 1/2. [This type of problem is why in my own code base
I switched to a segment representation much like that of oriented bounding
boxes: center point, unit-length direction, and extent (half-length). My
distance calculator returns distance exactly zero and closest point
the origin.]

As another poster indicated, it is also possible that the numerical value
for D0 is a small negative number. Theoretically, D0 >= 0 because
it is the squared length of the cross product of the direction vectors.
Your code has an explicit test for D0 < 0, but it does not appear to
recognize that this condition is incompatible with the theory. Regardless,
the other poster's point was that numerical round-off errors might cause
you to misclassify a situation. Two segments are deemed parallel when
they are not. Or two segments are deemed not parallel when they are.
An interesting way to see how robust your code: When the input
produces a value such as a determinant, and that value is close to the
"switch" you use to toggle between the "then" and "else" clauses, use a
debugger and set the program counter to point to the clause opposite
the one you would take under normal execution. If you get a drastically
different result, you have some more designing to do.

At any rate, I have seen examples of problems like this when using
'float'. Use of 'double' tends to hide the problems and makes
counterexamples more difficult to locate. Enough diligence will lead
to them.

> We had already a similar discussion about intersections
> of Bézier curves, represented by straight line segments.

Yes, we have. And some of us keep pointing out that mathematics
is not enough. Floating-point issues must be considered in any
implementation of such algorithms.

--
Dave Eberly
http://www.geometrictools.com


.

Relevant Pages

  • Re: Query on distance between 2D line segments
    ... > direction of the first segment be U1 and the unit-length direction of the ... > second segment be U2. ... > close to the true value of zero. ... You may always rely on the uncertainty in the input data set. ...
    (comp.graphics.algorithms)
  • Re: [crash, bisected] Re: [PATCH 3/4] x86_64: Fold pda into per cpu area
    ... to initialize the kernel segments. ... The zero-based PDA mechanism requires the introduction of a new ELF ... Why are we doing a zero based pda? ... offset in the segment register to allow for everything to work as ...
    (Linux-Kernel)
  • Re: Efficency and the standard library
    ... that she or he could evaluate the code, and either find bugs or verify ... The problem with the Nilges string replace program is largely that the design ... Check for a previous segment: ... See if it has a zero or negative length ...
    (comp.lang.c)
  • Re: Storage of Variables
    ... code at entry to a function zero out the stack. ... In C, there is really no such concept as a heap, though you could ... the bss segment. ... And to initialize the data segment, ...
    (comp.os.linux.misc)
  • Re: Is a line segment composed of points?
    ... piece if of zero length, then it is not a result of cutting at ... You end up with an infinite number of pieces, ... segment is a matter of how you define line segment. ... first segment is not reached by a division process. ...
    (sci.math)