home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / algebrai.c next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  12.6 KB  |  598 lines
  1. /*
  2.  * algebraic.c -- handle higher order surfaces.
  3.  * 
  4.  * Copyright (c) by Han-Wen Nienhuys, 1993
  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.  
  25. extern struct methods my_methods;
  26.  
  27.  
  28. /* the polystack with polynomials */
  29. PRIVATE spoly  *polystack = NULL;
  30. PRIVATE double *dpolystack = NULL;
  31.  
  32. PRIVATE int     polystackp,    /* and it's pointer */
  33.                 maxpolystack;
  34.  
  35. PRIVATE vector  raypos,        /* current ray.pos */
  36.                 raydir;        /* and ray.dir */
  37.  
  38. PRIVATE void    do_poly_op(struct poly_instruction *instr);
  39. PRIVATE void    do_poly_diff_op(struct poly_instruction *instr);
  40.  
  41. PRIVATE void    print_polystack(void);
  42. PRIVATE void    print_instr(struct poly_instruction *instr);
  43. PRIVATE void    print_diffstack(void);
  44.  
  45. /* execute an instruction */
  46. PRIVATE void
  47. exec(struct poly_instruction *instr)
  48. {
  49.     switch (instr->type) {
  50.     case POLY_CONS:        /* push a constant on the stack */
  51.     polystack[polystackp].coef[0] =
  52.         instr->data.factor;
  53.     polystack[polystackp].ord = 0;
  54.     polystackp++;
  55.     break;
  56.     case POLY_OP:        /* evaluate operator. */
  57.     do_poly_op(instr);
  58.     break;
  59.     case POLY_X:        /* push X on stack */
  60.     polystack[polystackp].coef[0] =
  61.         raypos.x;
  62.     if ((polystack[polystackp].coef[1] = raydir.x))
  63.         polystack[polystackp].ord = 1;
  64.     else
  65.         polystack[polystackp].ord = 0;
  66.  
  67.     polystackp++;
  68.     break;
  69.     case POLY_Y:        /* push  Y on stack */
  70.     polystack[polystackp].coef[0] = raypos.y;
  71.     if ((polystack[polystackp].coef[1] = raydir.y))
  72.         polystack[polystackp].ord = 1;
  73.     else
  74.         polystack[polystackp].ord = 0;
  75.     polystackp++;
  76.     break;
  77.     case POLY_Z:        /* and Z */
  78.     polystack[polystackp].coef[0] = raypos.z;
  79.     if ((polystack[polystackp].coef[1] = raydir.z))
  80.         polystack[polystackp].ord = 1;
  81.     else
  82.         polystack[polystackp].ord = 0;
  83.     polystackp++;
  84.     break;
  85.     default:
  86.     assert(FALSE);
  87.  
  88.     }
  89. }
  90.  
  91.  
  92. /* do an operation on the stack */
  93. PRIVATE void
  94. do_poly_op(struct poly_instruction *instr)
  95. {
  96.     spoly           temp;
  97.     int             exponent;
  98.  
  99.     switch (instr->data.op) {
  100.     case 'n':            /* negate */
  101.     poly_neg(polystack + polystackp - 1,
  102.          polystack + polystackp - 1);
  103.     break;
  104.     case '-':            /* substract */
  105.     poly_sub(polystack + polystackp - 2,
  106.          polystack + polystackp - 2, polystack + polystackp - 1);
  107.     polystackp--;
  108.  
  109.     break;
  110.     case '+':            /* add */
  111.     poly_add(polystack + polystackp - 2, polystack + polystackp - 2,
  112.          polystack + polystackp - 1);
  113.     polystackp--;
  114.     break;
  115.     case '*':            /* multiply */
  116.     poly_multiply(&temp, polystack + polystackp - 2,
  117.               polystack + polystackp - 1);
  118.     poly_copy(polystack + polystackp - 2, &temp);
  119.     polystackp--;
  120.     break;
  121.     case '^':            /* power */
  122.     exponent = (int) polystack[polystackp - 1].coef[0];
  123.     poly_power(&temp, exponent, polystack + polystackp - 2);
  124.     poly_copy(polystack + polystackp - 2, &temp);
  125.     polystackp--;
  126.     break;
  127.     }
  128. }
  129.  
  130. /* calculate a polynomial from its code */
  131. PRIVATE void
  132. calc_poly(spoly *dst, struct poly_instruction *instr, int
  133.       len)
  134. {
  135.     int             i;
  136.  
  137.     /* init stack */
  138.     polystackp = 0;
  139.  
  140.     /* execute poly instructions to form polynomial */
  141.     for (i = 0; i < len; i++) {
  142.     exec(instr + i);
  143.     }
  144.  
  145.     assert(polystackp == 1);
  146.  
  147.     poly_copy(dst, polystack);
  148.  
  149. }
  150.  
  151. /* zapp */
  152. PRIVATE void
  153. free_algebraic(object *o)
  154. {
  155.     struct algebraic_data *a = o->data.algebraic;
  156.  
  157.     free(a->code);
  158.     free((void *) a);
  159.     o->type = NOSHAPE;
  160. }
  161.  
  162. /* evaluate the polynomial in X */
  163. PRIVATE bool
  164. inside_algebraic(object *o, vector x)
  165. {
  166.     spoly           tmp;
  167.  
  168.     raypos = x;
  169.     setvector(raydir, 0, 0, 0);
  170.  
  171.     calc_poly(&tmp, o->data.algebraic->code, o->data.algebraic->codesize);
  172.     poly_clean(&tmp);
  173.  
  174.     assert(tmp.ord == 0);
  175.  
  176.     return (tmp.coef[0] > 0) ^ o->inverted;
  177. }
  178.  
  179. /* debugging stuff */
  180. PRIVATE void
  181. print_algebraic(object *o)
  182. {
  183.  
  184. #ifdef DEBUG
  185.     int             i;
  186.     struct algebraic_data *a;
  187.  
  188.     a = o->data.algebraic;
  189.  
  190.  
  191.     printf("\t%s%s", (a->isclosed) ? "closed " : "", (a->usesturm) ? "sturm solver" : "");
  192.     printf("code: \n[");
  193.     for (i = 0; i < a->codesize; i++)
  194.     print_instr(a->code + i);
  195.     printf("]\n");
  196. #endif
  197. }
  198.  
  199. /* just the intersections. */
  200. PRIVATE int
  201. intersect_algebraic(struct algebraic_data *t, struct ray r,
  202.             double *rhits, bool chkall)
  203. {
  204.     spoly           tosolve;
  205.     int             n;
  206.  
  207.  
  208.     my_methods.test++;
  209.  
  210.     /* back to object coords */
  211.     raypos = r.pos;
  212.     raydir = r.dir;
  213.  
  214.     /* construct polynomial */
  215.     calc_poly(&tosolve, t->code, t->codesize);
  216.  
  217.     /* and solve it */
  218.     if (t->usesturm)
  219.     n = sturm_solve_rt_poly(&tosolve, rhits, chkall);
  220.     else
  221.     n = solve_rt_poly(&tosolve, rhits, chkall);
  222.  
  223.     if (n)
  224.     my_methods.hit++;
  225.  
  226.     return n;
  227. }
  228.  
  229. /*
  230.  * find each and every intersection with a polysurface, and stuff them
  231.  * into a queue
  232.  */
  233. PRIVATE bool
  234. all_algebraic_intersections(dqueue * q, object *o, struct ray *r,
  235.                 int flags, bool *isinside)
  236. {
  237.     int             n,
  238.                     j;
  239.     double          inter[MAX_ORDER];
  240.     dqueue          q_ent;
  241.     struct ray      localray;
  242.     bool            loc_ins;
  243.  
  244.  
  245.     localray = *r;
  246.     transform_ray(&localray, o);
  247.  
  248.     n = intersect_algebraic(o->data.algebraic, localray, inter, flags & CHKALL);
  249.  
  250.  
  251.     if (flags & CHKINSIDE) {
  252.     if (o->data.algebraic->isclosed)    /* closed surface */
  253.         loc_ins = *isinside = (n % 2) ^ o->inverted;
  254.     else {
  255.         /* not closed, then use brute force */
  256.         vector          sample;
  257.  
  258.         svproduct(sample, 2 * tolerance, r->dir);
  259.         vadd(sample, sample, r->pos);
  260.         loc_ins = *isinside = inside_algebraic(o, sample);
  261.     }
  262.     }
  263.     if (!n)
  264.     return FALSE;
  265.  
  266.  
  267.     if (flags & CHKALL && n > 1)/* check everything? */
  268.     for (j = 0; j < n; j++) {
  269.         q_ent.t = inter[j];
  270.         q_ent.obj = o;
  271.         q_ent.entering = (loc_ins = !loc_ins);
  272.         add_to_queue(q, q_ent);
  273.     } else {
  274.  
  275.     /* record only one hit. */
  276.  
  277.     q_ent.t = inter[0];
  278.     q_ent.obj = o;
  279.     q_ent.entering = (loc_ins = !loc_ins);
  280.     add_to_queue(q, q_ent);
  281.     }
  282.  
  283.     return TRUE;
  284. }
  285.  
  286.  
  287.  
  288. PRIVATE vector  normal_loc;
  289. PRIVATE int     diff_to;
  290.  
  291.  
  292. /* handle instructions, for derivative */
  293. PRIVATE void
  294. exec_diff(struct poly_instruction *instr)
  295. {
  296.     switch (instr->type) {
  297.     case POLY_OP:
  298.     do_poly_diff_op(instr);
  299.     break;
  300.     case POLY_X:
  301.     dpolystack[polystackp] = normal_loc.x;
  302.     dpolystack[polystackp + 1] = (diff_to == POLY_X);
  303.     polystackp += 2;
  304.     break;
  305.     case POLY_Y:
  306.     dpolystack[polystackp] = normal_loc.y;
  307.     dpolystack[polystackp + 1] = (diff_to == POLY_Y);
  308.     polystackp += 2;
  309.     break;
  310.     case POLY_Z:
  311.     dpolystack[polystackp] = normal_loc.z;
  312.     dpolystack[polystackp + 1] = (diff_to == POLY_Z);
  313.     polystackp += 2;
  314.     break;
  315.     case POLY_CONS:
  316.     dpolystack[polystackp] = instr->data.factor;
  317.     dpolystack[polystackp + 1] = 0;
  318.     polystackp += 2;
  319.     break;
  320.     default:
  321.     assert(FALSE);
  322.     }
  323. }
  324.  
  325.  
  326. /*
  327.  * calculate the resulting polynomial for p1 OP p2, and it's derivative.
  328.  */
  329. PRIVATE void
  330. do_poly_diff_op(struct poly_instruction *instr)
  331. {
  332.     double          p,
  333.                     Dp;
  334.  
  335.     switch (instr->data.op) {
  336.     case '+':
  337.     polystackp -= 2;
  338.     dpolystack[polystackp - 2] = dpolystack[polystackp - 2]
  339.         + dpolystack[polystackp];
  340.  
  341.     /* (p1+p2)' = p1' + p2' */
  342.     dpolystack[polystackp - 1] = dpolystack[polystackp - 1]
  343.         + dpolystack[polystackp + 1];
  344.     break;
  345.  
  346.     case '-':
  347.     polystackp -= 2;
  348.     dpolystack[polystackp - 2] = dpolystack[polystackp - 2]
  349.         - dpolystack[polystackp];
  350.  
  351.     /* (p1-p2)' = p1' - p2' */
  352.     dpolystack[polystackp - 1] = dpolystack[polystackp - 1]
  353.         - dpolystack[polystackp + 1];
  354.     break;
  355.  
  356.     case '*':
  357.  
  358.     polystackp -= 2;
  359.     p = dpolystack[polystackp - 2] * dpolystack[polystackp];
  360.  
  361.     /* (p1*p2)' = p1'*p2 + p2'*p1 */
  362.     Dp = dpolystack[polystackp - 2] * dpolystack[polystackp + 1] +
  363.         dpolystack[polystackp - 1] * dpolystack[polystackp];
  364.     dpolystack[polystackp - 2] = p;
  365.     dpolystack[polystackp - 1] = Dp;
  366.     break;
  367.     case '^':
  368.  
  369.     polystackp -= 2;
  370.     p = pow(dpolystack[polystackp - 2], dpolystack[polystackp]);
  371.  
  372.     /* (p^a)' = a * p^(a-1) p' */
  373.     Dp = pow(dpolystack[polystackp - 2], dpolystack[polystackp] - 1.0)
  374.         * dpolystack[polystackp - 1] * dpolystack[polystackp];
  375.     dpolystack[polystackp - 2] = p;
  376.     dpolystack[polystackp - 1] = Dp;
  377.     break;
  378.     case 'n':
  379.  
  380.     dpolystack[polystackp - 2] = -dpolystack[polystackp - 2];
  381.  
  382.     /* (-p)' = - p' */
  383.     dpolystack[polystackp - 1] = -dpolystack[polystackp - 1];
  384.     break;
  385.     default:
  386.     assert(FALSE);
  387.     }
  388. }
  389.  
  390.  
  391. /*
  392.  * calculate normal. calculate the gradient for each component
  393.  */
  394. PRIVATE double
  395. calc_normal(struct poly_instruction *p, int size)
  396. {
  397.     int             i;
  398.     double          retval;
  399.  
  400.     polystackp = 0;
  401.  
  402.     /* calcalute the gradient. */
  403.     for (i = 0; i < size; i++) {
  404.     exec_diff(p + i);
  405.     }
  406.     assert(polystackp <= 2);
  407.  
  408.  
  409.     retval = dpolystack[1];
  410.  
  411.     return retval;
  412. }
  413.  
  414. /*
  415.  * returns the normal to the algebraic
  416.  */
  417. PRIVATE vector
  418. algebraic_normal(struct intersect i, vector loc)
  419. {
  420.     vector          n;
  421.     struct algebraic_data *a;
  422.  
  423.     a = i.q_ent.obj->data.algebraic;
  424.  
  425.     /* back to object coords */
  426.     normal_loc = mvproduct(*i.q_ent.obj->inv_trans, loc);
  427.  
  428.     diff_to = POLY_X;
  429.     n.x = calc_normal(a->code, a->codesize);
  430.     diff_to = POLY_Y;
  431.     n.y = calc_normal(a->code, a->codesize);
  432.     diff_to = POLY_Z;
  433.     n.z = calc_normal(a->code, a->codesize);
  434.  
  435.     n = transform_normal(*i.q_ent.obj->inv_trans, n);
  436.     norm(n, n);
  437.  
  438.     return n;
  439. }
  440.  
  441. /* initialize a algebraic_data struct */
  442. PRIVATE void
  443. init_algebraic(struct algebraic_data *t)
  444. {
  445.     t->code = NULL;
  446.     t->codesize = 0;
  447.  
  448.     t->usesturm = FALSE;
  449.     t->isclosed = FALSE;
  450. }
  451.  
  452. /* alloc and init algebraic */
  453. PRIVATE struct algebraic_data *
  454. get_new_algebraic(void)
  455. {
  456.     struct algebraic_data *a;
  457.  
  458.     a = ALLOC(struct algebraic_data);
  459.  
  460.     CHECK_MEM(a, my_methods.name);
  461.     init_algebraic(a);
  462.  
  463.     return a;
  464. }
  465.  
  466. /* copy the shape data */
  467. PRIVATE void
  468. copy_algebraic(object *dst, object *src)
  469. {
  470.     struct algebraic_data *s,
  471.                    *d;
  472.  
  473.     assert(dst != NULL && src != NULL);
  474.  
  475.     /* alloc a new one */
  476.     if (dst->type != ALGEBRAIC) {
  477.     dst->data.algebraic = get_new_algebraic();
  478.     dst->type = ALGEBRAIC;
  479.     }
  480.     d = dst->data.algebraic;
  481.     s = src->data.algebraic;
  482.     *d = *s;
  483.  
  484.  
  485.     /* alloc new code */
  486.     d->code = malloc(sizeof(struct poly_instruction) * s->codesize);
  487.  
  488.     CHECK_MEM(d->code, "polycode");
  489.     memcpy(d->code, s->code, s->codesize * sizeof(struct poly_instruction));
  490. }
  491.  
  492. PRIVATE void
  493. precompute_algebraic(object *o)
  494. {
  495.     struct algebraic_data *a;
  496.  
  497.     a = o->data.algebraic;
  498.  
  499.     if (a->codesize > maxpolystack) {
  500.     maxpolystack = a->codesize + 1;
  501.     polystack = realloc(polystack, maxpolystack * sizeof(spoly));
  502.  
  503.     CHECK_MEM(polystack, "polystack");
  504.     dpolystack = realloc(dpolystack, 2 * maxpolystack * sizeof(double));
  505.  
  506.     CHECK_MEM(dpolystack, "dpolystack");
  507.     }
  508.     my_methods.howmuch++;
  509.     global_stats.prims++;
  510. }
  511.  
  512.  
  513. PRIVATE struct methods my_methods =
  514. {
  515.     all_algebraic_intersections,
  516.     algebraic_normal,
  517.     copy_algebraic,
  518.     inside_algebraic,
  519.     generic_rotate_object,
  520.     generic_translate_object,
  521.     generic_scale_object,
  522.     free_algebraic,
  523.     print_algebraic,
  524.     precompute_algebraic,
  525.     "algebraic",
  526. };
  527.  
  528. /* alloc/init */
  529. PUBLIC object  *
  530. get_new_algebraic_object(void)
  531. {
  532.     object         *o;
  533.  
  534.     o = get_new_object();
  535.  
  536.     o->data.algebraic = get_new_algebraic();
  537.     o->methods = &my_methods;
  538.     o->type = ALGEBRAIC;
  539.  
  540.     return o;
  541. }
  542.  
  543. #ifdef DEBUG
  544.  
  545. /* dump stack */
  546. PRIVATE void
  547. print_polystack(void)
  548. {
  549.     int             i;
  550.  
  551.     printf("\n----\n");
  552.     for (i = 0; i < polystackp; i++) {
  553.     printf("%d : ", i);
  554.     print_poly(polystack + i);
  555.     }
  556. }
  557.  
  558. /* print out the stack for calculating the derivative */
  559. PRIVATE void
  560. print_diffstack(void)
  561. {
  562.     int             i;
  563.  
  564.     printf("\n----\n");
  565.     for (i = 0; i < polystackp; i += 2) {
  566.     printf("%d : %f, %f\n", i / 2, dpolystack[i], dpolystack[i + 1]);
  567.     }
  568. }
  569.  
  570. /* print instruction */
  571. PRIVATE void
  572. print_instr(struct poly_instruction *instr)
  573. {
  574.     switch (instr->type) {
  575.     case POLY_CONS:        /* push a constant on the stack */
  576.     printf("c %f ", instr->data.factor);
  577.     break;
  578.     case POLY_OP:        /* evaluate operator. */
  579.     printf("op %c ", instr->data.op);
  580.     break;
  581.     case POLY_X:        /* push X on stack */
  582.     printf("X ");
  583.     break;
  584.     case POLY_Y:        /* push  Y on stack */
  585.     printf("Y ");
  586.     break;
  587.     case POLY_Z:        /* and Z */
  588.     printf("Z ");
  589.     break;
  590.     default:
  591.     assert(FALSE);
  592.     }
  593. }
  594.  
  595. #endif
  596.  
  597. /* eof */
  598.