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