By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
435,425 Members | 2,549 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 435,425 IT Pros & Developers. It's quick & easy.

Point on Line Segment in 2D. Which code is faster ? Can you improve it ?

P: n/a
Hi,

I needed a method to determine if a point was on a line segment in 2D. So I
googled for some help and so far I have evaluated two methods.

The first method was only a formula, the second method was a piece of C code
which turned out to be incorrect and incomplete but by modifieing it would
still be usuable.

The first method was this piece of text:

"
if x is between StartX and EndX or if y is between StartY and EndY, then the
point (x,y) is on the line.

Another way is the parametric form of a line,ie:
X = StartX + t0 * (EndX - StartX)
Y = StartY + t1 * (EndY - StartY)
calculate t0 and t1, if t0 and t1 are equal, then the point is on the line
if the parameter (t0 = t1) is between 0 and 1, the the point is
on the segment.
"

The second method was this piece of text:

"
Very simple,

First test if the x-value of the point is within the interval StartX,EndX.
If this is not true the point is not on the line.

int on_line(x, y, StartX, StartY, EndX, EndY)
{ int dx, dy, dx1, dy1;

if (x<StartX || x>EndX) return 0;
dx = EndX-StartX;
dy = EndY-StartY;
dx1 = x-StartX;
dy1 = y-StartY;

return (dx*dy1 == dy*dx1);
}
"

So after reading the first method I needed to be able to calculate T0 and T1
since all other variables are known the formula's can be rewritten as to
calculate T0 and T1.

Changing the formula's to get the final formula:

Step1 original:
X = StartX + t0 * (EndX - StartX)
Y = StartY + t1 * (EndY - StartY)

Step2 delta's introduced:
DeltaX = EndX - StartX
DeltaY = EndY - Starty
X = StartX + t0 * (DeltaX)
Y = StartY + t1 * (DeltaY)

Step3 switch
StartX + t0 * (DeltaX) = X
StartY + t1 * (DeltaY) = Y

Step4 start to isolate t0 and t1 by moving StartX and StartY to the
otherside
(thanks to my math teacher for those usefull lessons =D)
t0 * (DeltaX) = X - StartX
t1 * (DeltaY) = Y - StartY

Step5 isolate to and t1 further by moving delta's to the other side
(thanks to my math teacher again for those usefull lessons =D)
t0 = (X - StartX)/DeltaX
t1 = (Y - StartY)/DeltaY

If you look at the original formula you will notice how t0 and t1 are sort
of distances on the line segment itself ranging from 0.0 to 1.0 (like a
percentage).

The second method looks a little bit similair but works totally different.

It multiplies the delta's with the line segment with the delta's with the
start and point. Which kinda looks like a dot product.. I dont know what
that's all about ;) Maybe this is just another way of calculating how much
distance each component has. (There are two components x and y, to be on the
line they both need to have the same distance from the start point for
example ;) ofcourse this could mean a circle... but since you working with
delta's there is only one possible direction)

First of all the C method is incorrect unless the coordinates are pre-sorted
from left to right. If the coordinates are not sorted, for example
(100,10)-(10,50) (going from right to left) the method will simply exit
because of the first if statement which is just wrong. So this is easily
fixed by removing that if statement.

However I also consider the C method to be incomplete. It's fine if one
wants to calculate if a point is on an infinite line. It's impossible to
tell if the original poster wanted a line segment or just a line... he does
mention an start and end point... For my purposes however I need to be able
to tell if a point is on a line segment. So this C method needs to be
modified/expanded to accomadate for it.

The fascination starts with the C method looking very fast. Just a few lines
of code and a multiplication.
It also looks a bit dangerous because of the multiplication which could
overflow... but all in all not to bad/shabby and definetly worth a try.

The first method looks kinda big and slower because of more variables/code
and divisions etc, but it does look complete and stable and robust.

So now for a reality check.

I have implemented both versions in pascal/delphi which I believe to produce
correct results everytime.

The question is which method is faster. I would place my money on the first
method. Simply because the code is shorter and it has to do less floating
point comparisions and less jumps than the second method.

The C method (second method) actually needs a bounding box to be constructed
to be usefull (for line segment tests). The bounding box is necessary to
determine if the point is inside the line segment. Therefore I believe the C
method in this case to be slower. It requires more instructions, more
comparisions and more jumps.

I think I have done a fair job of optimizing both implementations by trying
to reduce the number of comparisions, jumps and optimizing for branch
prediction (most likely case fall through etc).

I spent lot's of time optimizing method1 to reduce the number of compares by
first checking one
delta before checking the next etc.. that simply saved one or two compares
and jumps compared
to previous versions.

I also spent some time optimizing method2 but probably a little bit less
time. I inlined the bounding box calculation and I also postponed it until
after the if statement etc as to not execute if the if statement is false.

Finally the point needs to be checked to be inside or outside the bounding
box. Checking for outside is easier/shorter for multiple bounding so I coded
it like that. However in this case there is only one bounding box so in this
case it doesn't matter and it could also be checked to be inside. The
generated assembler would however be the same so it doesn't matter.

For my purposes the diagonal case will be the most likely and occuring case
so it's most interesting to see which method would be faster.

The first method requires 53 instructions for the diagonal case (with point
on line)

The second method requires 81 instructions for the diagonal case (with point
on line)

Both methods require the same ammounts of pushes onto the stack so those
pushes have not been counted.

The current implementations look like this:

v3 is the first method
v5 is the second method

While writing this post I discovered a flaw in v3, fixed it, and even
optimized it.
I think v3 should now be much faster for case horizontal and case vertical
than v5.

So the main focus is on the diagonal case.

My question to you the delphi, electronics and computer architecture
community is:

Which code do you think is faster ?

Can you explain it purely based on theory without doing any performance
measurements ? ;)

For example: number of instructions, instruction cycles, branch prediction,
pipe line optimizations, pairing, etc.

I shall provide you with the current pascal/delphi and the assembler code:

If you think you really are an expert in the optimization
field/instruction/hardware/assembler field then next question will be
ofcourse:

Can you optimize these methods or come up with a faster method ?

Quite frankly I am obsessed with this code. It's such a trivial code that I
dont want my PC to be spending to much time on it... currently this code is
not being used in any algorithms. I could use alternative information from a
line segment intersection test algorithm which can already determine if one
of the 4 points is on a line segment yes/no etc... but that would make the
routine's header much more complex and maybe I dont want that at first at
least.... so I might choose to use a routine like below to keep things a bit
more simple and pragmatic... etc.

So after each line segment intersection test this routine would also need to
be called, so it would be called a lot, so I dont want it to be a lot of
overhead if you know what I mean ;)

And even if this code will never be used it's still quite fascinating and
educational to see/learn how floating point optimizations work etc... that's
the main reason why I am investigating this and what you to have a look at
this because I am puzzled =D

Ok now the subjective stuff:

My gutt feeling says v3 is faster than v5 because v3 uses less compares and
less jumps which should be good for branch prediction etc.. v3 also has less
instructions.

However v3 uses two divisions (fdiv) which might be a lot slower than
multiplications (fmul) which version v5 uses two of. However my gutt feeling
says v5 has the bounding box thing going on etc so it needs to do more
stuff/instructions which could make it slow ;)

Add the unknown factors like branch prediction, pipe lining, pairing and god
knows what else and you start to see the doubt hehehehehe.

Can you sort it out ? ;) =D

Or do you need performance measurements like me ?

And let's be honest about that... How are we going to do accurate
performance measurements ?

Lot's of context switching, multi threading, harddisk activity going on, on
windows xp, so I am not too sure if that is reliable but ok :)

Here is the pascal/delphi code:

function point_on_line_segment_in_2d_v3( StartX, StartY : double;
EndX, EndY : double;
PointX, PointY : double ) : boolean;
var
vDeltaX : double;
vDeltaY : double;

vDistanceXcomponent : double;
vDistanceYcomponent : double;
begin
// default is false
result := false;

vDeltaX := EndX - StartX;
vDeltaY := EndY - StartY;

// if vDeltaX is not zero then that means the line is horizontal or
diagonal
if (vDeltaX <> 0) then
begin
// if vDeltaY is not zero then that means the line is diagonal
if (vDeltaY <> 0) then
begin
// line is diagonal

// calculate both components and look if they are the same
// if yes then see if they both lie between 0 and 1.
vDistanceXcomponent := (PointX - StartX) / vDeltaX;
vDistanceYcomponent := (PointY - StartY) / vDeltaY;

if (vDistanceXcomponent = vDistanceYcomponent) then
begin
if (vDistanceXcomponent>0) and (vDistanceXcomponent<1) then
begin
result := true;
end;
end;

end else
// else it means it is horizontal
begin

// line is horizontal

// check if the point is at the same y position as the horizontal line
if (PointY = StartY) then
begin

if StartX < EndX then
begin
if (PointX > StartX) and (PointX < EndX) then
begin
result := true;
end;
end else
begin
if (PointX > EndX) and (PointX < StartX) then
begin
result := true;
end;
end;

// not necessary simply check if between startx and endx or
// vica versa depending on if endx is greater then startx etc
// since fdiv seems to be slow let's do it with some compares
// see above ;)
{
// calculate and look at the horizontal component if it lies between 0
and 1.
vDistanceXcomponent := (PointX - StartX) / vDeltaX;

if (vDistanceXComponent>0) and (vDistanceXComponent<1) then
begin
result := true;
end;
}
end;

end;

end else
// else that means the line is vertical or a point
// if vDeltaY is not zero then that means the line is vertical
if (vDeltaY <> 0) then
begin
// line is vertical

// check if point is at same x position as the vertical line.
if (PointX = StartX) then
begin

if StartY < EndY then
begin
if (PointY > StartY) and (PointY < EndY) then
begin
result := true;
end;
end else
begin
if (PointY > EndY) and (PointY < StartY) then
begin
result := true;
end;
end;

{
// calculate and see if the vertical component lies between 0 and 1
vDistanceYcomponent := (PointY - StartY) / vDeltaY;

if (vDistanceYComponent>0) and (vDistanceYComponent<1) then
begin
result := true;
end;
}
end;

// not necessary to do ith with divisions etc.

end;
// else the line is not a line but a point so exit
end;

function point_on_line_segment_in_2d_v5( StartX, StartY : double;
EndX, EndY : double;
PointX, PointY : double ) : boolean;
var
vDeltaX : double;
vDeltaY : double;

vDistanceX : double;
vDistanceY : double;

vBoundingBoxX1 : double;
vBoundingBoxY1 : double;
vBoundingBoxX2 : double;
vBoundingBoxY2 : double;
begin
result := false;

vDeltaX := EndX - StartX;
vDeltaY := EndY - StartY;

vDistanceX := PointX - StartX;
vDistanceY := PointY - StartY;

if (vDeltaX * vDistanceY) = (vDeltaY * vDistanceX) then
begin

if (StartX<EndX) then
begin
vBoundingBoxX1 := StartX;
vBoundingBoxX2 := EndX;
end else
begin
vBoundingBoxX1 := EndX;
vBoundingBoxX2 := StartX;
end;

if (StartY<EndY) then
begin
vBoundingBoxY1 := StartY;
vBoundingBoxY2 := EndY;
end else
begin
vBoundingBoxY1 := EndY;
vBoundingBoxY2 := StartY;
end;

// I'll check if the point is outside to determine the overlap etc.
// which is consistent which how it's done in the line segment
intersection/overlap
// test where two bounding boxes need to be checked for overlap.
// determining if they are outside is much easier and shorter then
checking if each
// point is inside ;)
// so let's not be stupid and stick with it ;)
// determining overlap should be done this way:
// checking for outside or on it if not then assume it's inside.
if (PointX <= vBoundingBoxX1) or (PointX >= vBoundingBoxX2) or
(PointY <= vBoundingBoxY1) or (PointY >= vBoundingBoxY2) then exit;

result := true;

// which is the same as:
// assembler generated is exactly the same if I am not mistaken
// this code looks cleaner tough... ;)
// if (PointX > vBoundingBoxX1) and (PointX < vBoundingBoxX2) and
// (PointY > vBoundingBoxY1) and (PointY < vBoundingBoxY2) then
// begin
// result := true;
// end;
end;
end;

And here is some simple test code:

var
a1x, a1y,
a2x, a2y,
px, py : double;
begin
a1x := 10;
a1y := 10;

a2x := 100;
a2y := 100;

px := 50;
py := 50;

