Hello,
I am using Triangle intersection algorithm of Tomas Moller in my win32 project to check intersection between two triangles which works perfectly. However I have problem to find intersection points between two given triangles.
I’ve tried to debug its algorithm but couldn’t find a way to get an answer for my question.
It would be appreciated if you can help me to have algorithm to find coordinates of each points?
// Triangle intersection. Implementation of Tomas Moller, 1997.
// See article "A Fast Triangle-Triangle Intersection Test", Journal of Graphics Tools, 2(2), 1997.
class TriangleIntersection {
private:
// Sort in descending order
void Sort(CoordPoints2D v) {
if (v.x > v.y) {
float c;
c = v.x;
v.x = v.y;
v.y = c;
}
}
// Edge to edge test is based on Franlin Antonio's gem: "Faster Line Segment Intersection", in Graphics Gems III, pp. 199-202
bool EdgeToEdgeTest(vec3 v0, vec3 v1, vec3 u0, vec3 u1, int i0, int i1) {
float Ax, Ay, Bx, By, Cx, Cy, e, d, f;
Ax = v1[i0] - v0[i0];
Ay = v1[i1] - v0[i1];
Bx = u0[i0] - u1[i0];
By = u0[i1] - u1[i1];
Cx = v0[i0] - u0[i0];
Cy = v0[i1] - u0[i1];
f = Ay * Bx - Ax * By;
d = By * Cx - Bx * Cy;
if ((f > 0 && d >= 0 && d <= f) || (f < 0 && d <= 0 && d >= f)) {
e = Ax * Cy - Ay * Cx;
if (f > 0) {
if (e >= 0 && e <= f)
return true;
}
else {
if (e <= 0 && e >= f)
return true;
}
}
return false;
}
// Edge to triangle edges test
bool EdgeAgainstTriangleEdges(vec3 v0, vec3 v1, vec3 u0, vec3 u1, vec3 u2, short i0, short i1) {
// Test edge u0, u1 against v0, v1
if (EdgeToEdgeTest(v0, v1, u0, u1, i0, i1))
return true;
// Test edge u1, u2 against v0, v1
if (EdgeToEdgeTest(v0, v1, u1, u2, i0, i1))
return true;
// Test edge u2, u1 against v0, v1
if (EdgeToEdgeTest(v0, v1, u2, u0, i0, i1))
return true;
return false;
}
// Check whether the point is in triangle
bool PointInTriangle(vec3 v0, vec3 u0, vec3 u1, vec3 u2, short i0, short i1) {
float a, b, c, d0, d1, d2;
// Is triangle 1 completly inside triangle 2? let's Check if v0 is inside triangle (u0, u1, u2)
a = u1[i1] - u0[i1];
b = -(u1[i0] - u0[i0]);
c = -a * u0[i0] - b * u0[i1];
d0 = a * v0[i0] + b * v0[i1] + c;
a = u2[i1] - u1[i1];
b = -(u2[i0] - u1[i0]);
c = -a * u1[i0] - b * u1[i1];
d1 = a * v0[i0] + b * v0[i1] + c;
a = u0[i1] - u2[i1];
b = -(u0[i0] - u2[i0]);
c = -a * u2[i0] - b * u2[i1];
d2 = a * v0[i0] + b * v0[i1] + c;
if (d0 * d1 > 0.0f && d0 * d2 > 0.0f)
return true;
return false;
}
bool TriTriCoplanar(vec3 N, vec3 v0, vec3 v1, vec3 v2, vec3 u0, vec3 u1, vec3 u2) {
float *A = new float[3];
short i0, i1;
// First project onto an axis-aligned plane, that maximizes the area of the triangles, compute indices: i0, i1
A[0] = fabs(N[0]);
A[1] = fabs(N[1]);
A[2] = fabs(N[2]);
if (A[0] > A[1]) {
// A[0] is greatest
if (A[0] > A[2]) {
i0 = 1;
i1 = 2;
}
// A[2] is greatest
else {
i0 = 0;
i1 = 1;
}
}
else {
// A[2] is greatest
if (A[2] > A[1]) {
i0 = 0;
i1 = 1;
}
// A[1] is greatest
else {
i0 = 0;
i1 = 2;
}
}
// Test 1st edge of triangle 1 against all the edges of triangle 2
if (EdgeAgainstTriangleEdges(v0, v1, u0, u1, u2, i0, i1))
return true;
// Test 2nd edge of triangle 1 against all the edges of triangle 2
if (EdgeAgainstTriangleEdges(v1, v2, u0, u1, u2, i0, i1))
return true;
// Test 3rd edge of triangle 1 against all the edges of triangle 2
if (EdgeAgainstTriangleEdges(v2, v0, u0, u1, u2, i0, i1))
return true;
// Is triangle 1 is completely inside in triangle 2?
if (PointInTriangle(v0, u0, u1, u2, i0, i1))
return true;
// Is triangle 2 is completely inside in triangle 1?
if (PointInTriangle(u0, v0, v1, v2, i0, i1))
return true;
return false;
}
bool ComputeIntervals(float VV0, float VV1, float VV2, float D0, float D1, float D2, float D0D1, float D0D2, float &A, float &B, float &C, float &X0, float &X1) {
// We know that D0D2 <= 0.0 that is D0, D1 are on the same side, D2 on the other or on the plane
if (D0D1 > 0.0f) {
A = VV2;
B = (VV0 - VV2) * D2;
C = (VV1 - VV2) * D2;
X0 = D2 - D0;
X1 = D2 - D1;
}
// We know that d0d1 <= 0.0
else if (D0D2 > 0.0f) {
A = VV1;
B = (VV0 - VV1) * D1;
C = (VV2 - VV1) * D1;
X0 = D1 - D0;
X1 = D1 - D2;
}
// We know that d0d1 <= 0.0 or that D0 != 0.0
else if (D1 * D2 > 0.0f || D0 != 0.0f) {
A = VV0;
B = (VV1 - VV0) * D0;
C = (VV2 - VV0) * D0;
X0 = D0 - D1;
X1 = D0 - D2;
}
else if (D1 != 0.0f) {
A = VV1;
B = (VV0 - VV1) * D1;
C = (VV2 - VV1) * D1;
X0 = D1 - D0;
X1 = D1 - D2;
}
else if (D2 != 0.0f) {
A = VV2;
B = (VV0 - VV2) * D2;
C = (VV1 - VV2) * D2;
X0 = D2 - D0;
X1 = D2 - D1;
}
else
return true;
return false;
}
public:
// Checks if the triangle V(v0, v1, v2) intersects the triangle U(u0, u1, u2).
// Returns true if V intersects U, otherwise false
bool TriangleToTriangleIntersect(vec3 v0, vec3 v1, vec3 v2, vec3 u0, vec3 u1, vec3 u2) {
vec3 edge1, edge2;
vec3 n1, n2;
vec3 dd;
CoordPoints2D isect1, isect2;
float du0, du1, du2, dv0, dv1, dv2, d1, d2;
float du0du1, du0du2, dv0dv1, dv0dv2;
float vp0, vp1, vp2;
float up0, up1, up2;
float bb, cc, max;
short index;
// Compute plane equation of triangle(v0, v1, v2)
edge1 = v1 - v0;
edge2 = v2 - v0;
n1 = cross(edge1, edge2);
d1 = -dot(n1, v0);
// Plane equation 1: N1.X +d1 = 0, let's put u0, u1, u2 into plane equation 1 to compute signed distances to the plane
du0 = dot(n1, u0) + d1;
du1 = dot(n1, u1) + d1;
du2 = dot(n1, u2) + d1;
// Coplanarity robustness check
if (fabs(du0) < Epsilon)
du0 = 0.0f;
if (fabs(du1) < Epsilon)
du1 = 0.0f;
if (fabs(du2) < Epsilon)
du2 = 0.0f;
du0du1 = du0 * du1;
du0du2 = du0 * du2;
// Same sign on all of them + not equal 0 ? No intersection occurs
if (du0du1 > 0.0f && du0du2 > 0.0f)
return false;
// Compute plane of triangle (u0, u1, u2)
edge1 = u1 - u0;
edge2 = u2 - u0;
n2 = cross(edge1, edge2);
d2 = -dot(n2, u0);
// Plane equation 2: n2.X + d2 = 0, let's put v0,v1,v2 into plane equation 2
dv0 = dot(n2, v0) + d2;
dv1 = dot(n2, v1) + d2;
dv2 = dot(n2, v2) + d2;
if (fabs(dv0) < Epsilon)
dv0 = 0.0f;
if (fabs(dv1) < Epsilon)
dv1 = 0.0f;
if (fabs(dv2) < Epsilon)
dv2 = 0.0f;
dv0dv1 = dv0 * dv1;
dv0dv2 = dv0 * dv2;
// Same sign on all of them + not equal 0 ?, No intersection occurs
if (dv0dv1 > 0.0f && dv0dv2 > 0.0f)
return false;
// Compute direction of intersection line
dd = cross(n1, n2);
// Compute and index to the largest component of D
max = (float)fabs(dd[0]);
index = 0;
bb = (float)fabs(dd[1]);
cc = (float)fabs(dd[2]);
if (bb > max) {
max = bb;
index = 1;
}
if (cc > max) {
max = cc;
index = 2;
}
// This is the simplified projection onto L
vp0 = v0[index];
vp1 = v1[index];
vp2 = v2[index];
up0 = u0[index];
up1 = u1[index];
up2 = u2[index];
// Compute interval for triangle 1
float a = 0, b = 0, c = 0, x0 = 0, x1 = 0;
if (ComputeIntervals(vp0, vp1, vp2, dv0, dv1, dv2, dv0dv1, dv0dv2, a, b, c, x0, x1))
return TriTriCoplanar(n1, v0, v1, v2, u0, u1, u2);
// Compute interval for triangle 2
float d = 0, e = 0, f = 0, y0 = 0, y1 = 0;
if (ComputeIntervals(up0, up1, up2, du0, du1, du2, du0du1, du0du2, d, e, f, y0, y1))
return TriTriCoplanar(n2, v0, v1, v2, u0, u1, u2);
float xx, yy, xxyy, tmp;
xx = x0 * x1;
yy = y0 * y1;
xxyy = xx * yy;
tmp = a * xxyy;
isect1.x = tmp + b * x1 * yy;
isect1.y = tmp + c * x0 * yy;
tmp = d * xxyy;
isect2.x = tmp + e * xx * y1;
isect2.y = tmp + f * xx * y0;
Sort(isect1);
Sort(isect2);
return !(isect1.y < isect2.x || isect2.y < isect1.x);
}
};