RayTriangleIntersection.h

#pragma once
 
#include <glm/geometric.hpp>
#include <glm/vec3.hpp>
 
namespace Cogs
{
  namespace Geometry
  {
    inline bool intersect(
      const glm::vec3 & v0,
      const glm::vec3 & v1,
      const glm::vec3 & v2,
      const glm::vec3 & orig,
      const glm::vec3 & dir,
      glm::vec3 & intersection,
      glm::vec3 & barycentric,
      bool & front)
    {
 
      glm::vec3 edge1 = v1 - v0;
      glm::vec3 edge2 = v2 - v0;
 
      glm::vec3 pvec = glm::cross(dir, edge2);
 
      // If determinant is near zero, ray lies in plane of triangle
      float det = glm::dot(edge1, pvec);
 
      if (std::fabs(det) < std::numeric_limits<float>::epsilon()) return false;
 
      // Does ray hit front or back of triangle
      if (det > 0.0) front = true;
      else front = false;
 
      // Create some more intuitive barycentric coordinate names
      float u, v, w;
      float inv_det = 1.0f / det;
 
      // Calculate distance from v0 to ray origin
      glm::vec3 tvec = orig - v0;
 
      // Calculate U parameter and test bounds
      u = glm::dot(tvec, pvec) * inv_det;
 
      if (u < 0.0f || u > 1.0f)
        return false;
 
      // Prepare to test V parameter
      glm::vec3 qvec = glm::cross(tvec, edge1);
 
      // Calculate V parameter and test bounds
      v = glm::dot(dir, qvec) * inv_det;
 
      if (v < 0.0f || u + v > 1.0f)
        return false;
 
      // Third barycentric coordinate
      w = 1.0f - u - v;
 
      // Calculate t and intersection point
      float t = glm::dot(edge2, qvec) * inv_det;
 
      glm::vec3 itmp = orig + t * dir;
 
      intersection = itmp;
 
      // Set the barycentric coordinates before returning
      barycentric[0] = static_cast<float>(w);
      barycentric[1] = static_cast<float>(u);
      barycentric[2] = static_cast<float>(v);
 
      return true;
    }
 
    static double distanceToQuad(
      const double xmin, const double ymin,
      const double xmax, const double ymax,
      const double x, const double y,
      double & cx, double & cy)
    {
      if (x < xmin) {
        if (y < ymin) {
          cx = xmin;
          cy = ymin;
          return (x - xmin)*(x - xmin) + (y - ymin)*(y - ymin);
        } else if (y > ymax) {
          cx = xmin;
          cy = ymax;
          return (x - xmin)*(x - xmin) + (y - ymax)*(y - ymax);
        } else {
          cx = xmin;
          cy = y;
          return (x - xmin)*(x - xmin);
        }
      } else if (x > xmax) {
        if (y < ymin) {
          cx = xmax;
          cy = ymin;
          return (x - xmax)*(x - xmax) + (y - ymin) * (y - ymin);
        } else if (y > ymax) {
          cx = xmax;
          cy = ymax;
          return (x - xmax)*(x - xmax) + (y - ymax)*(y - ymax);
        } else {
          cx = xmax;
          cy = y;
          return (x - xmax)*(x - xmax);
        }
      } else {
        if (y < ymin) {
          cx = x;
          cy = ymin;
          return (y - ymin)*(y - ymin);
        } else if (y > ymax) {
          cx = x;
          cy = ymax;
          return (y - ymax)*(y - ymax);
        } else {
          // inside rectangle
          cx = x;
          cy = y;
          return -1.0;
        }
      }
    }
  }
}