if point_on_line_segment_in_2d_v3( a1x, a1y,
a2x, a2y,
px, py ) then
begin
ShowMessage('yes');
end else
begin
ShowMessage('no');
end;

if point_on_line_segment_in_2d_v5( a1x, a1y,
a2x, a2y,
px, py ) then
begin
ShowMessage('yes');
end else
begin
ShowMessage('no');
end;
Here is the assembler code. (since delphi's cpu window does not allow easy
copy/paste this code has been taking from w32dasm disassembler, which
shouldn't matter. Actually the disassembly is better since it shows exactly
where jumps jump too ;)

(Taken from the "release" version without debugging information... I think
the assembler is the same for the debugging version so it doesnt matter
(?)):

// assembler for v3:

// the call:

004523DE FF742404 push [esp+04]
:004523E2 FF742404 push [esp+04]
:004523E6 FF742414 push [esp+14]
:004523EA FF742414 push [esp+14]
:004523EE FF742424 push [esp+24]
:004523F2 FF742424 push [esp+24]
:004523F6 FF742434 push [esp+34]
:004523FA FF742434 push [esp+34]
:004523FE FF742444 push [esp+44]
:00452402 FF742444 push [esp+44]
:00452406 FF742454 push [esp+54]
:0045240A FF742454 push [esp+54]
:0045240E E839FDFFFF call 0045214C
:00452413 84C0 test al, al
:00452415 740C je 00452423

// the routine:

:0045214C 55 push ebp
:0045214D 8BEC mov ebp, esp
:0045214F 83C4E0 add esp, FFFFFFE0
:00452152 33D2 xor edx, edx
:00452154 DD4520 fld qword ptr [ebp+20]
:00452157 DC6530 fsub qword ptr [ebp+30]
:0045215A DD5DF8 fstp qword ptr [ebp-08]
:0045215D 9B wait
:0045215E DD4518 fld qword ptr [ebp+18]
:00452161 DC6528 fsub qword ptr [ebp+28]
:00452164 DD5DF0 fstp qword ptr [ebp-10]
:00452167 9B wait
:00452168 DD45F8 fld qword ptr [ebp-08]
:0045216B D81D88224500 fcomp dword ptr [00452288]
:00452171 DFE0 fstsw ax
:00452173 9E sahf
:00452174 0F84B0000000 je 0045222A
:0045217A DD45F0 fld qword ptr [ebp-10]
:0045217D D81D88224500 fcomp dword ptr [00452288]
:00452183 DFE0 fstsw ax

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:00452114(C)
|
:00452185 9E sahf
:00452186 7454 je 004521DC

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:00452112(C)
|
:00452188 DD4510 fld qword ptr [ebp+10]
:0045218B DC6530 fsub qword ptr [ebp+30]
:0045218E DC75F8 fdiv qword ptr [ebp-08]
:00452191 DD5DE8 fstp qword ptr [ebp-18]
:00452194 9B wait
:00452195 DD4508 fld qword ptr [ebp+08]
:00452198 DC6528 fsub qword ptr [ebp+28]
:0045219B DC75F0 fdiv qword ptr [ebp-10]
:0045219E DD5DE0 fstp qword ptr [ebp-20]
:004521A1 9B wait
:004521A2 DD45E8 fld qword ptr [ebp-18]
:004521A5 DC5DE0 fcomp qword ptr [ebp-20]
:004521A8 DFE0 fstsw ax
:004521AA 9E sahf
:004521AB 0F85CF000000 jne 00452280
:004521B1 DD45E8 fld qword ptr [ebp-18]
:004521B4 D81D88224500 fcomp dword ptr [00452288]
:004521BA DFE0 fstsw ax
:004521BC 9E sahf
:004521BD 0F86BD000000 jbe 00452280
:004521C3 DD45E8 fld qword ptr [ebp-18]
:004521C6 D81D8C224500 fcomp dword ptr [0045228C]
:004521CC DFE0 fstsw ax
:004521CE 9E sahf
:004521CF 0F83AB000000 jnb 00452280
:004521D5 B201 mov dl, 01
:004521D7 E9A4000000 jmp 00452280

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:00452186(C)
|
:004521DC DD4508 fld qword ptr [ebp+08]
:004521DF DC5D28 fcomp qword ptr [ebp+28]
:004521E2 DFE0 fstsw ax
:004521E4 9E sahf
:004521E5 0F8595000000 jne 00452280
:004521EB DD4530 fld qword ptr [ebp+30]
:004521EE DC5D20 fcomp qword ptr [ebp+20]
:004521F1 DFE0 fstsw ax
:004521F3 9E sahf
:004521F4 731A jnb 00452210
:004521F6 DD4510 fld qword ptr [ebp+10]
:004521F9 DC5D30 fcomp qword ptr [ebp+30]
:004521FC DFE0 fstsw ax
:004521FE 9E sahf
:004521FF 767F jbe 00452280
:00452201 DD4510 fld qword ptr [ebp+10]
:00452204 DC5D20 fcomp qword ptr [ebp+20]
:00452207 DFE0 fstsw ax
:00452209 9E sahf
:0045220A 7374 jnb 00452280
:0045220C B201 mov dl, 01
:0045220E EB70 jmp 00452280

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:004521F4(C)
|
:00452210 DD4510 fld qword ptr [ebp+10]
:00452213 DC5D20 fcomp qword ptr [ebp+20]
:00452216 DFE0 fstsw ax
:00452218 9E sahf
:00452219 7665 jbe 00452280
:0045221B DD4510 fld qword ptr [ebp+10]
:0045221E DC5D30 fcomp qword ptr [ebp+30]
:00452221 DFE0 fstsw ax
:00452223 9E sahf
:00452224 735A jnb 00452280
:00452226 B201 mov dl, 01
:00452228 EB56 jmp 00452280

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:00452174(C)
|
:0045222A DD45F0 fld qword ptr [ebp-10]
:0045222D D81D88224500 fcomp dword ptr [00452288]
:00452233 DFE0 fstsw ax
:00452235 9E sahf
:00452236 7448 je 00452280
:00452238 DD4510 fld qword ptr [ebp+10]
:0045223B DC5D30 fcomp qword ptr [ebp+30]
:0045223E DFE0 fstsw ax
:00452240 9E sahf
:00452241 753D jne 00452280
:00452243 DD4528 fld qword ptr [ebp+28]
:00452246 DC5D18 fcomp qword ptr [ebp+18]
:00452249 DFE0 fstsw ax
:0045224B 9E sahf
:0045224C 731A jnb 00452268
:0045224E DD4508 fld qword ptr [ebp+08]
:00452251 DC5D28 fcomp qword ptr [ebp+28]
:00452254 DFE0 fstsw ax
:00452256 9E sahf
:00452257 7627 jbe 00452280
:00452259 DD4508 fld qword ptr [ebp+08]
:0045225C DC5D18 fcomp qword ptr [ebp+18]
:0045225F DFE0 fstsw ax
:00452261 9E sahf
:00452262 731C jnb 00452280
:00452264 B201 mov dl, 01
:00452266 EB18 jmp 00452280

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:0045224C(C)
|
:00452268 DD4508 fld qword ptr [ebp+08]
:0045226B DC5D18 fcomp qword ptr [ebp+18]
:0045226E DFE0 fstsw ax
:00452270 9E sahf
:00452271 760D jbe 00452280
:00452273 DD4508 fld qword ptr [ebp+08]
:00452276 DC5D28 fcomp qword ptr [ebp+28]
:00452279 DFE0 fstsw ax
:0045227B 9E sahf
:0045227C 7302 jnb 00452280
:0045227E B201 mov dl, 01

* Referenced by a (U)nconditional or (C)onditional Jump at Addresses:
|:004521AB(C), :004521BD(C), :004521CF(C), :004521D7(U), :004521E5(C)
|:004521FF(C), :0045220A(C), :0045220E(U), :00452219(C), :00452224(C)
|:00452228(U), :00452236(C), :00452241(C), :00452257(C), :00452262(C)
|:00452266(U), :00452271(C), :0045227C(C)
|
:00452280 8BC2 mov eax, edx
:00452282 8BE5 mov esp, ebp
:00452284 5D pop ebp
:00452285 C23000 ret 0030

// assembler for v5:

// the call:

:0045242D FF742404 push [esp+04]
:00452431 FF742404 push [esp+04]
:00452435 FF742414 push [esp+14]
:00452439 FF742414 push [esp+14]
:0045243D FF742424 push [esp+24]
:00452441 FF742424 push [esp+24]
:00452445 FF742434 push [esp+34]
:00452449 FF742434 push [esp+34]
:0045244D FF742444 push [esp+44]
:00452451 FF742444 push [esp+44]
:00452455 FF742454 push [esp+54]
:00452459 FF742454 push [esp+54]
:0045245D E82EFEFFFF call 00452290
:00452462 84C0 test al, al
:00452464 740C je 00452472
// the routine:

* Referenced by a CALL at Address:
|:0045245D
|
:00452290 55 push ebp
:00452291 8BEC mov ebp, esp
:00452293 83C4C0 add esp, FFFFFFC0
:00452296 33D2 xor edx, edx
:00452298 DD4520 fld qword ptr [ebp+20]
:0045229B DC6530 fsub qword ptr [ebp+30]
:0045229E DD5DF8 fstp qword ptr [ebp-08]
:004522A1 9B wait
:004522A2 DD4518 fld qword ptr [ebp+18]
:004522A5 DC6528 fsub qword ptr [ebp+28]
:004522A8 DD5DF0 fstp qword ptr [ebp-10]
:004522AB 9B wait
:004522AC DD4510 fld qword ptr [ebp+10]
:004522AF DC6530 fsub qword ptr [ebp+30]
:004522B2 DD5DE8 fstp qword ptr [ebp-18]
:004522B5 9B wait
:004522B6 DD4508 fld qword ptr [ebp+08]
:004522B9 DC6528 fsub qword ptr [ebp+28]
:004522BC DD5DE0 fstp qword ptr [ebp-20]
:004522BF 9B wait
:004522C0 DD45F8 fld qword ptr [ebp-08]
:004522C3 DC4DE0 fmul qword ptr [ebp-20]
:004522C6 DD45F0 fld qword ptr [ebp-10]
:004522C9 DC4DE8 fmul qword ptr [ebp-18]
:004522CC DED9 fcompp
:004522CE DFE0 fstsw ax
:004522D0 9E sahf
:004522D1 0F85A8000000 jne 0045237F
:004522D7 DD4530 fld qword ptr [ebp+30]
:004522DA DC5D20 fcomp qword ptr [ebp+20]
:004522DD DFE0 fstsw ax
:004522DF 9E sahf
:004522E0 731A jnb 004522FC
:004522E2 8B4530 mov eax, dword ptr [ebp+30]
:004522E5 8945D8 mov dword ptr [ebp-28], eax
:004522E8 8B4534 mov eax, dword ptr [ebp+34]
:004522EB 8945DC mov dword ptr [ebp-24], eax
:004522EE 8B4520 mov eax, dword ptr [ebp+20]
:004522F1 8945C8 mov dword ptr [ebp-38], eax
:004522F4 8B4524 mov eax, dword ptr [ebp+24]
:004522F7 8945CC mov dword ptr [ebp-34], eax
:004522FA EB18 jmp 00452314

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:004522E0(C)
|
:004522FC 8B4520 mov eax, dword ptr [ebp+20]
:004522FF 8945D8 mov dword ptr [ebp-28], eax
:00452302 8B4524 mov eax, dword ptr [ebp+24]
:00452305 8945DC mov dword ptr [ebp-24], eax
:00452308 8B4530 mov eax, dword ptr [ebp+30]
:0045230B 8945C8 mov dword ptr [ebp-38], eax
:0045230E 8B4534 mov eax, dword ptr [ebp+34]
:00452311 8945CC mov dword ptr [ebp-34], eax

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:004522FA(U)
|
:00452314 DD4528 fld qword ptr [ebp+28]
:00452317 DC5D18 fcomp qword ptr [ebp+18]
:0045231A DFE0 fstsw ax
:0045231C 9E sahf
:0045231D 731A jnb 00452339
:0045231F 8B4528 mov eax, dword ptr [ebp+28]
:00452322 8945D0 mov dword ptr [ebp-30], eax
:00452325 8B452C mov eax, dword ptr [ebp+2C]
:00452328 8945D4 mov dword ptr [ebp-2C], eax
:0045232B 8B4518 mov eax, dword ptr [ebp+18]
:0045232E 8945C0 mov dword ptr [ebp-40], eax
:00452331 8B451C mov eax, dword ptr [ebp+1C]
:00452334 8945C4 mov dword ptr [ebp-3C], eax
:00452337 EB18 jmp 00452351

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:0045231D(C)
|
:00452339 8B4518 mov eax, dword ptr [ebp+18]
:0045233C 8945D0 mov dword ptr [ebp-30], eax
:0045233F 8B451C mov eax, dword ptr [ebp+1C]
:00452342 8945D4 mov dword ptr [ebp-2C], eax
:00452345 8B4528 mov eax, dword ptr [ebp+28]
:00452348 8945C0 mov dword ptr [ebp-40], eax
:0045234B 8B452C mov eax, dword ptr [ebp+2C]
:0045234E 8945C4 mov dword ptr [ebp-3C], eax

