Saturday, January 25, 2014

Accurate point in triangle test

Many resources on the web deal with the 2D point in triangle test. This post will summarize the most famous methods to solve it and will show why in some cases they are not accurate and can lead to errors. As a conclusion, we will expose a new algorithm.

At first, we need to remember the problem. We are in 2D. We have a triangle T defined by 3 points p1(x1, y1), p2(x2, y2), p3(x3, y3) and a single point p(x, y). Does this single point p lies inside the triangle T ?



1st method : barycentric coordinate system

Barycentric coordinate allows to express new p coordinates as a linear combination of p1, p2, p3. More precisely, it defines 3 scalars a, b, c such that :

x = a * x1 + b * x2  + c * x3
y = a * y1 + b * y2 + c * y3
a + b + c = 1

The way to compute a, b, c is not difficult :

a = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3)) / ((y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3))
b = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3)) / ((y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3))
c = 1 - a - b

Then we just need to apply the interesting following property :

p lies in T if and only if 0 <= a <= 1 and 0 <= b <= 1 and 0 <= c <= 1

Code sample :

function pointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var denominator:Number = ((y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3));
var a:Number = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3)) / denominator;
var b:Number = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3)) / denominator;
var c:Number = 1 - a - b;

return 0 <= a && a <= 1 && 0 <= b && b <= 1 && 0 <= c && c <= 1;
}

2nd method : parametric equations system

Here the idea is to consider the parametric expressions of the 2 edges [p1, p2] and [p1, p3] in T :

x(t1) = t1*(x2 - x1)
y(t1) = t1*(y2 - y1)

x(t2) = t2*(x3 - x1)
y(t2) = t2*(y3 - y1)

Then express p(x, y) as a linear combination of them :

x = x1 + x(t1) + x(t2)
y = y1 + y(t1) + y(t2)

Solving the system, it gives to us :

t1 = (x*(y3 - y1) + y*(x1 - x3) - x1*y3 + y1*x3) / (x1*(y2 - y3) + y1*(x3 - x2) + x2*y3 - y2*x3)
t2 = (x*(y2 - y1) + y*(x1 - x2) - x1*y2 + y1*x2) / -(x1*(y2 - y3) + y1*(x3 - x2) + x2*y3 - y2*x3)

Finally, we just need to apply the interesting following property :

p lies in T if and only if 0 <= t1 <= 1 and 0 <= t2 <= 1 and t1 + t2 <= 1

Code sample :

function pointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var denominator:Number = (x1*(y2 - y3) + y1*(x3 - x2) + x2*y3 - y2*x3);
var t1:Number = (x*(y3 - y1) + y*(x1 - x3) - x1*y3 + y1*x3) / denominator;
var t2:Number = (x*(y2 - y1) + y*(x1 - x2) - x1*y2 + y1*x2) / -denominator;
var s:Number = t1 + t2;

return 0 <= t1 && t1 <= 1 && 0 <= t2 && t2 <= 1 && s <= 1;
}

3rd method : check sides with dot product

Maybe the most famous method, based on dot product. We assume that p1, p2, p3 are ordered in counterclockwise. Then we can check if p lies at left of the 3 oriented edges [p1, p2], [p2, p3] and [p3, p1].

For that, first we need to consider  the 3 vectors v1, v2 and v3 that are respectively left-orthogonal to [p1, p2][p2, p3] and [p3, p1] :

v1 = <y2 - y1, -x2 + x1>
v2 = <y3 - y2, -x3 + x2>
v3 = <y1 - y3, -x1 + x3>

Then we get the 3 following vectors :

v1' = <x - x1, y - y1>
v2' = <x - x2, y - y2>
v3' = <x - x3, y - y3>

At last, we compute the 3 dot products :
dot1 = v1 . v1' = (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1)
dot2 = v1 . v2' = (y3 - y2)*(x - x2) + (-x3 + x2)*(y - y2)
dot3 = v3 . v3' = (y1 - y3)*(x - x3) + (-x1 + x3)*(y - y3)

Finally, we can apply the interesting property :

p lies in T if and only if 0 <= dot1 and 0 <= dot2 and 0 <= dot3

Code sample :

function side(x1, y1, x2, y2, x, y:Number):Number
{
return (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1);
}

function pointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean

