| Serial | SSE/NEON | AVX2 | AVX512 |
 |
 |
 |
 |
qTriangle is a personal study to design a quick way to detect if a point is within a Triangle by means of vectorization
The domain of the Point-In-Triangle problem is determining if a cartesian coordinate happens to land upon the interior of a triangle. In this case the 2D case of triangles will be examined and will require some kind of surface area for a point to land on so a case in which all three points are collinear(which is the worst case of having a very slim triangle) are ruled out.
This problem comes up a lot in the domain of computer graphics and gameplay programming at times. Sometimes it's testing if a single point lands upon a polygon(made up of triangles) and sometimes it's testing if thousands of points happen to land on a triangle or not(such as when rendering a vector triangle against a regular grid during rasterization).
There are several methods to test if a point happens to land within a triangle in 2D space, each have their own pros and cons and scaling properties.
The cross-product operation in vector algebra is a binary operation that takes two vectors in 3D space and creates a new vector that is perpendicular to them both. You can think of the two vectors as describing some kind of plane in 2D space, and the cross product creates a new vector that is perpendicular to this plane. Though, there are two ways to create a vector perpendicular to a plane, by going "in" the plane and going "out" the plane.
The Cross-Product has a lot of useful properties but one of interest at the moment is the magnitude of the resulting vector of the cross product which will be the area of the parallelogram that the two original vectors create. Since the points are on the X-Y plane in 3D space this magnitude will always be the Z-component of the cross product since all vectors perpendicular to the orthogonal X-Y plane will take the form [ 0, 0, (some value)]. This becomes useful later on.
The particular numerical value of this area is not of importance either but rather the parity of the area is of interest(while a negative-surface-area does not make sense, this value tells us something about two directional vectors). Notice that whenever the direction to the moving point goes to the "left" of the black directional vector, that the area becomes negative, but when it is to the "right", the area is positive. This is due to the right hand rule where the direction of positive-rotation being clockwise or counter-clockwise causes the order of the original two directional vectors to determine the proper orientation of the cross product.
The parity of the cross-product-area depending on which side the direction of the "point" lands on solves the problem at the same tone of turning each edge into a linear inequality then testing if the point solves each of them at once but in a somewhat more optimal way.
The three positional vectors of the triangle must be in clockwise or counter-clockwise order so that three directional vectors can be determinately created.
EdgeDir0 = Vertex1 - Vertex0
EdgeDir1 = Vertex2 - Vertex1
EdgeDir2 = Vertex0 - Vertex2
Then, three additional directional vectors can be made that point from the triangle vertex position to the point that is being tested against
PointDir0 = Point - Vertex0
PointDir1 = Point - Vertex1
PointDir2 = Point - Vertex2
Now, finding out if a point lands within the triangle is determined by using three cross-products, and checking if each area is positive. If they are all positive. Then the point is to the "right" of all the edges. If any of them are negative, then it is not within the triangle.
| EdgeDir0 × PointDir0 | >= 0 &&
| EdgeDir1 × PointDir1 | >= 0 &&
| EdgeDir2 × PointDir2 | >= 0
Previously it was determined that all cross products against the X-Y plane will take the form [ 0, 0, (some value)]. With this the arithmetic behind the cross-product operation can be much more simplified. Since the cross product can be calculated using partial determinants of a 3x3 matrix then attention only has to be given to the calculations that determine the Z-component alone which is but a 2x2 determinant of the two input vectors.
This means that if I had two vectors A and B on the X-Y plane. The cross-product's magnitude is simply:
A.x * B.y - A.y * B.x;
which reduces the previous arithmetic to:
EdgeDir0.x * PointDir0.y - EdgeDir0.y * PointDir0.x >= 0 &&
EdgeDir1.x * PointDir1.y - EdgeDir1.y * PointDir1.x >= 0 &&
EdgeDir2.x * PointDir2.y - EdgeDir2.y * PointDir2.x >= 0
The full pseudo-code:
// Point - Position that is being tested
// Vertex0,1,2 - Vertices of the triangle in **clockwise order**
// Directional vertices along the edges of the triangle in clockwise order
EdgeDir0 = Vertex1 - Vertex0
EdgeDir1 = Vertex2 - Vertex1
EdgeDir2 = Vertex0 - Vertex2
// Directional vertices pointing from the triangle vertices to the point
PointDir0 = Point - Vertex0
PointDir1 = Point - Vertex1
PointDir2 = Point - Vertex2
// Test if each cross-product results in a positive area
if(
EdgeDir0.x * PointDir0.y - EdgeDir0.y * PointDir0.x >= 0 &&
EdgeDir1.x * PointDir1.y - EdgeDir1.y * PointDir1.x >= 0 &&
EdgeDir2.x * PointDir2.y - EdgeDir2.y * PointDir2.x >= 0
)
{
// CurPoint is in triangle!
}If I was to throw thousands of points at a triangle in a for-loop using this algorithm then not all variables have to be re-calculated for each point.
The vectors EdgeDir0, EdgeDir1, EdgeDir2 only have to be calculated once. For each point the vectors PointDir0, PointDir1, PointDir2 have to be recreated.
EdgeDir0 = Vertex1 - Vertex0
EdgeDir1 = Vertex2 - Vertex1
EdgeDir2 = Vertex0 - Vertex2
foreach(CurPoint in LotsOfPoints)
{
PointDir0 = Point - Vertex0
PointDir1 = Point - Vertex1
PointDir2 = Point - Vertex2
if(
EdgeDir0.x * PointDir0.y - EdgeDir0.y * PointDir0.x >= 0 &&
EdgeDir1.x * PointDir1.y - EdgeDir1.y * PointDir1.x >= 0 &&
EdgeDir2.x * PointDir2.y - EdgeDir2.y * PointDir2.x >= 0
)
{
// CurPoint is in triangle!
}
}
Which results in the total overhead for each point being
| Subtractions | Multiplications | Comparisons |
|---|---|---|
| 6 | 6 | 3 |
With some barycentric coordinate trickery one can derive a coordiante system that allows a triangle to be described as a linear "mixture" of its three vertices, with some constraints on how they mix.
The Barycentric coordinates of a triangle involves it's three position-vectors p1, p2, p3.
Using these three points, any new point p' within this triangle can be generated by mixing the three vertex positions according to three scalar weights w1, w2, w3 such that:
p' = w1 * p1 + w2 * p2 + w3 * p3.
This is just a linear combination of three points, but p' isnt meant to just be any wild combination of three points. What you actually want is for the combination of three points to always land on the triangular surface that they all contain. This is called a convex combination where if you want to only reach every point within this triangle then you must bind these three weight weight to the conditions of being non-negative and also summing to the exact value 1.0. Making this the full barycentric equation for a triangle defined by three points:
In this case though. You already have p' and want to determine if it is within this triangle! Time to cleverly work backwards!
So given a p1, p2, p3 and p', you have to figure out the three positive unknowns w1, w2, w3 that balance these conditions. Some linear alg!
So starting with this:
How about expanding those vectors into their individual(2D) dimensions.
Now that looks like a matrix! Take out all those weights and put it into a vector!
So the solution is to invert that 3x3 matrix there and multiplying it by the point we are testing it against to get the resulting three weights. And once you get the three weights. All you have to do is test if they are positive(Note: If this smells somewhat like a close derivation of the cross-product method you're right!)
Matrix inverse... This is where it gets a little hairy
And with this, w1, w2, and w3 all reduce to relatively "simple" linear equations and all that would have to be checked is if they are greater than 0!
Here is the equivalent code in GLSL. Keep in mind that GLSL defines matrices and vectors as column-major. So going vec3(1) means a 1x3(width x height) arrangement of elements and mat3(vec(1),vec(2),vec(3)) means a 3x3 matrix with the columns(the vertical ones) being defined by {1,1,1},{2,2,2},{3,3,3}.
bool PointInTriangleBarycentric(
in vec2 Triangle[3],
in vec2 Point
)
{
mat3 Barycentric = inverse(
mat3(
Triangle[0], 1.0f,
Triangle[1], 1.0f,
Triangle[2], 1.0f
)
);
vec3 Weights = Barycentric * vec3( Point, 1.0 );
if(
// Weights.x >= 0.0f &&
// Weights.y >= 0.0f &&
// Weights.z >= 0.0f
all( greaterThanEqual( Weights, vec3(0.0f) ) )
)
{
return true;
}
return false;
}Though, this is a pretty naive approach. It still has its uses. If you wanted to test a million points against a single triangle. You'd have to do a single matrix inverse(pretty expensive!) and then the overhead for testing each individual point you want to test would be: a matrix-vector multiplication(a 3x3 matrix times an ℝ³ vector, which is pretty much just three dot-products) and three comparisons:
| Additions | Multiplications | Comparisons |
|---|---|---|
| 6 | 9 | 3 |
Though, a matrix inverse is pretty expensive operation to do. Especially in a context of having possibly thousands upon millions of triangles and a matrix inverse each involving almost dozens of multiplications, additions, and a division on top of it all. Some more clever observations can allow for some slim optimizations. Especially since all that matters is if the weights are positive or not.
The two directional vectors are derived from three points that describe a plane while the two scalars are typically denoted as U and V and determine how much these two directional vectors should mix together to create another point on this plane. Since a triangle has three points, these two vectors can be derived by picking any one point of the triangle and obtaining two directional vectors from this point to the two other points.
The U and V scalar values are typically normalized within the [0.0,1.0] range to easily translate these values into percentages that determine how much much the two vectors should contribute to the resulting vector.
Ex, if I had the two directional vectors ( 8, 2 ) and ( 7, 1 ) and the U V values 0.5 1.0 respectively. The resulting point is ( 8, 2 ) * 0.5 + ( 7, 1 ) * 1.0 = ( 11, 3 ) which in english would be something like I want 50% of the ( 8, 2 ) direction and 140% of ( 7, 1 ).
The actual triangle shape is made by constraining the U and V values so that rather than describing a plane very generally it instead will stay within the bounds of a triangle.
The first two constraints are U >= 0 and V >= 0 which guarentee that the U,V coordinates are always on the positive side of the two vectors and do not go backwards, off the triangle. The third constraint is U + V <= 1 is just the line y = 1 - x turned into an inequality such that all solutions to the inequality equate to a point within a triangle. The actual derivation of this involves some barycentric coordinate limbo.
Given a triangle, the two directional vectors are easy to calculate. Pick any of the three points of a triangle, get the vector direction from this point, to the two other points ( which is just a vector subtraction). After the two directional vectors are obtained, now all that has to be done to see if a point is within a triangle is projecting this point against the two directional vectors to get the U and V values to test against the three conditions.
Projecting a point against these two vectors to get the U and V values is but trivial dot-product arithmetic. First, another directional vector has to be created as a positional-vector and a directional-vector wouldn't make sense in this instance. The very same point that was selected to in step 1 to generate the first two directional vectors must be used once more to generate a third directional vector between this triangle vertex and the point being sampled against (another vector subtraction). This new vector will then be dot-product-ed against the two directional edges ( vertial vector multiplication and horizontal addition ) to finally determine the U and V values to test against U >= 0, V >= 0, and U + V <= 1.
