home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / triglcub.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-14  |  10.6 KB  |  288 lines
  1. #include <math.h>
  2.  
  3. #define SIGN3( A ) (((A).x<0)?4:0 | ((A).y<0)?2:0 | ((A).z<0)?1:0)
  4. #define CROSS( A, B, C ) { \
  5.   (C).x =  (A).y * (B).z - (A).z * (B).y; \
  6.   (C).y = -(A).x * (B).z + (A).z * (B).x; \
  7.   (C).z =  (A).x * (B).y - (A).y * (B).x; \
  8.    }
  9. #define SUB( A, B, C ) { \
  10.   (C).x =  (A).x - (B).x; \
  11.   (C).y =  (A).y - (B).y; \
  12.   (C).z =  (A).z - (B).z; \
  13.    }
  14. #define LERP( A, B, C) ((B)+(A)*((C)-(B)))
  15. #define MIN3(a,b,c) ((((a)<(b))&&((a)<(c))) ? (a) : (((b)<(c)) ? (b) : (c)))
  16. #define MAX3(a,b,c) ((((a)>(b))&&((a)>(c))) ? (a) : (((b)>(c)) ? (b) : (c)))
  17. #define INSIDE 0
  18. #define OUTSIDE 1
  19.  
  20. typedef struct {
  21.       float           x;
  22.       float           y;
  23.       float           z;
  24.       } Point3;
  25.  
  26. typedef struct{
  27.    Point3 v1;                 /* Vertex1 */
  28.    Point3 v2;                 /* Vertex2 */
  29.    Point3 v3;                 /* Vertex3 */
  30.    } Triangle3; 
  31.    
  32. /*___________________________________________________________________________*/
  33.  
  34. /* Which of the six face-plane(s) is point P outside of? */
  35.  
  36. long face_plane(Point3 p)
  37. {
  38. long outcode;
  39.  
  40.    outcode = 0;
  41.    if (p.x >  .5) outcode |= 0x01;
  42.    if (p.x < -.5) outcode |= 0x02;
  43.    if (p.y >  .5) outcode |= 0x04;
  44.    if (p.y < -.5) outcode |= 0x08;
  45.    if (p.z >  .5) outcode |= 0x10;
  46.    if (p.z < -.5) outcode |= 0x20;
  47.    return(outcode);
  48. }
  49.  
  50. /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
  51.  
  52. /* Which of the twelve edge plane(s) is point P outside of? */
  53.  
  54. long bevel_2d(Point3 p)
  55. {
  56. long outcode;
  57.  
  58.    outcode = 0;
  59.    if ( p.x + p.y > 1.0) outcode |= 0x001;
  60.    if ( p.x - p.y > 1.0) outcode |= 0x002;
  61.    if (-p.x + p.y > 1.0) outcode |= 0x004;
  62.    if (-p.x - p.y > 1.0) outcode |= 0x008;
  63.    if ( p.x + p.z > 1.0) outcode |= 0x010;
  64.    if ( p.x - p.z > 1.0) outcode |= 0x020;
  65.    if (-p.x + p.z > 1.0) outcode |= 0x040;
  66.    if (-p.x - p.z > 1.0) outcode |= 0x080;
  67.    if ( p.y + p.z > 1.0) outcode |= 0x100;
  68.    if ( p.y - p.z > 1.0) outcode |= 0x200;
  69.    if (-p.y + p.z > 1.0) outcode |= 0x400;
  70.    if (-p.y - p.z > 1.0) outcode |= 0x800;
  71.    return(outcode);
  72. }
  73.  
  74. /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
  75.  
  76. /* Which of the eight corner plane(s) is point P outside of? */
  77.  
  78. long bevel_3d(Point3 p)
  79. {
  80. long outcode;
  81.  
  82.    outcode = 0;
  83.    if (( p.x + p.y + p.z) > 1.5) outcode |= 0x01;
  84.    if (( p.x + p.y - p.z) > 1.5) outcode |= 0x02;
  85.    if (( p.x - p.y + p.z) > 1.5) outcode |= 0x04;
  86.    if (( p.x - p.y - p.z) > 1.5) outcode |= 0x08;
  87.    if ((-p.x + p.y + p.z) > 1.5) outcode |= 0x10;
  88.    if ((-p.x + p.y - p.z) > 1.5) outcode |= 0x20;
  89.    if ((-p.x - p.y + p.z) > 1.5) outcode |= 0x40;
  90.    if ((-p.x - p.y - p.z) > 1.5) outcode |= 0x80;
  91.    return(outcode);
  92. }
  93.  
  94. /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
  95.  
  96. /* Test the point "alpha" of the way from P1 to P2 */
  97. /* See if it is on a face of the cube              */
  98. /* Consider only faces in "mask"                   */
  99.  
  100. long check_point(Point3 p1, Point3 p2, float alpha, long mask)
  101. {
  102. Point3 plane_point;
  103.  
  104.    plane_point.x = LERP(alpha, p1.x, p2.x);
  105.    plane_point.y = LERP(alpha, p1.y, p2.y);
  106.    plane_point.z = LERP(alpha, p1.z, p2.z);
  107.    return(face_plane(plane_point) & mask);
  108. }
  109.  
  110. /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
  111.  
  112. /* Compute intersection of P1 --> P2 line segment with face planes */
  113. /* Then test intersection point to see if it is on cube face       */
  114. /* Consider only face planes in "outcode_diff"                     */
  115. /* Note: Zero bits in "outcode_diff" means face line is outside of */
  116.  
  117. long check_line(Point3 p1, Point3 p2, long outcode_diff)
  118. {
  119.  
  120.    if ((0x01 & outcode_diff) != 0)
  121.       if (check_point(p1,p2,( .5-p2.x)/(p1.x-p2.x),0x3e) == INSIDE) return(INSIDE);
  122.    if ((0x02 & outcode_diff) != 0)
  123.       if (check_point(p1,p2,(-.5-p2.x)/(p1.x-p2.x),0x3d) == INSIDE) return(INSIDE);
  124.    if ((0x04 & outcode_diff) != 0) 
  125.       if (check_point(p1,p2,( .5-p2.y)/(p1.y-p2.y),0x3b) == INSIDE) return(INSIDE);
  126.    if ((0x08 & outcode_diff) != 0) 
  127.       if (check_point(p1,p2,(-.5-p2.y)/(p1.y-p2.y),0x37) == INSIDE) return(INSIDE);
  128.    if ((0x10 & outcode_diff) != 0) 
  129.       if (check_point(p1,p2,( .5-p2.z)/(p1.z-p2.z),0x2f) == INSIDE) return(INSIDE);
  130.    if ((0x20 & outcode_diff) != 0) 
  131.       if (check_point(p1,p2,(-.5-p2.z)/(p1.z-p2.z),0x1f) == INSIDE) return(INSIDE);
  132.    return(OUTSIDE);
  133. }
  134.  
  135. /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
  136.  
  137. /* Test if 3D point is inside 3D triangle */
  138.  
  139. long point_triangle_intersection(Point3 p, Triangle3 t)
  140. {
  141. long sign12,sign23,sign31;
  142. Point3 vect12,vect23,vect31,vect1h,vect2h,vect3h;
  143. Point3 cross12_1p,cross23_2p,cross31_3p;
  144.  
  145. /* First, a quick bounding-box test:                               */
  146. /* If P is outside triangle bbox, there cannot be an intersection. */
  147.  
  148.    if (p.x > MAX3(t.v1.x, t.v2.x, t.v3.x)) return(OUTSIDE);  
  149.    if (p.y > MAX3(t.v1.y, t.v2.y, t.v3.y)) return(OUTSIDE);
  150.    if (p.z > MAX3(t.v1.z, t.v2.z, t.v3.z)) return(OUTSIDE);
  151.    if (p.x < MIN3(t.v1.x, t.v2.x, t.v3.x)) return(OUTSIDE);
  152.    if (p.y < MIN3(t.v1.y, t.v2.y, t.v3.y)) return(OUTSIDE);
  153.    if (p.z < MIN3(t.v1.z, t.v2.z, t.v3.z)) return(OUTSIDE);
  154.  
  155. /* For each triangle side, make a vector out of it by subtracting vertexes; */
  156. /* make another vector from one vertex to point P.                          */
  157. /* The crossproduct of these two vectors is orthogonal to both and the      */
  158. /* signs of its X,Y,Z components indicate whether P was to the inside or    */
  159. /* to the outside of this triangle side.                                    */
  160.  
  161.    SUB(t.v1, t.v2, vect12)
  162.    SUB(t.v1,    p, vect1h);
  163.    CROSS(vect12, vect1h, cross12_1p)
  164.    sign12 = SIGN3(cross12_1p);      /* Extract X,Y,Z signs as 0...7 integer */
  165.  
  166.    SUB(t.v2, t.v3, vect23)
  167.    SUB(t.v2,    p, vect2h);
  168.    CROSS(vect23, vect2h, cross23_2p)
  169.    sign23 = SIGN3(cross23_2p);
  170.  
  171.    SUB(t.v3, t.v1, vect31)
  172.    SUB(t.v3,    p, vect3h);
  173.    CROSS(vect31, vect3h, cross31_3p)
  174.    sign31 = SIGN3(cross31_3p);
  175.  
  176. /* If all three crossproduct vectors agree in their component signs,  */
  177. /* then the point must be inside all three.                           */
  178. /* P cannot be OUTSIDE all three sides simultaneously.                */
  179.  
  180.    if ((sign12 == sign23) && (sign23 == sign31))
  181.       return(INSIDE);
  182.    else
  183.       return(OUTSIDE);
  184. }
  185.  
  186. /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
  187.  
  188. /**********************************************/
  189. /* This is the main algorithm procedure.      */
  190. /* Triangle t is compared with a unit cube,   */
  191. /* centered on the origin.                    */
  192. /* It returns INSIDE (0) or OUTSIDE(1) if t   */
  193. /* intersects or does not intersect the cube. */
  194. /**********************************************/
  195.  
  196. long t_c_intersection(Triangle3 t)
  197. {
  198. long v1_test,v2_test,v3_test;
  199. float d;
  200. Point3 vect12,vect13,norm;
  201. Point3 hitpp,hitpn,hitnp,hitnn;
  202.  
  203. /* First compare all three vertexes with all six face-planes */
  204. /* If any vertex is inside the cube, return immediately!     */
  205.  
  206.    if ((v1_test = face_plane(t.v1)) == INSIDE) return(INSIDE);
  207.    if ((v2_test = face_plane(t.v2)) == INSIDE) return(INSIDE);
  208.    if ((v3_test = face_plane(t.v3)) == INSIDE) return(INSIDE);
  209.  
  210. /* If all three vertexes were outside of one or more face-planes, */
  211. /* return immediately with a trivial rejection!                   */
  212.  
  213.    if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);
  214.  
  215. /* Now do the same trivial rejection test for the 12 edge planes */
  216.  
  217.    v1_test |= bevel_2d(t.v1) << 8; 
  218.    v2_test |= bevel_2d(t.v2) << 8; 
  219.    v3_test |= bevel_2d(t.v3) << 8;
  220.    if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);  
  221.  
  222. /* Now do the same trivial rejection test for the 8 corner planes */
  223.  
  224.    v1_test |= bevel_3d(t.v1) << 24; 
  225.    v2_test |= bevel_3d(t.v2) << 24; 
  226.    v3_test |= bevel_3d(t.v3) << 24; 
  227.    if ((v1_test & v2_test & v3_test) != 0) return(OUTSIDE);   
  228.  
  229. /* If vertex 1 and 2, as a pair, cannot be trivially rejected */
  230. /* by the above tests, then see if the v1-->v2 triangle edge  */
  231. /* intersects the cube.  Do the same for v1-->v3 and v2-->v3. */
  232. /* Pass to the intersection algorithm the "OR" of the outcode */
  233. /* bits, so that only those cube faces which are spanned by   */
  234. /* each triangle edge need be tested.                         */
  235.  
  236.    if ((v1_test & v2_test) == 0)
  237.       if (check_line(t.v1,t.v2,v1_test|v2_test) == INSIDE) return(INSIDE);
  238.    if ((v1_test & v3_test) == 0)
  239.       if (check_line(t.v1,t.v3,v1_test|v3_test) == INSIDE) return(INSIDE);
  240.    if ((v2_test & v3_test) == 0)
  241.       if (check_line(t.v2,t.v3,v2_test|v3_test) == INSIDE) return(INSIDE);
  242.  
  243. /* By now, we know that the triangle is not off to any side,     */
  244. /* and that its sides do not penetrate the cube.  We must now    */
  245. /* test for the cube intersecting the interior of the triangle.  */
  246. /* We do this by looking for intersections between the cube      */
  247. /* diagonals and the triangle...first finding the intersection   */
  248. /* of the four diagonals with the plane of the triangle, and     */
  249. /* then if that intersection is inside the cube, pursuing        */
  250. /* whether the intersection point is inside the triangle itself. */
  251.  
  252. /* To find plane of the triangle, first perform crossproduct on  */
  253. /* two triangle side vectors to compute the normal vector.       */  
  254.                                 
  255.    SUB(t.v1,t.v2,vect12);
  256.    SUB(t.v1,t.v3,vect13);
  257.    CROSS(vect12,vect13,norm)
  258.  
  259. /* The normal vector "norm" X,Y,Z components are the coefficients */
  260. /* of the triangles AX + BY + CZ + D = 0 plane equation.  If we   */
  261. /* solve the plane equation for X=Y=Z (a diagonal), we get        */
  262. /* -D/(A+B+C) as a metric of the distance from cube center to the */
  263. /* diagonal/plane intersection.  If this is between -0.5 and 0.5, */
  264. /* the intersection is inside the cube.  If so, we continue by    */
  265. /* doing a point/triangle intersection.                           */
  266. /* Do this for all four diagonals.                                */
  267.  
  268.    d = norm.x * t.v1.x + norm.y * t.v1.y + norm.z * t.v1.z;
  269.    hitpp.x = hitpp.y = hitpp.z = d / (norm.x + norm.y + norm.z);
  270.    if (fabsf(hitpp.x) <= 0.5)
  271.       if (point_triangle_intersection(hitpp,t) == INSIDE) return(INSIDE);
  272.    hitpn.z = -(hitpn.x = hitpn.y = d / (norm.x + norm.y - norm.z));
  273.    if (fabsf(hitpn.x) <= 0.5)
  274.       if (point_triangle_intersection(hitpn,t) == INSIDE) return(INSIDE);
  275.    hitnp.y = -(hitnp.x = hitnp.z = d / (norm.x - norm.y + norm.z));
  276.    if (fabsf(hitnp.x) <= 0.5)
  277.       if (point_triangle_intersection(hitnp,t) == INSIDE) return(INSIDE);
  278.    hitnn.y = hitnn.z = -(hitnn.x = d / (norm.x - norm.y - norm.z));
  279.    if (fabsf(hitnn.x) <= 0.5)
  280.       if (point_triangle_intersection(hitnn,t) == INSIDE) return(INSIDE);
  281.  
  282. /* No edge touched the cube; no cube diagonal touched the triangle. */
  283. /* We're done...there was no intersection.                          */
  284.  
  285.    return(OUTSIDE);
  286.  
  287. }
  288.