added dev/MinGfx/

1 files changed, 313 insertions, 0 deletions

diff --git a/dev/MinGfx/src/ray.cc b/dev/MinGfx/src/ray.cc
new file mode 100644
index 0000000..5d406ef
--- /dev/null
+++ b/dev/MinGfx/src/ray.cc
@@ -0,0 +1,313 @@
+/*
+ Copyright (c) 2017,2018 Regents of the University of Minnesota.
+ All Rights Reserved.
+ See corresponding header file for details.
+ */
+
+#include "ray.h"
+
+namespace mingfx {
+    
+    Ray::Ray() : p_(Point3::Origin()), d_(-Vector3::UnitZ()) {
+    }
+    
+    Ray::Ray(const Point3 &origin, const Vector3 &direction) : p_(origin), d_(direction) {
+    }
+    
+    Ray::~Ray() {
+    }
+
+    
+    bool Ray::operator==(const Ray& other) const {
+        return (p_ == other.origin()) && (d_ == other.direction());
+    }
+
+    
+    bool Ray::operator!=(const Ray& other) const {
+        return (p_ != other.origin()) || (d_ != other.direction());
+    }
+
+    
+    float Ray::Length() const {
+        return d_.Length();
+    }
+
+    
+    Point3 Ray::origin() const {
+        return p_;
+    }
+
+    
+    Vector3 Ray::direction() const {
+        return d_;
+    }
+
+    
+    
+    bool Ray::IntersectPlane(const Point3 &planePt, const Vector3 &planeNormal,
+                             float *iTime, Point3 *iPoint) const
+    {
+        float dot = planeNormal.Dot(d_);
+        
+        // return false if we would hit the back face of the plane
+        if (dot > 0.0) {
+            return false;
+        }
+        
+        // return false if the ray and plane are parallel
+        if (std::fabs(dot) < MINGFX_MATH_EPSILON) {
+            return false;
+        }
+        
+        float denom = planeNormal.Dot(d_);
+        if (std::abs(denom) > MINGFX_MATH_EPSILON) {
+            *iTime = (planePt - p_).Dot(planeNormal) / denom;
+            if (*iTime >= 0) {
+                *iPoint = p_ + (*iTime)*d_;
+                return true;
+            }
+        }
+        return false;
+    }
+    
+    
+    bool Ray::IntersectTriangle(const Point3 &vertex0, const Point3 &vertex1, const Point3 &vertex2,
+                                float *iTime, Point3 *iPoint) const
+    {
+        // Implementation of the Möller–Trumbore intersection algorithm
+        // https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
+
+        Vector3 edge1, edge2, h, s, q;
+        float a,f,u,v;
+        edge1 = vertex1 - vertex0;
+        edge2 = vertex2 - vertex0;
+        h = d_.Cross(edge2);
+        a = edge1.Dot(h);
+        if (a > -MINGFX_MATH_EPSILON && a < MINGFX_MATH_EPSILON)
+            return false;
+        f = 1.f/a;
+        s = p_ - vertex0;
+        u = f * (s.Dot(h));
+        if (u < 0.0 || u > 1.f)
+            return false;
+        q = s.Cross(edge1);
+        v = f * d_.Dot(q);
+        if ((v < 0.0) || (u + v > 1.0f))
+            return false;
+        // At this stage we can compute t to find out where the intersection point is on the line.
+        *iTime = f * edge2.Dot(q);
+        if (*iTime > MINGFX_MATH_EPSILON) { // ray intersection
+            *iPoint = p_ + d_ * (*iTime);
+            return true;
+        }
+        else // This means that there is a line intersection but not a ray intersection.
+            return false;
+        
+        
+        /***
+        // A basic implementation
+         
+        // find the point of intersection of the ray with the plane of the triangle
+        Vector3 AB = B - A;
+        Vector3 AC = C - A;
+        Vector3 cross = AB.Cross(AC);
+        Vector3 N = cross.ToUnit();
+        if (!IntersectPlane(A, N, iTime, iPoint)) {
+            return false;
+        }
+        
+        // check to see if iPoint lies within the triangle
+        Vector3 edge1 = B - A;
+        Vector3 v1 = *iPoint - A;
+        Vector3 check1 = edge1.Cross(v1);
+        if (N.Dot(check1) < 0.0) {
+            return false;
+        }
+
+        Vector3 edge2 = C - B;
+        Vector3 v2 = *iPoint - B;
+        Vector3 check2 = edge2.Cross(v2);
+        if (N.Dot(check2) < 0.0) {
+            return false;
+        }
+
+        Vector3 edge3 = A - C;
+        Vector3 v3 = *iPoint - C;
+        Vector3 check3 = edge3.Cross(v3);
+        if (N.Dot(check3) < 0.0) {
+            return false;
+        }
+
+        return true;
+         ***/
+    }
+
+    
+    bool Ray::IntersectQuad(const Point3 &A, const Point3 &B, const Point3 &C, const Point3 &D,
+                            float *iTime, Point3 *iPoint) const
+    {
+        if (IntersectTriangle(A, B, C, iTime, iPoint)) {
+            return true;
+        }
+        else if (IntersectTriangle(A, C, D, iTime, iPoint)) {
+            return true;
+        }
+        else {
+            return false;
+        }
+    }
+
+    
+    bool Ray::IntersectSphere(const Point3 &center, float radius,
+                              float *iTime, Point3 *iPoint) const
+    {
+        Point3 P = p_ + (Point3::Origin() - center);
+        Vector3 D = d_;
+        
+        // A = (Dx^2 + Dy^2 + Dz^2)
+        const double A = ((double)D[0]*D[0] + (double)D[1]*D[1] + (double)D[2]*D[2]);
+        
+        // B = (Px * Dx + Py * Dy + Pz * Dz)
+        const double B = ((double)P[0]*D[0] + (double)P[1]*D[1] + (double)P[2]*D[2]);
+        
+        // C = (Px^2 + Py^2 + Pz^2 - (sphere radius)^2)
+        const double C = ((double)P[0]*P[0] + (double)P[1]*P[1] + (double)P[2]*P[2] - (double)radius*radius);
+        
+        // Discriminant of quadratic = B^2 - A * C
+        double discriminant =  B*B - A*C;
+        
+        if (discriminant < 0.0) {
+            return false;
+        }
+        else {
+            double discRoot = sqrt(discriminant);
+            double t1 = (-B - discRoot) / A;
+            double t2 = (-B + discRoot) / A;
+            bool hit1 = false;
+            bool hit2 = false;
+            if (t1 > MINGFX_MATH_EPSILON) {
+                hit1 = true;
+                *iTime = (float)t1;
+            }
+            if (t2 > MINGFX_MATH_EPSILON) {
+                hit2 = true;
+                *iTime = (float)t2;
+            }
+            if ((!hit1) && (!hit2)) {
+                return false;
+            }
+            if ((hit1) && (hit2)) {
+                if (t1 < t2) {
+                    *iTime = (float)t1;
+                }
+            }
+            
+            *iPoint = p_ + (*iTime)*d_;
+            return true;
+        }
+    }
+
+    
+    bool Ray::IntersectMesh(const Mesh &mesh,
+                            float *iTime, Point3 *iPoint, int *iTriangleID) const
+    {
+        *iTime = -1.0;
+        for (int i=0; i<mesh.num_triangles(); i++) {
+            Point3 p;
+            float t;
+            std::vector<unsigned int> indices = mesh.read_triangle_indices_data(i);
+            if (IntersectTriangle(mesh.read_vertex_data(indices[0]), mesh.read_vertex_data(indices[1]), mesh.read_vertex_data(indices[2]), &t, &p)) {
+                if ((*iTime < 0.0) || (t < *iTime)) {
+                    *iPoint = p;
+                    *iTime = t;
+                    *iTriangleID = i;
+                }
+            }
+        }
+        return (*iTime > 0.0);
+    }
+
+    bool Ray::FastIntersectMesh(Mesh *mesh, float *iTime,
+                                Point3 *iPoint, int *iTriangleID) const
+    {
+        std::vector<int> tri_ids = mesh->bvh_ptr()->IntersectAndReturnUserData(*this);
+        if (tri_ids.size()) {
+            *iTime = -1.0;
+            for (int i=0; i<tri_ids.size(); i++) {
+                Point3 p;
+                float t;
+                std::vector<unsigned int> indices = mesh->read_triangle_indices_data(tri_ids[i]);
+                if (IntersectTriangle(mesh->read_vertex_data(indices[0]), mesh->read_vertex_data(indices[1]), mesh->read_vertex_data(indices[2]), &t, &p)) {
+                    if ((*iTime < 0.0) || (t < *iTime)) {
+                        *iPoint = p;
+                        *iTime = t;
+                        *iTriangleID = i;
+                    }
+                }
+            }
+            return (*iTime > 0.0);
+        }
+        else {
+            return false;
+        }
+    }
+    
+    
+    bool Ray::IntersectAABB(const AABB &box, float *iTime) const {
+        // https://gamedev.stackexchange.com/questions/18436/most-efficient-aabb-vs-ray-collision-algorithms
+        
+        Point3 origin = p_;
+        
+        Vector3 invdir = d_;
+        invdir[0] = 1.0f / invdir[0];
+        invdir[1] = 1.0f / invdir[1];
+        invdir[2] = 1.0f / invdir[2];
+        
+        float t1 = (box.min()[0] - origin[0])*invdir[0];
+        float t2 = (box.max()[0] - origin[0])*invdir[0];
+        float t3 = (box.min()[1] - origin[1])*invdir[1];
+        float t4 = (box.max()[1] - origin[1])*invdir[1];
+        float t5 = (box.min()[2] - origin[2])*invdir[2];
+        float t6 = (box.max()[2] - origin[2])*invdir[2];
+        
+        float tmin = std::max(std::max(std::min(t1, t2), std::min(t3, t4)), std::min(t5, t6));
+        float tmax = std::min(std::min(std::max(t1, t2), std::max(t3, t4)), std::max(t5, t6));
+        
+        // if tmax < 0, ray (line) is intersecting AABB, but the whole AABB is behind us
+        if (tmax < 0) {
+            *iTime = tmax;
+            return false;
+        }
+        
+        // if tmin > tmax, ray doesn't intersect AABB
+        if (tmin > tmax) {
+            *iTime = tmax;
+            return false;
+        }
+        
+        *iTime = tmin;
+        return true;
+    }
+    
+    
+    void Ray::set(Point3 newOrigin, Vector3 newDir) {
+        p_ = newOrigin;
+        d_ = newDir;
+    }
+
+    std::ostream & operator<< ( std::ostream &os, const Ray &r) {
+        return os << r.origin() << " " << r.direction();
+    }
+    
+    std::istream & operator>> ( std::istream &is, Ray &r) {
+        // format:  (x, y, z)
+        char dummy;
+        Point3 p;
+        Vector3 d;
+        is >> p >> dummy >> d;
+        r.set(p, d);
+        return is;
+    }
+    
+    
+} // end namespace