{
var checkSide1:Boolean = side(x1, y1, x2, y2, x, y) >= 0;
var checkSide2:Boolean = side(x2, y2, x3, y3, x, y) >= 0;
var checkSide3:Boolean = side(x3, y3, x1, y1, x, y) >= 0;
return checkSide1 && checkSide2 && checkSide3;
}


These 3 methods are quite good to solve the point in triangle test. Purely mathematically speaking, they must validate any point inside the triangle and even those lying exactly on the boundary (on any edge).


Accuracy problems

Despite the strong mathematical background of our methods, in some cases they can lead to a lack of accuracy because the floating-point number system have limited size and most of the time it deals with approximations. The problem occurs sometimes when the point p should exactly on one triangle's edge ; the approximations lead to fail the test.

Example 1*:
We consider the triangle T defined by 3 points p1(x1, y1), p2(x2, y2), p3(x3, y3) with values :

x1 = 1/10
y1 = 1/9
x2 = 100/8
y2 = 100/3
x3 = 100/4
y3 = 100/9

and a single point p(x, y) lying exactly on the segment [p1, p2] :
x = x1 + (3/7)*(x2 - x1)
y = y1 + (3/7)*(y2 - y1)



If we apply the barycentric method, we get the 3 following values for a, b, c :

a : 0.5714285714285715
b : 0.42857142857142855
c : -5.551115123125783e-17

Because c < 0, the test fails to validate the point inside the triangle. In many applications, this situation is not really a problem because a lack of accuracy is not a tragedy.

But in some situations, it can be really annoying.

Example 2*:
We consider 2 triangles and T' defined respectively by the points p1(x1, y1), p2(x2, y2), p3(x3, y3) and p1'(x1', y1'), p2'(x2', y2'), p3'(x3', y3'). Values are :

x1 = 1/10
y1 = 1/9
x2 = 100/8
y2 = 100/3
x3 = 100/4
y3 = 100/9

x1' = x1
y1' = y1
x2' = x2
y2' = y2
x3' = -100/8
y3' = 100/6

and a single point p(x, y) lying exactly on the segment [p1, p2] :
x = x1 + (3/7)*(x2 - x1)
y = y1 + (3/7)*(y2 - y1)



The situation is quite simple : 2 non-overlapping triangles sharing one edge and one single point lying on this edge. We know from the previous example that the barycentric method fails to validate p inside the triangle T :

a : 0.5714285714285715
b : 0.42857142857142855
c : -5.551115123125783e-17

But more surprisely, the method fails again when applied to the triangle T' :

a : 0.5714285714285715
b : 0.4285714285714285
c : -1.1102230246251565e-16

Mathematically speaking, the point p belongs to both triangles. In practice, we could tolerate that a lack of accuracy leads to validate only one triangle belonging. But here we face a complete invalidation and the point is detected as outside both triangles.

We could suppose that the barycentric-based method is intrinsically inaccurate and one other method will lead to satisfiable results. But in reality the problem remains.

Example 3*:
We consider the oriented edge E defined by 2 points p1(x1, y1) and p2(x2, y2) with values :

x1 = 1/10
y1 = 1/9
x2 = 100/8
y2 = 100/3

and a single point p(x, y) lying exactly on edge E :
x = x1 + (3/7)*(x2 - x1)
y = y1 + (3/7)*(y2 - y1)

If we apply the check side method with E and p, we get the following dot product :
dot = (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1)
dot = -2.842170943040401e-14

Then if we apply the same method with E' the reversed edge of E, we get the same result :
dot = (y1 - y2)*(x - x2) + (-x1 + x2)*(y - y2)
dot = -2.842170943040401e-14

It means that the check side test gives us 2 contradictory results. The point looks on the right sides of both the edge E and his reversed ; that is of course mathematically impossible. We face the same situation as using barycentric method : if the edge E is shared by 2 non-overlapping triangles T and T', despite the fact that p lies mathematically on E, our method will lead us to tragically conclude that p is outside both triangles T and T'.


Any accurate solution ?

The answer is yes. We can write a complete algorithm leading to a safe point in triangle test by combining many familiar algorithms. The core of the method is to assume a real thickness value for the triangle's edges and vertices ; it contrasts with the original purely mathematical situation where the triangle's edges and vertices have an infinitesimal thickness. We will call it epsilon, because in practice we will keep it very small (near 0.001).



The steps of the algorithm are:

1: use a bounding box as a fast pre-validation test
Because our new algorithm will involve more computations, it is convenient to add a pre-test to reject very quickly most of the false cases.

