home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / macraysh.sit / Code / Source / triangle.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-10  |  11.7 KB  |  474 lines
  1. /*
  2.  * triangle.c
  3.  *
  4.  * Copyright (C) 1989, 1991, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * $Id: triangle.c,v 4.0 91/07/17 14:39:38 kolb Exp Locker: kolb $
  17.  *
  18.  * $Log:    triangle.c,v $
  19.  * Revision 4.0  91/07/17  14:39:38  kolb
  20.  * Initial version.
  21.  * 
  22.  */
  23. #include "geom.h"
  24. #include "triangle.h"
  25.  
  26. Methods *iTriangleMethods = NULL;
  27. static char triName[] = "triangle";
  28.  
  29. unsigned long TriTests, TriHits;
  30.  
  31. static void TriangleSetdPdUV();
  32.  
  33. /*
  34.  * Create and return reference to a triangle.
  35.  */
  36. Triangle *
  37. TriangleCreate(type, p1, p2, p3, n1, n2, n3, u1, u2, u3, flipflag)
  38. int    type;
  39. Vector    *p1, *p2, *p3, *n1, *n2, *n3;
  40. Vec2d    *u1, *u2, *u3;
  41. int flipflag;
  42. {
  43.     Triangle *triangle;
  44.     Vector ptmp, anorm;
  45.     Float d;
  46.  
  47.     /*
  48.      * Allocate new triangle and primitive to point to it.
  49.      */
  50.     triangle = (Triangle *)share_malloc(sizeof(Triangle));
  51.     triangle->type = type;    /* so inttri can tell the difference */
  52.  
  53.     VecSub(*p2, *p1, &triangle->e[0]);
  54.     VecSub(*p3, *p2, &triangle->e[1]);
  55.     VecSub(*p1, *p3, &triangle->e[2]);
  56.  
  57.     /* Find plane normal. */
  58.     VecCross(&triangle->e[0], &triangle->e[1], &ptmp);
  59.     triangle->nrm = ptmp;
  60.     if (VecNormalize(&triangle->nrm) == 0.) {
  61.         RLerror(RL_ADVISE, "Degenerate triangle.\n",0,0,0);
  62.         return (Triangle *)NULL;
  63.     }
  64.  
  65.     if (flipflag)
  66.         VecScale(-1, triangle->nrm, &triangle->nrm);
  67.  
  68.     triangle->d = dotp(&triangle->nrm, p1);
  69.  
  70.     triangle->p[0] = *p1;
  71.     triangle->p[1] = *p2;
  72.     triangle->p[2] = *p3;
  73.  
  74.     triangle->vnorm = NULL;
  75.     if (type == PHONGTRI) {
  76.         if (VecNormalize(n1) == 0. || VecNormalize(n2) == 0. ||
  77.             VecNormalize(n3) == 0.) {
  78.             RLerror(RL_WARN, "Degenerate vertex normal.\n",0,0,0);
  79.             return (Triangle *)NULL;
  80.         }
  81.         triangle->vnorm = (Vector *)Malloc(3 * sizeof(Vector));
  82.         triangle->vnorm[0] = *n1;
  83.         triangle->vnorm[1] = *n2;
  84.         triangle->vnorm[2] = *n3;
  85.         if (flipflag) {
  86.             /* Flip the vertex normals */
  87.             VecScale(-1, triangle->vnorm[0],
  88.                 &triangle->vnorm[0]);
  89.             VecScale(-1, triangle->vnorm[1],
  90.                 &triangle->vnorm[1]);
  91.             VecScale(-1, triangle->vnorm[2],
  92.                 &triangle->vnorm[2]);
  93.         } else if (dotp(&triangle->vnorm[0], &triangle->nrm) < 0.) {
  94.             /*
  95.              * Reverse direction of surface normal on Phong
  96.              * triangle if the surface normal points "away"
  97.              * from the first vertex normal.
  98.              * Note that this means that we trust the vertex
  99.              * normals rather than trust that the user gave the
  100.              * vertices in the correct order.
  101.              */
  102.             RLerror(RL_ADVISE, "Inconsistant triangle normals.\n",0,0,0);
  103.             VecScale(-1., triangle->nrm, &triangle->nrm);
  104.             VecScale(-1., ptmp, &ptmp);
  105.             triangle->d = -triangle->d;
  106.             VecScale(-1., triangle->e[0], &triangle->e[0]);
  107.             VecScale(-1., triangle->e[1], &triangle->e[1]);
  108.             VecScale(-1., triangle->e[2], &triangle->e[2]);
  109.         }
  110.     }
  111.  
  112.     /*
  113.      * If UV coordinates are given for the vertices, allocate and
  114.      * store them.
  115.      */
  116.     triangle->uv = NULL;
  117.     if (u1 && u2 && u3) {
  118.         triangle->uv = (Vec2d *)Malloc(3 * sizeof(Vec2d));
  119.         triangle->uv[0] = *u1;
  120.         triangle->uv[1] = *u2;
  121.         triangle->uv[2] = *u3;
  122.         /* Calculate dpdu and dpdv vectors */
  123.         triangle->dpdu = (Vector *)Malloc(sizeof(Vector));
  124.         triangle->dpdv = (Vector *)Malloc(sizeof(Vector));
  125.         TriangleSetdPdUV(triangle->p, triangle->uv,
  126.                  triangle->dpdu, triangle->dpdv);
  127.     } else {
  128.         triangle->uv = (Vec2d *)NULL;
  129.     }
  130.  
  131.     /*
  132.      * Find "dominant" part of normal vector.
  133.      */
  134.     anorm.x = fabs(ptmp.x);
  135.     anorm.y = fabs(ptmp.y);
  136.     anorm.z = fabs(ptmp.z);
  137.  
  138.     /*
  139.      * Scale edges by dominant part of normal.  This makes intersection
  140.      * testing a bit faster.
  141.      */ 
  142.     if (anorm.x > anorm.y && anorm.x > anorm.z) {
  143.         triangle->index = XNORMAL;
  144.         d = 1. / ptmp.x;
  145.     } else if (anorm.y > anorm.z) {
  146.         triangle->index = YNORMAL;
  147.         d = 1. / ptmp.y;
  148.     } else {
  149.         triangle->index = ZNORMAL;
  150.         d = 1. /ptmp.z;
  151.     }
  152.  
  153.     VecScale(d, triangle->e[0], &triangle->e[0]);
  154.     VecScale(d, triangle->e[1], &triangle->e[1]);
  155.     VecScale(d, triangle->e[2], &triangle->e[2]);
  156.  
  157.     return triangle;
  158. }
  159.  
  160. Methods *
  161. TriangleMethods()
  162. {
  163.     if (iTriangleMethods == (Methods *)NULL) {
  164.         iTriangleMethods = MethodsCreate();
  165.         iTriangleMethods->create = (GeomCreateFunc *)TriangleCreate;
  166.         iTriangleMethods->methods = TriangleMethods;
  167.         iTriangleMethods->name = TriangleName;
  168.         iTriangleMethods->intersect = TriangleIntersect;
  169.         iTriangleMethods->normal = TriangleNormal;
  170.         iTriangleMethods->uv = TriangleUV;
  171.         iTriangleMethods->bounds = TriangleBounds;
  172.         iTriangleMethods->stats = TriangleStats;
  173.         iTriangleMethods->checkbounds = TRUE;
  174.         iTriangleMethods->closed = FALSE;
  175.     }
  176.     return iTriangleMethods;
  177. }
  178.  
  179. /*
  180.  * Intersect ray with triangle.  This is an optimized version of the
  181.  * intersection routine from Snyder and Barr's '87 SIGGRAPH paper.
  182.  */
  183. int
  184. TriangleIntersect(tri, ray, mindist, maxdist)
  185. Triangle *tri;
  186. Ray *ray;
  187. Float mindist, *maxdist;
  188. {
  189.     Float qi1, qi2, s, k, b0, b1, b2;
  190.     Vector pos, dir;
  191.  
  192.     TriTests++;
  193.     pos = ray->pos;
  194.     dir = ray->dir;
  195.     /*
  196.      * Plane intersection.
  197.      */
  198.     k = dotp(&tri->nrm, &dir);
  199.     if (fabs(k) < EPSILON)
  200.         return FALSE;
  201.     s = (tri->d - dotp(&tri->nrm, &pos)) / k;
  202.     if (s < mindist || s > *maxdist)
  203.         return FALSE;
  204.     
  205.     if (tri->index == XNORMAL) {
  206.         qi1 = pos.y + s * dir.y;
  207.         qi2 = pos.z + s * dir.z;
  208.         b0 = tri->e[1].y * (qi2 - tri->p[1].z) -
  209.                 tri->e[1].z * (qi1 - tri->p[1].y);
  210.         if (b0 < 0. || b0 > 1.)
  211.             return FALSE;
  212.         b1 = tri->e[2].y * (qi2 - tri->p[2].z) -
  213.                 tri->e[2].z * (qi1 - tri->p[2].y);
  214.         if (b1 < 0. || b1 > 1.)
  215.             return FALSE;
  216.         b2 = tri->e[0].y * (qi2 - tri->p[0].z) -
  217.                 tri->e[0].z * (qi1 - tri->p[0].y);
  218.         if (b2 < 0. || b2 > 1.)
  219.             return FALSE;
  220.     } else if (tri->index == YNORMAL) {
  221.         qi1 = pos.x + s * dir.x;
  222.         qi2 = pos.z + s * dir.z;
  223.         b0 = tri->e[1].z * (qi1 - tri->p[1].x) -
  224.             tri->e[1].x * (qi2 - tri->p[1].z);
  225.         if (b0 < 0. || b0 > 1.)
  226.             return FALSE;
  227.         b1 = tri->e[2].z * (qi1 - tri->p[2].x) -
  228.             tri->e[2].x * (qi2 - tri->p[2].z);
  229.         if (b1 < 0. || b1 > 1.)
  230.             return FALSE;
  231.         b2 = tri->e[0].z * (qi1 - tri->p[0].x) -
  232.             tri->e[0].x * (qi2 - tri->p[0].z);
  233.         if (b2 < 0. || b2 > 1.)
  234.             return FALSE;
  235.     } else {
  236.         qi1 = pos.x + s * dir.x;
  237.         qi2 = pos.y + s * dir.y;
  238.         b0 = tri->e[1].x * (qi2 - tri->p[1].y) -
  239.             tri->e[1].y * (qi1 - tri->p[1].x);
  240.         if (b0 < 0. || b0 > 1.)
  241.             return FALSE;
  242.         b1 = tri->e[2].x * (qi2 - tri->p[2].y) -
  243.                 tri->e[2].y * (qi1 - tri->p[2].x);
  244.         if (b1 < 0. || b1 > 1.)
  245.             return FALSE;
  246.         b2 = tri->e[0].x * (qi2 - tri->p[0].y) -
  247.                 tri->e[0].y * (qi1 - tri->p[0].x);
  248.         if (b2 < 0. || b2 > 1.)
  249.             return FALSE;
  250.     }
  251.  
  252.     tri->b[0] = b0;
  253.     tri->b[1] = b1;
  254.     tri->b[2] = b2;
  255.  
  256.     TriHits++;
  257.     *maxdist = s;
  258.     return TRUE;
  259. }
  260.  
  261. int
  262. TriangleNormal(tri, pos, nrm, gnrm)
  263. Triangle *tri;
  264. Vector *pos, *nrm, *gnrm;
  265. {
  266.     *gnrm = tri->nrm;
  267.  
  268.     if (tri->type == FLATTRI) {
  269.         *nrm = tri->nrm;
  270.         return FALSE;
  271.     }
  272.  
  273.     /*
  274.      * Interpolate normals of Phong-shaded triangles.
  275.      */
  276.     nrm->x = tri->b[0]*tri->vnorm[0].x+tri->b[1]*tri->vnorm[1].x+
  277.         tri->b[2]*tri->vnorm[2].x;
  278.     nrm->y = tri->b[0]*tri->vnorm[0].y+tri->b[1]*tri->vnorm[1].y+
  279.         tri->b[2]*tri->vnorm[2].y;
  280.     nrm->z = tri->b[0]*tri->vnorm[0].z+tri->b[1]*tri->vnorm[1].z+
  281.         tri->b[2]*tri->vnorm[2].z;
  282.     (void)VecNormalize(nrm);
  283.     return TRUE;
  284. }
  285.  
  286. /*ARGSUSED*/
  287. void
  288. TriangleUV(tri, pos, norm, uv, dpdu, dpdv)
  289. Triangle *tri;
  290. Vector *pos, *norm, *dpdu, *dpdv;
  291. Vec2d *uv;
  292. {
  293.     Float d;
  294.  
  295.     /*
  296.      * Normalize barycentric coordinates.
  297.      */
  298.     d = tri->b[0]+tri->b[1]+tri->b[2];
  299.  
  300.     tri->b[0] /= d;
  301.     tri->b[1] /= d; 
  302.     tri->b[2] /= d;
  303.  
  304.     if (dpdu) {
  305.         if (tri->uv == (Vec2d *)NULL) {
  306.             *dpdu = tri->e[0];
  307.             (void)VecNormalize(dpdu);
  308.             VecSub(tri->p[0], *pos, dpdv);
  309.             (void)VecNormalize(dpdv);
  310.         } else {
  311.             *dpdu = *tri->dpdu;
  312.             *dpdv = *tri->dpdv;
  313.         }
  314.     }
  315.  
  316.     if (tri->uv == (Vec2d *)NULL) {
  317.         uv->v = tri->b[2];
  318.         if (equal(uv->v, 1.))
  319.             uv->u = 0.;
  320.         else
  321.             uv->u = tri->b[1] / (tri->b[0] + tri->b[1]);
  322.     } else {
  323.         /*
  324.          * Compute UV by taking weighted sum of UV coordinates.
  325.          */
  326.         uv->u = tri->b[0]*tri->uv[0].u + tri->b[1]*tri->uv[1].u +
  327.             tri->b[2]*tri->uv[2].u;
  328.         uv->v = tri->b[0]*tri->uv[0].v + tri->b[1]*tri->uv[1].v +
  329.             tri->b[2]*tri->uv[2].v;
  330.     }
  331. }
  332.  
  333. void
  334. TriangleBounds(tri, bounds)
  335. Triangle *tri;
  336. Float bounds[2][3];
  337. {
  338.     bounds[LOW][X] = bounds[HIGH][X] = tri->p[0].x;
  339.     bounds[LOW][Y] = bounds[HIGH][Y] = tri->p[0].y;
  340.     bounds[LOW][Z] = bounds[HIGH][Z] = tri->p[0].z;
  341.  
  342.     if (tri->p[1].x < bounds[LOW][X]) bounds[LOW][X] = tri->p[1].x;
  343.     if (tri->p[1].x > bounds[HIGH][X]) bounds[HIGH][X] = tri->p[1].x;
  344.     if (tri->p[2].x < bounds[LOW][X]) bounds[LOW][X] = tri->p[2].x;
  345.     if (tri->p[2].x > bounds[HIGH][X]) bounds[HIGH][X] = tri->p[2].x;
  346.  
  347.     if (tri->p[1].y < bounds[LOW][Y]) bounds[LOW][Y] = tri->p[1].y;
  348.     if (tri->p[1].y > bounds[HIGH][Y]) bounds[HIGH][Y] = tri->p[1].y;
  349.     if (tri->p[2].y < bounds[LOW][Y]) bounds[LOW][Y] = tri->p[2].y;
  350.     if (tri->p[2].y > bounds[HIGH][Y]) bounds[HIGH][Y] = tri->p[2].y;
  351.  
  352.     if (tri->p[1].z < bounds[LOW][Z]) bounds[LOW][Z] = tri->p[1].z;
  353.     if (tri->p[1].z > bounds[HIGH][Z]) bounds[HIGH][Z] = tri->p[1].z;
  354.     if (tri->p[2].z < bounds[LOW][Z]) bounds[LOW][Z] = tri->p[2].z;
  355.     if (tri->p[2].z > bounds[HIGH][Z]) bounds[HIGH][Z] = tri->p[2].z;
  356. }
  357.  
  358. char *
  359. TriangleName()
  360. {
  361.     return triName;
  362. }
  363.  
  364. void
  365. TriangleStats(tests, hits)
  366. unsigned long *tests, *hits;
  367. {
  368.     *tests = TriTests;
  369.     *hits = TriHits;
  370. }
  371.  
  372. /*
  373.  * Given three vertices of a triangle and the uv coordinates associated
  374.  * with each, compute directions of u and v axes.
  375.  */
  376. static void
  377. TriangleSetdPdUV(p, t, dpdu, dpdv)
  378. Vector p[3];            /* Triangle vertices */
  379. Vec2d t[3];            /* uv coordinates for each vertex */
  380. Vector *dpdu, *dpdv;        /* u and v axes (return values) */
  381. {
  382.     Float scale;
  383.     int hi, mid, lo;
  384.     Vector base;
  385.  
  386.     /* sort u coordinates */
  387.     if (t[2].u > t[1].u) {
  388.         if (t[1].u > t[0].u) {
  389.             hi = 2; mid = 1; lo = 0;
  390.         } else if (t[2].u > t[0].u) {
  391.             hi = 2; mid = 0; lo = 1;
  392.         } else {
  393.             hi = 0; mid = 2; lo = 1;
  394.         }
  395.     } else {
  396.         if (t[2].u > t[0].u) {
  397.             hi = 1; mid = 2; lo = 0;
  398.         } else if (t[1].u > t[0].u) {
  399.             hi = 1; mid = 0; lo = 2;
  400.         } else {
  401.             hi = 0; mid = 1; lo = 2;
  402.         }
  403.     }
  404.     if (t[hi].u - t[lo].u == 0.) {
  405.         /* degenerate axis */
  406.         dpdv->x = dpdv->y = dpdv->z = 0.;
  407.     } else {
  408.         /*
  409.          * Given u coordinates of vertices forming the
  410.          * 'long' edge, find where 'middle'
  411.          * vertex falls on that edge given its u coordinate.
  412.          */
  413.         scale = (t[mid].u - t[lo].u) / (t[hi].u - t[lo].u);
  414.         VecComb(1.0 - scale, p[lo], scale, p[hi], &base);
  415.         /*
  416.          * v axis extends from computed basepoint to
  417.          * middle vertex -- but in which direction?
  418.          */
  419.         if (t[mid].v < ((1.0 - scale)*t[lo].v + scale*t[hi].v))
  420.             VecSub(base, p[mid], dpdv);
  421.         else
  422.             VecSub(p[mid], base, dpdv);
  423.         (void)VecNormalize(dpdv);
  424.     }
  425.  
  426.     /* sort v coordinates */
  427.     if (t[2].v > t[1].v) {
  428.         if (t[1].v > t[0].v) {
  429.             hi = 2; mid = 1; lo = 0;
  430.         } else if (t[2].v > t[0].v) {
  431.             hi = 2; mid = 0; lo = 1;
  432.         } else {
  433.             hi = 0; mid = 2; lo = 1;
  434.         }
  435.     } else {
  436.         if (t[2].v > t[0].v) {
  437.             hi = 1; mid = 2; lo = 0;
  438.         } else if (t[1].v > t[0].v) {
  439.             hi = 1; mid = 0; lo = 2;
  440.         } else {
  441.             hi = 0; mid = 1; lo = 2;
  442.         }
  443.     }
  444.     if (t[hi].v - t[lo].v == 0.) {
  445.         /* degenerate axis */
  446.         dpdu->x = dpdu->y = dpdu->z = 0.;
  447.     } else {
  448.         /*
  449.          * Given v coordinates of vertices forming the
  450.          * 'long' edge, find where 'middle'
  451.          * vertex falls on that edge given its v coordinate.
  452.          */
  453.         scale = (t[mid].v - t[lo].v) / (t[hi].v - t[lo].v);
  454.         VecComb(1.0 - scale, p[lo], scale, p[hi], &base);
  455.         /*
  456.          * u axis extends from computed basepoint to
  457.          * middle vertex -- but in which direction?
  458.          */
  459.         if (t[mid].u < ((1.0 - scale)*t[lo].u + scale*t[hi].u))
  460.             VecSub(base, p[mid], dpdu);
  461.         else
  462.             VecSub(p[mid], base, dpdu);
  463.         (void)VecNormalize(dpdu);
  464.     }
  465. }
  466.  
  467. void
  468. TriangleMethodRegister(meth)
  469. UserMethodType meth;
  470. {
  471.     if (iTriangleMethods)
  472.         iTriangleMethods->user = meth;
  473. }
  474.