$search
00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include <string.h> 00004 #include <assert.h> 00005 #include <math.h> 00006 00063 #include "raytri.h" 00064 00065 namespace ConvexDecomposition 00066 { 00067 00068 00069 /* a = b - c */ 00070 #define vector(a,b,c) \ 00071 (a)[0] = (b)[0] - (c)[0]; \ 00072 (a)[1] = (b)[1] - (c)[1]; \ 00073 (a)[2] = (b)[2] - (c)[2]; 00074 00075 00076 00077 #define innerProduct(v,q) \ 00078 ((v)[0] * (q)[0] + \ 00079 (v)[1] * (q)[1] + \ 00080 (v)[2] * (q)[2]) 00081 00082 #define crossProduct(a,b,c) \ 00083 (a)[0] = (b)[1] * (c)[2] - (c)[1] * (b)[2]; \ 00084 (a)[1] = (b)[2] * (c)[0] - (c)[2] * (b)[0]; \ 00085 (a)[2] = (b)[0] * (c)[1] - (c)[0] * (b)[1]; 00086 00087 bool rayIntersectsTriangle(const double *p,const double *d,const double *v0,const double *v1,const double *v2,double &t) 00088 { 00089 00090 double e1[3],e2[3],h[3],s[3],q[3]; 00091 double a,f,u,v; 00092 00093 vector(e1,v1,v0); 00094 vector(e2,v2,v0); 00095 crossProduct(h,d,e2); 00096 a = innerProduct(e1,h); 00097 00098 if (a > -0.00001 && a < 0.00001) 00099 return(false); 00100 00101 f = 1/a; 00102 vector(s,p,v0); 00103 u = f * (innerProduct(s,h)); 00104 00105 if (u < 0.0 || u > 1.0) 00106 return(false); 00107 00108 crossProduct(q,s,e1); 00109 v = f * innerProduct(d,q); 00110 if (v < 0.0 || u + v > 1.0) 00111 return(false); 00112 // at this stage we can compute t to find out where 00113 // the intersection point is on the line 00114 t = f * innerProduct(e2,q); 00115 if (t > 0) // ray intersection 00116 return(true); 00117 else // this means that there is a line intersection 00118 // but not a ray intersection 00119 return (false); 00120 } 00121 00122 00123 bool lineIntersectsTriangle(const double *rayStart,const double *rayEnd,const double *p1,const double *p2,const double *p3,double *sect) 00124 { 00125 double dir[3]; 00126 00127 dir[0] = rayEnd[0] - rayStart[0]; 00128 dir[1] = rayEnd[1] - rayStart[1]; 00129 dir[2] = rayEnd[2] - rayStart[2]; 00130 00131 double d = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]); 00132 double r = 1.0f / d; 00133 00134 dir[0]*=r; 00135 dir[1]*=r; 00136 dir[2]*=r; 00137 00138 00139 double t; 00140 00141 bool ret = rayIntersectsTriangle(rayStart, dir, p1, p2, p3, t ); 00142 00143 if ( ret ) 00144 { 00145 if ( t > d ) 00146 { 00147 sect[0] = rayStart[0] + dir[0]*t; 00148 sect[1] = rayStart[1] + dir[1]*t; 00149 sect[2] = rayStart[2] + dir[2]*t; 00150 } 00151 else 00152 { 00153 ret = false; 00154 } 00155 } 00156 00157 return ret; 00158 } 00159 00160 }; // end of namespace