* Referenced by a (U)nconditional or (C)onditional Jump at Address:
|:00452337(U)
|
:00452351 DD4510 fld qword ptr [ebp+10]
:00452354 DC5DD8 fcomp qword ptr [ebp-28]
:00452357 DFE0 fstsw ax
:00452359 9E sahf
:0045235A 7623 jbe 0045237F
:0045235C DD4510 fld qword ptr [ebp+10]
:0045235F DC5DC8 fcomp qword ptr [ebp-38]
:00452362 DFE0 fstsw ax
:00452364 9E sahf
:00452365 7318 jnb 0045237F
:00452367 DD4508 fld qword ptr [ebp+08]
:0045236A DC5DD0 fcomp qword ptr [ebp-30]
:0045236D DFE0 fstsw ax
:0045236F 9E sahf
:00452370 760D jbe 0045237F
:00452372 DD4508 fld qword ptr [ebp+08]
:00452375 DC5DC0 fcomp qword ptr [ebp-40]
:00452378 DFE0 fstsw ax
:0045237A 9E sahf
:0045237B 7302 jnb 0045237F
:0045237D B201 mov dl, 01

* Referenced by a (U)nconditional or (C)onditional Jump at Addresses:
|:004522D1(C), :0045235A(C), :00452365(C), :00452370(C), :0045237B(C)
|
:0045237F 8BC2 mov eax, edx
:00452381 8BE5 mov esp, ebp
:00452383 5D pop ebp
:00452384 C23000 ret 0030

The assembler is provided in case you dont have any pascal/delphi compilers
at hand ;) and allows one to easily see what is going on.

I think I have provided enough assembler listing by providing the call and
the routine. There was some extra stuff at the beginning and the end like
dup(0) stuff and nop etc... but I dont know why it's there etc.. it's
probably not needed to understand what is going on with these routines, how
they work and how fast they will work ;)

I am curious if anybody is able to tell anything from these listings ;)

If not then that's ok I can still fall back on performance
testing/measurements =D though looking at it from a theoretical point of
view is very interesting ;)

P.S.: Ok I am obsessed with this code lol I admit =D

Bye,
Skybuck.
Nov 15 '05 #1
Share this Question
Share on Google+
65 Replies


P: n/a
Skybuck Flying wrote:
I needed a method to determine if a point was on a line segment in 2D. So I
googled for some help and so far I have evaluated two methods.


function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
begin
result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001;
end;

Cheers,
Nicholas Sherlock
Nov 15 '05 #2

P: n/a

"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote:
I needed a method to determine if a point was on a line segment in 2D. So I googled for some help and so far I have evaluated two methods.


function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
begin
result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001;
end;


Congratulations, you just reposted the C/second method ;)

Though the < 0.0000001 is new which I dont quite understand ;)

Bye,
Skybuck.
Nov 15 '05 #3

P: n/a

"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote:
I needed a method to determine if a point was on a line segment in 2D. So I googled for some help and so far I have evaluated two methods.


function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
begin
result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001;
end;


Besides it says OnLine ?

I need OnLineSegment...

See the difference ?

Bye,
Skybuck.
Nov 15 '05 #4

P: n/a
Nicholas Sherlock wrote:
Skybuck Flying wrote:
I needed a method to determine if a point was on a line segment in 2D.
So I
googled for some help and so far I have evaluated two methods.


function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
begin
result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001;
end;


Sorry, I missed some of this code from my app:

function PointOnLine2D(x1, y1, x2, y2, x3, y3: double): boolean;
begin
result :=
(((x3 >= x1) and (x3 <= x2)) or ((x3 >= x2) and (x3 <= x1))) and
(((y3 >= y1) and (y3 <= y2)) or ((y3 >= y2) and (y3 <= y1))) and
(abs(((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1))) < 0.00001);
end;

I had to rip it out of some other bit of code, so check it carefully
against your data.

Cheers,
Nicholas Sherlock
Nov 15 '05 #5

P: n/a
Skybuck Flying wrote:
Besides it says OnLine ?

I need OnLineSegment...

See the difference ?


The other code I posted does segment.

Cheers,
Nicholas Sherlock
Nov 15 '05 #6

P: n/a
Skybuck Flying wrote:
Though the < 0.0000001 is new which I dont quite understand ;)


Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding a
fudge factor.

Cheers,
Nicholas Sherlock
Nov 15 '05 #7

P: n/a

"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote:
Besides it says OnLine ?

I need OnLineSegment...

See the difference ?


The other code I posted does segment.


Actually I have a confesion to make as well ;)

The second method (the C method is flawed).

I am writing a test/measure program, I still need to input some more
verification data ;)

Tomorrow or so it will be done and then I will post it so everybody can
verify and test there goddamn methods lol =D

Bye,
Skybuck.
Nov 15 '05 #8

P: n/a
* Removed from followup: comp.lang.asm sci.electronics.design

In article <dg**********@lust.ihug.co.nz>, Nicholas Sherlock
<n_********@hotmail.com> wrote:

|| Nicholas Sherlock wrote:
|| > Skybuck Flying wrote:
|| >
|| >> I needed a method to determine if a point was on a line segment in 2D.
|| >> So I
|| >> googled for some help and so far I have evaluated two methods.
|| >
|| > function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
|| > begin
|| > result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001;
|| > end;
||
|| Sorry, I missed some of this code from my app:
||
|| function PointOnLine2D(x1, y1, x2, y2, x3, y3: double): boolean;
|| begin
|| result :=
|| (((x3 >= x1) and (x3 <= x2)) or ((x3 >= x2) and (x3 <= x1))) and
|| (((y3 >= y1) and (y3 <= y2)) or ((y3 >= y2) and (y3 <= y1))) and
|| (abs(((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1))) < 0.00001);
|| end;
||
|| I had to rip it out of some other bit of code, so check it carefully
|| against your data.
||
|| Cheers,
|| Nicholas Sherlock

Let's leave the two bounds out for clarity:
((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < epsilon
Distribute:
+x2y3 -x2y1 -x1y3 +x1y1 -x3y2 +x3y1 +x1y2 -x1y1 < epsilon
Cancel:
+x2y3 -x2y1 -x1y3 -x3y2 +x3y1 +x1y2 < epsilon
Reorder:
-x1y3 +x1y2 -x2y1 +x3y1 +x2y3 -x3y2 < epsilon
Inverse distribute:
(y2-y3)*x1 + (x3-x2)*y1 + x2*y3-x3*y2 < epsilon
Name constants:
a*x1 + b*y1 + c < epsilon

a,b,c are constant. For a given line, you can precompute:

a = y2-y3
b = x3-x2
c = x2*y3 - x3*y2

This helps when you want to check a lot of points against the same
line, because it involves only two multiplications, two additions and
a comparison against epsilon. The original formula had 5 additions
(or subtractions) instead.

It corresponds to the general formula for a line in 2d:

ax + by + c = 0

Incidentally, vector (a,b) is a normal vector of the line: perpendicular
to the direction of the line.

For checking the bounds of the line segment, just precomputing the
bounding box and checking against that will do.

Ciao. Vincent.
--
Vincent Zweije <zw****@xs4all.nl> | "If you're flamed in a group you
<http://www.xs4all.nl/~zweije/> | don't read, does anybody get burnt?"
[Xhost should be taken out and shot] | -- Paul Tomblin on a.s.r.
Nov 15 '05 #9

P: n/a

"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote:
Though the < 0.0000001 is new which I dont quite understand ;)


Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding a
fudge factor.


Yes I figured that much. So this would mean the function returns true if
abs(blabla) is near zero ?

Why use 0.00001 why not 0.0000001 or 0.00000000000001 ?

Personally I don't like epsilons like this... it leaves room for
error/malfunction ?

function PointOnLine2D(x1, y1, x2, y2, x3, y3: double): boolean;
begin
result :=
(((x3 >= x1) and (x3 <= x2)) or ((x3 >= x2) and (x3 <= x1))) and
(((y3 >= y1) and (y3 <= y2)) or ((y3 >= y2) and (y3 <= y1))) and
(abs(((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1))) < 0.00001);
end;

I have asked a question on the math newsgroup and also algorithms newsgroup
if the points of the line segment are to be considered part of the line
segment.

This code considers the points to be part of the line segment seeing the >=
and <= instead of > and <.

I myself choose to exclude the points from the line segment. Since I will
probably have to treat overlapping points as a special case, so those will
have to be specially detected etc by another routine to keep things simple
;) As somebody else already pointed out it's free to interpretation ? I
wonder if others agree with that, so far I have only see one reply to it.

Anyway this construction is kinda interesting since it contains
controversial/oppositing conditions so not all conditions will be checked
since some will jump out/forward prematurely etc... I wonder what the worst
case would be... Maybe 8 compares ? In that case it would be worse than both
other methods which only have 7 to 5 compares. (7 for the C method which is
like this one... and only 5 for my version ^^ =D )

Bye,
Skybuck.
Nov 15 '05 #10

P: n/a

"Vincent Zweije" <vz*****@sense.znet> wrote in message
news:dg**********@love.zweije.nl.eu.org...
* Removed from followup: comp.lang.asm sci.electronics.design

In article <dg**********@lust.ihug.co.nz>, Nicholas Sherlock
<n_********@hotmail.com> wrote:

|| Nicholas Sherlock wrote:
|| > Skybuck Flying wrote:
|| >
|| >> I needed a method to determine if a point was on a line segment in 2D. || >> So I
|| >> googled for some help and so far I have evaluated two methods.
|| >
|| > function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
|| > begin
|| > result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001; || > end;
||
|| Sorry, I missed some of this code from my app:
||
|| function PointOnLine2D(x1, y1, x2, y2, x3, y3: double): boolean;
|| begin
|| result :=
|| (((x3 >= x1) and (x3 <= x2)) or ((x3 >= x2) and (x3 <= x1))) and
|| (((y3 >= y1) and (y3 <= y2)) or ((y3 >= y2) and (y3 <= y1))) and
|| (abs(((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1))) < 0.00001); || end;
||
|| I had to rip it out of some other bit of code, so check it carefully
|| against your data.
||
|| Cheers,
|| Nicholas Sherlock

Let's leave the two bounds out for clarity:
((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < epsilon
Distribute:
+x2y3 -x2y1 -x1y3 +x1y1 -x3y2 +x3y1 +x1y2 -x1y1 < epsilon
Cancel:
+x2y3 -x2y1 -x1y3 -x3y2 +x3y1 +x1y2 < epsilon
Reorder:
-x1y3 +x1y2 -x2y1 +x3y1 +x2y3 -x3y2 < epsilon
Inverse distribute:
(y2-y3)*x1 + (x3-x2)*y1 + x2*y3-x3*y2 < epsilon
Name constants:
a*x1 + b*y1 + c < epsilon

a,b,c are constant. For a given line, you can precompute:

a = y2-y3
b = x3-x2
c = x2*y3 - x3*y2

This helps when you want to check a lot of points against the same
line, because it involves only two multiplications, two additions and
a comparison against epsilon. The original formula had 5 additions
(or subtractions) instead.

It corresponds to the general formula for a line in 2d:

ax + by + c = 0

Incidentally, vector (a,b) is a normal vector of the line: perpendicular
to the direction of the line.

For checking the bounds of the line segment, just precomputing the
bounding box and checking against that will do.
Ok, while you make some fine points/tips, only theory and stuff like that
won't do.

For example, I consider the ammount of comparisions also very important for
performance reasons.

Also the way it is implemented is important for branch prediction.

So the only persons getting a cooky is those persons that provide a full
correctly functioning implementation.

Tomorrow I will post a test program to test for all kinds of special cases
etc.

( I have done some prematurely testing and I think method 2 is flawed etc...
maybe some slight modifications can fix it... but more on that tomorrow ;) )

It would be cool if you or anybody else could supply different
implementations of this routine etc... in c, delphi, java, or any other
language... as long as it can hopefully be converted back to pascal/delphi
;)

