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

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

--
	Sander

+++ Out of cheese error +++

(snip)

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

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_length0 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_length0 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.

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

I don't see the relevance of the case.

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

No on the contrary ;) See above lol.

Bye, Skybuck.

etc...

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.

Excellent idea... done that ;)

Bye, Skybuck.

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

are

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.

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.

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/

we

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/

Moron.

Nicholas Sherlock

[...]

Yes, you do have a fixed point implementation. It has a constant delta value of 1. It's called integer artithmetic.

In a fixed point numbering system, numbers are represented as integer multiples of some base value, which is sometimes called the "delta". For example, a fixed-point representation with a delta of 0.01 might be useful for representing money, with dollar amounts being represented as a whole number of cents.

Some languages have direct support for such types. In others, you can always simulate it by using integer arithmetic and keeping track of the delta yourself.

Well, not quite; you've swapped the coordinates. You mean (299.1,153.6), not (153.6,299.1).

If your end points are (0,0) and (997,512), and you're using integer arithmetic (with a delta of 1), there are no representable points on the line segment excluding the end points. (299.1,153.6) is not representable.

Equivalently, if you're using a fixed-point delta of 0.01, and your end points are (0,0) and (9.97,5.12) there are still no representable points on the line segment excluding the end points. (2.991,1.536) is not representable.

Now if you're willing to have your is_this_point_on_this_line_segment() function always return false for some line segments, that's ok. If you expect it to return true for some significant number of points, either you can use a numeric representation that supports arbitrary precision (integer, fixed-point, and floating-point won't do; rational numbers with unlimited range on the numerator and denominator probably will), or you can define some concept of "close enough".

What might turn out to be more useful is a function that, given the end points of a line segment an arbitrary point, tells you the distance between the line segment and the point. Defining the "distance" when the point is beyond the end points is left as an exercise.

--
Keith Thompson (The_Other_Keith) kst-u@mib.org  
San Diego Supercomputer Center               
We must do something.  This is something.  Therefore, we must do this.

However

why I

meaning

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.

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

Bye, Skybuck =D

fixed

Ok, fixed point implementations therefore have the same rounding problems as floating point implementations ;) =D Since both work with a "point" =D

that

No, the original poster has swapped the coordinates ;)

In another post the original poster mentioned this line segment:

" 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. "

Yes "close enough" would be the epsilon method... However using multiple floating point operations will eventually create numerical drift so at some point even the epsilon method will probably fail because the rounding errors get to large !? ;)

Neh, this solution faces the same rounding problems.

Bye, Skybuck.

Read Knuth. Google brounding. Get educated. It's all quite simple. Generalize to n dimensions, two is boring. You have not yet defined the problem you wish to solve. See if you can do that. Bet you can't. "Is this point on this line segment" is not a problem definition. If you are talking mathematics the solution is simple and known, if you provide the rest of the problem definition. If you are talking computers and programming, you didn't pose any interesting problem. The answer is always "yes and no" without further definition of what you mean.

Go for it!

--

  ... Hank

http://home.earthlink.net/~horedson
http://home.earthlink.net/~w0rli

Hmmm? I google'd "brounding" and it seems to have something to do with electrical grounding (but appears to have an inflection I didn't quite catch... a brazed on grounding point, maybe??)

I don't see what brounding has to do with the topic at hand?

--
These .signatures are sold by volume, and not by weight.