The bounding box is simply the min/max of the x/y values among the 3 triangle's vertices, slightly inflated by the epsilon value.



Code sample :

const EPSILON:Number = 0.001;

function pointInTriangleBoundingBox(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var xMin:Number = Math.min(x1, Math.min(x2, x3)) - EPSILON;
var xMax:Number = Math.max(x1, Math.max(x2, x3)) + EPSILON;
var yMin:Number = Math.min(y1, Math.min(y2, y3)) - EPSILON;
var yMax:Number = Math.max(y1, Math.max(y2, y3)) + EPSILON;

if ( x < xMin || xMax < x || y < yMin || yMax < y )
return false;
else
return true;
}

2: use any method studied prevously (barycentric, parametric or dot product)
If the test is positive, we can trust it and stop the algorithm immediatly. But if the test is negative, maybe we face the situation with the point lying on one triangle's edge : then we need more investigations, involving methods using the epsilon value.

3: use the point to segment distance
For every edge of the triangle, we compute the shortest distance between the edge and the point to evaluate. If the distance is shorter that epsilon, we can validate definitely the test.

In detail, given 3 points  p, p1 and p2, a very tricky use of the dot product allows us to check efficiently the relative position of the orthogonal projection p' of p on the infinite line passing through p1 and p2. If the projection lies between p1 and p2 then we compute the distance p and p'. Otherwise, we compute the distance between p and the nearest among p1 and p2.


Finally, because our algorithm use distances for comparisons only, we will restrict our computations to square distances only (faster because we can omit the square root).

Code sample:

function distanceSquarePointToSegment(x1, y1, x2, y2, x, y:Number):Number
{
var p1_p2_squareLength:Number = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
var dotProduct:Number = ((x - x1)*(x2 - x1) + (y - y1)*(y2 - y1)) / p1_p2_squareLength;
if ( dotProduct < 0 )
{
return (x - x1)*(x - x1) + (y - y1)*(y - y1);
}
else if ( dotProduct <= 1 )
{
var p_p1_squareLength:Number = (x1 - x)*(x1 - x) + (y1 - y)*(y1 - y);
return p_p1_squareLength - dotProduct * dotProduct * p1_p2_squareLength;
}
else
{
return (x - x2)*(x - x2) + (y - y2)*(y - y2);
}
}


Final sample code

Here is the code illustrating the steps we described in the previous section.

const EPSILON:Number = 0.001;
const EPSILON_SQUARE:Number = EPSILON*EPSILON;

function side(x1, y1, x2, y2, x, y:Number):Number
{
return (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1);
}

function naivePointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var checkSide1:Boolean = side(x1, y1, x2, y2, x, y) >= 0;
var checkSide2:Boolean = side(x2, y2, x3, y3, x, y) >= 0;
var checkSide3:Boolean = side(x3, y3, x1, y1, x, y) >= 0;
return checkSide1 && checkSide2 && checkSide3;
}

function pointInTriangleBoundingBox(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var xMin:Number = Math.min(x1, Math.min(x2, x3)) - EPSILON;
var xMax:Number = Math.max(x1, Math.max(x2, x3)) + EPSILON;
var yMin:Number = Math.min(y1, Math.min(y2, y3)) - EPSILON;
var yMax:Number = Math.max(y1, Math.max(y2, y3)) + EPSILON;

if ( x < xMin || xMax < x || y < yMin || yMax < y )
return false;
else
return true;
}

function distanceSquarePointToSegment(x1, y1, x2, y2, x, y:Number):Number
{
var p1_p2_squareLength:Number = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
var dotProduct:Number = ((x - x1)*(x2 - x1) + (y - y1)*(y2 - y1)) / p1_p2_squareLength;
if ( dotProduct < 0 )
{
return (x - x1)*(x - x1) + (y - y1)*(y - y1);
}
else if ( dotProduct <= 1 )
{
var p_p1_squareLength:Number = (x1 - x)*(x1 - x) + (y1 - y)*(y1 - y);
return p_p1_squareLength - dotProduct * dotProduct * p1_p2_squareLength;
}
else
{
return (x - x2)*(x - x2) + (y - y2)*(y - y2);
}
}

function accuratePointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
if (! pointInTriangleBoundingBox(x1, y1, x2, y2, x3, y3, x, y))
return false;

if (naivePointInTriangle(x1, y1, x2, y2, x3, y3, x, y))
return true;