Ofcourse everybody should make an attempt to optimize it as good as
possible.

And yes I must give you some good credits because you did mention some good
tricks/optimizations... like precomputing the lines and bounding boxes
etc....

Though on the other hand this would require the routine to be integrated
with external algorithms etc...

Soooooooooooooooooo that's not really ideal etc... since then everything has
to be changed, and re-checked and re-tested etc.

Soooo I could consider it to be a little bit cheating ;)

Yes I am going to vote it to be cheating for now.... no pre-emptive look ups
etc...

The routine should be a compact/atomic/reusable entity/routine which can be
used for any algorithm straight out of the box and any data etc.

Ofcourse you are free to implement caching techniques in case you do want to
try and do look-ups etc..

Those are ofcourse not cheating etc... but everything has to be inside a
simple routine like this:

function ( x1,y1, x2,y2, x3,y3 : double ) : boolean;
begin

<anything inside here goes> ;)

end;

For your sake and optimizations sake I will allow object oriented routines
as well:

For example

type
Tgeometry public
<anything inside here is allowed>
end;

function Tgeometry.PointOnLineSegment( x1,y1, x2,y2, x3,y3 : double ) :
boolean;
begin

<anything inside here is allowed>

end;

In case of pascal I will also allow unit variables which are scope limited
as follows

unit

interface

< header>

implementation

var
<anything inside here is allowed>

<implementation>

end.

But ofcourse everything then has to be in a seperate unit ;)
So this means your lookup technique is possible and a simple example would
be:
if (X1=StoredX1) and (Y1=StoredY1) and
(X2=StoredX2) and (Y2=StoredY2) then
begin
// use whatever it is you cached...

// or
if (PX = StoredPX) and (PY=StoredPY) then
begin
Return StoredSolution;
end;
end;

Etc....

Ofcourse this would require all kinds of lookups and assignments... so I
dont think a cache like this will speed things up or be helpfull etc...
ofcourse you are free to proof me wrong.

In the end all methods/routines/implementations will be verified and
measured performance/time wise ;)

Bye,
Skybuck.

Ciao. Vincent.
--
Vincent Zweije <zw****@xs4all.nl> | "If you're flamed in a group you
<http://www.xs4all.nl/~zweije/> | don't read, does anybody get burnt?" [Xhost should be taken out and shot] | -- Paul Tomblin on

a.s.r.
Nov 15 '05 #11

P: n/a
"Skybuck Flying" <no****@hotmail.com> wrote in message
news:dg**********@news6.zwoll1.ov.home.nl...
"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote:
Though the < 0.0000001 is new which I dont quite understand ;)


Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding
a fudge factor.


Yes I figured that much. So this would mean the function returns true
if abs(blabla) is near zero ?

Why use 0.00001 why not 0.0000001 or 0.00000000000001 ?

Personally I don't like epsilons like this... it leaves room for
error/malfunction ?


Some real life logic: once in a graphics rendering lab assignment, we
were instructed to approximate Bezier curves. This is an iterative
process; the stop condition was for the error to become less then
half a pixel. For visualisation on a raster display, that is exactly
what's required.

Groetjes,
Maarten Wiltink
Nov 15 '05 #12

P: n/a

"Maarten Wiltink" <ma*****@kittensandcats.net> wrote in message
news:43***********************@news.xs4all.nl...
"Skybuck Flying" <no****@hotmail.com> wrote in message
news:dg**********@news6.zwoll1.ov.home.nl...
"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote: Though the < 0.0000001 is new which I dont quite understand ;)

Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding
a fudge factor.


Yes I figured that much. So this would mean the function returns true
if abs(blabla) is near zero ?

Why use 0.00001 why not 0.0000001 or 0.00000000000001 ?

Personally I don't like epsilons like this... it leaves room for
error/malfunction ?


Some real life logic: once in a graphics rendering lab assignment, we
were instructed to approximate Bezier curves. This is an iterative
process; the stop condition was for the error to become less then
half a pixel. For visualisation on a raster display, that is exactly
what's required.


That certainly won't do in this case etc.

It should be as exact/accurate as possible.

Bye,
Skybuck.
Nov 15 '05 #13

P: n/a
Skybuck Flying wrote:
"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
(snip)
Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding a
fudge factor.

Yes I figured that much. So this would mean the function returns true if
abs(blabla) is near zero ? Why use 0.00001 why not 0.0000001 or 0.00000000000001 ? Personally I don't like epsilons like this... it leaves room for
error/malfunction ?


Then don't use floating point math.

As you don't indicate the origin of the problem, we can't help any
more than that. In some cases it can be done in fixed point,
especially if the segment ends and point being checked are fixed point
values.

Otherwise, if you select two end points somewhat randomly there is a
high probability that no floating point values lie on the segment
between them.

-- glen

Nov 15 '05 #14

P: n/a
Skybuck Flying wrote:

(snip)
That certainly won't do in this case etc. It should be as exact/accurate as possible.


You really need to tell us what "this case" is.

You expect everyone to help you, but don't seem
interested in helping much yourself.

-- glen

Nov 15 '05 #15

P: n/a
Skybuck Flying wrote:
"Hans-Bernhard Broeker" <br*****@physik.rwth-aachen.de> wrote in message
news:3p************@news.dfncis.de...
In comp.graphics.algorithms Skybuck Flying <no****@hotmail.com> wrote:

[Massive crossposting reduced to single group. Please, Skybuck, try
to follow netiquette a bit better in the future. If you feel a need
to X-post, at be a bit more selective to where --- this had *no*
business being posted to comp.arch or sci.electronics. And *always*
select a single group for the follow-ups.]
No, I like to have input from all those newsgroups since this is about
optimizations and those newsgroup are related to software and hardware
optimizations. C/Delphi/Asm for software/cpu optimizations and
comp.arch/asm/eletronics about hardware optimizations and algorithms
possibly about algorithm optimization. Those newsgroups are likely to
contain people who know something about this.


I'm with Hans on this one. Skybuck, by your logic you should have
posted to ALL newsgroups, since there's probably someone in every
newsgroup who will be somehow "related to optimization". :-(

Comp.graphics.algorithms was the right place. Comp.lang.c is hard to
justify, since your question has nothing to do with C. But there's no
way to justify posting to sci.electronics.design, comp.arch, or
alt.comp.lang.borland-delphi. Nor comp.lang.asm, unless you want a
response written in assembly language.

That said, other appropriate forums might have included comp.programming
or *maybe* comp.theory.

....

You need to define the task a good deal more cleanly before an
implementation can reasonably be suggested. Missing information:

1) What's a line segment, in this case? It may come as a surprise to
you, but ideas about the details of this actually differ. So: how do
you actually specify this line segment?

A 2D line segment defined by (x1,y1) - (x2,y2)


Integer or floating point? In either case, how big are your points and
how fat is the line? Is it acceptable to return a false positive or
negative? If so, how often? If it's not acceptable, then you have a
problem. On a any machine that lacks infinite numerical precision, you
can never know whether an infinitely small point lies on an infinitely
thin line segment.

In the end, you're going to have to do something like this:

"Draw two lines, one from each point in the line segment to the
candidate point. If these lines' slopes are equal but differ in sign
(e.g. X and -X), then your point lies on the segment."

Of course, what does it mean for slope values to be "equal" on a digital
computer? Unless you can represent the problem entirely symbolically,
"equal" will be necessarily limited by the representation of the data
and the accuracy of the math. In the real world, epsilon is a necessary
evil.
2) When exactly do you want the test to return "true"? I.e.: how
close to the ideal segment does a point have to be to still be
considered "on" it?

As precise as possible. (Floating point precision, no rounding)


This is an oxymoron. Floating point is inescapably about rounding.
Digital representations of floating point are approximations, as are the
math operations that manipulate them. The programmer has to manage
numerical imprecision on computers or live with wrong answers.

For some perspective, read David Goldberg's "What Every Computer
Scientist Should Know About Floating-Point Arithmetic":

http://docs.sun.com/source/806-3568/ncg_goldberg.html

Randy
Nov 15 '05 #16

P: n/a
[some off-topic newsgroups trimmed]

Skybuck Flying wrote:
"Maarten Wiltink" <ma*****@kittensandcats.net> wrote:
"Skybuck Flying" <no****@hotmail.com> wrote:
"Nicholas Sherlock" <n_********@hotmail.com> wrote:
Skybuck Flying wrote:

>Though the < 0.0000001 is new which I dont quite understand ;)

Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding
a fudge factor.

Yes I figured that much. So this would mean the function returns true
if abs(blabla) is near zero ?

Why use 0.00001 why not 0.0000001 or 0.00000000000001 ?

Personally I don't like epsilons like this... it leaves room for
error/malfunction ?


Some real life logic: once in a graphics rendering lab assignment, we
were instructed to approximate Bezier curves. This is an iterative
process; the stop condition was for the error to become less then
half a pixel. For visualisation on a raster display, that is exactly
what's required.


That certainly won't do in this case etc.

It should be as exact/accurate as possible.


Even trying to detect whether a point is on a line using floating point
arithmetic almost certainly indicates a specification error. It's impossible
to do that exactly.

--
David Hopwood <da******************@blueyonder.co.uk>
Nov 15 '05 #17

P: n/a
"Skybuck Flying" <no****@hotmail.com> writes:
It should be as exact/accurate as possible.


Then use rational numbers?

-k
--
If I haven't seen further, it is by standing in the footprints of giants
Nov 15 '05 #18

P: n/a

"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:gv********************@comcast.com...
Skybuck Flying wrote:

(snip)
That certainly won't do in this case etc.

It should be as exact/accurate as possible.


You really need to tell us what "this case" is.


You ll see quite soon the test program is done ;)

Everybody should simply do his/her best at "as accurate" as possible.

Though it should also be fast.

So you can use single or double floating point format.

Whatever floats your boat ;)

Though I would suggest doubles since those can handle higher/better
precision.. bigger/smaller numbers etc.

Using epsilon's and stuff like should probably be prevented as much as
possible since those are common forms of errors etc and limitations etc...
;)

Bye,
Skybuck.
Nov 15 '05 #19

P: n/a

"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:dP********************@comcast.com...
Skybuck Flying wrote:
"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...


(snip)
Floating point math is always inaccurate as most numbers cannot be
exactly represented. This bit basically takes care of that by adding a
fudge factor.

Yes I figured that much. So this would mean the function returns true if
abs(blabla) is near zero ?

Why use 0.00001 why not 0.0000001 or 0.00000000000001 ?

Personally I don't like epsilons like this... it leaves room for
error/malfunction ?


Then don't use floating point math.

As you don't indicate the origin of the problem, we can't help any
more than that. In some cases it can be done in fixed point,
especially if the segment ends and point being checked are fixed point
values.

Otherwise, if you select two end points somewhat randomly there is a
high probability that no floating point values lie on the segment
between them.


You are more then welcome to prove this very soon.
I'll shall do or attempt to do two things:

"allow dll plugins" for the test program and
"allow data/verification" files for the test program.

As to provide a binary and test files for those people who don't have a
pascal compiler, or can't program or don't want to program or just want to
supply some verification data =D

It's gonna be quite cool lol.

Bye,
Skybuck.
Nov 15 '05 #20

P: n/a
Skybuck Flying wrote:
Using epsilon's and stuff like should probably be prevented as much as
possible since those are common forms of errors etc and limitations
etc...


.... otherwise you might wind up learning how to use
relevant features of the C standard library.

--
pete
Nov 15 '05 #21

P: n/a
"Skybuck Flying" <no****@hotmail.com> wrote:
Using epsilon's and stuff like should probably be prevented as much as
possible since those are common forms of errors etc and limitations etc...


You have a gross and dangerous misunderstanding of principles of doing
numerical processing with the computer. I suggest you study the
fundmentals carefully first before wasting your time writing what will
undoubtedly be some ridicolously naive and rather useless code.

Nov 15 '05 #22

P: n/a
Skybuck Flying wrote:
Hi,

I needed a method to determine if a point was on a line segment in 2D. So I
googled for some help and so far I have evaluated two methods.


Here's another one to play with. It goes from the assumption that if
the slope from p1 to p2 is the same as the slope from p1 to p3, then p3
is colinear with p1 and p2. Then it checks x coordinates to see if p3
is on the segment between p1 and p2.

It expresses slope as rational numbers, so the arithmetic is integral.
It probably needs to be bulletproofed.

