How to determine if a point lies OVER a triangle in 3D

Question

I need an example of fast algorithm allowing to calculate if a point lies over a triangle in 3D. I mean if the projection of this point on a plane containing given triangle is inside of this triangle.

I need to calculate distance between a point and a triangle (between a point and the face of this triangle if its projection lies inside the triangle or between a point and an edge of a triangle if its projection lays outside the triangle).

I hope I made it clear enough. I found some examples for 2D using barycentric coordinates but can't find any for 3D. Is there a faster way than calculating projection of a point, projecting this projected point and a given triangle to 2D and solving standard "point in triangle" problem?

Answer 1

If the triangle's vertices are A, B, C and the point is P, then begin by finding the triangle's normal N. For this just compute N = (B-A) X (C-A), where X is the vector cross product.

For the moment, assume P lies on the same side of ABC as its normal.

Consider the 3d pyramid with faces ABC, ABP, BCP, CAP. The projection of P onto ABC is inside it if and only if the dihedral angles between ABC and each of the other 3 triangles are all less than 90 degrees. In turn, these angles are equal to the angle between N and the respective outward-facing triangle normal! So our algorithm is this:

Let N = (B-A) X (C-A), N1 = (B-A) X (P-A), N2 = (C-B) X (P-B), N3 = (A-C) X (P-C)
return N1 * N >= 0 and N2 * N >= 0 and N3 * N >= 0;

The stars are dot products.

We still need to consider the case where P lies on the opposite side of ABC as its normal. Interestingly, in this case the vectors N1, N2, N3 now point into the pyramid, where in the above case they point outward. This cancels the opposing normal, and the algorithm above still provides the right answer. (Don't you love it when that happens?)

Cross products in 3d each require 6 multiplies and 3 subtractions. Dot products are 3 multiplies and 2 additions. On average (considering e.g. N2 and N3 need not be calculated if N1 * N < 0), the algorithm needs 2.5 cross products and 1.5 dot products. So this ought to be pretty fast.

If the triangles can be poorly formed, then you might want to use Newell's algorithm in place of the arbitrarily chosen cross products.

Note that edge cases where any triangle turns out to be degenerate (a line or point) are not handled here. You'd have to do this with special case code, which is not so bad because the zero normal says much about the geometry of ABC and P.

Here is C code, which uses a simple identity to reuse operands better than the math above:

#include <stdio.h>

void diff(double *r, double *a, double *b) {
  r[0] = a[0] - b[0];
  r[1] = a[1] - b[1];
  r[2] = a[2] - b[2];
}

void cross(double *r, double *a, double *b) {
  r[0] = a[1] * b[2] - a[2] * b[1];
  r[1] = a[2] * b[0] - a[0] * b[2];
  r[2] = a[0] * b[1] - a[1] * b[0];
}

double dot(double *a, double *b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

int point_over_triangle(double *a, double *b, double *c, double *p) {
  double ba[3], cb[3], ac[3], px[3], n[3], nx[3];

  diff(ba, b, a);
  diff(cb, c, b);
  diff(ac, a, c);

  cross(n, ac, ba);  // Same as n = ba X ca

  diff(px, p, a);
  cross(nx, ba, px);
  if (dot(nx, n) < 0) return 0;

  diff(px, p, b);
  cross(nx, cb, px);
  if (dot(nx, n) < 0) return 0;

  diff(px, p, c);
  cross(nx, ac, px);
  if (dot(nx, n) < 0) return 0;

  return 1;
}

int main(void) {
  double a[] = { 1, 1, 0 };
  double b[] = { 0, 1, 1 };
  double c[] = { 1, 0, 1 };
  double p[] = { 0, 0, 0 };

  printf("%s\n", point_over_triangle(a, b, c, p) ? "over" : "not over");

  return 0;
}

I've tested it lightly and it seems to be working fine.

Answer 2

Let's assume that the vertices of the triangle are v, w, and the origin 0. Let's call the point p.

For the benefit of other readers, here's the barycentric approach for 2D point-in-triangle, to which you alluded. We solve the following system in variables beta:

[v.x w.x] [beta.v]   [p.x]
[v.y w.y] [beta.w] = [p.y] .

Test whether 0 <= beta.v && 0 <= beta.w && beta.v + beta.w <= 1.

For 3D projected-point-in-triangle, we have a similar but overdetermined system:

[v.x w.x] [beta.v]   [p.x]
[v.y w.y] [beta.w] = [p.y] .
[v.z w.z]            [p.z]

The linear least squares solution gives coefficients beta for the point closest to p on the plane spanned by v and w, i.e., the projection. For your application, a solution via the following normal equations likely will suffice:

[v.x v.y v.z] [v.x w.x] [beta.v]   [v.x v.y v.z] [p.x]
[w.x w.y w.z] [v.y w.y] [beta.w] = [w.x w.y w.z] [p.y] ,
              [v.z w.z]                          [p.z]

from which we can reduce the problem to the 2D case using five dot products. This should be comparable in complexity to the method that Nico suggested but without the singularity.