home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / quadric.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  8.2 KB  |  374 lines
  1. /*
  2.  * quadric.c -- standard quadratic shaperoutines.
  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.  
  21. #include "ray.h"
  22. #include "proto.h"
  23. #include "extern.h"
  24.  
  25. extern struct methods my_methods;
  26.  
  27.  
  28. PRIVATE bool    inside_quadric(object *o, vector x);
  29. PUBLIC struct quadric_data *get_new_quadric(void);
  30. PRIVATE int     intersect_quadric(struct quadric_data *qp, struct ray *r, double *roots);
  31.  
  32. PRIVATE bool
  33. all_quadric_intersections(dqueue * q, object *o, struct ray *r,
  34.               int flags, bool *isinside)
  35. {
  36.     int             no_int;
  37.     dqueue          q_ent;
  38.     double          roots[2];
  39.     vector          sample;
  40.  
  41.     no_int = intersect_quadric(o->data.quadric, r, roots);
  42.     if (flags & CHKINSIDE)
  43.     if (o->data.quadric->isclosed) {
  44.         *isinside = (no_int % 2) ^ o->inverted;
  45.     } else {
  46.         svproduct(sample, 2 * tolerance, r->dir);
  47.         vadd(sample, sample, r->pos);
  48.         *isinside = inside_quadric(o, sample);
  49.     }
  50.     if (!no_int)
  51.     return FALSE;
  52.  
  53.     q_ent.t = roots[0];
  54.     q_ent.obj = o;
  55.     q_ent.entering = !(*isinside);
  56.     add_to_queue(q, q_ent);
  57.  
  58.     if (no_int > 1 && flags & CHKALL) {
  59.     q_ent.t = roots[1];
  60.     q_ent.obj = o;
  61.     q_ent.entering = *isinside;
  62.     add_to_queue(q, q_ent);
  63.     }
  64.     return TRUE;
  65. }
  66.  
  67.  
  68. /* intersect ray with quadratic form. */
  69. PRIVATE int
  70. intersect_quadric(struct quadric_data *qp, struct ray *r, double *roots)
  71. {
  72.     double          a1,
  73.                     a2,
  74.                     a3,
  75.                     b1,
  76.                     b2,
  77.                     b3,
  78.                     xx,
  79.                     yy,
  80.                     zz,
  81.                     xy,
  82.                     yz,
  83.                     xz,
  84.                     x0,
  85.                     y0,
  86.                     z0;
  87.     spoly           quad_poly;
  88.     int             sols;
  89.  
  90.     my_methods.test++;
  91.  
  92.     a1 = r->dir.x;
  93.     a2 = r->dir.y;
  94.     a3 = r->dir.z;
  95.     b1 = r->pos.x;
  96.     b2 = r->pos.y;
  97.     b3 = r->pos.z;
  98.  
  99.     xx = qp->xx;
  100.     yy = qp->yy;
  101.     zz = qp->zz;
  102.     xy = qp->xy;
  103.     yz = qp->yz;
  104.     xz = qp->xz;
  105.     x0 = qp->x0;
  106.     y0 = qp->y0;
  107.     z0 = qp->z0;
  108.  
  109.     /* coefficients of t^2, t^1 and t^0 */
  110.     quad_poly.coef[2] = xx * sqr(a1) + yy * sqr(a2) + zz * sqr(a3)
  111.     + xy * a1 * a2 + yz * a2 * a3 + xz * a3 * a1;
  112.  
  113.     quad_poly.coef[1] = 2 * (xx * a1 * b1 + yy * a2 * b2 + zz * a3 * b3)
  114.     + xy * (a1 * b2 + b1 * a2)
  115.     + yz * (a2 * b3 + a3 * b2)
  116.     + xz * (a1 * b3 + a3 * b1)
  117.     + x0 * a1 + y0 * a2 + z0 * a3;
  118.  
  119.     quad_poly.coef[0] = xx * sqr(b1) + yy * sqr(b2) + zz * sqr(b3)
  120.     + xy * b1 * b2 + yz * b2 * b3 + xz * b1 * b3 +
  121.     x0 * b1 + y0 * b2 + z0 * b3 + qp->a;
  122.  
  123.  
  124.     quad_poly.ord = 2;
  125.     sols = solve_rt_poly(&quad_poly, roots, TRUE);    /* solve polynomial. */
  126.  
  127.     if (sols)
  128.     my_methods.hit++;
  129.  
  130.     return sols;
  131. }
  132.  
  133.  
  134. /* normal to a quadric. It's just the gradient at the desired point */
  135. PRIVATE vector
  136. quadric_normal(struct intersect i, vector x)
  137. {
  138.     vector          n;
  139.     object         *o = i.q_ent.obj;
  140.     struct quadric_data *qp;
  141.  
  142.     qp = o->data.quadric;
  143.  
  144.  
  145.     /* calculate gradient: it's svproduct((d/dx, d/dy, d/dz) , function) */
  146.     n.x = 2 * qp->xx * x.x + qp->xy * x.y + qp->xz * x.z + qp->x0;
  147.     n.y = 2 * qp->yy * x.y + qp->yz * x.z + qp->xy * x.x + qp->y0;
  148.     n.z = 2 * qp->zz * x.z + qp->xz * x.x + qp->yz * x.y + qp->z0;
  149.  
  150.     if (quickveclen(n) < EPSILON) {    /* null vector */
  151.     return n;        /* just return, don't normalize */
  152.     }
  153.     norm(n, n);
  154.     return n;
  155. }
  156.  
  157. PRIVATE bool
  158. inside_quadric(object *o, vector x)
  159. {
  160.     struct quadric_data *q = o->data.quadric;
  161.  
  162.     if (q->xx * sqr(x.x) + q->yy * sqr(x.y) + q->zz * sqr(x.z)
  163.     + q->xy * x.x * x.y + q->yz * x.z * x.y + q->xz * x.x * x.z
  164.     + q->x0 * x.x + q->y0 * x.y + q->z0 * x.z + q->a < 0)
  165.     return !o->inverted;
  166.     else
  167.     return o->inverted;
  168. }
  169.  
  170. PRIVATE void
  171. translate_quadric(object *o, vector t)
  172. {
  173.     struct quadric_data quad,
  174.                    *q;
  175.  
  176.     q = o->data.quadric;
  177.  
  178.     quad = *q;
  179.  
  180.     /*
  181.      * substitute x:=x-t.x, y:=y-t.y, z:=z-t.z in the quadratic equation.
  182.      * Is there a smarter way of doing this?
  183.      */
  184.  
  185.     q->a += -(quad.x0 * t.x + quad.y0 * t.y + quad.z0 * t.z) +
  186.     quad.xy * t.x * t.y + quad.yz * t.y * t.z + quad.xz * t.z * t.x +
  187.     quad.xx * sqr(t.x) + quad.yy * sqr(t.y) + quad.zz * sqr(t.z);
  188.  
  189.     q->x0 += -2 * t.x * quad.xx + -t.y * quad.xy - t.z * quad.xz;
  190.     q->y0 += -2 * t.y * quad.yy + -t.z * quad.yz - t.x * quad.xy;
  191.     q->z0 += -2 * t.z * quad.zz + -t.x * quad.xz - t.y * quad.yz;
  192. }
  193.  
  194.  
  195. /*
  196.  * rotate quadric: a quadric is determined by it's equation Q: (Ax,x) +
  197.  * (b,x) + c = 0. (A matrix, b vector, c constant) so if x \in RQ (R is a
  198.  * rotation, then R^{-1}x \in Q, so (AR^-1x, x) + (b, R^-1x) + c = 0, or
  199.  * (RAR^Tx,x) + (Rx,b)  +c = 0.
  200.  */
  201.  
  202. PRIVATE void
  203. rotate_quadric(object *o, matrix rotmat)
  204. {
  205.     matrix          quad,
  206.                     prod,
  207.                     transm;
  208.     vector          lin;
  209.     struct quadric_data *q;
  210.  
  211.     q = o->data.quadric;
  212.  
  213.     if (o == NULL || q == NULL)
  214.     assert(FALSE);
  215.  
  216.  
  217.     null_matrix(quad);
  218.  
  219.     quad[0][0] = q->xx;
  220.     quad[1][1] = q->yy;
  221.     quad[2][2] = q->zz;
  222.     quad[0][1] = quad[1][0] = q->xy / 2.0;
  223.     quad[1][2] = quad[2][1] = q->yz / 2.0;
  224.     quad[0][2] = quad[2][0] = q->xz / 2.0;
  225.  
  226.     transpose_matrix(transm, rotmat);
  227.  
  228.     /* calc AR^T, store result in prod, quad is no longer needed */
  229.     mmproduct(prod, quad, transm);
  230.  
  231.     /* calc RAR^T, store result in quad */
  232.     mmproduct(quad, rotmat, prod);
  233.  
  234.     q->xx = quad[0][0];
  235.     q->yy = quad[1][1];
  236.     q->zz = quad[2][2];
  237.  
  238.     q->xy = quad[0][1] + quad[1][0];
  239.     q->yz = quad[1][2] + quad[2][1];
  240.     q->xz = quad[0][2] + quad[2][0];
  241.  
  242.     lin.x = q->x0;
  243.     lin.y = q->y0;
  244.     lin.z = q->z0;
  245.     rotate_vector(&lin, rotmat);
  246.     q->x0 = lin.x;
  247.     q->y0 = lin.y;
  248.     q->z0 = lin.z;
  249. }
  250.  
  251. PRIVATE void
  252. scale_quadric(object *o, vector s)
  253. {
  254.     double          a,
  255.                     b,
  256.                     c;
  257.     struct quadric_data *q;
  258.  
  259.     q = o->data.quadric;
  260.     if (q == NULL)
  261.     assert(FALSE);
  262.  
  263.     a = 1 / s.x;
  264.     b = 1 / s.y;
  265.     c = 1 / s.z;
  266.  
  267.     q->xx *= sqr(a);
  268.     q->yy *= sqr(b);
  269.     q->zz *= sqr(c);
  270.  
  271.     q->xy *= a * b;
  272.     q->yz *= b * c;
  273.     q->xz *= c * a;
  274.  
  275.     q->x0 *= a;
  276.     q->y0 *= b;
  277.     q->z0 *= c;
  278. }
  279.  
  280.  
  281. PRIVATE void
  282. copy_quadric(object *dst, object *src)
  283. {
  284.     if (src == NULL || dst == NULL)
  285.     assert(FALSE);
  286.     if (dst->type != QUADRIC)
  287.     dst->data.quadric = get_new_quadric();
  288.     *dst->data.quadric = *src->data.quadric;
  289.     dst->type = src->type;
  290. }
  291.  
  292.  
  293. PRIVATE void
  294. print_quadric(object *o)
  295. {
  296. #ifdef DEBUG
  297.     struct quadric_data *q;
  298.  
  299.     q = o->data.quadric;
  300.  
  301.     printf("%5lfx^2 + %5lfy^2 + %5lfz^2 + %5lfxy\n"
  302.        "\t+ %5lfyz +%5lfxz %5lfx\n"
  303.        "\t+ %5lfy + %5lfz + %5lf = 0\n",
  304.        q->xx, q->yy, q->zz,
  305.        q->xy, q->yz, q->xz,
  306.        q->x0, q->y0, q->z0,
  307.        q->a);
  308.  
  309. #endif
  310. }
  311.  
  312. PRIVATE void
  313. free_quadric(object *o)
  314. {
  315.     free((void *) o->data.quadric);
  316.     o->type = NOSHAPE;
  317. }
  318.  
  319. PUBLIC void
  320. init_quadric(struct quadric_data *p)
  321. {
  322.     p->xx = p->yy = p->zz = p->xy = p->yz = p->xz = p->x0 = p->y0 = p->z0 = p->a = 0.0;
  323.     p->isclosed = FALSE;
  324. }
  325.  
  326. PUBLIC struct quadric_data *
  327. get_new_quadric(void)
  328. {
  329.     struct quadric_data *p;
  330.     p = ALLOC(struct quadric_data);
  331.  
  332.     CHECK_MEM(p, my_methods.name);
  333.  
  334.     init_quadric(p);
  335.  
  336.     return p;
  337. }
  338.  
  339. PRIVATE void
  340. precompute_quadric(object *o)
  341. {
  342.     my_methods.howmuch++;
  343.     global_stats.prims++;
  344. }
  345.  
  346. PRIVATE struct methods my_methods =
  347. {
  348.     all_quadric_intersections,
  349.     quadric_normal,
  350.     copy_quadric,
  351.     inside_quadric,
  352.     rotate_quadric,
  353.     translate_quadric,
  354.     scale_quadric,
  355.     free_quadric,
  356.     print_quadric,
  357.  
  358.     precompute_quadric,
  359.     "quadric",
  360. };
  361.  
  362. PUBLIC object  *
  363. get_new_quadric_object(void)
  364. {
  365.     object         *o;
  366.  
  367.     o = get_new_object();
  368.     o->data.quadric = get_new_quadric();
  369.     o->methods = &my_methods;
  370.     o->type = QUADRIC;
  371.  
  372.     return o;
  373. }
  374.