计算机图形学里最重要的是检测基本几何体之间的相交测试（检测两个基本几何体是否相交），以及求出交点在空间中的具体位置。一旦解决了基本的求交问题，那么对于复杂的几何体，如果视其为许多基本几何体的组合，问题也就迎刃而解了。
上国外网站淘了一些源码，自己把它整理了一下，变成了直接就能用的独立函数，贴在这里。原作者网址已经有一定年头了，貌似是Sun公司做CPU那会的年代......所以并不知道里面所做的一些优化对于现代的CPU有没有什么帮助，没准会有副作用。
但至少是一个robust的实现。
一、射线——三角形相交：
1 #pragma once 2 /* Ray-Triangle Intersection Test Routines */ 3 /* Different optimizations of my and Ben Trumbore's */ 4 /* code from journals of graphics tools (JGT) */ 5 /* http://www.acm.org/jgt/ */ 6 /* by Tomas Moller, May 2000 */ 7 8 #include <math.h> 9 10 #define EPSILON 0.000001 11 #define CROSS(dest,v1,v2) \ 12 dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \ 13 dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \ 14 dest[2]=v1[0]*v2[1]-v1[1]*v2[0]; 15 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) 16 #define SUB(dest,v1,v2) \ 17 dest[0]=v1[0]-v2[0]; \ 18 dest[1]=v1[1]-v2[1]; \ 19 dest[2]=v1[2]-v2[2]; 20 21 /* the original jgt code */ 22 int intersect_triangle(double orig[3], double dir[3], 23 double vert0[3], double vert1[3], double vert2[3], 24 double *t, double *u, double *v) 25 { 26 double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3]; 27 double det, inv_det; 28 29 /* find vectors for two edges sharing vert0 */ 30 SUB(edge1, vert1, vert0); 31 SUB(edge2, vert2, vert0); 32 33 /* begin calculating determinant - also used to calculate U parameter */ 34 CROSS(pvec, dir, edge2); 35 36 /* if determinant is near zero, ray lies in plane of triangle */ 37 det = DOT(edge1, pvec); 38 39 if (det > -EPSILON && det < EPSILON) 40 return 0; 41 inv_det = 1.0 / det; 42 43 /* calculate distance from vert0 to ray origin */ 44 SUB(tvec, orig, vert0); 45 46 /* calculate U parameter and test bounds */ 47 *u = DOT(tvec, pvec) * inv_det; 48 if (*u < 0.0 || *u > 1.0) 49 return 0; 50 51 /* prepare to test V parameter */ 52 CROSS(qvec, tvec, edge1); 53 54 /* calculate V parameter and test bounds */ 55 *v = DOT(dir, qvec) * inv_det; 56 if (*v < 0.0 || *u + *v > 1.0) 57 return 0; 58 59 /* calculate t, ray intersects triangle */ 60 *t = DOT(edge2, qvec) * inv_det; 61 62 return 1; 63 } 64 65 66 /* code rewritten to do tests on the sign of the determinant */ 67 /* the division is at the end in the code */ 68 int intersect_triangle1(double orig[3], double dir[3], 69 double vert0[3], double vert1[3], double vert2[3], 70 double *t, double *u, double *v) 71 { 72 double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3]; 73 double det, inv_det; 74 75 /* find vectors for two edges sharing vert0 */ 76 SUB(edge1, vert1, vert0); 77 SUB(edge2, vert2, vert0); 78 79 /* begin calculating determinant - also used to calculate U parameter */ 80 CROSS(pvec, dir, edge2); 81 82 /* if determinant is near zero, ray lies in plane of triangle */ 83 det = DOT(edge1, pvec); 84 85 if (det > EPSILON) 86 { 87 /* calculate distance from vert0 to ray origin */ 88 SUB(tvec, orig, vert0); 89 90 /* calculate U parameter and test bounds */ 91 *u = DOT(tvec, pvec); 92 if (*u < 0.0 || *u > det) 93 return 0; 94 95 /* prepare to test V parameter */ 96 CROSS(qvec, tvec, edge1); 97 98 /* calculate V parameter and test bounds */ 99 *v = DOT(dir, qvec); 100 if (*v < 0.0 || *u + *v > det) 101 return 0; 102 103 } 104 else if (det < -EPSILON) 105 { 106 /* calculate distance from vert0 to ray origin */ 107 SUB(tvec, orig, vert0); 108 109 /* calculate U parameter and test bounds */ 110 *u = DOT(tvec, pvec); 111 /* printf("*u=%f\n",(float)*u); */ 112 /* printf("det=%f\n",det); */ 113 if (*u > 0.0 || *u < det) 114 return 0; 115 116 /* prepare to test V parameter */ 117 CROSS(qvec, tvec, edge1); 118 119 /* calculate V parameter and test bounds */ 120 *v = DOT(dir, qvec); 121 if (*v > 0.0 || *u + *v < det) 122 return 0; 123 } 124 else return 0; /* ray is parallell to the plane of the triangle */ 125 126 127 inv_det = 1.0 / det; 128 129 /* calculate t, ray intersects triangle */ 130 *t = DOT(edge2, qvec) * inv_det; 131 (*u) *= inv_det; 132 (*v) *= inv_det; 133 134 return 1; 135 } 136 137 /* code rewritten to do tests on the sign of the determinant */ 138 /* the division is before the test of the sign of the det */ 139 int intersect_triangle2(double orig[3], double dir[3], 140 double vert0[3], double vert1[3], double vert2[3], 141 double *t, double *u, double *v) 142 { 143 double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3]; 144 double det, inv_det; 145 146 /* find vectors for two edges sharing vert0 */ 147 SUB(edge1, vert1, vert0); 148 SUB(edge2, vert2, vert0); 149 150 /* begin calculating determinant - also used to calculate U parameter */ 151 CROSS(pvec, dir, edge2); 152 153 /* if determinant is near zero, ray lies in plane of triangle */ 154 det = DOT(edge1, pvec); 155 156 /* calculate distance from vert0 to ray origin */ 157 SUB(tvec, orig, vert0); 158 inv_det = 1.0 / det; 159 160 if (det > EPSILON) 161 { 162 /* calculate U parameter and test bounds */ 163 *u = DOT(tvec, pvec); 164 if (*u < 0.0 || *u > det) 165 return 0; 166 167 /* prepare to test V parameter */ 168 CROSS(qvec, tvec, edge1); 169 170 /* calculate V parameter and test bounds */ 171 *v = DOT(dir, qvec); 172 if (*v < 0.0 || *u + *v > det) 173 return 0; 174 175 } 176 else if (det < -EPSILON) 177 { 178 /* calculate U parameter and test bounds */ 179 *u = DOT(tvec, pvec); 180 if (*u > 0.0 || *u < det) 181 return 0; 182 183 /* prepare to test V parameter */ 184 CROSS(qvec, tvec, edge1); 185 186 /* calculate V parameter and test bounds */ 187 *v = DOT(dir, qvec); 188 if (*v > 0.0 || *u + *v < det) 189 return 0; 190 } 191 else return 0; /* ray is parallell to the plane of the triangle */ 192 193 /* calculate t, ray intersects triangle */ 194 *t = DOT(edge2, qvec) * inv_det; 195 (*u) *= inv_det; 196 (*v) *= inv_det; 197 198 return 1; 199 } 200 201 /* code rewritten to do tests on the sign of the determinant */ 202 /* the division is before the test of the sign of the det */ 203 /* and one CROSS has been moved out from the if-else if-else */ 204 int intersect_triangle3(double orig[3], double dir[3], 205 double vert0[3], double vert1[3], double vert2[3], 206 double *t, double *u, double *v) 207 { 208 double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3]; 209 double det, inv_det; 210 211 /* find vectors for two edges sharing vert0 */ 212 SUB(edge1, vert1, vert0); 213 SUB(edge2, vert2, vert0); 214 215 /* begin calculating determinant - also used to calculate U parameter */ 216 CROSS(pvec, dir, edge2); 217 218 /* if determinant is near zero, ray lies in plane of triangle */ 219 det = DOT(edge1, pvec); 220 221 /* calculate distance from vert0 to ray origin */ 222 SUB(tvec, orig, vert0); 223 inv_det = 1.0 / det; 224 225 CROSS(qvec, tvec, edge1); 226 227 if (det > EPSILON) 228 { 229 *u = DOT(tvec, pvec); 230 if (*u < 0.0 || *u > det) 231 return 0; 232 233 /* calculate V parameter and test bounds */ 234 *v = DOT(dir, qvec); 235 if (*v < 0.0 || *u + *v > det) 236 return 0; 237 238 } 239 else if (det < -EPSILON) 240 { 241 /* calculate U parameter and test bounds */ 242 *u = DOT(tvec, pvec); 243 if (*u > 0.0 || *u < det) 244 return 0; 245 246 /* calculate V parameter and test bounds */ 247 *v = DOT(dir, qvec); 248 if (*v > 0.0 || *u + *v < det) 249 return 0; 250 } 251 else return 0; /* ray is parallell to the plane of the triangle */ 252 253 *t = DOT(edge2, qvec) * inv_det; 254 (*u) *= inv_det; 255 (*v) *= inv_det; 256 257 return 1; 258 }
二、两个三角形相交：
1 #pragma once 2 /* Triangle/triangle intersection test routine, 3 * by Tomas Moller, 1997. 4 * See article "A Fast Triangle-Triangle Intersection Test", 5 * Journal of Graphics Tools, 2(2), 1997 6 * 7 * Updated June 1999: removed the divisions -- a little faster now! 8 * Updated October 1999: added {} to CROSS and SUB macros 9 * 10 * int NoDivTriTriIsect(float V0[3],float V1[3],float V2[3], 11 * float U0[3],float U1[3],float U2[3]) 12 * 13 * parameters: vertices of triangle 1: V0,V1,V2 14 * vertices of triangle 2: U0,U1,U2 15 * result : returns 1 if the triangles intersect, otherwise 0 16 * 17 */ 18 19 #include <math.h> 20 #define FABS(x) (float(fabs(x))) /* implement as is fastest on your machine */ 21 22 /* if USE_EPSILON_TEST is true then we do a check: 23 if |dv|<EPSILON then dv=0.0; 24 else no check is done (which is less robust) 25 */ 26 #define USE_EPSILON_TEST TRUE 27 #define EPSILON 0.000001 28 29 30 /* some macros */ 31 #define CROSS(dest,v1,v2){ \ 32 dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \ 33 dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \ 34 dest[2]=v1[0]*v2[1]-v1[1]*v2[0];} 35 36 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) 37 38 #define SUB(dest,v1,v2){ \ 39 dest[0]=v1[0]-v2[0]; \ 40 dest[1]=v1[1]-v2[1]; \ 41 dest[2]=v1[2]-v2[2];} 42 43 /* sort so that a<=b */ 44 #define SORT(a,b) \ 45 if(a>b) \ 46 { \ 47 float c; \ 48 c=a; \ 49 a=b; \ 50 b=c; \ 51 } 52 53 54 /* this edge to edge test is based on Franlin Antonio's gem: 55 "Faster Line Segment Intersection", in Graphics Gems III, 56 pp. 199-202 */ 57 #define EDGE_EDGE_TEST(V0,U0,U1) \ 58 Bx=U0[i0]-U1[i0]; \ 59 By=U0[i1]-U1[i1]; \ 60 Cx=V0[i0]-U0[i0]; \ 61 Cy=V0[i1]-U0[i1]; \ 62 f=Ay*Bx-Ax*By; \ 63 d=By*Cx-Bx*Cy; \ 64 if((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f)) \ 65 { \ 66 e=Ax*Cy-Ay*Cx; \ 67 if(f>0) \ 68 { \ 69 if(e>=0 && e<=f) return 1; \ 70 } \ 71 else \ 72 { \ 73 if(e<=0 && e>=f) return 1; \ 74 } \ 75 } 76 77 #define EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2) \ 78 { \ 79 float Ax,Ay,Bx,By,Cx,Cy,e,d,f; \ 80 Ax=V1[i0]-V0[i0]; \ 81 Ay=V1[i1]-V0[i1]; \ 82 /* test edge U0,U1 against V0,V1 */ \ 83 EDGE_EDGE_TEST(V0,U0,U1); \ 84 /* test edge U1,U2 against V0,V1 */ \ 85 EDGE_EDGE_TEST(V0,U1,U2); \ 86 /* test edge U2,U1 against V0,V1 */ \ 87 EDGE_EDGE_TEST(V0,U2,U0); \ 88 } 89 90 #define POINT_IN_TRI(V0,U0,U1,U2) \ 91 { \ 92 float a,b,c,d0,d1,d2; \ 93 /* is T1 completly inside T2? */ \ 94 /* check if V0 is inside tri(U0,U1,U2) */ \ 95 a=U1[i1]-U0[i1]; \ 96 b=-(U1[i0]-U0[i0]); \ 97 c=-a*U0[i0]-b*U0[i1]; \ 98 d0=a*V0[i0]+b*V0[i1]+c; \ 99 \ 100 a=U2[i1]-U1[i1]; \ 101 b=-(U2[i0]-U1[i0]); \ 102 c=-a*U1[i0]-b*U1[i1]; \ 103 d1=a*V0[i0]+b*V0[i1]+c; \ 104 \ 105 a=U0[i1]-U2[i1]; \ 106 b=-(U0[i0]-U2[i0]); \ 107 c=-a*U2[i0]-b*U2[i1]; \ 108 d2=a*V0[i0]+b*V0[i1]+c; \ 109 if(d0*d1>0.0) \ 110 { \ 111 if(d0*d2>0.0) return 1; \ 112 } \ 113 } 114 115 int coplanar_tri_tri(float N[3], float V0[3], float V1[3], float V2[3], 116 float U0[3], float U1[3], float U2[3]) 117 { 118 float A[3]; 119 short i0, i1; 120 /* first project onto an axis-aligned plane, that maximizes the area */ 121 /* of the triangles, compute indices: i0,i1. */ 122 A[0] = FABS(N[0]); 123 A[1] = FABS(N[1]); 124 A[2] = FABS(N[2]); 125 if (A[0]>A[1]) 126 { 127 if (A[0]>A[2]) 128 { 129 i0 = 1; /* A[0] is greatest */ 130 i1 = 2; 131 } 132 else 133 { 134 i0 = 0; /* A[2] is greatest */ 135 i1 = 1; 136 } 137 } 138 else /* A[0]<=A[1] */ 139 { 140 if (A[2]>A[1]) 141 { 142 i0 = 0; /* A[2] is greatest */ 143 i1 = 1; 144 } 145 else 146 { 147 i0 = 0; /* A[1] is greatest */ 148 i1 = 2; 149 } 150 } 151 152 /* test all edges of triangle 1 against the edges of triangle 2 */ 153 EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2); 154 EDGE_AGAINST_TRI_EDGES(V1, V2, U0, U1, U2); 155 EDGE_AGAINST_TRI_EDGES(V2, V0, U0, U1, U2); 156 157 /* finally, test if tri1 is totally contained in tri2 or vice versa */ 158 POINT_IN_TRI(V0, U0, U1, U2); 159 POINT_IN_TRI(U0, V0, V1, V2); 160 161 return 0; 162 } 163 164 165 166 #define NEWCOMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,A,B,C,X0,X1) \ 167 { \ 168 if(D0D1>0.0f) \ 169 { \ 170 /* here we know that D0D2<=0.0 */ \ 171 /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \ 172 A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \ 173 } \ 174 else if(D0D2>0.0f)\ 175 { \ 176 /* here we know that d0d1<=0.0 */ \ 177 A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \ 178 } \ 179 else if(D1*D2>0.0f || D0!=0.0f) \ 180 { \ 181 /* here we know that d0d1<=0.0 or that D0!=0.0 */ \ 182 A=VV0; B=(VV1-VV0)*D0; C=(VV2-VV0)*D0; X0=D0-D1; X1=D0-D2; \ 183 } \ 184 else if(D1!=0.0f) \ 185 { \ 186 A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \ 187 } \ 188 else if(D2!=0.0f) \ 189 { \ 190 A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \ 191 } \ 192 else \ 193 { \ 194 /* triangles are coplanar */ \ 195 return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2); \ 196 } \ 197 } 198 199 200 201 int NoDivTriTriIsect(float V0[3], float V1[3], float V2[3], 202 float U0[3], float U1[3], float U2[3]) 203 { 204 float E1[3], E2[3]; 205 float N1[3], N2[3], d1, d2; 206 float du0, du1, du2, dv0, dv1, dv2; 207 float D[3]; 208 float isect1[2], isect2[2]; 209 float du0du1, du0du2, dv0dv1, dv0dv2; 210 short index; 211 float vp0, vp1, vp2; 212 float up0, up1, up2; 213 float bb, cc, max; 214 215 /* compute plane equation of triangle(V0,V1,V2) */ 216 SUB(E1, V1, V0); 217 SUB(E2, V2, V0); 218 CROSS(N1, E1, E2); 219 d1 = -DOT(N1, V0); 220 /* plane equation 1: N1.X+d1=0 */ 221 222 /* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/ 223 du0 = DOT(N1, U0) + d1; 224 du1 = DOT(N1, U1) + d1; 225 du2 = DOT(N1, U2) + d1; 226 227 /* coplanarity robustness check */ 228 #if USE_EPSILON_TEST==TRUE 229 if (FABS(du0)<EPSILON) du0 = 0.0; 230 if (FABS(du1)<EPSILON) du1 = 0.0; 231 if (FABS(du2)<EPSILON) du2 = 0.0; 232 #endif 233 du0du1 = du0*du1; 234 du0du2 = du0*du2; 235 236 if (du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */ 237 return 0; /* no intersection occurs */ 238 239 /* compute plane of triangle (U0,U1,U2) */ 240 SUB(E1, U1, U0); 241 SUB(E2, U2, U0); 242 CROSS(N2, E1, E2); 243 d2 = -DOT(N2, U0); 244 /* plane equation 2: N2.X+d2=0 */ 245 246 /* put V0,V1,V2 into plane equation 2 */ 247 dv0 = DOT(N2, V0) + d2; 248 dv1 = DOT(N2, V1) + d2; 249 dv2 = DOT(N2, V2) + d2; 250 251 #if USE_EPSILON_TEST==TRUE 252 if (FABS(dv0)<EPSILON) dv0 = 0.0; 253 if (FABS(dv1)<EPSILON) dv1 = 0.0; 254 if (FABS(dv2)<EPSILON) dv2 = 0.0; 255 #endif 256 257 dv0dv1 = dv0*dv1; 258 dv0dv2 = dv0*dv2; 259 260 if (dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */ 261 return 0; /* no intersection occurs */ 262 263 /* compute direction of intersection line */ 264 CROSS(D, N1, N2); 265 266 /* compute and index to the largest component of D */ 267 max = (float)FABS(D[0]); 268 index = 0; 269 bb = (float)FABS(D[1]); 270 cc = (float)FABS(D[2]); 271 if (bb>max) max = bb, index = 1; 272 if (cc>max) max = cc, index = 2; 273 274 /* this is the simplified projection onto L*/ 275 vp0 = V0[index]; 276 vp1 = V1[index]; 277 vp2 = V2[index]; 278 279 up0 = U0[index]; 280 up1 = U1[index]; 281 up2 = U2[index]; 282 283 /* compute interval for triangle 1 */ 284 float a, b, c, x0, x1; 285 NEWCOMPUTE_INTERVALS(vp0, vp1, vp2, dv0, dv1, dv2, dv0dv1, dv0dv2, a, b, c, x0, x1); 286 287 /* compute interval for triangle 2 */ 288 float d, e, f, y0, y1; 289 NEWCOMPUTE_INTERVALS(up0, up1, up2, du0, du1, du2, du0du1, du0du2, d, e, f, y0, y1); 290 291 float xx, yy, xxyy, tmp; 292 xx = x0*x1; 293 yy = y0*y1; 294 xxyy = xx*yy; 295 296 tmp = a*xxyy; 297 isect1[0] = tmp + b*x1*yy; 298 isect1[1] = tmp + c*x0*yy; 299 300 tmp = d*xxyy; 301 isect2[0] = tmp + e*xx*y1; 302 isect2[1] = tmp + f*xx*y0; 303 304 SORT(isect1[0], isect1[1]); 305 SORT(isect2[0], isect2[1]); 306 307 if (isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0; 308 return 1; 309 }
三、三角形与轴对齐盒体相交：
1 #pragma once 2 3 /********************************************************/ 4 /* AABB-triangle overlap test code */ 5 /* by Tomas Akenine-Möller */ 6 /* Function: int triBoxOverlap(float boxcenter[3], */ 7 /* float boxhalfsize[3],float triverts[3][3]); */ 8 /* History: */ 9 /* 2001-03-05: released the code in its first version */ 10 /* 2001-06-18: changed the order of the tests, faster */ 11 /* */ 12 /* Acknowledgement: Many thanks to Pierre Terdiman for */ 13 /* suggestions and discussions on how to optimize code. */ 14 /* Thanks to David Hunt for finding a ">="-bug! */ 15 /********************************************************/ 16 17 #include <math.h> 18 #include <stdio.h> 19 20 #define X 0 21 #define Y 1 22 #define Z 2 23 24 #define CROSS(dest,v1,v2) \ 25 dest[0] = v1[1] * v2[2] - v1[2] * v2[1]; \ 26 dest[1] = v1[2] * v2[0] - v1[0] * v2[2]; \ 27 dest[2] = v1[0] * v2[1] - v1[1] * v2[0]; 28 29 #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]) 30 31 #define SUB(dest,v1,v2) \ 32 dest[0] = v1[0] - v2[0]; \ 33 dest[1] = v1[1] - v2[1]; \ 34 dest[2] = v1[2] - v2[2]; 35 36 #define FINDMINMAX(x0,x1,x2,min,max) \ 37 min = max = x0; \ 38 if (x1<min) min = x1; \ 39 if (x1>max) max = x1; \ 40 if (x2<min) min = x2; \ 41 if (x2>max) max = x2; 42 43 int planeBoxOverlap(float normal[3], float vert[3], float maxbox[3]) // -NJMP- 44 { 45 int q; 46 float vmin[3], vmax[3], v; 47 for (q = X; q <= Z; q++) 48 { 49 v = vert[q]; // -NJMP- 50 if (normal[q]>0.0f) 51 { 52 vmin[q] = -maxbox[q] - v; // -NJMP- 53 vmax[q] = maxbox[q] - v; // -NJMP- 54 } 55 else 56 { 57 vmin[q] = maxbox[q] - v; // -NJMP- 58 vmax[q] = -maxbox[q] - v; // -NJMP- 59 } 60 } 61 if (DOT(normal, vmin)>0.0f) return 0; // -NJMP- 62 if (DOT(normal, vmax) >= 0.0f) return 1; // -NJMP- 63 return 0; 64 } 65 66 /*======================== X-tests ========================*/ 67 68 #define AXISTEST_X01(a, b, fa, fb) \ 69 p0 = a*v0[Y] - b*v0[Z]; \ 70 p2 = a*v2[Y] - b*v2[Z]; \ 71 if (p0<p2) { min = p0; max = p2; }\ 72 else { min = p2; max = p0; } \ 73 rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \ 74 if (min>rad || max<-rad) return 0; 75 76 #define AXISTEST_X2(a, b, fa, fb) \ 77 p0 = a*v0[Y] - b*v0[Z]; \ 78 p1 = a*v1[Y] - b*v1[Z]; \ 79 if (p0<p1) { min = p0; max = p1; }\ 80 else { min = p1; max = p0; } \ 81 rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \ 82 if (min>rad || max<-rad) return 0; 83 84 /*======================== Y-tests ========================*/ 85 86 #define AXISTEST_Y02(a, b, fa, fb) \ 87 p0 = -a*v0[X] + b*v0[Z]; \ 88 p2 = -a*v2[X] + b*v2[Z]; \ 89 if (p0<p2) { min = p0; max = p2; }\ 90 else { min = p2; max = p0; } \ 91 rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \ 92 if (min>rad || max<-rad) return 0; 93 94 #define AXISTEST_Y1(a, b, fa, fb) \ 95 p0 = -a*v0[X] + b*v0[Z]; \ 96 p1 = -a*v1[X] + b*v1[Z]; \ 97 if (p0<p1) { min = p0; max = p1; }\ 98 else { min = p1; max = p0; } \ 99 rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \ 100 if (min>rad || max<-rad) return 0; 101 102 /*======================== Z-tests ========================*/ 103 104 #define AXISTEST_Z12(a, b, fa, fb) \ 105 p1 = a*v1[X] - b*v1[Y]; \ 106 p2 = a*v2[X] - b*v2[Y]; \ 107 if (p2<p1) { min = p2; max = p1; }\ 108 else { min = p1; max = p2; } \ 109 rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \ 110 if (min>rad || max<-rad) return 0; 111 112 #define AXISTEST_Z0(a, b, fa, fb) \ 113 p0 = a*v0[X] - b*v0[Y]; \ 114 p1 = a*v1[X] - b*v1[Y]; \ 115 if (p0<p1) { min = p0; max = p1; }\ 116 else { min = p1; max = p0; } \ 117 rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \ 118 if (min>rad || max<-rad) return 0; 119 120 int triBoxOverlap(float* vert1, float* vert2, float* vert3, float* boxcenter, float* boxhalfsize) 121 { 122 /* use separating axis theorem to test overlap between triangle and box */ 123 /* need to test for overlap in these directions: */ 124 /* 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle */ 125 /* we do not even need to test these) */ 126 /* 2) normal of the triangle */ 127 /* 3) crossproduct(edge from tri, {x,y,z}-directin) */ 128 /* this gives 3x3=9 more tests */ 129 float v0[3], v1[3], v2[3]; 130 float min, max, p0, p1, p2, rad, fex, fey, fez; // -NJMP- "d" local variable removed 131 float normal[3], e0[3], e1[3], e2[3]; 132 133 /* This is the fastest branch on Sun */ 134 /* move everything so that the boxcenter is in (0,0,0) */ 135 136 SUB(v0, vert1, boxcenter); 137 SUB(v1, vert2, boxcenter); 138 SUB(v2, vert3, boxcenter); 139 140 /* compute triangle edges */ 141 SUB(e0, v1, v0); /* tri edge 0 */ 142 SUB(e1, v2, v1); /* tri edge 1 */ 143 SUB(e2, v0, v2); /* tri edge 2 */ 144 145 /* Bullet 3: */ 146 /* test the 9 tests first (this was faster) */ 147 fex = fabsf(e0[X]); 148 fey = fabsf(e0[Y]); 149 fez = fabsf(e0[Z]); 150 AXISTEST_X01(e0[Z], e0[Y], fez, fey); 151 AXISTEST_Y02(e0[Z], e0[X], fez, fex); 152 AXISTEST_Z12(e0[Y], e0[X], fey, fex); 153 154 fex = fabsf(e1[X]); 155 fey = fabsf(e1[Y]); 156 fez = fabsf(e1[Z]); 157 AXISTEST_X01(e1[Z], e1[Y], fez, fey); 158 AXISTEST_Y02(e1[Z], e1[X], fez, fex); 159 AXISTEST_Z0(e1[Y], e1[X], fey, fex); 160 161 fex = fabsf(e2[X]); 162 fey = fabsf(e2[Y]); 163 fez = fabsf(e2[Z]); 164 AXISTEST_X2(e2[Z], e2[Y], fez, fey); 165 AXISTEST_Y1(e2[Z], e2[X], fez, fex); 166 AXISTEST_Z12(e2[Y], e2[X], fey, fex); 167 168 /* Bullet 1: */ 169 /* first test overlap in the {x,y,z}-directions */ 170 /* find min, max of the triangle each direction, and test for overlap in */ 171 /* that direction -- this is equivalent to testing a minimal AABB around */ 172 /* the triangle against the AABB */ 173 174 /* test in X-direction */ 175 FINDMINMAX(v0[X], v1[X], v2[X], min, max); 176 if (min>boxhalfsize[X] || max<-boxhalfsize[X]) return 0; 177 178 /* test in Y-direction */ 179 FINDMINMAX(v0[Y], v1[Y], v2[Y], min, max); 180 if (min>boxhalfsize[Y] || max<-boxhalfsize[Y]) return 0; 181 182 /* test in Z-direction */ 183 FINDMINMAX(v0[Z], v1[Z], v2[Z], min, max); 184 if (min>boxhalfsize[Z] || max<-boxhalfsize[Z]) return 0; 185 186 /* Bullet 2: */ 187 /* test if the box intersects the plane of the triangle */ 188 /* compute plane equation of triangle: normal*x+d=0 */ 189 CROSS(normal, e0, e1); 190 if (!planeBoxOverlap(normal, v0, boxhalfsize)) return 0; // -NJMP- 191 192 return 1; /* box and triangle overlaps */ 193 }
7140
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
