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         VectorCopy( origin, ray->origin );
  30         VectorCopy( direction, ray->direction );
  31 }
  32
  33 void ray_transform( ray_t *ray, const m4x4_t matrix ){
  34         m4x4_transform_point( matrix, ray->origin );
  35         m4x4_transform_normal( matrix, ray->direction );
  36 }
  37
  38 vec_t ray_intersect_point( const ray_t *ray, const vec3_t point, vec_t epsilon, vec_t divergence ){
  39         vec3_t displacement;
  40         vec_t depth;
  41
  42         // calc displacement of test point from ray origin
  43         VectorSubtract( point, ray->origin, displacement );
  44         // calc length of displacement vector along ray direction
  45         depth = DotProduct( displacement, ray->direction );
  46         if ( depth < 0.0f ) {
  47                 return (vec_t)VEC_MAX;
  48         }
  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 ) {
  55                 return (vec_t)VEC_MAX;
  56         }
  57         return depth;
  58 }
  59
  60 // Tomas Moller and Ben Trumbore. Fast, minimum storage ray-triangle intersection. Journal of graphics tools, 2(1):21-28, 1997
  61
  62 #define EPSILON 0.000001
  63
  64 vec_t ray_intersect_triangle( const ray_t *ray, qboolean bCullBack, const vec3_t vert0, const vec3_t vert1, const vec3_t vert2 ){
  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                 if ( det < EPSILON ) {
  82                         return depth;
  83                 }
  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
  94                 // prepare to test V parameter
  95                 CrossProduct( tvec, edge1, qvec );
  96
  97                 // calculate V parameter and test bounds
  98                 v = DotProduct( ray->direction, qvec );
  99                 if ( v < 0.0 || u + v > det ) {
 100                         return depth;
 101                 }
 102
 103                 // calculate t, scale parameters, ray intersects triangle
 104                 depth = DotProduct( edge2, qvec );
 105                 inv_det = 1.0f / det;
 106                 depth *= inv_det;
 107                 //u *= inv_det;
 108                 //v *= inv_det;
 109         }
 110         else
 111         {
 112                 /* the non-culling branch */
 113                 if ( det > -EPSILON && det < EPSILON ) {
 114                         return depth;
 115                 }
 116                 inv_det = 1.0f / det;
 117
 118                 /* calculate distance from vert0 to ray origin */
 119                 VectorSubtract( ray->origin, vert0, tvec );
 120
 121                 /* calculate U parameter and test bounds */
 122                 u = DotProduct( tvec, pvec ) * inv_det;
 123                 if ( u < 0.0 || u > 1.0 ) {
 124                         return depth;
 125                 }
 126
 127                 /* prepare to test V parameter */
 128                 CrossProduct( tvec, edge1, qvec );
 129
 130                 /* calculate V parameter and test bounds */
 131                 v = DotProduct( ray->direction, qvec ) * inv_det;
 132                 if ( v < 0.0 || u + v > 1.0 ) {
 133                         return depth;
 134                 }
 135
 136                 /* calculate t, ray intersects triangle */
 137                 depth = DotProduct( edge2, qvec ) * inv_det;
 138         }
 139         return depth;
 140 }