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