8

I am trying to verify when 3 points (double) are collinear in 2-D. I have found different Pascal functions that return true if this is verified; those functions use integer to specify X and Y coordinates. I need a more precise calculation at least to the first 3 digits of the decimal part of X and Y expressed as double type. Who can help me with this?

I found this function:

function Collinear(x1, y1, x2, y2, x3, y3: Double): Boolean;
begin
  Result := (((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) = 0);
end;

But I guess the calculation would never be 0. Should I use something like that?

function Collinear(x1, y1, x2, y2, x3, y3: Double): Boolean;
begin
  Result := (((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) < 0.01);
end;
Rob Kennedy
  • 161,384
  • 21
  • 275
  • 467
marcostT
  • 573
  • 1
  • 10
  • 20
  • 1
    analytic geometry can help, 2 points constitutes a straight line, so there be system of 3 linear equations, and if it is consistent - points are collinear. Re: `0.01` - have to precalculate "machine epsilon" due to discrete arithmetics. – Free Consulting Feb 15 '11 at 23:09
  • @Worm machine epsilon isn't particularly relevant to this. You also need to account for the scale of the input values and "how" co-linear you want the points to be. If you use `eps` as a tolerance then you'll most likely never find any points to be co-linear unless two of the points are identical! – David Heffernan Feb 16 '11 at 06:30

2 Answers2

9

The scalar product equation you calculate is zero if and only if the three points are co-linear. However, on a finite precision machine you don't want to test for equality to zero but instead you test for zero up to some small tolerance.

Since the equation can be negative as well as positive your test isn't going to work. It will return false positives when the equation evaluates to a large negative value. Thus you need to test that the absolute value is small:

function Collinear(const x1, y1, x2, y2, x3, y3: Double): Boolean;
const
  tolerance = 0.01;//need a rationale for this magic number
begin
  Result := abs((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) < tolerance;
end;

Exactly how to choose tolerance depends on information you haven't provided. Where do the values come from? Are they dimensional?

David Heffernan
  • 601,492
  • 42
  • 1,072
  • 1,490
  • 3
    Performance tip: if you add `inline;` after `Boolean;`, this function becomes about 10x faster (on this machine). – Wouter van Nifterick Feb 15 '11 at 22:46
  • Agreed, this is a great function to inline. – Andreas Rejbrand Feb 15 '11 at 22:51
  • My analytic geometry is rusty, but how can you do dot product with 3 vectors? – arthurprs Feb 16 '11 at 02:25
  • Nice, It's not dot (scalar) product, but a simplification of line-point distance which should be 0, http://mathworld.wolfram.com/Collinear.html – arthurprs Feb 16 '11 at 02:33
  • @arthurps It most certainly is scalar product (or can be thought of that way). If v1=(x1,y1), v2=(x2,y2) and v3=(x3,y3) are the three points then you form the direction vectors d1=v2-v1 and d2=v3-v1. Then find the vector perpendicular to d2, d2_perp which can be written down by inspection trivially. Finally the co-linearity condition amounts to d1.d2_perp=0. – David Heffernan Feb 16 '11 at 06:29
  • Hi David. Yes got it. But it seems that the function is true only when the line is vertical or horizontal. Even if I increase the tolerance. - Exactly how to choose tolerance depends on information you haven't provided - What information you need? – marcostT Feb 16 '11 at 09:18
  • @MarcosT What is the scale of the values, e.g. typical values. Do you get values like 0.001, 1, 1000, 1e20 etc. – David Heffernan Feb 16 '11 at 09:37
  • @Wouter I don't believe a 10x speedup for inlining that: 6 QWORD reads, 4 adds, 2 mults and a compare does not come to 1/10th of the cost of a function call in my experience. I think your benchmark is wrong. – David Heffernan Feb 16 '11 at 09:38
3

David's code will work, but you should get the tolerance as a function of the parameters, like so:

function Collinear(const x1, y1, x2, y2, x3, y3: Double): Boolean; inline;   
var  
  tolerance: double;  
begin  
  tolerance := abs(max(x1,x2,x3,y1,y2,y3)) * 0.000001;  
  Result := abs((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) < tolerance;  
end;

If you don't do that and use a constant instead you can run into weird errors with large values for x1..y3.

Johan
  • 74,508
  • 24
  • 191
  • 319
  • The tolerance calculation will also need an abs call Consider the case where all coordinates are negative but equal. The left side of the "... < tolerance" expression will be zero (because everything is equal), but tolerance will be negative. Therefore `tolerance := abs(max(x1,x2,x3,y1,y2,y3)) * 0.000001;` – OldPeculier Jun 24 '13 at 22:19