Count points inside triangle fast

Question

I need to find the number of points in a given list that are inside a triangle. The problem here is, there can be up to a million points.

I tried a simple approach: if the area of the triangle is equal to the sum of the areas of 3 triangles formed by taking 2 of the triangle's points at a time and the point to check, its inside. This doesn't have any precision errors, since I don't divide by two to find the area.

But I need something faster. The goal is speed. Can it be made faster by some sort of preprocessing, ignoring some points based on some criteria or something similar?

EDIT: Forgot to add a critical detail. The points once given, are fixed. Then the points are static, and need to be checked for up to a million triangles...

EDIT 2: Turns out that a good(maybe optimal too) way to do it was using a line sweep. Still, thanks for your answers!

Answer 1

According to the computational geometry weenies the fastest way to do this is by a barycentric coordinates transformation. In a case where you have a fixed triangle and many test points, then this approach will be especially fast, because once you have computed the barycentric coordinates of the 3 points of the triangle you have done most of the work. Here is the full algorithm where ABC are the triangle and P is the point under test:

// Compute vectors        
v0 = C - A
v1 = B - A
v2 = P - A

// Compute dot products
dot00 = dot(v0, v0)
dot01 = dot(v0, v1)
dot02 = dot(v0, v2)
dot11 = dot(v1, v1)
dot12 = dot(v1, v2)

// Compute barycentric coordinates
invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
u = (dot11 * dot02 - dot01 * dot12) * invDenom
v = (dot00 * dot12 - dot01 * dot02) * invDenom

// Check if point is in triangle
return (u >= 0) && (v >= 0) && (u + v < 1)

Here the barycentric coordinates are computed with respect to A, but B or C would have worked as well.

To test additional points you only need to recompute v2, dot02, dot12, u and v. Quantities such as invDenom stay the same.

Answer 2

A point is inside a triangle if it's to the left (right) of every side. You can calculate the cross products (just one component of it, actually) of a vector constructed of a point to be tested and one of the triangle vertices and the 3 vectors lying on the sides of the triangle (all in clock-wise or all in counter-clock-wise direction). See if the computed component of all 3 has the same sign (all 3 negative or all 3 positive). That'll tell you the point is in. Fast, no problems with precision, at least, if you use integers for the task.

You may stop further calculations for every point once you see it's to the wrong side of one of the triangle sides.

Sample code in C:

#include <stdio.h>

#define SCREEN_HEIGHT 22
#define SCREEN_WIDTH  78

// Simulated frame buffer
char Screen[SCREEN_HEIGHT][SCREEN_WIDTH];

void SetPixel(int x, int y, char color)
{
  if ((x < 0) || (x >= SCREEN_WIDTH) ||
      (y < 0) || (y >= SCREEN_HEIGHT))
    return;

  Screen[y][x] = color;
}

void Visualize(void)
{
  int x, y;

  for (y = 0; y < SCREEN_HEIGHT; y++)
  {
    for (x = 0; x < SCREEN_WIDTH; x++)
      printf("%c", Screen[y][x]);

    printf("\n");
  }
}

typedef struct
{
  int x, y;
} Point2D;

int main(void)
{
  // triangle vertices
  Point2D vx0 = { SCREEN_WIDTH / 2, SCREEN_HEIGHT / 7 };
  Point2D vx1 = { SCREEN_WIDTH * 6 / 7, SCREEN_HEIGHT * 2 / 3 };
  Point2D vx2 = { SCREEN_WIDTH / 7, SCREEN_HEIGHT * 6 / 7 };
  // vectors lying on triangle sides
  Point2D v0, v1, v2;
  // current point coordinates
  int x, y;

  // calculate side vectors

  v0.x = vx1.x - vx0.x;
  v0.y = vx1.y - vx0.y;

  v1.x = vx2.x - vx1.x;
  v1.y = vx2.y - vx1.y;

  v2.x = vx0.x - vx2.x;
  v2.y = vx0.y - vx2.y;

  // process all points

  for (y = 0; y < SCREEN_HEIGHT; y++)
    for (x = 0; x < SCREEN_WIDTH; x++)
    {
      int z1 = (x - vx0.x) * v0.y - (y - vx0.y) * v0.x;
      int z2 = (x - vx1.x) * v1.y - (y - vx1.y) * v1.x;
      int z3 = (x - vx2.x) * v2.y - (y - vx2.y) * v2.x;

      if ((z1 * z2 > 0) && (z1 * z3 > 0))
        SetPixel(x, y, '+'); // point is to the right (left) of all vectors
      else
        SetPixel(x, y, '-');
    }

  // draw triangle vertices

  SetPixel(vx0.x, vx0.y, '0');
  SetPixel(vx1.x, vx1.y, '1');
  SetPixel(vx2.x, vx2.y, '2');

  // visualize the result

  Visualize();

  return 0;
}

Output (ideone):

------------------------------------------------------------------------------
------------------------------------------------------------------------------
------------------------------------------------------------------------------
---------------------------------------0--------------------------------------
--------------------------------------++++------------------------------------
------------------------------------++++++++----------------------------------
----------------------------------+++++++++++++-------------------------------
--------------------------------+++++++++++++++++-----------------------------
------------------------------++++++++++++++++++++++--------------------------
----------------------------++++++++++++++++++++++++++------------------------
--------------------------+++++++++++++++++++++++++++++++---------------------
-------------------------++++++++++++++++++++++++++++++++++-------------------
-----------------------+++++++++++++++++++++++++++++++++++++++----------------
---------------------+++++++++++++++++++++++++++++++++++++++++++--------------
-------------------+++++++++++++++++++++++++++++++++++++++++++++++1-----------
-----------------++++++++++++++++++++++++++++++++++++-------------------------
---------------++++++++++++++++++++++++---------------------------------------
-------------++++++++++++-----------------------------------------------------
-----------2------------------------------------------------------------------
------------------------------------------------------------------------------
------------------------------------------------------------------------------
------------------------------------------------------------------------------

Answer 3

Simple pre-filter is to eliminate any points whose coordinates obviously lie outside of the bounds of the triangle's corners. e.g.

  a
+
|\
| \ b
|c \
+---+  d

A and D are obviously outside. A's Y coord is far above the triangle's maximum Y, and D is obviously beyond the triangle's maximum X.

That leaves B and C for testing.

Answer 4

You can also use a quadtree to accelerate the computation.

Compute a quadtree for the triangle (stopping at arbitrary resolution) and for each node (square), store a flag saying whether the node is completely inside, completely outside, or partially inside the triangle. A node that's partially inside the triangle can potentially have children (depending on depth)

For each point, traverse the quadtree. If we visit a node that's completely outside or inside the triangle, we're all set. If we're unsure about if we're in the triangle (node is partially inside the triangle) and the current node has children, we test against its children recursively. If we hit a leaf node that's partially inside the triangle, do an analytical point - triangle containment check.