home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / sphere.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  8.1 KB  |  376 lines
  1. /*
  2.  * sphere.c -- standard sphere routines.
  3.  * 
  4.  * (c) 1993, 1994 by Han-Wen Nienhuys <hanwen@stack.urc.tue.nl>
  5.  * 
  6.  * This is based on work by George Kyriazis.
  7.  * 
  8.  * This program is free software; you can redistribute it and/or modify it
  9.  * under the terms of the GNU General Public License as published by the
  10.  * Free Software Foundation;
  11.  * 
  12.  * This program is distributed in the hope that it will be useful, but
  13.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  14.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15.  * General Public License for more details.
  16.  * 
  17.  * You should have received a copy of the GNU General Public License along
  18.  * with this program; if not, write to the Free Software Foundation, Inc.,
  19.  * 675 Mass Ave, Cambridge, MA 02139, USA.
  20.  */
  21.  
  22. #include "ray.h"
  23. #include "proto.h"
  24. #include "extern.h"
  25.  
  26. extern struct methods my_methods;
  27.  
  28.  
  29. PRIVATE bool
  30. all_sphere_intersections(dqueue * q, object *o, struct ray *r,
  31.              int flags,    /* give all intersections */
  32.              bool *isinside)
  33. {                /* for use with CSG: isinside ==
  34.                  * inside_sphere() */
  35.     dqueue          q_ent;
  36.     int             no_int;
  37.     double          d[2];
  38.  
  39.  
  40.     if (o->inv_trans != NULL) {
  41.     struct ray      localray;
  42.  
  43.     localray = *r;
  44.     transform_ray(&localray, o);
  45.     no_int = intersect_sphere(o->data.sphere->center, o->data.sphere->rsq, &localray, d);
  46.     } else
  47.     no_int = intersect_sphere(o->data.sphere->center, o->data.sphere->rsq, r, d);
  48.  
  49.     *isinside = (no_int % 2) ^ o->inverted;
  50.  
  51.     if (!no_int)
  52.     return FALSE;
  53.  
  54.  
  55.     q_ent.t = d[0];
  56.     q_ent.obj = o;
  57.     q_ent.entering = !(*isinside);
  58.     add_to_queue(q, q_ent);
  59.  
  60.     if ((flags & CHKALL) && no_int == 2) {
  61.     q_ent.t = d[1];
  62.     q_ent.entering = (*isinside);
  63.     add_to_queue(q, q_ent);
  64.     }
  65.     return TRUE;
  66. }
  67.  
  68. /* Intersect ray with sphere */
  69. PUBLIC int
  70. intersect_sphere(vector center, double rsq, struct ray *r, double *hits)
  71. {
  72.     vector          center_direction;
  73.     double          b,
  74.                     c,
  75.                     d,
  76.                     dir_projection,
  77.                     dir_len,
  78.                     cd_lensq;
  79.     double          solmax,
  80.                     solmin;
  81.  
  82.  
  83.     my_methods.test++;
  84.  
  85.     /*
  86.      * solve: x = (pos+t*dir) and |x-center|^2  = radius^2
  87.      */
  88.     vsub(center_direction, r->pos, center);    /* vector from center to
  89.                          * r.pos */
  90.     dir_projection = vdot(r->dir, center_direction);
  91.     cd_lensq = vdot(center_direction, center_direction);    /* square of length of
  92.                                  * center_direction */
  93.  
  94.     if (cd_lensq >= rsq &&    /* we're not inside the sphere */
  95.     dir_projection > 0) {    /* and we're not looking at the sphere */
  96.     return 0;
  97.     }
  98.     b = 2 * dir_projection;
  99.     c = cd_lensq - rsq;
  100.  
  101.     /* with the advent of extrusions, |r.dir| isn't always 1 */
  102.     dir_len = vdot(r->dir, r->dir);
  103.     d = b * b - 4 * dir_len * c;
  104.  
  105.     if (d < 0)
  106.     return 0;
  107.  
  108.     dir_len *= 2;
  109.     d = sqrt(d);
  110.     solmax = (-b + d) / (dir_len);
  111.     solmin = (-b - d) / (dir_len);
  112.  
  113.     if (solmax < tolerance)
  114.     return 0;        /* maximum is  behind us. Shake it. */
  115.  
  116.     /* This is an intersection! */
  117.     my_methods.hit++;
  118.     if (solmin < tolerance) {
  119.     *hits = solmax;
  120.     return 1;
  121.     } else {
  122.     *hits++ = solmin;
  123.     *hits = solmax;
  124.     return 2;
  125.     }
  126.  
  127. }
  128.  
  129. PRIVATE vector
  130. sphere_normal(struct intersect i, vector x)
  131. {
  132.     vector          n;
  133.     object         *o = i.q_ent.obj;
  134.  
  135.  
  136.     if (o->inv_trans != NULL)
  137.     x = mvproduct(*o->inv_trans, x);
  138.  
  139.  
  140.     /* calculate the normal.  It is just the direction of the radius */
  141.     /* norm(i.x - cntr) */
  142.     vsub(n, x, o->data.sphere->center);
  143.  
  144.     if (o->inv_trans != NULL)
  145.     n = transform_normal(*o->inv_trans, n);
  146.  
  147.     norm(n, n);
  148.  
  149.     return n;
  150. }
  151.  
  152. PRIVATE void
  153. precompute_sphere(object *o)
  154. {
  155.     struct sphere_data *sph = o->data.sphere;
  156.  
  157.     sph->rsq = sqr(sph->radius);
  158.     sph->radius = ABS(sph->radius);
  159.  
  160.     /* bounding volume calculation */
  161.  
  162.     /*
  163.      * The bounding volume of a sphere is found.  The algorithm performs
  164.      * three passes, one for each dimension.  In each pass the extrema of
  165.      * the surface are found as values of the parameters u and v.  These
  166.      * values are then used to calculate the position of the two points in
  167.      * the current dimension. Finally, these points are sorted to find the
  168.      * minimum and maximum extents. The canonical sphere has a radius of
  169.      * 1.0 and is centered at the origin.
  170.      */
  171.  
  172.  
  173.  
  174.     {
  175.     double          min[3],
  176.                     max[3];    /* returned minimum and maximum of extent */
  177.     matrix          ctm,
  178.                     tempmatrix;    /* cumulative transformation
  179.                      * matrix */
  180.  
  181.     int             i;
  182.     float           u1,
  183.                     u2,
  184.                     v1,
  185.                     v2,
  186.                     tmp,
  187.                     denominator;
  188.     vector          scal;
  189.  
  190.     setvector(scal, sph->radius, sph->radius, sph->radius);
  191.     scale_matrix(ctm, scal);
  192.     translate_matrix(tempmatrix, sph->center);
  193.     mmproduct(tempmatrix, tempmatrix, ctm);
  194.  
  195.     if (o->inv_trans) {
  196.         invert_trans(tempmatrix, *o->inv_trans);
  197.         mmproduct(tempmatrix, tempmatrix, ctm);
  198.     }
  199.     transpose_matrix(ctm, tempmatrix);
  200.  
  201.     for (i = 0; i < 3; i++) {
  202.  
  203.         /* calculate first extremum. */
  204.         u1 = safe_arctangent(ctm[2][i], ctm[0][i]);
  205.         denominator = ctm[0][i] * cos(u1) + ctm[2][i] * sin(u1);
  206.         v1 = safe_arctangent(ctm[1][i], denominator);
  207.  
  208.         /* second extremum is +/- PI from u1, negative of v1 */
  209.         if (u1 <= 0)
  210.         u2 = u1 + M_PI;
  211.         else
  212.         u2 = u1 - M_PI;
  213.         v2 = -v1;
  214.  
  215.         /* find and sort extrema locations in this dimension */
  216.         min[i] =
  217.         ctm[0][i] * cos(u1) * cos(v1) +
  218.         ctm[1][i] * sin(v1) +
  219.         ctm[2][i] * sin(u1) * cos(v1) +
  220.         ctm[3][i];
  221.         max[i] =
  222.         ctm[0][i] * cos(u2) * cos(v2) +
  223.         ctm[1][i] * sin(v2) +
  224.         ctm[2][i] * sin(u2) * cos(v2) +
  225.         ctm[3][i];
  226.         if (min[i] > max[i]) {
  227.         tmp = max[i];
  228.         max[i] = min[i];
  229.         min[i] = tmp;
  230.         }
  231.     }
  232.     o->bmin.x = min[0];
  233.     o->bmin.y = min[1];
  234.     o->bmin.z = min[2];
  235.     o->bmax.x = max[0];
  236.     o->bmax.y = max[1];
  237.     o->bmax.z = max[2];
  238.     }
  239.  
  240.     my_methods.howmuch++;
  241.     global_stats.prims++;
  242. }
  243.  
  244. PRIVATE bool
  245. inside_sphere(object *o, vector x)
  246. {
  247.     if (o->inv_trans)
  248.     x = mvproduct(*o->inv_trans, x);
  249.     vsub(x, x, o->data.sphere->center);
  250.     if (vdot(x, x) < o->data.sphere->rsq)
  251.     return !o->inverted;
  252.     else
  253.     return o->inverted;
  254. }
  255.  
  256. PRIVATE void
  257. translate_sphere(object *o, vector t)
  258. {
  259.     struct sphere_data *s;
  260.  
  261.     s = o->data.sphere;
  262.     if (o->inv_trans)
  263.     generic_translate_object(o, t);
  264.     else
  265.     vadd(s->center, s->center, t);
  266. }
  267.  
  268. PRIVATE void
  269. scale_sphere(object *o, vector s)
  270. {
  271.     struct sphere_data *sph;
  272.     double          factor;
  273.  
  274.     sph = o->data.sphere;
  275.  
  276.  
  277.     if (ABS(s.x) != ABS(s.y) || ABS(s.y) != ABS(s.z)) {
  278.     /* scale unevenly */
  279.     generic_scale_object(o, s);
  280.     } else {
  281.     factor = ABS(s.x);
  282.     sph->radius *= factor;
  283.     vcproduct(sph->center, s, sph->center);
  284.     }
  285. }
  286.  
  287. PRIVATE void
  288. rotate_sphere(object *o, matrix m)
  289. {
  290.     struct sphere_data *sph;
  291.  
  292.     sph = o->data.sphere;
  293.     if (o->inv_trans != NULL)
  294.     generic_rotate_object(o, m);
  295.     else
  296.     rotate_vector(&sph->center, m);
  297. }
  298.  
  299. PRIVATE void
  300. init_sphere(struct sphere_data *s)
  301. {
  302.     setvector(s->center, 0, 0, 0);
  303.     s->radius = 1.0;
  304. }
  305.  
  306. PRIVATE struct sphere_data *
  307. get_new_sphere(void)
  308. {
  309.     struct sphere_data *s;
  310.  
  311.     s = ALLOC(struct sphere_data);
  312.  
  313.     CHECK_MEM(s, my_methods.name);
  314.  
  315.     init_sphere(s);
  316.     return s;
  317. }
  318.  
  319. PRIVATE void
  320. copy_sphere(object *dst, object *src)
  321. {
  322.     assert(dst != NULL && src != NULL);
  323.     if (dst->type != SPHERE)
  324.     dst->data.sphere = get_new_sphere();
  325.     *dst->data.sphere = *src->data.sphere;
  326.     dst->type = src->type;
  327. }
  328.  
  329. PRIVATE void
  330. free_sphere(object *s)
  331. {
  332.     free((void *) s->data.sphere);
  333.     s->type = NOSHAPE;
  334. }
  335.  
  336. PRIVATE void
  337. print_sphere(object *o)
  338. {
  339. #ifdef DEBUG
  340.     struct sphere_data *s = o->data.sphere;
  341.  
  342.     printf("radius %lf (%f) ", s->radius, s->rsq);
  343.     print_v("cntr", s->center);
  344. #endif
  345. }
  346.  
  347.  
  348. PRIVATE struct methods my_methods =
  349. {
  350.     all_sphere_intersections,
  351.     sphere_normal,
  352.     copy_sphere,
  353.     inside_sphere,
  354.     rotate_sphere,
  355.     translate_sphere,
  356.     scale_sphere,
  357.     free_sphere,
  358.     print_sphere,
  359.     precompute_sphere,
  360.     "sphere",
  361. };
  362.  
  363.  
  364. PUBLIC object  *
  365. get_new_sphere_object(void)
  366. {
  367.     object         *o;
  368.  
  369.     o = get_new_object();
  370.     o->data.sphere = get_new_sphere();
  371.     o->methods = &my_methods;
  372.     o->type = SPHERE;
  373.  
  374.     return o;
  375. }
  376.