Given a triangle ABC and a point P, let Q describe the point on ABC closest to P. One
way of obtaining Q is to rely on the fact that if P orthogonally projects inside ABC
the projection point is the closest point Q. If P projects outside ABC, the closest point
must instead lie on one of its edges. In this case, Q can be obtained by computing
the closest point to P for each of the line segments AB, BC, and CA and returning the
computed point closest to P. Although this works, it is not a very efficient approach.
A better solution is to compute which of the triangle’sVoronoi feature regions P is in.
Once determined, only the orthogonal projection of P onto the corresponding feature must be computed to obtain Q.
To see how P can be determined to be in a vertex Voronoi region, consider the
vertex Voronoi region of A. This region is determined as the intersection of the
negative halfspaces of two planes through A, one with a normal of B–A and the other
with the normal C–A (as illustrated in Figure 5.5).
Determining if P lies in one of the edgeVoronoi regions can be done in a number
of ways. It turns out that an efficient test is to effectively compute the barycentric
coordinates of the orthogonal projection R of P onto ABC. Recall from Section 3.4
that the barycentric coordinates of R are given as the ratios of (anything proportional to) the signed areas of triangles RAB, RBC, and RCA to the signed area of
ABC. Let n be the normal of ABC and let R = P − tn for some t. The barycentric coordinates (u, v, w) of R, R = uA + vB + wC, can then be computed from the
quantities
Vector n = Cross(b - a, c - a);
float rab = Dot(n, Cross(a - r, b - r)); // proportional to signed area of RAB
float rbc = Dot(n, Cross(b - r, c - r)); // proportional to signed area of RBC
float rca = Dot(n, Cross(c - r, a - r)); // proportional to signed area of RCA
float abc = rab + rbc + rca; // proportional to signed area of ABCas u = rbc/abc, v = rca/abc, and w = rab/abc. However, a little vector arithmetic
shows that these expressions simplify. For example, the expression for rab simplifies
as follows:
In other words, the barycentric coordinates of R can be obtained directly from P
without computing R.
For P to lie in an edge Voronoi region — for example, the one corresponding to
edge AB — P would have to lie outside or on AB, signified by rab ≤ 0, as well as within
the positive halfspaces of the planes (X − A) · (B − A) = 0 and (X − B) · (A − B) = 0.
Note that it is not sufficient just to test if P is outside AB, in that for a triangle with an
obtuse angle at A, P could be outside AB and actually be located in theVoronoi region
of edge CA (Figure 5.6). (Similarly, it is a common mistake to assume, for example,
that A is the closest point to P if P lies outside AB and (P − A) · (B − A) < 0.) If P is
not found to be in any of the vertex or edge Voronoi regions, Q must lie inside ABC and, in fact, be the orthogonal projection R, which can now be easily computed per the preceding. This information is now enough to produce a code solution.
Point ClosestPtPointTriangle(Point p, Point a, Point b, Point c)
{
Vector ab = b - a;
Vector ac = c - a;
Vector bc = c - b;
// Compute parametric position s for projection P’ of P on AB,
// P’ = A + s*AB, s = snom/(snom+sdenom)
float snom = Dot(p - a, ab), sdenom = Dot(p - b, a - b);
// Compute parametric position t for projection P’ of P on AC,
// P’ = A + t*AC, s = tnom/(tnom+tdenom)
float tnom = Dot(p - a, ac), tdenom = Dot(p - c, a - c);
if (snom <= 0.0f && tnom <= 0.0f) return a; // Vertex region early out
// Compute parametric position u for projection P’ of P on BC,
// P’ = B + u*BC, u = unom/(unom+udenom)
float unom = Dot(p - b, bc), udenom = Dot(p - c, b - c);
if (sdenom <= 0.0f && unom <= 0.0f) return b; // Vertex region early out
if (tdenom <= 0.0f && udenom <= 0.0f) return c; // Vertex region early out
// P is outside (or on) AB if the triple scalar product [N PA PB] <= 0
Vector n = Cross(b - a, c - a);
float vc = Dot(n, Cross(a - p, b - p));
// If P outside AB and within feature region of AB,
// return projection of P onto AB
if (vc <= 0.0f && snom >= 0.0f && sdenom >= 0.0f)
return a + snom / (snom + sdenom) * ab;
// P is outside (or on) BC if the triple scalar product [N PB PC] <= 0
float va = Dot(n, Cross(b - p, c - p));
// If P outside BC and within feature region of BC,
// return projection of P onto BC
if (va <= 0.0f && unom >= 0.0f && udenom >= 0.0f)
return b + unom / (unom + udenom) * bc;
// P is outside (or on) CA if the triple scalar product [N PC PA] <= 0
float vb = Dot(n, Cross(c - p, a - p));
// If P outside CA and within feature region of CA,
// return projection of P onto CA
if (vb <= 0.0f && tnom >= 0.0f && tdenom >= 0.0f)
return a + tnom / (tnom + tdenom) * ac;
// P must project inside face region. Compute Q using barycentric coordinates
float u = va / (va + vb + vc);
float v = vb / (va + vb + vc);
float w = 1.0f - u - v; // = vc / (va + vb + vc)
return u * a + v * b + w * c;
}
Comments