home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / torus.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  10.1 KB  |  474 lines
  1. /*
  2.  * torus.c -- this implements toruses.
  3.  * 
  4.  * (c) 1993, 1994 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 <math.h>
  21. #include <stdio.h>
  22.  
  23. #include "ray.h"
  24. #include "proto.h"
  25. #include "extern.h"
  26.  
  27.  
  28. extern struct methods my_methods;
  29.  
  30. /* give back memory */
  31. PRIVATE void
  32. free_torus(object *o)
  33. {
  34.     struct torus_data *t = o->data.torus;
  35.  
  36.     free((void *) t);
  37.     o->type = NOSHAPE;
  38. }
  39.  
  40. /*
  41.  * inside_xxx, this should be pretty optimized. It is used in every CSG
  42.  * test. (it could be better.. )
  43.  * 
  44.  * a torus eqn looks like factor^2 - stuff = 0.
  45.  */
  46.  
  47. PRIVATE bool
  48. inside_torus(object *o, vector i)
  49. {
  50.     struct torus_data *t = o->data.torus;
  51.     double          factor,
  52.                     stuff;
  53.     vector          loc;
  54.  
  55.     if (o->inv_trans)
  56.     loc = mvproduct(*o->inv_trans, i);
  57.  
  58.     factor = vdot(loc, loc) + sqr(t->Rad) - sqr(t->rad);
  59.     loc.y = 0;
  60.     stuff = 4 * sqr(t->Rad) * vdot(loc, loc);
  61.  
  62.     if (sqr(factor) - stuff < 0)
  63.     return !o->inverted;
  64.     else
  65.     return o->inverted;
  66. }
  67.  
  68. PRIVATE void
  69. print_torus(object *o)
  70. {
  71. #ifdef DEBUG
  72.     struct torus_data *p = o->data.torus;
  73.  
  74.     printf("Major radius %lf, Minor radius %lf\n", p->Rad, p->rad);
  75. #endif
  76. }
  77.  
  78.  
  79. /*
  80.  * the bounding was inspired on Graphics Gems II, elliptic torus solver.
  81.  */
  82. PRIVATE int
  83. intersect_torus(struct torus_data *tor, struct ray r, double *rhits, bool chkall)
  84. {
  85.     int             nhits;
  86.     double          yin,
  87.                     yout,
  88.                     brad;
  89.     spoly           tosolve;
  90.     double          sphere_hits[2];
  91.     vector          tmpvector;
  92.  
  93.  
  94.     my_methods.test++;
  95.  
  96.     /* Transform the intersection ray */
  97.     brad = tor->Rad + tor->rad + BOUNDFUDGE;
  98.     brad = sqr(brad);
  99.  
  100.     setvector(tmpvector, 0, 0, 0);
  101.  
  102.     /*
  103.      * Compute the intersection of the ray with a real minimum bounding
  104.      * sphere.
  105.      */
  106.     nhits = intersect_sphere(tmpvector, brad, &r, sphere_hits);
  107.  
  108.     if (nhits == 0)
  109.     return 0;
  110.  
  111.     if (nhits == 2) {
  112.     if (sphere_hits[0] > r.maxt)    /* something has been found
  113.                      * closer. Early abort. */
  114.         return 0;
  115.  
  116.     /*
  117.      * Bound the torus by two parallel planes
  118.      */
  119.     yin = r.pos.y + sphere_hits[0] * r.dir.y;
  120.     yout = r.pos.y + sphere_hits[1] * r.dir.y;
  121.  
  122.     if ((yin > tor->rad && yout > tor->rad) || (yin < -tor->rad && yout < -tor->rad))
  123.         return 0;
  124.  
  125.     /*
  126.      * translate the ray along it's direction. This should keep
  127.      * numerical inaccuracies a bit under control.
  128.      */
  129.  
  130.     /* does it,  actually? */
  131.     }
  132.     /* let the calculations begin! */
  133.     {
  134.     double          R,
  135.                     mr,
  136.                     quad1,
  137.                     lin1,
  138.                     cons1,
  139.                     quad0,
  140.                     lin0,
  141.                     cons0;
  142.     vector          pos,
  143.                     dir;
  144.  
  145.     /*
  146.      * svproduct(pos, baset, r.dir); vadd(pos, r.pos, pos);
  147.      */
  148.     pos = r.pos;
  149.  
  150.     dir = r.dir;
  151.  
  152.     mr = sqr(tor->rad);
  153.     R = sqr(tor->Rad);
  154.  
  155.     /*
  156.      * the equation is:
  157.      * 
  158.      * (z^2 + x^2 + y^2 + R - mr)^2 = 4R(x^2 + z^2)
  159.      * 
  160.      * note that R = Rad^2, and mr = rad^2
  161.      * 
  162.      * Calculate what happens if the ray is subst in (z^2 + x^2 + y^2 + R
  163.      * - r), which is vdot(X,X) + R-r.
  164.      */
  165.  
  166.     cons0 = vdot(pos, pos) + R - mr;
  167.     lin0 = 2 * vdot(pos, dir);
  168.     quad0 = vdot(dir, dir);
  169.  
  170.     /*
  171.      * now subst in right part:
  172.      * 
  173.      * 4R(x^2+y^2), which is 4R*vdot(X,X) with X.y = 0
  174.      */
  175.  
  176.     pos.y = dir.y = 0;
  177.     cons1 = R * 4 * vdot(pos, pos);
  178.     lin1 = R * 8 * vdot(pos, dir);
  179.     quad1 = R * 4 * vdot(dir, dir);
  180.  
  181.     /*
  182.      * now work out
  183.      * 
  184.      * (cons0 + t * lin0 + t^2 quad0) ^2 - (cons1 + t * lin1 + t^2 quad1)
  185.      */
  186.  
  187.     tosolve.coef[0] = sqr(cons0) - cons1;
  188.     tosolve.coef[1] = 2 * cons0 * lin0 - lin1;
  189.     tosolve.coef[2] = sqr(lin0) + 2 * cons0 * quad0 - quad1;
  190.     tosolve.coef[3] = 2 * lin0 * quad0;
  191.     tosolve.coef[4] = sqr(quad0);
  192.     tosolve.ord = 4;
  193.     }
  194.  
  195.     /* and solve it. */
  196.     if (tor->use_sturm)
  197.     nhits = sturm_solve_rt_poly(&tosolve, rhits, chkall);
  198.     else
  199.     nhits = solve_rt_poly(&tosolve, rhits, chkall);
  200.  
  201.     /*
  202.      * for (i = 0; i < nhits; i++) rhits[i] += ;
  203.      */
  204.  
  205.     if (nhits) {
  206.     my_methods.hit++;
  207.     }
  208.     return nhits;
  209. }
  210.  
  211.  
  212. PRIVATE bool
  213. all_torus_intersections(dqueue * q, object *o, struct ray *r,
  214.             int flags, bool *isinside)
  215. {
  216.     int             n,
  217.                     j;
  218.     double          inter[4];
  219.     dqueue          i;
  220.     bool            loc_ins;
  221.     struct ray      localray;
  222.  
  223.     localray = *r;
  224.     transform_ray(&localray, o);
  225.  
  226.     n = intersect_torus(o->data.torus, localray, inter, flags & CHKALL);
  227.     loc_ins = *isinside = (n % 2) ^ o->inverted;
  228.  
  229.     if (!n)
  230.     return FALSE;
  231.  
  232.  
  233.     if (flags & CHKALL && n > 1)/* check everything? */
  234.     for (j = 0; j < n; j++) {
  235.         i.t = inter[j];
  236.         i.obj = o;
  237.         i.entering = (loc_ins = !loc_ins);
  238.         add_to_queue(q, i);
  239.     } else {
  240.     /* record only one hit. */
  241.  
  242.     i.t = inter[0];
  243.     i.obj = o;
  244.     i.entering = (loc_ins = !loc_ins);
  245.     add_to_queue(q, i);
  246.     }
  247.  
  248.     return TRUE;
  249. }
  250.  
  251. /*
  252.  * returns the normal to the torus
  253.  * 
  254.  * The normal has the same direction as a vector from the loc's projection on
  255.  * the torus axis to the intersection point itself. Of course you could
  256.  * also take the gradient of the quartic formula, and then normalize it,
  257.  * but it's a load of work. This is simpler!
  258.  */
  259. PRIVATE vector
  260. torus_normal(struct intersect i, vector loc)
  261. {
  262.     struct torus_data *t;
  263.     vector          normal,
  264.                     heart_proj;
  265.     object         *o = i.q_ent.obj;
  266.  
  267.  
  268.     t = o->data.torus;
  269.  
  270.     if (o->inv_trans)
  271.     loc = mvproduct(*o->inv_trans, loc);
  272.  
  273.     heart_proj = loc;
  274.     heart_proj.y = 0;
  275.     norm(heart_proj, heart_proj);
  276.     svproduct(heart_proj, t->Rad, heart_proj);
  277.  
  278.     /* heart_proj now is the projection of loc on the tube axis. */
  279.     vsub(normal, loc, heart_proj);
  280.  
  281.     if (o->inv_trans)
  282.     normal = transform_normal(*o->inv_trans, normal);
  283.  
  284.     if (quickveclen(normal) < EPSILON && (debug_options & DEBUGRUNTIME))
  285.     warning("Weird torus!");
  286.     else
  287.     norm(normal, normal);
  288.  
  289.     return normal;
  290. }
  291.  
  292.  
  293. /* initialize a torus_data struct */
  294. PUBLIC void
  295. init_torus(struct torus_data *t)
  296. {
  297.     t->Rad = 1.0;
  298.     t->rad = 0.5;
  299.     t->use_sturm = FALSE;
  300. }
  301.  
  302. /* alloc and init torus */
  303. PUBLIC struct torus_data *
  304. get_new_torus(void)
  305. {
  306.     struct torus_data *p;
  307.     p = ALLOC(struct torus_data);
  308.  
  309.     CHECK_MEM(p, my_methods.name);
  310.     init_torus(p);
  311.     return p;
  312. }
  313.  
  314. /* copy the shape data */
  315. PRIVATE void
  316. copy_torus(object *dst, object *src)
  317. {
  318.     struct torus_data *s,
  319.                    *d;
  320.  
  321.     assert(dst != NULL && src != NULL);
  322.  
  323.     if (dst->type != TORUS)
  324.     dst->data.torus = get_new_torus();
  325.     *(d = dst->data.torus) = *(s = src->data.torus);
  326.     dst->type = src->type;
  327. }
  328.  
  329.  
  330. /* bounding volume calculation from Gems III */
  331. /*
  332.  * I did a little testing, it seems to work
  333.  * 
  334.  * 
  335.  * The bounding volume of a torus is found.  The algorithm performs three
  336.  * passes, one for each dimension.  In each pass the extrema of the
  337.  * surface are found as values of the parameters u and v.  These values
  338.  * are then used to calculate the position of the two points in the
  339.  * current dimension. Finally, these points are sorted to find the minimum
  340.  * and maximum extents. The torus is defined as the rotation of a circle
  341.  * that is perpendicular to the XZ plane.  This circle has radius q, and
  342.  * its center lies in the XZ plane and is rotated about the Y axis in a
  343.  * 
  344.  * circle of radius r.  The value of q must be less than that of r.
  345.  */
  346.  
  347. PRIVATE void
  348. precompute_torus(object *o)
  349. {
  350.     struct torus_data *t = o->data.torus;
  351.     double          min[3],
  352.                     max[3];    /* returned minimum and maximum of extent */
  353.     matrix          ctm,
  354.                     trans;    /* cumulative transformation matrix */
  355.     double          r;        /* major radius of torus */
  356.     double          q;        /* minor radius of torus */
  357.     int             i;
  358.     double          u1,
  359.                     u2,
  360.                     v1,
  361.                     v2,
  362.                     tmp,
  363.                     denominator;
  364.  
  365.  
  366.  
  367.     /* do bounding */
  368.     q = t->rad;
  369.     r = t->Rad;
  370.     if (o->inv_trans) {
  371.     invert_trans(trans, *o->inv_trans);
  372.  
  373.     /* Gems do everything transposed. */
  374.     transpose_matrix(ctm, trans);
  375.     } else
  376.     unit_matrix(ctm);
  377.  
  378.     for (i = 0; i < 3; i++) {
  379.  
  380.     /* calculate first extremum.  assure that -PI/2 <= v1 <= PI/2 */
  381.     u1 = safe_arctangent(ctm[2][i], ctm[0][i]);
  382.     denominator = ctm[0][i] * cos(u1) + ctm[2][i] * sin(u1);
  383.     v1 = safe_arctangent(ctm[1][i], denominator);
  384.  
  385.     /* second extremum is +/- PI from u1, negative of v1 */
  386.     if (u1 <= 0)
  387.         u2 = u1 + M_PI;
  388.     else
  389.         u2 = u1 - M_PI;
  390.     v2 = -v1;
  391.  
  392.     /* find and sort extrema locations in this dimension */
  393.     min[i] =
  394.         ctm[0][i] * cos(u1) * (r + q * cos(v1)) +
  395.         ctm[1][i] * sin(v1) * q +
  396.         ctm[2][i] * sin(u1) * (r + q * cos(v1)) +
  397.         ctm[3][i];
  398.     max[i] =
  399.         ctm[0][i] * cos(u2) * (r + q * cos(v2)) +
  400.         ctm[1][i] * sin(v2) * q +
  401.         ctm[2][i] * sin(u2) * (r + q * cos(v2)) +
  402.         ctm[3][i];
  403.     if (min[i] > max[i]) {
  404.         tmp = max[i];
  405.         max[i] = min[i];
  406.         min[i] = tmp;
  407.     }
  408.     }
  409.  
  410.     o->bmin.x = min[0];
  411.     o->bmin.y = min[1];
  412.     o->bmin.z = min[2];
  413.     o->bmax.x = max[0];
  414.     o->bmax.y = max[1];
  415.     o->bmax.z = max[2];
  416.  
  417.     my_methods.howmuch++;
  418.     global_stats.prims++;
  419. }
  420.  
  421. PRIVATE struct methods my_methods =
  422. {
  423.     all_torus_intersections,
  424.     torus_normal,
  425.     copy_torus,
  426.     inside_torus,
  427.     generic_rotate_object,
  428.     generic_translate_object,
  429.     generic_scale_object,
  430.     free_torus,
  431.     print_torus,
  432.     precompute_torus,
  433.     "torus",
  434. };
  435.  
  436.  
  437. /* alloc/init */
  438. PUBLIC object  *
  439. get_new_torus_object(void)
  440. {
  441.     object         *o;
  442.  
  443.     o = get_new_object();
  444.  
  445.     o->data.torus = get_new_torus();
  446.     o->methods = &my_methods;
  447.     o->type = TORUS;
  448.  
  449.     return o;
  450. }
  451.  
  452. /* eof */
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473. /* not really.... But who actually reads this crap? :) */
  474.