- template<class T> inline bool edgeIntersectTriangle3D(
- const Vector3<T>& v0, const Vector3<T>& v1, const Vector3<T>& v2,
- const Vector3<T>& edge0, const Vector3<T>& edge1)
- {
- // algorithm descrption see:
- // intersect ray with triangle plane
- Vector3<T> dir = edge1-edge0;
- T edgeLen = dir.length();
- dir /= edgeLen;
- Vector3<T> normal = (v1-v0) % (v2-v0);
- normal.normalize();
- T normalDotDir = normal*dir;
- if (abs(normalDotDir) < EPSILON) {
- // edge and triangle are parallel
- T d = -(v0*normal);
- if (abs(edge0*normal + d) < EPSILON || abs(edge1*normal + d) < EPSILON) {
- // edge and triangle are in the same plane
- // find and remove dominant axis
- Vector3<T> absolute = normal.abs();
- int dom = 0;
- if (absolute[1] > absolute[dom]) { dom = 1; }
- if (absolute[2] > absolute[dom]) { dom = 2; }
- int axis1 = (dom + 1) % 3;
- int axis2 = (dom + 2) % 3;
- // project to 2D space
- Vector2<T> v0p(v0[axis1], v0[axis2]);
- Vector2<T> v1p(v1[axis1], v1[axis2]);
- Vector2<T> v2p(v2[axis1], v2[axis2]);
- Vector2<T> edge0p(edge0[axis1], edge0[axis2]);
- Vector2<T> edge1p(edge1[axis1], edge1[axis2]);
- // 2D intersection
- if (pointInTriangle2D(v0p, v1p, v2p, edge0p) ||
- pointInTriangle2D(v0p, v1p, v2p, edge1p) ) {
- return true;
- }
- if (edgeIntersection2D(v0p, v1p-v0p, edge0p, edge1p-edge0p) ||
- edgeIntersection2D(v1p, v2p-v1p, edge0p, edge1p-edge0p) ||
- edgeIntersection2D(v2p, v0p-v2p, edge0p, edge1p-edge0p)) {
- return true;
- }
- return false;
- }
- else {
- return false;
- }
- }
- else {
- // edge and triangle are NOT parallel
- // compute intersection point
- T distance = (normal * (v0 -edge0)) / normalDotDir;
- if (distance < 0.0 || distance > edgeLen) { return false; }
- Vector3<T> point = edge0 + distance * dir;
- // find and remove dominant axis
- Vector3<T> absolute = normal.abs();
- int dom = 0;
- if (absolute[1] > absolute[dom]) { dom = 1; }
- if (absolute[2] > absolute[dom]) { dom = 2; }
- int axis1 = (dom + 1) % 3;
- int axis2 = (dom + 2) % 3;
- // compute barycentric coordinates
- Vector2<T> b = Vector2<T>(v1[axis1] - v0[axis1], v1[axis2] - v0[axis2]);
- Vector2<T> c = Vector2<T>(v2[axis1] - v0[axis1], v2[axis2] - v0[axis2]);
- Vector2<T> p = Vector2<T>(point[axis1] - v0[axis1], point[axis2] - v0[axis2]);
- T u = (p.y * c.x - p.x * c.y) / (b.y * c.x - b.x * c.y);
- if (u < 0.0) { return false; }
- T v = (p.y * b.x - p.x * b.y) / (c.y * b.x - c.x * b.y);
- if (v < 0.0 || u+v > 1.0) { return false; }
- return true;
- }
- }
- template<class T> inline bool edgeIntersection2D(const Vector2<T>& p1, const Vector2<T>& r1,
- const Vector2<T>& p2, const Vector2<T>& r2 )
- {
- T lenr1 = r1.length();
- T lenr2 = r2.length();
- Vector2<T> nr1(r1.x/lenr1,r1.y/lenr1);
- Vector2<T> nr2(r2.x/lenr2,r2.y/lenr2);
- Vector2<T> b = p2 - p1;
- T det = nr2.x*nr1.y - nr1.x*nr2.y;
- if (-EPSILON < det && det < EPSILON) {
- // edges parallel
- Vector2<T> br2 = b + r2;
- Vector2<T> normal(-nr1.y, nr1.x);
- if (abs(b * normal) > EPSILON && abs(br2 * normal) > EPSILON) {
- return false;
- }
- T a = b*nr1;
- T b = br2*nr1;
- return min(a,b) <= lenr1+EPSILON && max(a,b) >= -EPSILON;
- }
- else {
- T s = (nr2.x*b.y - b.x*nr2.y) / det;
- T t = (nr1.x*b.y - b.x*nr1.y) / det;
- return s > -EPSILON && s < lenr1+EPSILON && t > -EPSILON && t <= lenr2+EPSILON;
- }
- }
- /** checks if point P is inside triangle ABC. */
- template<class T> inline bool pointInTriangle2D(const Vector2<T>& A, const Vector2<T>& B,
- const Vector2<T>& C, const Vector2<T>& P)
- {
- Vector2<T> AB = B - A;
- Vector2<T> BC = C - B;
- Vector2<T> CA = A - C;
- bool inside = (AB % CA)*(AB % (P-A)) <= 0;
- inside &= (BC % AB)*(BC % (P-B)) <= 0;
- inside &= (CA % BC)*(CA % (P-C)) <= 0;
- return inside;
- }