if (distanceSquarePointToSegment(x1, y1, x2, y2, x, y) <= EPSILON_SQUARE)
return true;
if (distanceSquarePointToSegment(x2, y2, x3, y3, x, y) <= EPSILON_SQUARE)
return true;
if (distanceSquarePointToSegment(x3, y3, x1, y1, x, y) <= EPSILON_SQUARE)
return true;

return false;
}


I hope you appreciated.


* all the tests are done in Actionscript 3, typing values with the native Number type (IEEE-754 double-precision floating-point number).


35 comments:

  1. Thank you! Your algorythm for checking if point is in triangle is very cool!

    ReplyDelete
  2. thanks for your insights! though im still looking for how to determine if a point is on a triangle in three dimension space

    ReplyDelete
  3. On the dot product interiority test, I don't think you have to assume that p1, p2, p3 are in counter-clockwise order. Just continue with the rest, and if the dot products /are all the same sign/ then you're on the same side of all, and thus inside. (You could not be on the /outside/ of them all is why this works.)

    ReplyDelete
  4. Nice article!

    I am trying out the second one (parametric equations system). I am wondering however: in the calculations the denominator could be zero. Could this happen in practise? What would this mean geometrically if this was the case?

    ReplyDelete
    Replies
    1. Usually, 0 value for denominator happens for degenerate cases.

      For example, when you try to find the intersection of 2 lines using a parametric system, you face a 0 value when lines are parallel.

      In our case our system try to find the p coordinates as a linear combination of [p1, p2] and {p1, p3]. So what are our degenerate cases ? For example, when p2 = p3 ! Meaning your triangle is collapsed as a flat line !

      Delete
    2. so first you have to test if the triangle is in fact a triangle?

      Delete
    3. If you are not sure that your triangle is safe, better is to check it before perfoming calculation. Check if any 2 points don't overlap. But problems can still occur with points being too close from each other !

      That is why in my library (Daedalus Lib) I proceeded differently. I choosed to merge points that are too close each other. For that I checked if the distance is smaller than the epsilon value.

      Delete
  5. Hi,

    Can you tell how can i interpolate z using method 1. I also have a z coordinate for each point. I think there will be a d along with a, b, c but don't know how to calculate it.

    Thanks

    ReplyDelete
    Replies
    1. Assuming you have P(x, y, z) and a triangle defined by P1(x1, y1, z1), P2(x2, y2, z2), P3(x3, y3, z3), you can try to solve:

      x = a * x1 + b * x2 + c * x3
      y = a * y1 + b * y2 + c * y3
      z = a * z1 + b * z2 + c * z3
      a + b + c = 1

      But having 3 free variables and 4 constraints means you are surely in trouble.

      You must add a new free degree in your system. For example, you could consider to have a new axis pointing orthogonaly to the plane defined by the triangle ; let's call it Ph(xh, yh, zh). Then you consider the system:

      x = a * x1 + b * x2 + c * x3 + t * xh
      y = a * y1 + b * y2 + c * y3 + t * yh
      z = a * z1 + b * z2 + c * z3 + t * zh
      a + b + c = 1

      and you solve it for a, b, c, t.

      The t value can be:
      t=0 if P is on the plane defined by the triangle
      t>0 if P is above
      t<0 if P is below

      Delete
  6. This is a great article! Very informative.
    I have one question though. Wouldn't it also be possible to acchive accurate results with only the barycentric method if you use it taking EPSILON into account?
    Something like this:

    (0 - EPSILON) <= a <= (1 + EPSILON) and
    (0 - EPSILON) <= b <= (1 + EPSILON) and
    (0 - EPSILON) <= c <= (1 + EPSILON)

    Computation-wise that would be faster than the distance check on each segment. Wouldn't it be just as accurate?

    ReplyDelete
    Replies
    1. The problem here is you don't have the control of the acceptance distance. This distance will be variable and will depend on the size of your triangle. So you face the danger to have wrong results with large triangles.

      Delete
    2. Aha I see that makes sense!
      What I am trying to do is a is-point-on-quad-test. The quad might be distorted in whatever way the user drags the corners. So my idea was to divide the quad into two triangles and simply test if the point is on either of those. There might be an arbitrary amount of quads, so I am a little cautious about computation time.
      Since the trinagles share a border is it save to assume, that the barycentric test including Epsilon will not fail, if the point lies directly on the border?

      Delete
    3. If your purpose is to identify the quad the user has clicked on, you can use your method, I think you will not face problems.

      Delete
    4. Awesome thank you very much for the help!
      I will go ahead and implement it that way. Also thanks again for this great article! I will definitely come back to it, the next time I need a accurate triangle/point test!

      Delete
  7. Aren't the first and second methods the same thing?

    ReplyDelete
    Replies
    1. Yes from a computational point of view, they are very close. But the concepts involved in the reasonings are different, that's why I found interesting to expose both.

      Delete
  8. A graphical area is a semi-open interval - for consistency this is usually done with lower bound inclusive, upper bound exclusive. A point on the boundary of two regions should not be in both; it should be in the region where the point is inclusive and excluded from the region where the point is exclusive. Very few graphics packages handle this issue properly - most have a fencepost error on the boundary. It should be possible to draw two squares from (0,0) to (10,10) and from (10,0) to (20,10) and have them abut exactly, as the pixels should be drawn from 0..9 and 10..19 - in which case a point such as (10,5) would belong to the (10,0) to (20,10) square only.

    ReplyDelete
  9. Very good article. Thanks a lot it was very usefull to me.

    ReplyDelete
  10. Hello, what about the winding number method [1] for the step 2?

    thanks!

    [1] http://geomalgorithms.com/a03-_inclusion.html

    ReplyDelete
    Replies
    1. Winding number method is a powerful method for n-sides polygons. But applied and optimized to 3 sides polygons, it is equivalent to the "check sides" method I exposed ;)

      Delete
    2. Thanks for you reply!

      Delete
  11. Ok, suppose that we have a general convex polygon, N vertices.

    We seek to place or reject a test point, P within the polygon.

    We select as many test triangles as possible, N points 3 at a time, and apply the test on each triangle.

    The test point is in the polygon if it rests within at least one of the [N,3] of the triangles.

    Right?

    ReplyDelete
    Replies
    1. Take care how you build your list of triangles. Considering ALL possible triangles made from the vertices list is not efficient. On the other hand, trying to optimize your triangles list, you could miss one or several.

      You should consider to build a Delaunay triangulation from your n-vertices polygon.

      Or check for algorithms directly dedicated to n-sides polygons.

      Delete
    2. Well, for N<7 its not too bad, but the Delaunay procedure is certainly worth looking into. Can "CM" of the polygon, X=(sum Xi, i=1..N)/N, play a useful role in sorting vertices?

      Delete
  12. Hi! Your article is very helpful. I am a beginner and I would like to know how to use these codes using IDL (Interactive Data Language). Thank you and more power.

    ReplyDelete
    Replies
    1. Hi.

      You must "translate" my AS3 code into IDL code.

      It is not difficult because my code only use basic features (variable, function, if-else structure) and a few maths being common in all programming languages.

      At this point I have never worked with IDL so I can't help you more.

      Delete
  13. This comment has been removed by the author.

    ReplyDelete
  14. Hi!
    Excellent article!
    Will this work with lat,lng coordinates (of earth, assuming perfect sphere)?

    ReplyDelete
    Replies
    1. Unfortunately it doesn't work on a sphere. If you flatten your local geometry to fit an euclidean 2D space, your triangle will be made from curves, that makes things a bit more complicated.

      Delete
  15. Dear Sir/Madam,

    Hi, I am new to this topic, I have tried your code but I cant get the expected result, below is some of the information:

    It is given that, I have the following coordinates:
    MinX = (834186.3789)
    MinY = (820799.625)
    MaxX = (834910.8751)
    MaxY = (821412.6876)

    Point A = (MaxX, 0)
    Point B = (MinX, 0)
    Point C = (0, MaxY)

    Point P = (834620, 821059)

    And form a Triangle ABC

    You can find that,the range of Px is between Ax and Bx; And Py is less than MaxY, thus, I expect that Point P should be in the Triangle ABC.

    but with the work it done; The program on the above return FALSE; which is not as expected.

    Can you please help me to have a look, is there any I misunderstood or misused your code?

    Many Thanks.

    ReplyDelete
    Replies
    1. Oops, to be much easy to verify; I updated the coordinate of Point P as below:
      Point P = (MinX, MinY)

      Thus, P = (834186.3789, 820799.625)

      As long as P now is intercept point to hit edge AC, it should inside the triangle ABC.

      but still, with the logic mentioned on the above, it return a false.

      Delete