## Freesteel Blog » Finding points in triangles

## Finding points in triangles

Thursday, May 17th, 2007 at 2:51 pm

Just recently I discovered a flaw in my usual algorithm for choosing the triangle that a point is inside. This means the error is in everything I have written up till now. It’s a pretty rare case, and manifested at the university where they were trimming a line through the origin with a triangulation of a cylinder.

The problem resolves to the following question. Given a point in an area, and a triangulation of that area, which single triangle contains the point?

We want a single answer, calculated using all the reliability of floating point arithmetic, where there is no such thing as being perfectly on an edge. (If it is on an edge, which is practically impossible unless the edge is aligned with an axis or at 45 degrees, then it doesn’t matter which side)

**Point in triangle test**

The first part of the usual solution is to solve it for a single triangle. If you have a reliable test which tells you whether your point is in a given triangle, then all you need to do is loop through the triangles in the area, and return the one which satisfies the test.

We are given a point `p` and a triangle with vertices `a`, `b`, and `c`

You can find on-line some bad implementations of the point in triangle test, and I have even once witnessed an implementation of the “worst algorithm in the world for testing points” in a triangle, known as the angle summation test.

One common technique is to find the coordinates of the point in the triangle by solving the equation:

p = a + s(b - a) + t(c - a)

for the scalars `s` and `t`. The point is inside the triangle when `s > 0`, `t > 0`, and `s + t < 1`.

The actual way forward is to use dot-products. If you have a vector `v = (x, y)`, then the clockwise perpendicular `v ^{CT}` is

`(y, -x)`. We can say the point

`p`is to the right of the oriented line

`ab`if:

(p - a) .(b - a)^{CT}> 0

This is using the dot-product operator. The condition is determined by 5 additions and two multiplications (a subtraction is the same as an addition with the sign changed).

If `p` is to the right of the oriented lines `ab`, `bc`, and `ca`, then it’s inside the triangle (and the triangle has a clockwise orientation). If it’s to the left of all these lines, then it’s also inside the triangle (and the triangle has a anti-clockwise orientation).

So, we have a quick test for the sideness between a point and a line, and if the point is on the same side of the three lines joining three points, then it must be inside the defined triangle. There are obvious ways to optimize the function to quit early and reuse calculations.

**Point in triangulation test**

The trouble with these tests is the boundaries. Owing to the behavior of floating point arithmetic, a point that is very close to the line might give differing results depending on how you do the test. It is possible to get the contradictory result that `p` is to the right of `ab`, and also to the right of `ba`. The important asymmetry is the `(p - [a or b])` value, which will lose variable amounts of precision according to whether the coordinates of `a` or `b` have different enough exponents to vanish different parts of the mantissas of the coordinates of `p`. The fact that there are two coordinates means they can act in opposite directions to ensure that the results are truly unpredictable.

The result is, if you apply the point in triangle test on each triangle independently, you might wind up with the point being in no triangle or in two triangles when it is too close to an edge and the sidedness is decided inconsistently.

Until two days ago I thought all you needed to do to fix it was decide which side the point was in relation to each edge once, and stick to the answer. That way when the point was on an edge it would be sure to be in exactly one of the two adjacent triangles.

Unfortunately, I found the counter-example in the form of the point **b** in the picture above where it almost perfectly coincides with a vertex of the triangulation. Remember, the decision for a point that is in-line with an edge can go either way. If all the edges connecting to the vertex settle for sides that happen to all be in the same direction, then the algorithm will determine that the point is in none of them!

**Nasty fixes**

There are two quick fixes which you can apply if you can’t be bothered to get it right. The first is to add a sub-tolerant value that is too small for anyone to notice to the position of your point so it will be unlikely to coincide with any of the vertices. If there was such a thing as an inspector for auditing geometric modeling source code, they would know how to look for this bodge.

The second and most common fix is to bring in an epsilon value. It’s kind of the opposite attitude to dithering the point so it is safely away from the special cases, because you instead draw it into them. When the given point is with an epsilon distance of an edge, you treat it as though it were exactly on the edge, and if it’s within epsilon of a vertex, you make a special case out of that too.

This fix looks appealing because it appears to be covering all bases properly, subject to messing with the incoming data just as much as in the first fix, although not all the time. However, it does introduce discontinuities in behavior, and sometimes the special cases are so hard and rarely occurring that you are tempted to cut corners. The case I have in mind is when the triangulated plane is self-overlapping because it comes from the flat projection of a triangulated surface, and one of the points is surrounded by a shape that looks like a partly-folded umbrella.

I really don’t like epsilon value fixes, although they are rife throughout the geometric modeling world, and I am in a tiny minority for recognizing that they are a bad thing in the very long term. Why? If you set the value too small or too large, the system breaks down in various different ways. Therefore there is a range of values where things all work. As the system gets larger in terms of how much code has been added or how widely it is applied, the interval where it works can only get narrower. One day you might find that it disappears entirely, and there can be no value that makes the system reliable.

**Nice fixes**

I’ve not worked it out yet, but I believe there is a robust answer to the problem where you first cut the triangulated area with a horizontal line through the point so that you get left with a partition of intervals along that line. Then you locate the `x` value along that line and read off which interval and therefore triangle it is within.

The best answer is probably to make an analogue of my line cutting offset polygon function, and forgo any chance of making it work for an offset of zero. This is not the same as the an epsilon solution, since the range of values for which it works only has a lower bound. You can set it to take offsets that are large enough to close any gaps so that it can accept triangulations that are not watertight and avoid having to build a complex algorithm to “heal” the cracks.

There are several such healing algorithms in the world, and I imagine that they are unbelievably complicated. I don’t intend ever to write one, because I believe it is always be possible to reconsider what it is being used for and whether the watertightness condition is absolutely necessary.

## 2 Comments

1.SJ replies at 17th May 2007, 11:49 pm :Assuming this is a 2d problem and the triangulation is clean (with edges which meet at a common point doing so exactly and sharing a single vertex record), then you can build a robust ‘which triangle’ test from a single primitive – calculating the sign of the area of a triangle from its vertex coordinates (which is just a cross product or a 3×3 determinant of the homogeneous coordinates if you prefer, don’t bother to divide by two – the sign is all you need). So long as this gives consistent results for cyclic permutations of the vertex order (which you can arrange by sorting the 3 coordinates lexicographically on input) then there is no need to use infinite precision arithmetic to guarantee always getting a sensible answer – the answer will just be wrong if the primitive is ‘wrong’ when your test point lies very near a line or a vertex, and you can even cope if a triangle of the triangulation itself comes out the wrong way round (so long as you check them all for this). Your ‘partition of intervals’ idea is exactly how efficient and robust plane sweep operations like Fortune’s voronoi algorithm work but they do so exactly by using determinants to keep the lines meeting at a vertex of interest in sorted order without actually calculating where they cross the horizontal line of interest. For many test points pre-sort the vertices into lexicographic order; for many triangles cull out those whose bounding boxes cannot overlap the test point.

2.Freesteel » Blog Ar&hellip replies at 11th June 2007, 11:02 am :[…] Point inside triangle correct

June 11th, 2007 by Julian

Following from the incorrect method for finding which triangle of a mesh a given point is ins […]

## Leave a comment