#include <stdio.h>

typedef struct {
long x;
long y;
} point_t;

typedef struct {
long num;
long denom;
} rational_t;

#define REDUCE(frac) \
do { \
long e1=frac.num, \
e2=frac.denom, \
quot, rem = -1; \
long tmp; \
int done = 0; \
while(!done) \
{ \
if (e1 < e2) \
{ \
tmp = e1; \
e1 = e2; \
e2 = tmp; \
} \
quot = e1/e2; \
rem = e1 % e2; \
if (rem) \
e1 = rem; \
else \
done = 1; \
} \
frac.num /= e2; \
frac.denom /= e2; \
} while(0)

/**
* We are assuming that p1 and p2 have been ordered
* so that p1.x <= p2.x
*/
int isOnLine(point_t p1, point_t p2, point_t p3)
{
rational_t m1 = {p3.y - p1.y, p3.x - p1.x};
rational_t m2 = {p2.y - p1.y, p2.x - p1.x};

/**
* special case -- p1.x == p2.x
*/
if (p1.x == p2.x)
{
return p3.x == p1.x &&
(p1.y <= p3.y && p3.y <= p2.y ||
p1.y >= p3.y && p3.y >= p2.y);
}

/**
* Reduce both fractions before comparing. This is where
* any performance issues would be.
*/
REDUCE(m1);
REDUCE(m2);

if (m1.num == m2.num && m1.denom == m2.denom)
{
return p1.x <= p3.x && p3.x <= p2.x;
}

return 0;
}

Nov 15 '05 #23

P: n/a
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:dP********************@comcast.com...
(snip)
Otherwise, if you select two end points somewhat randomly there is a
high probability that no floating point values lie on the segment
between them.

You are more then welcome to prove this very soon.
I'll shall do or attempt to do two things:


I pick the points (0,0) and (512,997) slightly easier to see as
integers, but you can shift the binary point over and make them floating
point. Find any point on the segment between them.

-- glen

Nov 15 '05 #24

P: n/a
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote:
Skybuck Flying wrote:
It should be as exact/accurate as possible.
You really need to tell us what "this case" is.


You ll see quite soon the test program is done ;)


Ok, you are pissing people off with statements like this.

Here's the thing, the "epsilon thing" is far and away the most
practical way of doing this; but just as a matter of correctness, you
actually should not choose a constant epsilon. The correct epsilon
will depend on the magnitude of the coordinates for the original
segment points. For example it would be easy to construct an example
where the best correct epsilon is a million.

Ok, perhaps a better way to ask you for more information, is this.From what sort of *SOURCE* are your input points comming from? Are

they just chosen literally at random? (Using say, the Mersenne Twister
random number generator, which includes floating point random.) Or
(more likely) are they actually constructed to be *ON* the line with
some likelihood, but for some reason you loose track of this fact, and
wish to recalculate it?

This is important because, the *WAY* you construct the point (say, but
being given the x intercept, then computing the y that is supposed to
correspond to it) may actually give an exact matching algorithm without
the need for any epsilons. One could, for example, just try to
reproduce the procedure for finding the point, from one of the
coordinates, and see if the other matches exactly.

The problem with this is that starting with an x and computing a y, or
the other way around will not necessarily yield the same results.
Further more if you compute the point as (alpha * p_0 + (1-alpha)* p_1)
(which might be *WAY* harder to match -- intuitively it seems so to me,
but I don't have a 100% clear idea), you will again get different LSB
round offs.

So I hope you understand that these accuracy issues are pretty integral
to what it is that you want to do. Without more information from you,
I don't believe that anyone can suggest anything else more with any
expectation of it necessarily working better.

--
Paul Hsieh
http://www.pobox.com/~qed/
http://bstring.sf.net/

Nov 15 '05 #25

P: n/a

<we******@gmail.com> wrote in message
news:11**********************@g49g2000cwa.googlegr oups.com...
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote:
Skybuck Flying wrote:
> It should be as exact/accurate as possible.

You really need to tell us what "this case" is.


You ll see quite soon the test program is done ;)


Ok, you are pissing people off with statements like this.

Here's the thing, the "epsilon thing" is far and away the most
practical way of doing this; but just as a matter of correctness, you
actually should not choose a constant epsilon. The correct epsilon
will depend on the magnitude of the coordinates for the original
segment points. For example it would be easy to construct an example
where the best correct epsilon is a million.

Ok, perhaps a better way to ask you for more information, is this.
From what sort of *SOURCE* are your input points comming from? Are

they just chosen literally at random? (Using say, the Mersenne Twister
random number generator, which includes floating point random.) Or
(more likely) are they actually constructed to be *ON* the line with
some likelihood, but for some reason you loose track of this fact, and
wish to recalculate it?

This is important because, the *WAY* you construct the point (say, but
being given the x intercept, then computing the y that is supposed to
correspond to it) may actually give an exact matching algorithm without
the need for any epsilons. One could, for example, just try to
reproduce the procedure for finding the point, from one of the
coordinates, and see if the other matches exactly.

The problem with this is that starting with an x and computing a y, or
the other way around will not necessarily yield the same results.
Further more if you compute the point as (alpha * p_0 + (1-alpha)* p_1)
(which might be *WAY* harder to match -- intuitively it seems so to me,
but I don't have a 100% clear idea), you will again get different LSB
round offs.

So I hope you understand that these accuracy issues are pretty integral
to what it is that you want to do. Without more information from you,
I don't believe that anyone can suggest anything else more with any
expectation of it necessarily working better.


The test program is "done" well not really but a first version is available
on the net.

VerificationProgramVersion009.zip

https://sourceforge.net/project/show...group_id=53726

See the other usenet thread called "VerificationProgram009 etc" for more
details etc ;)

Bye,
Skybuck.

Nov 15 '05 #26

P: n/a

"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:0t********************@comcast.com...
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:dP********************@comcast.com...


(snip)
Otherwise, if you select two end points somewhat randomly there is a
high probability that no floating point values lie on the segment
between them.

You are more then welcome to prove this very soon.
I'll shall do or attempt to do two things:


I pick the points (0,0) and (512,997) slightly easier to see as
integers, but you can shift the binary point over and make them floating
point. Find any point on the segment between them.


The test program is "done" well not really but a first version is available
on the net.

Which also includes your line. "My" second implementation has problems with
this line. The second method uses some sort of dot product or cross product.
Which is kinda surprising because I think Nicholas's implementation also
works like that... oh well ;)

My first implementation however... which I do fully understand still works
perfectly =D

Can you find/think up any other possibly problems ? ;)

See this link for the source code, binaries, data, etc, etc, etc.

VerificationProgramVersion009.zip

https://sourceforge.net/project/show...group_id=53726

See the other usenet thread called "VerificationProgram009 etc" for more
details etc ;)

Bye,
Skybuck.
Nov 15 '05 #27

P: n/a

"John Bode" <jo*******@my-deja.com> wrote in message
news:11*********************@z14g2000cwz.googlegro ups.com...
Skybuck Flying wrote:
Hi,

I needed a method to determine if a point was on a line segment in 2D. So I googled for some help and so far I have evaluated two methods.


Here's another one to play with. It goes from the assumption that if
the slope from p1 to p2 is the same as the slope from p1 to p3, then p3
is colinear with p1 and p2. Then it checks x coordinates to see if p3
is on the segment between p1 and p2.

It expresses slope as rational numbers, so the arithmetic is integral.
It probably needs to be bulletproofed.

#include <stdio.h>

typedef struct {
long x;
long y;
} point_t;

typedef struct {
long num;
long denom;
} rational_t;

#define REDUCE(frac) \
do { \
long e1=frac.num, \
e2=frac.denom, \
quot, rem = -1; \
long tmp; \
int done = 0; \
while(!done) \
{ \
if (e1 < e2) \
{ \
tmp = e1; \
e1 = e2; \
e2 = tmp; \
} \
quot = e1/e2; \
rem = e1 % e2; \
if (rem) \
e1 = rem; \
else \
done = 1; \
} \
frac.num /= e2; \
frac.denom /= e2; \
} while(0)

/**
* We are assuming that p1 and p2 have been ordered
* so that p1.x <= p2.x
*/
int isOnLine(point_t p1, point_t p2, point_t p3)
{
rational_t m1 = {p3.y - p1.y, p3.x - p1.x};
rational_t m2 = {p2.y - p1.y, p2.x - p1.x};

/**
* special case -- p1.x == p2.x
*/
if (p1.x == p2.x)
{
return p3.x == p1.x &&
(p1.y <= p3.y && p3.y <= p2.y ||
p1.y >= p3.y && p3.y >= p2.y);
}

/**
* Reduce both fractions before comparing. This is where
* any performance issues would be.
*/
REDUCE(m1);
REDUCE(m2);

if (m1.num == m2.num && m1.denom == m2.denom)
{
return p1.x <= p3.x && p3.x <= p2.x;
}

return 0;
}


Hi, cool stuff. I haven't had the time yet to look into this.

It would help if someone could make it more suited for my test program and
maybe make a nice little dll ?

Anyway maybe somebody can convert this to delphi as well.

The test program is "done" well not really but a first version is available
on the net.

VerificationProgramVersion009.zip

https://sourceforge.net/project/show...group_id=53726

See the other usenet thread called "VerificationProgram009 etc" for more
details etc ;)

Bye,
Skybuck.
Nov 15 '05 #28

P: n/a
Why don't you make a region slightly larger than the line, and use
PtInRegion to check if the click is on the line.

var
hndRgn : hRGN;

const
LineSt : TPoint = (X:0; Y:0);
LineEnd : TPoint = (X:997; Y:512);
LineWidth : integer = 2;

procedure TForm1.FormCreate(Sender: TObject);
var
PointArray : array[0..3] of TPoint;
LW : integer;
begin
LW := LineWidth;
PointArray[0] := Point(LineSt.X - LW, LineSt.Y - LW);
PointArray[1] := Point(LineEnd.X + LW, LineEnd.Y - LW);
PointArray[2] := Point(LineEnd.X + LW, LineEnd.Y + LW);
PointArray[3] := Point(LineSt.X - LW, LineSt.Y + LW);
hndRgn := CreatePolygonRgn(PointArray, 4, WINDING);
end;

procedure TForm1.FormMouseUp(Sender: TObject; Button: TMouseButton;
Shift: TShiftState; X, Y: Integer);
var
OnLine : boolean;
begin
StartPeriod; // timer start
OnLine := PtInRegion(hndRgn, X, Y);
EndPeriod; // timer stop
{indicate if on line}
if OnLine then
Shape1.Brush.Color := clLime
else
Shape1.Brush.Color := clRed;
{display time to check in region}
Label1.Caption := Format('%d microsecs',
[trunc(ScaleTime(ElapsedTime, ttMicroSecs))]);
end;

This gives about 13 microsecs for the first click (2 microsecs for
later ones) at the start of the line, and 17 microsecs for the first
click (7 microsecs for later ones) at the end of the line. This is with
a 2.5GHz PC.

I can't believe that this is too slow for you <g>.

Alan Lloyd

Nov 15 '05 #29

P: n/a
In comp.arch glen herrmannsfeldt <ga*@ugcs.caltech.edu> wrote:
Skybuck Flying wrote:

(snip)
That certainly won't do in this case etc.
It should be as exact/accurate as possible.


You really need to tell us what "this case" is.

You expect everyone to help you, but don't seem
interested in helping much yourself.


Its rather clear he is interested in generating lots of usenet traffic
and not much else...

-- glen


--
Sander

+++ Out of cheese error +++
Nov 15 '05 #30

P: n/a

"John Bode" <jo*******@my-deja.com> wrote in message
news:11*********************@z14g2000cwz.googlegro ups.com...
Skybuck Flying wrote:
Hi,

I needed a method to determine if a point was on a line segment in 2D. So I googled for some help and so far I have evaluated two methods.


Here's another one to play with. It goes from the assumption that if
the slope from p1 to p2 is the same as the slope from p1 to p3, then p3
is colinear with p1 and p2. Then it checks x coordinates to see if p3
is on the segment between p1 and p2.

It expresses slope as rational numbers, so the arithmetic is integral.
It probably needs to be bulletproofed.

#include <stdio.h>

typedef struct {
long x;
long y;
} point_t;

typedef struct {
long num;
long denom;
} rational_t;

