X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=simgear%2Fmath%2FSGIntersect.hxx;h=46bdb427d0e027c388be42a3997858eeed1ec4ac;hb=914d3e6a2b323cf9f186cbef2aef7865ea07b309;hp=4d87d68844c19af03a4a1dcd707f8740a314653c;hpb=26cb8ec4f1db18084e4dcf1c1f9b5baee92d4e7e;p=simgear.git

diff --git a/simgear/math/SGIntersect.hxx b/simgear/math/SGIntersect.hxx
index 4d87d688..46bdb427 100644
--- a/simgear/math/SGIntersect.hxx
+++ b/simgear/math/SGIntersect.hxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+// Copyright (C) 2006-2009  Mathias Froehlich - Mathias.Froehlich@web.de
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Library General Public
@@ -18,42 +18,44 @@
 #ifndef SGIntersect_HXX
 #define SGIntersect_HXX
 
+#include <algorithm>
+
 template<typename T>
 inline bool
-intersects(const SGBox<T>& box, const SGSphere<T>& sphere)
+intersects(const SGSphere<T>& s1, const SGSphere<T>& s2)
 {
-  if (sphere.empty())
-    return false;
-  // Is more or less trivially included in the next tests
-  // if (box.empty())
-  //   return false;
-
-  if (sphere.getCenter().x() < box.getMin().x() - sphere.getRadius())
-    return false;
-  if (sphere.getCenter().y() < box.getMin().y() - sphere.getRadius())
+  if (s1.empty())
     return false;
-  if (sphere.getCenter().z() < box.getMin().z() - sphere.getRadius())
+  if (s2.empty())
     return false;
 
-  if (box.getMax().x() + sphere.getRadius() < sphere.getCenter().x())
-    return false;
-  if (box.getMax().y() + sphere.getRadius() < sphere.getCenter().y())
+  T dist = s1.getRadius() + s2.getRadius();
+  return distSqr(s1.getCenter(), s2.getCenter()) <= dist*dist;
+}
+
+
+template<typename T1, typename T2>
+inline bool
+intersects(const SGBox<T1>& box, const SGSphere<T2>& sphere)
+{
+  if (sphere.empty())
     return false;
-  if (box.getMax().z() + sphere.getRadius() < sphere.getCenter().z())
+  if (box.empty())
     return false;
 
-  return true;
+  SGVec3<T1> closest = box.getClosestPoint(sphere.getCenter());
+  return distSqr(closest, SGVec3<T1>(sphere.getCenter())) <= sphere.getRadius2();
 }
 // make it symmetric
-template<typename T>
+template<typename T1, typename T2>
 inline bool
-intersects(const SGSphere<T>& sphere, const SGBox<T>& box)
+intersects(const SGSphere<T1>& sphere, const SGBox<T2>& box)
 { return intersects(box, sphere); }
 
 
-template<typename T>
+template<typename T1, typename T2>
 inline bool
-intersects(const SGVec3<T>& v, const SGBox<T>& box)
+intersects(const SGVec3<T1>& v, const SGBox<T2>& box)
 {
   if (v[0] < box.getMin()[0])
     return false;
@@ -69,9 +71,9 @@ intersects(const SGVec3<T>& v, const SGBox<T>& box)
     return false;
   return true;
 }
-template<typename T>
+template<typename T1, typename T2>
 inline bool
-intersects(const SGBox<T>& box, const SGVec3<T>& v)
+intersects(const SGBox<T1>& box, const SGVec3<T2>& v)
 { return intersects(v, box); }
 
 
