home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / triangle.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  11.1 KB  |  541 lines
  1. /*
  2.  * triangle.c -- handle smooth triangles.
  3.  * 
  4.  * (c) 1993, 1994 by Han-Wen Nienhuys <hanwen@stack.urc.tue.nl>
  5.  * 
  6.  * This program is free software; you can redistribute it and/or modify it
  7.  * under the terms of the GNU General Public License as published by the
  8.  * Free Software Foundation;
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License along
  16.  * with this program; if not, write to the Free Software Foundation, Inc.,
  17.  * 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include "ray.h"
  21. #include "proto.h"
  22. #include "extern.h"
  23.  
  24. extern struct methods my_methods;
  25.  
  26. PRIVATE void
  27. free_triangle(object *o)
  28. {
  29.     free((void *) o->data.triangle);
  30.     if (o->data.triangle->normals != NULL)
  31.     free(o->data.triangle->normals);
  32.     o->type = NOSHAPE;
  33. }
  34.  
  35. PRIVATE bool
  36. inside_triangle(object *o, vector i)
  37. {
  38.     assert(FALSE);
  39. }
  40.  
  41. /* scale it. */
  42. PRIVATE void
  43. scale_triangle(object *o, vector s)
  44. {
  45.     struct triangle_data *t = o->data.triangle;
  46.     vector          invs;
  47.     double          len;
  48.     int             i;
  49.  
  50.     setvector(invs, 1 / s.x, 1 / s.y, 1 / s.z);
  51.     for (i = 0; i < 3; i++) {
  52.     vcproduct(t->vertices[i], t->vertices[i], s);
  53.     if (t->normals != NULL) {
  54.  
  55.         /* preserve length */
  56.         len = veclen(t->normals[i]);
  57.         vcproduct(t->normals[i], t->normals[i], invs);
  58.         norm(t->normals[i], t->normals[i]);
  59.         svproduct(t->normals[i], len, t->normals[i]);
  60.     }
  61.     }
  62. }
  63.  
  64. /* rotation */
  65. PRIVATE void
  66. rotate_triangle(object *o, matrix rotmat)
  67. {
  68.     struct triangle_data *t = o->data.triangle;
  69.     int             i;
  70.  
  71.  
  72.     for (i = 0; i < 3; i++) {
  73.     t->vertices[i] = mvproduct(rotmat, t->vertices[i]);
  74.     if (t->normals != NULL)
  75.         t->normals[i] = mvproduct(rotmat, t->normals[i]);
  76.     }
  77. }
  78.  
  79. /* and translation */
  80. PRIVATE void
  81. translate_triangle(object *o, vector tr)
  82. {
  83.     int             i;
  84.     struct triangle_data *t = o->data.triangle;
  85.  
  86.     for (i = 0; i < 3; i++) {
  87.     vadd(t->vertices[i], tr, t->vertices[i]);
  88.     }
  89. }
  90.  
  91. PRIVATE void
  92. print_triangle(object *o)
  93. {
  94. #ifdef DEBUG
  95.     int             i;
  96.     struct triangle_data *p = o->data.triangle;
  97.  
  98.     for (i = 0; i < 3; i++) {
  99.     print_v("vertex", p->vertices[i]);
  100.     if (p->normals != NULL)
  101.         print_v("normal", p->normals[i]);
  102.     }
  103.     printf("projecting along %c\n", p->project_to);
  104.     print_v("linfunc1", p->linfunc1);
  105.     print_v("linfunc2", p->linfunc2);
  106.     print_v("normal", p->n);
  107. #endif
  108. }
  109.  
  110. PRIVATE bool
  111. all_triangle_intersections(dqueue * q, object *o, struct ray *r, int flags, bool *isinside)
  112. {
  113.     dqueue          q_ent;
  114.     vector          x;
  115.     double          coef1,
  116.                     coef2;
  117.  
  118.     struct triangle_data *s = o->data.triangle;
  119.  
  120.     q_ent.obj = o;
  121.  
  122.  
  123.     my_methods.test++;
  124.  
  125.  
  126.     {                /* do plane intersection */
  127.     double          d;
  128.  
  129.     d = vdot(r->dir, s->n);
  130.     if (ISZERO(d))
  131.         return 0;
  132.  
  133.     q_ent.t = (s->mov - vdot(r->pos, s->n)) / d;
  134.     }
  135.  
  136.  
  137.     /* try to get out lazily */
  138.     if (q_ent.t < tolerance || q_ent.t > r->maxt)
  139.     return FALSE;
  140.  
  141.     /*
  142.      * write x-vertex2  as coef1 * (vertex0-vertex1)  + coef2 *
  143.      * (vertex0-vertex2)
  144.      */
  145.  
  146.     /*
  147.      * do the intersection. What really happens is the projected version
  148.      * of:
  149.      */
  150. #ifdef undefined
  151.  
  152.     vadd(x, x, r->pos);
  153.     svproduct(x, q_ent.t, r->dir);
  154.     vsub(x, x, s->vertices[2]);
  155.  
  156.     coef1 = vdot(x, s->linfunc1);
  157.  
  158.     if (coef1 < 0)
  159.     return FALSE;
  160.  
  161.     coef2 = vdot(x, s->linfunc2);
  162.     if (coef2 < 0)
  163.     return FALSE;
  164. #endif
  165.  
  166.  
  167.     switch (s->project_to) {
  168.     case 'x':
  169.     x.x = 0.0;
  170.     x.y = r->pos.y + q_ent.t * r->dir.y;
  171.     x.z = r->pos.z + q_ent.t * r->dir.z;
  172.  
  173.     /* bounding box test. Inspired on GGems */
  174.     if ((x.y < o->bmin.y) || (x.y > o->bmax.y) || (x.z < o->bmin.z) || (x.z > o->bmax.z))
  175.         return FALSE;
  176.  
  177.     x.y = x.y - s->vertices[2].y;
  178.     x.z = x.z - s->vertices[2].z;
  179.     coef1 = x.y * s->linfunc1.y + x.z * s->linfunc1.z;
  180.     if (coef1 < 0)
  181.         return FALSE;
  182.  
  183.     coef2 = x.y * s->linfunc2.y + x.z * s->linfunc2.z;
  184.  
  185.     if (coef2 < 0)
  186.         return FALSE;
  187.     break;
  188.     case 'y':
  189.  
  190.     x.x = r->pos.x + q_ent.t * r->dir.x;
  191.     x.y = 0.0;
  192.     x.z = r->pos.z + q_ent.t * r->dir.z;
  193.  
  194.     /* bounding box test. Inspired on GGems */
  195.     if ((x.x < o->bmin.x) || (x.x > o->bmax.x) || (x.z < o->bmin.z) || (x.z > o->bmax.z))
  196.         return FALSE;
  197.     x.x = x.x - s->vertices[2].x;
  198.  
  199.     x.z = x.z - s->vertices[2].z;
  200.     coef1 = x.x * s->linfunc1.x + x.z * s->linfunc1.z;
  201.     if (coef1 < 0)
  202.         return FALSE;
  203.  
  204.     coef2 = x.x * s->linfunc2.x + x.z * s->linfunc2.z;
  205.  
  206.     if (coef2 < 0)
  207.         return FALSE;
  208.     break;
  209.  
  210.     case 'z':
  211.     x.x = r->pos.x + q_ent.t * r->dir.x;
  212.     x.y = r->pos.y + q_ent.t * r->dir.y;
  213.     x.z = 0.0;
  214.  
  215.     /* bounding box test. Inspired on GGems */
  216.     if ((x.x <
  217.        o->bmin.x) || (x.x > o->bmax.x) || (x.y < o->bmin.y) || (x.y >
  218.                                   o->bmax.y))
  219.         return FALSE;
  220.  
  221.     x.x = x.x - s->vertices[2].x;
  222.     x.y = x.y - s->vertices[2].y;
  223.  
  224.     coef1 = x.y * s->linfunc1.y + x.x * s->linfunc1.x;
  225.     if (coef1 < 0)
  226.         return FALSE;
  227.  
  228.     coef2 = x.y * s->linfunc2.y + x.x * s->linfunc2.x;
  229.  
  230.     if (coef2 < 0)
  231.         return FALSE;
  232.     break;
  233.     default:
  234.     assert(FALSE);
  235.     }
  236.  
  237. #ifdef DEBUG            /* should be removed */
  238.     if (debug_options & DEBUGRUNTIME) {
  239.     vector          t1,
  240.                     t2,
  241.                     t3;
  242.  
  243.     svproduct(x, q_ent.t, r->dir);
  244.     vadd(x, x, r->pos);
  245.  
  246.     svproduct(t1, coef1, s->vertices[0]);
  247.     svproduct(t2, coef2, s->vertices[1]);
  248.     svproduct(t3, 1 - coef1 - coef2, s->vertices[2]);
  249.     vadd(t1, t1, t2);
  250.     vadd(t1, t1, t3);
  251.     vsub(t1, t1, x);
  252.     if (!ISZERO(veclen(t1)))
  253.         assert(FALSE);
  254.     }
  255. #endif
  256.  
  257.     if (coef1 + coef2 > 1)
  258.     return FALSE;
  259.     else {
  260.     add_to_queue(q, q_ent);
  261.     my_methods.hit++;
  262.     return TRUE;
  263.     }
  264. }
  265.  
  266. /*
  267.  * returns the normal to the triangle
  268.  */
  269. PRIVATE vector
  270. triangle_normal(struct intersect i, vector loc)
  271. {
  272.     struct triangle_data *s = i.q_ent.obj->data.triangle;
  273.  
  274.     if (s->normals != NULL) {
  275.     double          coef1,
  276.                     coef2,
  277.                     coef3;
  278.     vector          normal,
  279.                     normaltemp;
  280.  
  281.  
  282.     vsub(loc, loc, s->vertices[2]);
  283.     switch (s->project_to) {
  284.     case 'x':
  285.         loc.x = 0;
  286.         break;
  287.     case 'y':
  288.         loc.y = 0;
  289.         break;
  290.     case 'z':
  291.         loc.z = 0;
  292.         break;
  293.     default:
  294.         assert(FALSE);
  295.     }
  296.  
  297.  
  298.  
  299.     /* decompose into base formed by two edges */
  300.     coef1 = vdot(loc, s->linfunc1);
  301.     coef2 = vdot(loc, s->linfunc2);
  302.     coef3 = 1 - coef1 - coef2;
  303.     svproduct(normal, coef1, s->normals[0]);
  304.     svproduct(normaltemp, coef2, s->normals[1]);
  305.     vadd(normal, normaltemp, normal);
  306.     svproduct(normaltemp, coef3, s->normals[2]);
  307.     vadd(normal, normaltemp, normal);
  308.     norm(normal, normal);
  309.     return normal;
  310.     } else {
  311.     return s->n;
  312.     }
  313. }
  314.  
  315. /* initialize a triangle_data struct */
  316. PRIVATE void
  317. init_triangle(struct triangle_data *t)
  318. {
  319.     setvector(t->vertices[0], 0, 0, 0);
  320.     setvector(t->vertices[1], 0, 1, 0);
  321.     setvector(t->vertices[2], 1, 0, 0);
  322.     setvector(t->linfunc1, 0, 0, 0);
  323.     setvector(t->linfunc2, 0, 0, 0);
  324.     t->project_to = ' ';
  325. }
  326.  
  327. PRIVATE void
  328. init_smooth_triangle(struct triangle_data *t)
  329. {
  330.     setvector(t->vertices[0], 0, 0, 0);
  331.     setvector(t->vertices[1], 0, 1, 0);
  332.     setvector(t->vertices[2], 1, 0, 0);
  333.     setvector(t->normals[0], 0, 0, 0);
  334.     setvector(t->normals[1], 0, 1, 0);
  335.     setvector(t->normals[2], 1, 0, 0);
  336.     setvector(t->linfunc1, 0, 0, 0);
  337.     setvector(t->linfunc2, 0, 0, 0);
  338.     t->project_to = ' ';
  339. }
  340.  
  341. PRIVATE struct triangle_data *
  342. get_new_triangle(void)
  343. {
  344.     struct triangle_data *p;
  345.     p = ALLOC(struct triangle_data);
  346.  
  347.     CHECK_MEM(p, my_methods.name);
  348.     init_triangle(p);
  349.     return p;
  350. }
  351.  
  352. /* copy the shape data */
  353. PRIVATE void
  354. copy_triangle(object *dst, object *src)
  355. {
  356.     assert(dst != NULL && src != NULL);
  357.  
  358.     if (dst->type != TRIANGLE)
  359.     dst->data.triangle = get_new_triangle();
  360.     *dst->data.triangle = *src->data.triangle;
  361.     dst->type = src->type;
  362. }
  363.  
  364. PRIVATE void
  365. precompute_triangle(object *o)
  366. {
  367.     vector          bas1,
  368.                     bas2;
  369.     struct triangle_data *s = o->data.triangle;
  370.  
  371.     vector          orth,
  372.                     bas1_norm;
  373.     double          alpha,
  374.                     normlen;
  375.     int             i;
  376.  
  377.  
  378.     vsub(bas1, s->vertices[0], s->vertices[2]);
  379.     vsub(bas2, s->vertices[1], s->vertices[2]);
  380.  
  381.     vcross(s->n, bas1, bas2);
  382.     normlen = veclen(s->n);
  383.     if (normlen < EPSILON)
  384.     warning("degenerate triangle");
  385.     else
  386.     svproduct(s->n, 1 / normlen, s->n);
  387.  
  388.     s->mov = vdot(s->n, s->vertices[2]);
  389.     /* project along which axis ? */
  390.  
  391.     {
  392.     double          max;
  393.  
  394.     if (ABS(s->n.x) > ABS(s->n.y)) {
  395.         max = ABS(s->n.x);
  396.         s->project_to = 'x';
  397.     } else {
  398.         max = ABS(s->n.y);
  399.         s->project_to = 'y';
  400.     }
  401.     if (ABS(s->n.z) > max) {
  402.         s->project_to = 'z';
  403.     }
  404.     }
  405.  
  406.     /*
  407.      * we want to express x in
  408.      * 
  409.      * a_1 x_1 + a_2 x_2 + a_3 x_3,
  410.      * 
  411.      * with a_1 + a_2+ a_3 = 1 this is equivalent to (x - x_3)  = a_1 ( x_1 -
  412.      * x_3) + a_2 ( x_2 -x_3) So we have to express x-x_3 in the basis {
  413.      * (x_1-x_3), (x_2-x_3) }
  414.      * 
  415.      * The theory of finite dimensional vectorspaces states, that the
  416.      * coefficients a_1 and a_2 can be found by taking an improduct of
  417.      * suitable linfunc1 and linfunc2 with x - x_3. These linfunc[12] are
  418.      * computed here:
  419.      * 
  420.      * since this is a 2D problem, we might as well project the intersection
  421.      * and the vertices onto one of the coordinate planes.
  422.      * 
  423.      * first we compute orth, a vector orthogonal to bas1 = x1 - x3
  424.      */
  425.  
  426.  
  427.     switch (s->project_to) {
  428.     case 'x':
  429.     bas1.x = bas2.x = 0;
  430.     break;
  431.     case 'y':
  432.     bas1.y = bas2.y = 0;
  433.     break;
  434.     case 'z':
  435.     bas1.z = bas2.z = 0;
  436.     break;
  437.     default:
  438.     assert(FALSE);
  439.     }
  440.  
  441.     alpha = -vdot(bas1, bas2) / vdot(bas1, bas1);
  442.     svproduct(orth, alpha, bas1);
  443.     vadd(orth, orth, bas2);
  444.  
  445.     svproduct(s->linfunc2, 1 / vdot(orth, orth), orth);
  446.     svproduct(s->linfunc1, alpha, s->linfunc2);
  447.     svproduct(bas1_norm, 1 / vdot(bas1, bas1), bas1);
  448.     vadd(s->linfunc1, bas1_norm, s->linfunc1);
  449.  
  450.  
  451.     for (i = 0; i < 3; i++)
  452.     update_min_and_max(&o->bmin, &o->bmax, s->vertices[i]);
  453.  
  454.     my_methods.howmuch++;
  455.     global_stats.prims++;
  456. }
  457.  
  458. PRIVATE struct methods my_methods =
  459. {
  460.     all_triangle_intersections,
  461.     triangle_normal,
  462.     copy_triangle,
  463.     inside_triangle,
  464.     rotate_triangle,
  465.     translate_triangle,
  466.     scale_triangle,
  467.     free_triangle,
  468.     print_triangle,
  469.     precompute_triangle,
  470.     "triangle",
  471. };
  472.  
  473.  
  474. /* alloc/init */
  475. PUBLIC object  *
  476. get_new_triangle_object(void)
  477. {
  478.     object         *o;
  479.  
  480.     o = get_new_object();
  481.  
  482.     o->data.triangle = ALLOC(struct triangle_data);
  483.  
  484.     CHECK_MEM(o->data.triangle, my_methods.name);
  485.  
  486.     init_triangle(o->data.triangle);
  487.     o->methods = &my_methods;
  488.     o->type = TRIANGLE;
  489.  
  490.     return o;
  491. }
  492.  
  493. /* alloc/init */
  494. PUBLIC object  *
  495. get_new_smooth_triangle_object(void)
  496. {
  497.     object         *o;
  498.     struct triangle_data *p;
  499.  
  500.     o = get_new_object();
  501.  
  502.     p = ALLOC(struct triangle_data);
  503.  
  504.     CHECK_MEM(p, my_methods.name);
  505.     p->normals = malloc(3 * sizeof(vector));
  506.  
  507.     CHECK_MEM(p->normals, "triangle normals");
  508.     o->data.triangle = p;
  509.  
  510.     o->methods = &my_methods;
  511.     o->type = TRIANGLE;
  512.     init_smooth_triangle(o->data.triangle);
  513.  
  514.     return o;
  515. }
  516.  
  517.  
  518.  
  519. /* eof */
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540. /* not really.... But who actually reads this crap? :) */
  541.