#define REDUCE(frac) \
do { \
long e1=frac.num, \
e2=frac.denom, \
quot, rem = -1; \
long tmp; \
int done = 0; \
while(!done) \
{ \
if (e1 < e2) \
{ \
tmp = e1; \
e1 = e2; \
e2 = tmp; \
} \
quot = e1/e2; \
rem = e1 % e2; \
if (rem) \
e1 = rem; \
else \
done = 1; \
} \
frac.num /= e2; \
frac.denom /= e2; \
} while(0)

/**
* We are assuming that p1 and p2 have been ordered
* so that p1.x <= p2.x
*/
int isOnLine(point_t p1, point_t p2, point_t p3)
{
rational_t m1 = {p3.y - p1.y, p3.x - p1.x};
rational_t m2 = {p2.y - p1.y, p2.x - p1.x};

/**
* special case -- p1.x == p2.x
*/
if (p1.x == p2.x)
{
return p3.x == p1.x &&
(p1.y <= p3.y && p3.y <= p2.y ||
p1.y >= p3.y && p3.y >= p2.y);
}

/**
* Reduce both fractions before comparing. This is where
* any performance issues would be.
*/
REDUCE(m1);
REDUCE(m2);

if (m1.num == m2.num && m1.denom == m2.denom)
{
return p1.x <= p3.x && p3.x <= p2.x;
}

return 0;
}


Ok, your method has been added to the verification program.

However during verification the algorithm hangs at this data:

Verifieing:
Type index: 1
Type description: Diagonal line with point inside line segment
Data index: 2
Data description: positive area, third crossed

Which is this one:

AddGeneratedVerificationData( 'positive area, third crossed', 200,2000,
3000,100, 0.3, 0.0 );

Which means:
X1 := 200;
Y1 := 2000;
X2 := 3000;
Y2 := 100;
PX := 1040; // generated/calculated
PY := 1430; // generated/calculated

I put your source code in a library/DLL which looks like this:

JOHNBODERATIONALMETHOD_API BOOL __stdcall PointOnLineSegmentIn2D(
double StartX, double StartY,
double EndX, double EndY,
double PointX, double PointY )
{
point_t p1;
point_t p2;
point_t p3;
if (StartX <= EndX)
{
p1.x = long(StartX); // typecast to prevent loss of data warning
p1.y = long(StartY); // typecast to prevent loss of data warning
p2.x = long(EndX); // typecast to prevent loss of data warning
p2.y = long(EndY); // typecast to prevent loss of data warning
} else
{
p1.x = long(EndX); // typecast to prevent loss of data warning
p1.y = long(EndY); // typecast to prevent loss of data warning
p2.x = long(StartX); // typecast to prevent loss of data warning
p2.y = long(StartY); // typecast to prevent loss of data warning
}
p3.x = long(PointX); // typecast to prevent loss of data warning
p3.y = long(PointY); // typecast to prevent loss of data warning
return isOnLine(p1,p2,p3);
}

It might be hanging because Y2 is smaller than Y1 ? or maybe there is
another problem ? ;)

Bye,
Skybuck.
Nov 15 '05 #31

P: n/a
Hi John, your posted method has now also been added to the
program/project/etc in
VerificationProgram0.10 for PointOnLineSegmentIn2D =D
(Though the dll has been disabled by simply renaming it until the hang/bug
is fixed ;) )

Hello Peeps.

Here is what's new ;)

version 0.10 created on 18 and 19 september 2005 by Skybuck Flying

+ Some code moved to other/new units etc.
+ New source and binary implementation added: Nils from www.simdesign.nl
+ New source and binary implementation added: Skybuck Flying method 3
+ New binary implementation added: JohnBodeRationalMethod (source in visual
cpp).
+ ProjectGroups added
+ Each verification is displayed in case a binary plug in hangs.. then at
least
it's possible to see at which verification it hangs.

Get the source and binaries from:

https://sourceforge.net/project/show...group_id=53726

=D

Bye,
Skybuck.
"John Bode" <jo*******@my-deja.com> wrote in message
news:11*********************@z14g2000cwz.googlegro ups.com...
Skybuck Flying wrote:
Hi,

I needed a method to determine if a point was on a line segment in 2D. So I googled for some help and so far I have evaluated two methods.


Here's another one to play with. It goes from the assumption that if
the slope from p1 to p2 is the same as the slope from p1 to p3, then p3
is colinear with p1 and p2. Then it checks x coordinates to see if p3
is on the segment between p1 and p2.

It expresses slope as rational numbers, so the arithmetic is integral.
It probably needs to be bulletproofed.

#include <stdio.h>

typedef struct {
long x;
long y;
} point_t;

typedef struct {
long num;
long denom;
} rational_t;

#define REDUCE(frac) \
do { \
long e1=frac.num, \
e2=frac.denom, \
quot, rem = -1; \
long tmp; \
int done = 0; \
while(!done) \
{ \
if (e1 < e2) \
{ \
tmp = e1; \
e1 = e2; \
e2 = tmp; \
} \
quot = e1/e2; \
rem = e1 % e2; \
if (rem) \
e1 = rem; \
else \
done = 1; \
} \
frac.num /= e2; \
frac.denom /= e2; \
} while(0)

/**
* We are assuming that p1 and p2 have been ordered
* so that p1.x <= p2.x
*/
int isOnLine(point_t p1, point_t p2, point_t p3)
{
rational_t m1 = {p3.y - p1.y, p3.x - p1.x};
rational_t m2 = {p2.y - p1.y, p2.x - p1.x};

/**
* special case -- p1.x == p2.x
*/
if (p1.x == p2.x)
{
return p3.x == p1.x &&
(p1.y <= p3.y && p3.y <= p2.y ||
p1.y >= p3.y && p3.y >= p2.y);
}

/**
* Reduce both fractions before comparing. This is where
* any performance issues would be.
*/
REDUCE(m1);
REDUCE(m2);

if (m1.num == m2.num && m1.denom == m2.denom)
{
return p1.x <= p3.x && p3.x <= p2.x;
}

return 0;
}

Nov 15 '05 #32

P: n/a
<al********@aol.com> wrote in message
news:11*********************@f14g2000cwb.googlegro ups.com...
Why don't you make a region slightly larger than the line, and use
PtInRegion to check if the click is on the line.
Not too shabby ;)

This would/could be pretty good for gui's where the user needs to pick a
line etc.

That's just one application of this method.

Another application is for massive calculations of all kinds of geometry
algorithms etc.

Your method is kinda funny since it uses windows to do it's thing. Tomorrow
or so I will add it as well to see how it performance, speed wise, compared
to all the other methods =D

That's going to be much fun ;)

I also wonder how stable it would be... probably pretty stable though ;) I
also wonder how large or small the values could be.

So this method works with integers etc... small values would get rounded to
zero or so or one like
0.05 and 0.0002 etc.. and point 0.06 whatever... then this method would
always return true or something like that ;) could be fun and maybe even
handy as well =D

Bye,
Skybuck.

var
hndRgn : hRGN;

const
LineSt : TPoint = (X:0; Y:0);
LineEnd : TPoint = (X:997; Y:512);
LineWidth : integer = 2;

procedure TForm1.FormCreate(Sender: TObject);
var
PointArray : array[0..3] of TPoint;
LW : integer;
begin
LW := LineWidth;
PointArray[0] := Point(LineSt.X - LW, LineSt.Y - LW);
PointArray[1] := Point(LineEnd.X + LW, LineEnd.Y - LW);
PointArray[2] := Point(LineEnd.X + LW, LineEnd.Y + LW);
PointArray[3] := Point(LineSt.X - LW, LineSt.Y + LW);
hndRgn := CreatePolygonRgn(PointArray, 4, WINDING);
end;

procedure TForm1.FormMouseUp(Sender: TObject; Button: TMouseButton;
Shift: TShiftState; X, Y: Integer);
var
OnLine : boolean;
begin
StartPeriod; // timer start
OnLine := PtInRegion(hndRgn, X, Y);
EndPeriod; // timer stop
{indicate if on line}
if OnLine then
Shape1.Brush.Color := clLime
else
Shape1.Brush.Color := clRed;
{display time to check in region}
Label1.Caption := Format('%d microsecs',
[trunc(ScaleTime(ElapsedTime, ttMicroSecs))]);
end;

This gives about 13 microsecs for the first click (2 microsecs for
later ones) at the start of the line, and 17 microsecs for the first
click (7 microsecs for later ones) at the end of the line. This is with
a 2.5GHz PC.
Nice one ;)

I can't believe that this is too slow for you <g>.


For one click no :)

For massive calculations I dont know lol =D

Bye,
Skybuck.
Nov 15 '05 #33

P: n/a
we******@gmail.com wrote
(in article
<11**********************@g49g2000cwa.googlegroups .com>):
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote:
Skybuck Flying wrote:
It should be as exact/accurate as possible.

You really need to tell us what "this case" is.


You ll see quite soon the test program is done ;)


Ok, you are pissing people off with statements like this.


Shouldn't be a surprise. If you search newsgroup archives for
"Skybuck" you'll find his been doing that sort of thing in a lot
of different groups for a long time.
--
Randy Howard (2reply remove FOOBAR)

Nov 15 '05 #34

P: n/a
On Mon, 19 Sep 2005 08:05:41 GMT, Randy Howard
<ra*********@FOOverizonBAR.net> wrote:
we******@gmail.com wrote
(in article
<11**********************@g49g2000cwa.googlegroup s.com>):
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote:
Skybuck Flying wrote:
> It should be as exact/accurate as possible.

You really need to tell us what "this case" is.

You ll see quite soon the test program is done ;)


Ok, you are pissing people off with statements like this.


Shouldn't be a surprise. If you search newsgroup archives for
"Skybuck" you'll find his been doing that sort of thing in a lot
of different groups for a long time.


Yep. I changed his filter from specific newsgroup to global a long
time ago.
--
Al Balmer
Balmer Consulting
re************************@att.net
Nov 15 '05 #35

P: n/a

"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:0t********************@comcast.com...
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:dP********************@comcast.com...


(snip)
Otherwise, if you select two end points somewhat randomly there is a
high probability that no floating point values lie on the segment
between them.

You are more then welcome to prove this very soon.
I'll shall do or attempt to do two things:


I pick the points (0,0) and (512,997) slightly easier to see as
integers, but you can shift the binary point over and make them floating
point. Find any point on the segment between them.


I wrote a little program which uses rational numbers and floating point
numbers and shows the difference.

The rational number version is able to perfectly calculate the point.
The floating point version fails miserably at it ;)

The rational number version is about 20 times slower than the floating point
version ;)

program started

(Rational) PointX in rational notation: 768/5
(Rational) PointY in rational notation: 2991/10

(Rational) PointX in real notation: 153.6000000000000000
(Rational) PointY in real notation: 299.1000000000000000
(Rational) ticks: 282
(Rational) time (in seconds): 0.0000787809623849

(Floating point) PointX: 153.5999999999999940
(Floating point) PointY: 299.0999999999999660
(Floating point) ticks: 14
(Floating point) time (in seconds): 0.0000039111116078

press enter to exit

Here is the source code. Thanks to delphi fundamentals for the Trational
class etc =D

Now for the final question ;)

Can you still find a point which would in theory be on the line segment but
would still be mis-calculated by the rational number version ? ;)

program TestRationalNumbers;

{$APPTYPE CONSOLE}

uses
SysUtils,
cUtils in 'Utils\cUtils.pas',
cStrings in 'Utils\cStrings.pas',
cRational in 'Maths\cRational.pas',
cMaths in 'Maths\cMaths.pas',
Windows;
// rational number versions:

procedure calc_perpendicular( var x, y : Trational );
var
temp_x, temp_y : Trational;
begin
temp_x := Trational.Create;
temp_y := Trational.Create;

// temp_x := -y;
temp_x.Assign( y );
temp_x.Negate;

// temp_y := x;
temp_y.Assign( x );

// x := temp_x;
// y := temp_y;
x.Assign( temp_x );
y.Assign( temp_y );

temp_x.Destroy;
temp_y.Destroy;
end;

procedure calc_normalize( var x, y : Trational );
var
temp_squared_x : Trational;
temp_squared_y : Trational;
temp_length : Trational;

begin
temp_squared_x := Trational.Create;
temp_squared_y := Trational.Create;
temp_length := Trational.Create;

// temp_length := sqrt( (x*x) + (y*y) );
temp_squared_x.Assign( x );
temp_squared_y.Assign( y );

temp_squared_x.Sqr;
temp_squared_y.Sqr;

temp_length.Assign( temp_squared_x );
temp_length.Add( temp_squared_y );

temp_length.Sqrt;

