libs/mathlib/ray.c

   1 /*
   2 Copyright (C) 2001-2006, William Joseph.
   3 All Rights Reserved.
   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 #include <float.h>
  24
  25 vec3_t identity = { 0,0,0 };
  26
  27 void ray_construct_for_vec3(ray_t *ray, const vec3_t origin, const vec3_t direction)
  28 {
  29   VectorCopy(origin, ray->origin);
  30   VectorCopy(direction, ray->direction);
  31 }
  32
  33 void ray_transform(ray_t *ray, const m4x4_t matrix)
  34 {
  35   m4x4_transform_point(matrix, ray->origin);
  36   m4x4_transform_normal(matrix, ray->direction);
  37 }
  38
  39 vec_t ray_intersect_point(const ray_t *ray, const vec3_t point, vec_t epsilon, vec_t divergence)
  40 {
  41   vec3_t displacement;
  42   vec_t depth;
  43
  44   // calc displacement of test point from ray origin
  45   VectorSubtract(point, ray->origin, displacement);
  46   // calc length of displacement vector along ray direction
  47         depth = DotProduct(displacement, ray->direction);
  48   if(depth < 0.0f) return (vec_t)FLT_MAX;
  49   // calc position of closest point on ray to test point
  50   VectorMA (ray->origin, depth, ray->direction, displacement);
  51   // calc displacement of test point from closest point
  52         VectorSubtract(point, displacement, displacement);
  53   // calc length of displacement, subtract depth-dependant epsilon
  54   if (VectorLength(displacement) - (epsilon + (depth * divergence)) > 0.0f) return (vec_t)FLT_MAX;
  55   return depth;
  56 }
  57
  58 // Tomas Moller and Ben Trumbore. Fast, minimum storage ray-triangle intersection. Journal of graphics tools, 2(1):21-28, 1997
  59
  60 #define EPSILON 0.000001
  61
  62 vec_t ray_intersect_triangle(const ray_t *ray, qboolean bCullBack, const vec3_t vert0, const vec3_t vert1, const vec3_t vert2)
  63 {
  64   float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
  65   float det,inv_det;
  66   float u, v;
  67   vec_t depth = (vec_t)FLT_MAX;
  68
  69   /* find vectors for two edges sharing vert0 */
  70   VectorSubtract(vert1, vert0, edge1);
  71   VectorSubtract(vert2, vert0, edge2);
  72
  73   /* begin calculating determinant - also used to calculate U parameter */
  74   CrossProduct(ray->direction, edge2, pvec);
  75
  76   /* if determinant is near zero, ray lies in plane of triangle */
  77   det = DotProduct(edge1, pvec);
  78
  79   if (bCullBack == qtrue)
  80   {
  81     if (det < EPSILON)
  82       return depth;
  83
  84     // calculate distance from vert0 to ray origin
  85     VectorSubtract(ray->origin, vert0, tvec);
  86
  87     // calculate U parameter and test bounds
  88     u = DotProduct(tvec, pvec);
  89     if (u < 0.0 || u > det)
  90       return depth;
  91
  92     // prepare to test V parameter
  93     CrossProduct(tvec, edge1, qvec);
  94
  95     // calculate V parameter and test bounds
  96     v = DotProduct(ray->direction, qvec);
  97     if (v < 0.0 || u + v > det)
  98       return depth;
  99
 100     // calculate t, scale parameters, ray intersects triangle
 101     depth = DotProduct(edge2, qvec);
 102     inv_det = 1.0f / det;
 103     depth *= inv_det;
 104     //u *= inv_det;
 105     //v *= inv_det;
 106   }
 107   else
 108   {
 109     /* the non-culling branch */
 110     if (det > -EPSILON && det < EPSILON)
 111       return depth;
 112     inv_det = 1.0f / det;
 113
 114     /* calculate distance from vert0 to ray origin */
 115     VectorSubtract(ray->origin, vert0, tvec);
 116
 117     /* calculate U parameter and test bounds */
 118     u = DotProduct(tvec, pvec) * inv_det;
 119     if (u < 0.0 || u > 1.0)
 120       return depth;
 121
 122     /* prepare to test V parameter */
 123     CrossProduct(tvec, edge1, qvec);
 124
 125     /* calculate V parameter and test bounds */
 126     v = DotProduct(ray->direction, qvec) * inv_det;
 127     if (v < 0.0 || u + v > 1.0)
 128       return depth;
 129
 130     /* calculate t, ray intersects triangle */
 131     depth = DotProduct(edge2, qvec) * inv_det;
 132   }
 133   return depth;
 134 }
 135
 136 vec_t ray_intersect_plane(const ray_t* ray, const vec3_t normal, vec_t dist)
 137 {
 138   return -(DotProduct(normal, ray->origin) - dist) / DotProduct(ray->direction, normal);
 139 }
 140