@@ -251,11 +253,11 @@ intersects(SGVec3<T>& dst, const SGLineSegment<T>& lineSegment, const SGPlane<T>
 
   // The negative numerator for the \alpha expression
   T num = plane.getPositiveDist();
-  num -= dot(plane.getNormal(), lineSegment.getOrigin());
+  num -= dot(plane.getNormal(), lineSegment.getStart());
   
   // If the numerator is zero, we have the lines origin included in the plane
   if (fabs(num) <= SGLimits<T>::min()) {
-    dst = lineSegment.getOrigin();
+    dst = lineSegment.getStart();
     return true;
   }
 
@@ -277,7 +279,7 @@ intersects(SGVec3<T>& dst, const SGLineSegment<T>& lineSegment, const SGPlane<T>
   if (1 < alpha)
     return false;
 
-  dst = lineSegment.getOrigin() + alpha*lineSegment.getDirection();
+  dst = lineSegment.getStart() + alpha*lineSegment.getDirection();
   return true;
 }
 // make it symmetric
@@ -287,6 +289,38 @@ intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGLineSegment<T>& line
 { return intersects(dst, lineSegment, plane); }
 
 
+// Distance of a line segment to a point
+template<typename T>
+inline T
+distSqr(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
+{
+  SGVec3<T> ps = p - lineSeg.getStart();
+
+  T psdotdir = dot(ps, lineSeg.getDirection());
+  if (psdotdir <= 0)
+    return dot(ps, ps);
+
+  SGVec3<T> pe = p - lineSeg.getEnd();
+  if (0 <= dot(pe, lineSeg.getDirection()))
+    return dot(pe, pe);
+ 
+  return dot(ps, ps) - psdotdir*psdotdir/dot(lineSeg.getDirection(), lineSeg.getDirection());
+}
+// make it symmetric
+template<typename T>
+inline T
+distSqr(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
+{ return distSqr(lineSeg, p); }
+// with sqrt
+template<typename T>
+inline T
+dist(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
+{ return sqrt(distSqr(lineSeg, p)); }
+template<typename T>
+inline T
+dist(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
+{ return sqrt(distSqr(lineSeg, p)); }
+
 template<typename T>
 inline bool
 intersects(const SGRay<T>& ray, const SGSphere<T>& sphere)
@@ -497,13 +531,427 @@ template<typename T>
 inline bool
 intersects(const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
 {
-  // FIXME: for now just wrap the othr method. When that has prooven
+  // FIXME: for now just wrap the other method. When that has prooven
   // well optimized, implement that special case
   SGVec3<T> dummy;
   return intersects(dummy, tri, lineSegment, eps);
 }
 
 
+template<typename T>
+inline SGVec3<T>
+closestPoint(const SGTriangle<T>& tri, const SGVec3<T>& p)
+{
+  // This method minimizes the distance function Q(u, v) = || p - x ||
+  // where x is a point in the trialgle given by the vertices v_i
+  //  x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
+  // The theoretical analysis is somehow too long for a comment.
+  // May be it is sufficient to see that this code passes all the tests.
+
+  SGVec3<T> off = tri.getBaseVertex() - p;
+  T a = dot(tri.getEdge(0), tri.getEdge(0));
+  T b = dot(tri.getEdge(0), tri.getEdge(1));
+  T c = dot(tri.getEdge(1), tri.getEdge(1));
+  T d = dot(tri.getEdge(0), off);
+  T e = dot(tri.getEdge(1), off);
+  
+  T det = a*c - b*b;
+
+  T u = b*e - c*d;
+  T v = b*d - a*e;
+
+/*
+  // Regions
+  // \2|
+  //  \|
+  //   |\ 
+  // 3 |0\ 1
+  //----------
+  // 4 | 5 \ 6
+*/
+
+  if (u + v <= det) {
+    if (u < 0) {
+      if (v < 0) {
+        // region 4
+        if (d < 0) {
+          if (a <= -d) {
+            // u = 1;
+            // v = 0;
+            return tri.getBaseVertex() + tri.getEdge(0);
+          } else {
+            u = -d/a;
+            // v = 0;
+            return tri.getBaseVertex() + u*tri.getEdge(0);
+          }
+        } else {
+          if (0 < e) {
+            // u = 0;
+            // v = 0;
+            return tri.getBaseVertex();
+          } else if (c <= -e) {
+            // u = 0;
+            // v = 1;
+            return tri.getBaseVertex() + tri.getEdge(1);
+          } else {
+            // u = 0;
+            v = -e/c;
+            return tri.getBaseVertex() + v*tri.getEdge(1);
+          }
+        }
+      } else {
+        // region 3
+        if (0 <= e) {
+          // u = 0;
+          // v = 0;
+          return tri.getBaseVertex();
+        } else if (c <= -e) {
+          // u = 0;
+          // v = 1;
+          return tri.getBaseVertex() + tri.getEdge(1);
+        } else {
+          // u = 0;
+          v = -e/c;
+          return tri.getBaseVertex() + v*tri.getEdge(1);
+        }
+      }
+    } else if (v < 0) {
+      // region 5
+      if (0 <= d) {
+        // u = 0;
+        // v = 0;
+        return tri.getBaseVertex();
+      } else if (a <= -d) {
+        // u = 1;
+        // v = 0;
+        return tri.getBaseVertex() + tri.getEdge(0);
+      } else {
+        u = -d/a;
+        // v = 0;
+        return tri.getBaseVertex() + u*tri.getEdge(0);
+      }
+    } else {
+      // region 0
+      if (det <= SGLimits<T>::min()) {
+        u = 0;
+        v = 0;
+        return tri.getBaseVertex();
+      } else {
+        T invDet = 1/det;
+        u *= invDet;
+        v *= invDet;
+        return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
+      }
+    }
+  } else {
+    if (u < 0) {
+      // region 2
+      T tmp0 = b + d;
+      T tmp1 = c + e;
+      if (tmp0 < tmp1) {
+        T numer = tmp1 - tmp0;
+        T denom = a - 2*b + c;
+        if (denom <= numer) {
+          // u = 1;
+          // v = 0;
+          return tri.getBaseVertex() + tri.getEdge(0);
+        } else {
+          u = numer/denom;
+          v = 1 - u;
+          return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
+        }
+      } else {
+        if (tmp1 <= 0) {
+          // u = 0;
+          // v = 1;
+          return tri.getBaseVertex() + tri.getEdge(1);
+        } else if (0 <= e) {
+          // u = 0;
+          // v = 0;
+          return tri.getBaseVertex();
+        } else {
+          // u = 0;
+          v = -e/c;
+          return tri.getBaseVertex() + v*tri.getEdge(1);
+        }
+      }
+    } else if (v < 0) {
+      // region 6
+      T tmp0 = b + e;
+      T tmp1 = a + d;
+      if (tmp0 < tmp1) {
+        T numer = tmp1 - tmp0;
+        T denom = a - 2*b + c;
+        if (denom <= numer) {
+          // u = 0;
+          // v = 1;
+          return tri.getBaseVertex() + tri.getEdge(1);
+        } else {
+          v = numer/denom;
+          u = 1 - v;
+          return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
+        }
+      } else {
+        if (tmp1 < 0) {
+          // u = 1;
+          // v = 0;
+          return tri.getBaseVertex() + tri.getEdge(0);
+        } else if (0 <= d) {
+          // u = 0;
+          // v = 0;
+          return tri.getBaseVertex();
+        } else {
+          u = -d/a;
+          // v = 0;
+          return tri.getBaseVertex() + u*tri.getEdge(0);
+        }
+      }
+    } else {
+      // region 1
+      T numer = c + e - b - d;
+      if (numer <= 0) {
+        // u = 0;
+        // v = 1;
+        return tri.getVertex(2);
+      } else {
+        T denom = a - 2*b + c;
+        if (denom <= numer) {
+          // u = 1;
+          // v = 0;
+          return tri.getBaseVertex() + tri.getEdge(0);
+        } else {
+          u = numer/denom;
+          v = 1 - u;
+          return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
+        }
+      }
+    }
+  }
+}
+template<typename T>
+inline SGVec3<T>
+closestPoint(const SGVec3<T>& p, const SGTriangle<T>& tri)
+{ return closestPoint(tri, p); }
+
+template<typename T, typename T2>
+inline bool
+intersects(const SGTriangle<T>& tri, const SGSphere<T2>& sphere)
+{
+  // This method minimizes the distance function Q(u, v) = || p - x ||
+  // where x is a point in the trialgle given by the vertices v_i
+  //  x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
+  // The theoretical analysis is somehow too long for a comment.
+  // May be it is sufficient to see that this code passes all the tests.
+
+  SGVec3<T> off = tri.getBaseVertex() - SGVec3<T>(sphere.getCenter());
+  T baseDist2 = dot(off, off);
+  T a = dot(tri.getEdge(0), tri.getEdge(0));
+  T b = dot(tri.getEdge(0), tri.getEdge(1));
+  T c = dot(tri.getEdge(1), tri.getEdge(1));
+  T d = dot(tri.getEdge(0), off);
+  T e = dot(tri.getEdge(1), off);
+  
+  T det = a*c - b*b;
+
+  T u = b*e - c*d;
+  T v = b*d - a*e;
+
+/*
+  // Regions
+  // \2|
+  //  \|
+  //   |\ 
+  // 3 |0\ 1
+  //----------
+  // 4 | 5 \ 6
+*/
+
+  if (u + v <= det) {
+    if (u < 0) {
+      if (v < 0) {
+        // region 4
+        if (d < 0) {
+          if (a <= -d) {
+            // u = 1;
+            // v = 0;
+            T sqrDist = a + 2*d + baseDist2;
+            return sqrDist <= sphere.getRadius2();
+          } else {
+            u = -d/a;
+            // v = 0;
+            T sqrDist = d*u + baseDist2;
+            return sqrDist <= sphere.getRadius2();
+          }
+        } else {
+          if (0 < e) {
+            // u = 0;
+            // v = 0;
+            return baseDist2 <= sphere.getRadius2();
+          } else if (c <= -e) {
+            // u = 0;
+            // v = 1;
+            T sqrDist = c + 2*e + baseDist2;
+            return sqrDist <= sphere.getRadius2();
+          } else {
+            // u = 0;
+            v = -e/c;
+            T sqrDist = e*v + baseDist2;
+            return sqrDist <= sphere.getRadius2();
+          }
+        }
+      } else {
+        // region 3
+        if (0 <= e) {
+          // u = 0;
+          // v = 0;
+          return baseDist2 <= sphere.getRadius2();
+        } else if (c <= -e) {
+          // u = 0;
+          // v = 1;
+          T sqrDist = c + 2*e + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        } else {
+          // u = 0;
+          v = -e/c;
+          T sqrDist = e*v + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        }
+      }
+    } else if (v < 0) {
+      // region 5
+      if (0 <= d) {
+        // u = 0;
+        // v = 0;
+        return baseDist2 <= sphere.getRadius2();
+      } else if (a <= -d) {
+        // u = 1;
+        // v = 0;
+        T sqrDist = a + 2*d + baseDist2;
+        return sqrDist <= sphere.getRadius2();
+      } else {
+        u = -d/a;
+        // v = 0;
+        T sqrDist = d*u + baseDist2;
+        return sqrDist <= sphere.getRadius2();
+      }
+    } else {
+      // region 0
+      if (det <= SGLimits<T>::min()) {
+        // sqrDist = baseDist2;
+        u = 0;
+        v = 0;
+        return baseDist2 <= sphere.getRadius2();
+      } else {
+        T invDet = 1/det;
+        u *= invDet;
+        v *= invDet;
+        T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
+        return sqrDist <= sphere.getRadius2();
+      }
+    }
+  } else {
+    if (u < 0) {
+      // region 2
+      T tmp0 = b + d;
+      T tmp1 = c + e;
+      if (tmp0 < tmp1) {
+        T numer = tmp1 - tmp0;
+        T denom = a - 2*b + c;
+        if (denom <= numer) {
+          // u = 1;
+          // v = 0;
+          T sqrDist = a + 2*d + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        } else {
+          u = numer/denom;
+          v = 1 - u;
+          T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        }
+      } else {
+        if (tmp1 <= 0) {
+          // u = 0;
+          // v = 1;
+          T sqrDist = c + 2*e + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        } else if (0 <= e) {
+          // u = 0;
+          // v = 0;
+          return baseDist2 <= sphere.getRadius2();
+        } else {
+          // u = 0;
+          v = -e/c;
+          T sqrDist = e*v + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        }
+      }
+    } else if (v < 0) {
+      // region 6
+      T tmp0 = b + e;
+      T tmp1 = a + d;
+      if (tmp0 < tmp1) {
+        T numer = tmp1 - tmp0;
+        T denom = a - 2*b + c;
+        if (denom <= numer) {
+          // u = 0;
+          // v = 1;
+          T sqrDist = c + 2*e + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        } else {
+          v = numer/denom;
+          u = 1 - v;
+          T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e)+baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        }
+      } else {
+        if (tmp1 < 0) {
+          // u = 1;
+          // v = 0;
+          T sqrDist = a + 2*d + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        } else if (0 <= d) {
+          // sqrDist = baseDist2;
+          // u = 0;
+          // v = 0;
+          return baseDist2 <= sphere.getRadius2();
+        } else {
+          u = -d/a;
+          // v = 0;
+          T sqrDist = d*u + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        }
+      }
+    } else {
+      // region 1
+      T numer = c + e - b - d;
+      if (numer <= 0) {
+        // u = 0;
+        // v = 1;
+        T sqrDist = c + 2*e + baseDist2;
+        return sqrDist <= sphere.getRadius2();
+      } else {
+        T denom = a - 2*b + c;
+        if (denom <= numer) {
+          // u = 1;
+          // v = 0;
+          T sqrDist = a + 2*d + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        } else {
+          u = numer/denom;
+          v = 1 - u;
+          T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
+          return sqrDist <= sphere.getRadius2();
+        }
+      }
+    }
+  }
+}
+template<typename T1, typename T2>
+inline bool
+intersects(const SGSphere<T1>& sphere, const SGTriangle<T2>& tri)
+{ return intersects(tri, sphere); }
+
+
 template<typename T>
 inline bool
 intersects(const SGVec3<T>& v, const SGSphere<T>& sphere)
@@ -525,9 +973,9 @@ intersects(const SGBox<T>& box, const SGLineSegment<T>& lineSegment)
   // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
 
   SGVec3<T> c = lineSegment.getCenter() - box.getCenter();
-  SGVec3<T> w = 0.5*lineSegment.getDirection();
+  SGVec3<T> w = T(0.5)*lineSegment.getDirection();
   SGVec3<T> v(fabs(w.x()), fabs(w.y()), fabs(w.z()));
-  SGVec3<T> h = 0.5*box.getSize();
+  SGVec3<T> h = T(0.5)*box.getSize();
 
   if (fabs(c[0]) > v[0] + h[0])
     return false;
@@ -592,4 +1040,26 @@ inline bool
 intersects(const SGRay<T>& ray, const SGBox<T>& box)
 { return intersects(box, ray); }
 
+template<typename T1, typename T2>
+inline bool
+intersects(const SGBox<T1>& box1, const SGBox<T2>& box2)
+{
+  if (box2.getMax()[0] < box1.getMin()[0])
+    return false;
+  if (box1.getMax()[0] < box2.getMin()[0])
+    return false;
+
+  if (box2.getMax()[1] < box1.getMin()[1])
+    return false;
+  if (box1.getMax()[1] < box2.getMin()[1])
+    return false;
+
+  if (box2.getMax()[2] < box1.getMin()[2])
+    return false;
+  if (box1.getMax()[2] < box2.getMin()[2])
+    return false;
+
+  return true;
+}
+
 #endif