// if temp_length<>0 then
if not temp_length.IsZero then
begin
// x := x / temp_length;
// y := y / temp_length;

x.Divide( temp_length );
y.Divide( temp_length );
end;

temp_squared_x.Destroy;
temp_squared_y.Destroy;
temp_length.Destroy;
end;

procedure calc_normal( _x1,_y1, _x2,_y2 : Trational; var _normal_x,
_normal_y : Trational );
begin
// _normal_x := _x2-_x1;
// _normal_y := _y2-_y1;

_normal_x.Assign( _x2 );
_normal_y.Assign( _y2 );

_normal_x.Subtract( _x1 );
_normal_y.Subtract( _y1 );

// calculate perpendicular "vector".
calc_perpendicular( _normal_x, _normal_y );

// divide vector by it's own length there by creating a "normal" :)
calc_normalize( _normal_x, _normal_y );
end;

// original double versions
// using flawed floating point shit ;)

// helper procedure to generate points on or near line segments
procedure GeneratePointForLineSegment( StartX, StartY,
EndX, EndY : Trational;
DistanceFromStart : Trational; // distance from start
// < 0 is before start
// 0 is on start
// 0.5 is on center
// 1 is one end
// > 1 is after end
DistanceFromLine : Trational; // closest distance from line

// follows the normal
// positive is left side of line
// zero is on line
// negative is right side of line
var px, py : Trational );
var
NX, NY : Trational;

begin
NX := Trational.Create;
NY := Trational.Create;

// theory of method1 ;)
// PX := StartX + DistanceFromStart * (EndX - StartX);
// PY := StartY + DistanceFromStart * (EndY - StartY);

PX.Assign( EndX );
PX.Subtract( StartX );
PX.Multiply( DistanceFromStart );
PX.Add( StartX );

PY.Assign( EndY );
PY.Subtract( StartY );
PY.Multiply( DistanceFromStart );
PY.Add( StartY );

// create a perpendicular point... we do this by calculating a normal and
multiplieing
// this with the distance from the line and adding this to px and py ;)

// calculate normilized normal for line
calc_normal( StartX,StartY, EndX,EndY, NX,NY );

// offset the point along the normal with distance from line
// PX := PX + NX * DistanceFromLine;
// PY := PY + NY * DistanceFromLine;

NX.Multiply( DistanceFromLine );
NY.Multiply( DistanceFromLine );

PX.Add( NX );
PY.Add( NY );

NX.Destroy;
NY.Destroy;
end;

procedure org_calc_perpendicular( var x, y : double );
var
temp_x, temp_y : double;
begin
temp_x := -y;
temp_y := x;
x := temp_x;
y := temp_y;
end;

procedure org_calc_normalize( var x, y : double );
var
temp_length : double;
begin
temp_length := sqrt( (x*x) + (y*y) );
if temp_length<>0 then
begin
x := x / temp_length;
y := y / temp_length;
end;
end;

procedure org_calc_normal( _x1,_y1, _x2,_y2 : double; var _normal_x,
_normal_y : double );
begin
_normal_x := _x2-_x1;
_normal_y := _y2-_y1;

// calculate perpendicular "vector".
org_calc_perpendicular( _normal_x, _normal_y );

// divide vector by it's own length there by creating a "normal" :)
org_calc_normalize( _normal_x, _normal_y );
end;

// helper procedure to generate points on or near line segments
procedure org_GeneratePointForLineSegment( StartX, StartY,
EndX, EndY : double;
DistanceFromStart : double; // distance from start
// < 0 is before start
// 0 is on start
// 0.5 is on center
// 1 is one end
// > 1 is after end
DistanceFromLine : double; // closest distance from line

// follows the normal
// positive is left side of line
// zero is on line
// negative is right side of line
var px, py : double );
var
NX, NY : double;
begin
// theory of method1 ;)
PX := StartX + DistanceFromStart * (EndX - StartX);
PY := StartY + DistanceFromStart * (EndY - StartY);

// create a perpendicular point... we do this by calculating a normal and
multiplieing
// this with the distance from the line and adding this to px and py ;)

// calculate normilized normal for line
org_calc_normal( StartX,StartY, EndX,EndY, NX,NY );

// offset the point along the normal with distance from line
PX := PX + NX * DistanceFromLine;
PY := PY + NY * DistanceFromLine;
end;

var
vHighPerformanceCounterFrequency : int64;
vCounter1 : int64;
vCounter2 : int64;

vStartX,
vStartY,
vEndX,
vEndY,
vPointX,
vPointY,
vDistanceFromStart,
vDistanceFromLine : Trational;

vOrgStartX,
vOrgStartY,
vOrgEndX,
vOrgEndY,
vOrgPointX,
vOrgPointY,
vOrgDistanceFromStart,
vOrgDistanceFromLine : double;

vMeasurementsEnabled : boolean;

vRationalTicks : int64;
vRationalTime : double;

vFloatingPointTicks : int64;
vFloatingPointTime : double;

begin
writeln('program started');

vMeasurementsEnabled := true;;

vRationalTicks := 0;
vRationalTime := 0;

vFloatingPointTicks := 0;
vFloatingPointTime := 0;

if vMeasurementsEnabled then
begin
if not QueryPerformanceFrequency( vHighPerformanceCounterFrequency ) then
begin
writeln('High performance counter not present, speed measurements
disabled');
vMeasurementsEnabled := false;
end;
end;

vStartX := Trational.Create;
vStartY := Trational.Create;
vEndX := Trational.Create;
vEndY := Trational.Create;
vPointX := Trational.Create;
vPointY := Trational.Create;
vDistanceFromStart := Trational.Create;
vDistanceFromLine := Trational.Create;

vStartX.Assign( 0, 1 );
vStartY.Assign( 0, 1 );

vEndX.Assign( 512, 1 );
vEndY.Assign( 997, 1 );

vDistanceFromStart.Assign( 3, 10 );
vDistanceFromLine.Assign( 0, 1 );

if vMeasurementsEnabled then
begin
QueryPerformanceCounter( vCounter1 );
end;

GeneratePointForLineSegment(
vStartX, vStartY,
vEndX, vEndY,
vDistanceFromStart,
vDistanceFromLine,
vPointX, vPointY );

if vMeasurementsEnabled then
begin
QueryPerformanceCounter( vCounter2 );

vRationalTicks := vCounter2-vCounter1;
vRationalTime := ( (vCounter2-vCounter1)/
vHighPerformanceCounterFrequency );
end;

writeln;
writeln( '(Rational) PointX in rational notation: ', vPointX.AsString );
writeln( '(Rational) PointY in rational notation: ', vPointY.AsString );

writeln;
writeln( '(Rational) PointX in real notation: ', vPointX.AsReal : 16 :
16 );
writeln( '(Rational) PointY in real notation: ', vPointY.AsReal : 16 :
16 );

if vMeasurementsEnabled then
begin
writeln('(Rational) ticks: ', vRationalTicks );
writeln('(Rational) time (in seconds): ', vRationalTime : 16 : 16 );
end;

vOrgStartX := 0.0;
vOrgStartY := 0.0;
vOrgEndX := 512.0;
vOrgEndY := 997.0;
vOrgDistanceFromStart := 0.3;
vOrgDistanceFromLine := 0.0;

if vMeasurementsEnabled then
begin
QueryPerformanceCounter( vCounter1 );
end;

org_GeneratePointForLineSegment(
vOrgStartX, vOrgStartY,
vOrgEndX, vOrgEndY,
vOrgDistanceFromStart,
vOrgDistanceFromLine,
vOrgPointX, vOrgPointY );

if vMeasurementsEnabled then
begin
QueryPerformanceCounter( vCounter2 );

vFloatingPointTicks := vCounter2-vCounter1;
vFloatingPointTime := ( (vCounter2-vCounter1)/
vHighPerformanceCounterFrequency );
end;

writeln;
writeln( '(Floating point) PointX: ', vOrgPointX : 16 : 16 );
writeln( '(Floating point) PointY: ', vOrgPointY : 16 : 16 );

if vMeasurementsEnabled then
begin
writeln('(Floating point) ticks: ', vFloatingPointTicks );
writeln('(Floating point) time (in seconds): ', vFloatingPointTime : 16 :
16 );
end;

vStartX.Destroy;
vStartY.Destroy;
vEndX.Destroy;
vEndY.Destroy;
vPointX.Destroy;
vPointY.Destroy;
vDistanceFromStart.Destroy;
vDistanceFromLine.Destroy;

writeln;
writeln('press enter to exit');
readln;

writeln('program finished');
end.

Bye,
Skybuck.
Nov 15 '05 #36

P: n/a

"Randy" <jo*@burgershack.com> wrote in message
news:dg**********@joe.rice.edu...
Skybuck Flying wrote:
"Hans-Bernhard Broeker" <br*****@physik.rwth-aachen.de> wrote in message
news:3p************@news.dfncis.de...
In comp.graphics.algorithms Skybuck Flying <no****@hotmail.com> wrote:

[Massive crossposting reduced to single group. Please, Skybuck, try
to follow netiquette a bit better in the future. If you feel a need
to X-post, at be a bit more selective to where --- this had *no*
business being posted to comp.arch or sci.electronics. And *always*
select a single group for the follow-ups.]


No, I like to have input from all those newsgroups since this is about
optimizations and those newsgroup are related to software and hardware
optimizations. C/Delphi/Asm for software/cpu optimizations and
comp.arch/asm/eletronics about hardware optimizations and algorithms
possibly about algorithm optimization. Those newsgroups are likely to
contain people who know something about this.


I'm with Hans on this one. Skybuck, by your logic you should have
posted to ALL newsgroups, since there's probably someone in every
newsgroup who will be somehow "related to optimization". :-(

Comp.graphics.algorithms was the right place. Comp.lang.c is hard to
justify, since your question has nothing to do with C. But there's no
way to justify posting to sci.electronics.design, comp.arch, or
alt.comp.lang.borland-delphi. Nor comp.lang.asm, unless you want a
response written in assembly language.

That said, other appropriate forums might have included comp.programming
or *maybe* comp.theory.

...

You need to define the task a good deal more cleanly before an
implementation can reasonably be suggested. Missing information:

1) What's a line segment, in this case? It may come as a surprise to
you, but ideas about the details of this actually differ. So: how do
you actually specify this line segment?

A 2D line segment defined by (x1,y1) - (x2,y2)


Integer or floating point? In either case, how big are your points and
how fat is the line? Is it acceptable to return a false positive or
negative? If so, how often? If it's not acceptable, then you have a
problem. On a any machine that lacks infinite numerical precision, you
can never know whether an infinitely small point lies on an infinitely
thin line segment.

In the end, you're going to have to do something like this:

"Draw two lines, one from each point in the line segment to the
candidate point. If these lines' slopes are equal but differ in sign
(e.g. X and -X), then your point lies on the segment."

Of course, what does it mean for slope values to be "equal" on a digital
computer? Unless you can represent the problem entirely symbolically,
"equal" will be necessarily limited by the representation of the data
and the accuracy of the math. In the real world, epsilon is a necessary
evil.
2) When exactly do you want the test to return "true"? I.e.: how
close to the ideal segment does a point have to be to still be
considered "on" it?

As precise as possible. (Floating point precision, no rounding)


This is an oxymoron. Floating point is inescapably about rounding.
Digital representations of floating point are approximations, as are the
math operations that manipulate them. The programmer has to manage
numerical imprecision on computers or live with wrong answers.

For some perspective, read David Goldberg's "What Every Computer
Scientist Should Know About Floating-Point Arithmetic":

http://docs.sun.com/source/806-3568/ncg_goldberg.html


Yes I am sure you all right about this... But instead of wasting my time
reading through this highly theoretical shit I rather find practical
solutions to this potential floating point rounding problem. And using
epsilon's is probably not it since it still leaves room for errors etc.

So far I have found one nice solution... Trational class from delphi
fundamentals.

Another solution might be a fixed point library which might be faster =D

Bye,
Skybuck.
Nov 15 '05 #37

P: n/a

"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Nicholas Sherlock wrote:
Skybuck Flying wrote:
I needed a method to determine if a point was on a line segment in 2D.
So I
googled for some help and so far I have evaluated two methods.

>
function PointOnLine2D(x1,y1, x2,y2, x3,y3:double):boolean;
begin
result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.0000001;
end;


Sorry, I missed some of this code from my app:

