libs/mathlib/ray.c

   1 /*
   2 Copyright (C) 1999-2007 id Software, Inc. and contributors.
   3 For a list of contributors, see the accompanying CONTRIBUTORS file.
   4
   5 This file is part of GtkRadiant.
   6
   7 GtkRadiant is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2 of the License, or
  10 (at your option) any later version.
  11
  12 GtkRadiant is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with GtkRadiant; if not, write to the Free Software
  19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 */
  21
  22 #include "mathlib.h"
  23 /*! for memcpy */
  24 #include <memory.h>
  25
  26 vec3_t identity = { 0,0,0 };
  27
  28 void ray_construct_for_vec3(ray_t *ray, const vec3_t origin, const vec3_t direction)
  29 {
  30   VectorCopy(origin, ray->origin);
  31   VectorCopy(direction, ray->direction);
  32 }
  33
  34 void ray_transform(ray_t *ray, const m4x4_t matrix)
  35 {
  36   m4x4_transform_point(matrix, ray->origin);
  37   m4x4_transform_normal(matrix, ray->direction);
  38 }
  39
  40 vec_t ray_intersect_point(const ray_t *ray, const vec3_t point, vec_t epsilon, vec_t divergence)
  41 {
  42   vec3_t displacement;
  43   vec_t depth;
  44
  45   // calc displacement of test point from ray origin
  46   VectorSubtract(point, ray->origin, displacement);
  47   // calc length of displacement vector along ray direction
  48         depth = DotProduct(displacement, ray->direction);
  49   if(depth < 0.0f) return (vec_t)VEC_MAX;
  50   // calc position of closest point on ray to test point
  51   VectorMA (ray->origin, depth, ray->direction, displacement);
  52   // calc displacement of test point from closest point
  53         VectorSubtract(point, displacement, displacement);
  54   // calc length of displacement, subtract depth-dependant epsilon
  55   if (VectorLength(displacement) - (epsilon + (depth * divergence)) > 0.0f) return (vec_t)VEC_MAX;
  56   return depth;
  57 }
  58
  59 // Tomas Moller and Ben Trumbore. Fast, minimum storage ray-triangle intersection. Journal of graphics tools, 2(1):21-28, 1997
  60
  61 #define EPSILON 0.000001
  62
  63 vec_t ray_intersect_triangle(const ray_t *ray, qboolean bCullBack, const vec3_t vert0, const vec3_t vert1, const vec3_t vert2)
  64 {
  65   float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
  66   float det,inv_det;
  67   float u, v;
  68   vec_t depth = (vec_t)VEC_MAX;
  69
  70   /* find vectors for two edges sharing vert0 */
  71   VectorSubtract(vert1, vert0, edge1);
  72   VectorSubtract(vert2, vert0, edge2);
  73
  74   /* begin calculating determinant - also used to calculate U parameter */
  75   CrossProduct(ray->direction, edge2, pvec);
  76
  77   /* if determinant is near zero, ray lies in plane of triangle */
  78   det = DotProduct(edge1, pvec);
  79
  80   if (bCullBack == qtrue)
  81   {
  82     if (det < EPSILON)
  83       return depth;
  84
  85     // calculate distance from vert0 to ray origin
  86     VectorSubtract(ray->origin, vert0, tvec);
  87
  88     // calculate U parameter and test bounds
  89     u = DotProduct(tvec, pvec);
  90     if (u < 0.0 || u > det)
  91       return depth;
  92
  93     // prepare to test V parameter
  94     CrossProduct(tvec, edge1, qvec);
  95
  96     // calculate V parameter and test bounds
  97     v = DotProduct(ray->direction, qvec);
  98     if (v < 0.0 || u + v > det)
  99       return depth;
 100
 101     // calculate t, scale parameters, ray intersects triangle
 102     depth = DotProduct(edge2, qvec);
 103     inv_det = 1.0f / det;
 104     depth *= inv_det;
 105     //u *= inv_det;
 106     //v *= inv_det;
 107   }
 108   else
 109   {
 110     /* the non-culling branch */
 111     if (det > -EPSILON && det < EPSILON)
 112       return depth;
 113     inv_det = 1.0f / det;
 114
 115     /* calculate distance from vert0 to ray origin */
 116     VectorSubtract(ray->origin, vert0, tvec);
 117
 118     /* calculate U parameter and test bounds */
 119     u = DotProduct(tvec, pvec) * inv_det;
 120     if (u < 0.0 || u > 1.0)
 121       return depth;
 122
 123     /* prepare to test V parameter */
 124     CrossProduct(tvec, edge1, qvec);
 125
 126     /* calculate V parameter and test bounds */
 127     v = DotProduct(ray->direction, qvec) * inv_det;
 128     if (v < 0.0 || u + v > 1.0)
 129       return depth;
 130
 131     /* calculate t, ray intersects triangle */
 132     depth = DotProduct(edge2, qvec) * inv_det;
 133   }
 134   return depth;
 135 }