function PointOnLine2D(x1, y1, x2, y2, x3, y3: double): boolean;
begin
result :=
(((x3 >= x1) and (x3 <= x2)) or ((x3 >= x2) and (x3 <= x1))) and
(((y3 >= y1) and (y3 <= y2)) or ((y3 >= y2) and (y3 <= y1))) and
(abs(((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1))) < 0.00001);
end;


Ok I first I didn't understand this code, but not I get it.

The other C example was like this example.

The C example compares the first multiplication result with the second
multiplication result.

If they are the same then it's true else false.

Because floating point rounding errors could occur the C example could fail.
Which was demonstrated by my method 2 which was the ported C version.

So your/the fastgeo version solves this rounding error by using an
epsilon...

I still think it's a shabby solution since it could horribly fail if the
points are small enough.

Maybe something like:
0.00000003245
0.000000003245
0.000000003456
0.000000000123

etc...

OUCH

;)

Bye,
Skybuck =D
Nov 15 '05 #38

P: n/a

"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:gv********************@comcast.com...
Skybuck Flying wrote:

(snip)
That certainly won't do in this case etc.
It should be as exact/accurate as possible.


You really need to tell us what "this case" is.


I don't see the relevance of the case.

The math is clear and simple and is case independant =D

You expect everyone to help you, but don't seem
interested in helping much yourself.


No on the contrary ;) See above lol.

Bye,
Skybuck.
Nov 15 '05 #39

P: n/a

"nobody" <no****@nowhere.com> wrote in message
news:43********************************@4ax.com...
"Skybuck Flying" <no****@hotmail.com> wrote:
Using epsilon's and stuff like should probably be prevented as much as
possible since those are common forms of errors etc and limitations
etc...
You have a gross and dangerous misunderstanding of principles of doing
numerical processing with the computer. I suggest you study the
fundmentals carefully first before wasting your time writing what will
undoubtedly be some ridicolously naive and rather useless code.


I don't need to understand fully how floating point works. Sometimes I
forget about rounding errors, so this thread was a nice refreshment. However
in the back of my mind I know floating points are inaccurate, that is why I
requested the methods to be "as accurate as possible". "as possible" meaning
as possibly as *you* can get it without being to fricking slow lol =D

The epsilon method is flawed in my mind since it assumes the points are
quite large but in fact could be quite small thereby completely failing.

Maybe you have a gross and dangerous misunderstanding of using a smudge
factor/epsilon ;)

I also don't like to study flawed/worthless shit like floating point format
etc... I rather spent my time finding solutions to this shit lol =D. So far
I have found a better alternative which is a rational number library... it
still uses floating points inside which are used to represent the rational
numbers but since it uses fractions it should be a lot more accurate than
just using floating points directly ;)

Another solution might be a fixed point library =D

Bye,
Skybuck.
Nov 15 '05 #40

P: n/a
Skybuck Flying wrote:
"nobody" <no****@nowhere.com> wrote in message
news:43********************************@4ax.com...
"Skybuck Flying" <no****@hotmail.com> wrote:
Using epsilon's and stuff like should probably be prevented as much as
possible since those are common forms of errors etc and limitations

You have a gross and dangerous misunderstanding of principles of doing
numerical processing with the computer


The epsilon method is flawed in my mind since it assumes the points are
quite large but in fact could be quite small thereby completely failing.


Moron.

Nicholas Sherlock
Nov 15 '05 #41

P: n/a

<we******@gmail.com> wrote in message
news:11**********************@g49g2000cwa.googlegr oups.com...
Skybuck Flying wrote:
"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote:
Skybuck Flying wrote:
> It should be as exact/accurate as possible.

You really need to tell us what "this case" is.
You ll see quite soon the test program is done ;)


Ok, you are pissing people off with statements like this.


Lol, don't get carried away ;) Relax, take a break, it's never been my
intention to piss anybody off.

Though I can understand the frustration some people might be having because
the requirements seem
to be impossible to unite with floating points because of floating point
rounding errors.

Those people should really try and find some nice rational number or fixed
point library and poof your frustration will turn into pure relaxation lol
=D
Here's the thing, the "epsilon thing" is far and away the most
practical way of doing this; but just as a matter of correctness, you
actually should not choose a constant epsilon. The correct epsilon
will depend on the magnitude of the coordinates for the original
segment points. For example it would be easy to construct an example
where the best correct epsilon is a million.
And how would such a dynamic epsilon be determined or constructed ? ;)

Besides from that... the code will have to be modified to isolate the
rounding error... which might be quite difficult sometimes (?) or maybe
not.. but just extra boring work which also slows down the code ;)
Ok, perhaps a better way to ask you for more information, is this.
From what sort of *SOURCE* are your input points comming from? Are they just chosen literally at random? (Using say, the Mersenne Twister
random number generator, which includes floating point random.) Or
(more likely) are they actually constructed to be *ON* the line with
some likelihood, but for some reason you loose track of this fact, and
wish to recalculate it?

This is important because, the *WAY* you construct the point (say, but
being given the x intercept, then computing the y that is supposed to
correspond to it) may actually give an exact matching algorithm without
the need for any epsilons. One could, for example, just try to
reproduce the procedure for finding the point, from one of the
coordinates, and see if the other matches exactly.


This is a good point you make. The source that constructs the points could
introduce rounding errors...

Then the method which checks the points could also introduce the same
rounding errors thereby cancelling the rounding errors against each other
and still returning success while in reality the point is slighty
wrong/besides the line... other methods might thus correctly detect this.
The verification program would then verify this correct method as being
flawed which is not the case. It's exactly the opposite =D

Thus the verification program needs to be improved to use more accurate
point construction.. for example by using rational numbers or fixed point
numbers ;)
The problem with this is that starting with an x and computing a y, or
the other way around will not necessarily yield the same results.
Further more if you compute the point as (alpha * p_0 + (1-alpha)* p_1)
(which might be *WAY* harder to match -- intuitively it seems so to me,
but I don't have a 100% clear idea), you will again get different LSB
round offs.
Yeah, I understand what you mean =D
So I hope you understand that these accuracy issues are pretty integral
to what it is that you want to do. Without more information from you,
I don't believe that anyone can suggest anything else more with any
expectation of it necessarily working better.
At this point I disagree... math is math... it shouldn't matter what the
source of the information/variables are.

You all have been helpfull at pointing out the inaccuracy problems =D

Bye,
Skybuck.

--
Paul Hsieh
http://www.pobox.com/~qed/
http://bstring.sf.net/

Nov 15 '05 #42

P: n/a
Skybuck Flying wrote:

(snip)
I don't need to understand fully how floating point works. Sometimes I
forget about rounding errors, so this thread was a nice refreshment. However
in the back of my mind I know floating points are inaccurate, that is why I
requested the methods to be "as accurate as possible". "as possible" meaning
as possibly as *you* can get it without being to fricking slow lol =D The epsilon method is flawed in my mind since it assumes the points are
quite large but in fact could be quite small thereby completely failing.


It is not only floating point, but that usually makes it worse.

My example of (0,0) and (997,512) is fixed point. There are no points
on the line segment between them in fixed point.

-- glen

Nov 15 '05 #43

P: n/a

"David Hopwood" <da******************@blueyonder.co.uk> wrote in message
news:VP******************@fe1.news.blueyonder.co.u k...
[some off-topic newsgroups trimmed]

Skybuck Flying wrote:
"Maarten Wiltink" <ma*****@kittensandcats.net> wrote:
"Skybuck Flying" <no****@hotmail.com> wrote:
"Nicholas Sherlock" <n_********@hotmail.com> wrote:
>Skybuck Flying wrote:

>>Though the < 0.0000001 is new which I dont quite understand ;)
>
>Floating point math is always inaccurate as most numbers cannot be
>exactly represented. This bit basically takes care of that by adding
>a fudge factor.

Yes I figured that much. So this would mean the function returns true
if abs(blabla) is near zero ?

Why use 0.00001 why not 0.0000001 or 0.00000000000001 ?

Personally I don't like epsilons like this... it leaves room for
error/malfunction ?

Some real life logic: once in a graphics rendering lab assignment, we
were instructed to approximate Bezier curves. This is an iterative
process; the stop condition was for the error to become less then
half a pixel. For visualisation on a raster display, that is exactly
what's required.
That certainly won't do in this case etc.

It should be as exact/accurate as possible.


Even trying to detect whether a point is on a line using floating point
arithmetic almost certainly indicates a specification error. It's

impossible to do that exactly.
Not quite... it depends on the input/data and how the floating points are
used.

The rational number library uses floating points but since it uses them only
to represent fractions like
2/3 which could also be represented with integers... but it probably uses
floating points for more convenience and internal conversions etc.. it will
be a lot more accurate... but at some point it might again get inaccurate ;)

Bye,
Skybuck.

--
David Hopwood <da******************@blueyonder.co.uk>

Nov 15 '05 #44

P: n/a

"Ketil Malde" <ke********@ii.uib.no> wrote in message
news:eg************@polarvier.ii.uib.no...
"Skybuck Flying" <no****@hotmail.com> writes:
It should be as exact/accurate as possible.
Then use rational numbers?


Excellent idea... done that ;)

Bye,
Skybuck.

-k
--
If I haven't seen further, it is by standing in the footprints of giants

Nov 15 '05 #45

P: n/a
"Skybuck Flying" <no****@hotmail.com> writes:
Can you still find a point which would in theory be on the line segment but
would still be mis-calculated by the rational number version ? ;)


One that overflows your rational number implementation, perhaps? (If
it uses bignums, you may have to exhaust the heap).

-k
--
If I haven't seen further, it is by standing in the footprints of giants
Nov 15 '05 #46

P: n/a
Maarten Wiltink wrote:
Some real life logic: once in a graphics rendering lab assignment, we
were instructed to approximate Bezier curves. This is an iterative
process; the stop condition was for the error to become less then
half a pixel. For visualisation on a raster display, that is exactly
what's required.


And then the formula is

result:=((x2 - x1) * (y3 - y1)) - ((x3 - x1) * (y2 - y1)) < 0.25;

or < width*width*0.25, if you have a defined line width ("is on a line" begs
the question "line width", anyway).

--
Bernd Paysan
"If you want it done right, you have to do it yourself"
http://www.jwdt.com/~paysan/
Nov 15 '05 #47

P: n/a
"Skybuck Flying" <no****@hotmail.com> wrote:
The epsilon method is flawed in my mind since it assumes the points are
quite large but in fact could be quite small thereby completely failing.


And that's why you need to study some fundementals. Start with
absolute versus relative tolerance (GIYF) and decide which (or some
other measure) is better in your case.

Nov 15 '05 #48

P: n/a

"glen herrmannsfeldt" <ga*@ugcs.caltech.edu> wrote in message
news:4r********************@comcast.com...
Skybuck Flying wrote:

(snip)
I don't need to understand fully how floating point works. Sometimes I
forget about rounding errors, so this thread was a nice refreshment. However in the back of my mind I know floating points are inaccurate, that is why I requested the methods to be "as accurate as possible". "as possible" meaning as possibly as *you* can get it without being to fricking slow lol =D

The epsilon method is flawed in my mind since it assumes the points are
quite large but in fact could be quite small thereby completely failing.


It is not only floating point, but that usually makes it worse.

My example of (0,0) and (997,512) is fixed point. There are no points
on the line segment between them in fixed point.


What do you mean with fixed point ? I guess you are using some kind of fixed
point implementation, which I dont have ofcourse ;) ?

At least in the rational number implementation I have already proven that
there is a point on the line segment.

(Rational) PointX in real notation: 153.6000000000000000
(Rational) PointY in real notation: 299.1000000000000000

Quite easily actually.

It's just that the floating point method can't calculate it accurately.

Bye,
Skybuck.
Nov 15 '05 #49

P: n/a

"Nicholas Sherlock" <n_********@hotmail.com> wrote in message
news:dg**********@lust.ihug.co.nz...
Skybuck Flying wrote:
"nobody" <no****@nowhere.com> wrote in message
news:43********************************@4ax.com...
"Skybuck Flying" <no****@hotmail.com> wrote:
Using epsilon's and stuff like should probably be prevented as much as
possible since those are common forms of errors etc and limitations
You have a gross and dangerous misunderstanding of principles of doing
numerical processing with the computer


The epsilon method is flawed in my mind since it assumes the points are
quite large but in fact could be quite small thereby completely failing.


Moron.


Lol, Am I a moron because I am right and you were lazy and copied a flawed
method from some library ? ;)

Bye,
Skybuck =D

Nicholas Sherlock

Nov 15 '05 #50

65 Replies

This discussion thread is closed

Replies have been disabled for this discussion.