home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2213 / quadric.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-12-28  |  3.1 KB  |  179 lines
  1.  
  2. /*
  3.  * quadric.c - This module conatins all of the code that relates to
  4.  * quadratics.
  5.  * 
  6.  * Copyright (C) 1990, Kory Hamzeh
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include <malloc.h>
  11. #include <math.h>
  12.  
  13. #include "rt.h"
  14. #include "externs.h"
  15.  
  16.  
  17. int             Quadric_intersect(), Quadric_normal();
  18.  
  19. /*
  20.  * Build_quadric()
  21.  * 
  22.  * Given some info on a quadric, build the entire object structure.
  23.  */
  24.  
  25. Build_quadric(q)
  26. QUADRIC        *q;
  27. {
  28.     OBJECT         *o;
  29.  
  30.     if (nobjects == MAX_PRIMS)
  31.     {
  32.         fprintf(stderr, "%s: too many objects specified\n", my_name);
  33.         exit(1);
  34.     }
  35.  
  36.     if ((o = (OBJECT *) malloc(sizeof(OBJECT))) == NULL)
  37.     {
  38.         fprintf(stderr, "%s: malloc failed\n", my_name);
  39.         exit(1);
  40.     }
  41.  
  42.  
  43.     o->type = T_QUADRIC;
  44.     o->obj = q;
  45.     o->surf = cur_surface;
  46.     o->inter = Quadric_intersect;
  47.     o->normal = Quadric_normal;
  48.  
  49.     objects[nobjects++] = o;
  50.  
  51.     /*
  52.      * Calculate some constants that we will need in the intersect
  53.      * routine.
  54.      */
  55.  
  56.     q->a2 = q->a * 2.0;
  57.     q->b2 = q->b * 2.0;
  58.     q->c2 = q->c * 2.0;
  59.     q->d2 = q->d * 2.0;
  60.     q->f2 = q->f * 2.0;
  61.     q->g2 = q->g * 2.0;
  62.     q->i2 = q->i * 2.0;
  63.  
  64.     /*
  65.      * Setup the min and the max values for the bouding box. For now,
  66.      * just use what the user specifies.
  67.      */
  68.  
  69.     o->b_min = q->min;
  70.     o->b_max = q->max;
  71.  
  72. }
  73.  
  74.  
  75.  
  76. /*
  77.  * Quadric_intersect()
  78.  * 
  79.  * Check given quadratic for intersection with given ray. Return TRUE if an
  80.  * intersection takes place.
  81.  */
  82.  
  83. Quadric_intersect(obj, ray, inter)
  84. OBJECT         *obj;
  85. RAY            *ray;
  86. INTERSECT      *inter;
  87. {
  88.     QUADRIC        *q;
  89.     VECTOR          rd, rp;
  90.     double          aq, nbq, cq;
  91.     double          t, disc;
  92.     double          ka, kb;
  93.  
  94.     q = (QUADRIC *) obj->obj;
  95.  
  96.     rd = ray->dir;
  97.     rp = ray->pos;
  98.  
  99.     /*
  100.      * Compute Aq, Bq, Cq.
  101.      */
  102.  
  103.     aq = rd.x * (q->a * rd.x + q->b2 * rd.y + q->c2 * rd.z) +
  104.         rd.y * (q->e * rd.y + q->f2 * rd.z) +
  105.         q->h * rd.z * rd.z;
  106.  
  107.     nbq = rd.x * (q->a * rp.x + q->b * rp.y + q->c * rp.z + q->d) +
  108.         rd.y * (q->b * rp.x + q->e * rp.y + q->f * rp.z + q->g) +
  109.         rd.z * (q->c * rp.x + q->f * rp.y + q->h * rp.z + q->i);
  110.  
  111.     cq = rp.x * (q->a * rp.x + q->b2 * rp.y + q->c2 * rp.z + q->d2) +
  112.         rp.y * (q->e * rp.y + q->f2 * rp.z + q->g2) +
  113.         rp.z * (q->h * rp.z + q->i2) + q->j;
  114.  
  115.     if (fabs(aq) < MIN_T)
  116.     {
  117.         t = -cq / (2 * nbq);
  118.         if (t < MIN_T)
  119.             return (0);    /* no hit */
  120.  
  121.         inter->obj = obj;
  122.         inter->t = t;
  123.         inter->inside = 0;
  124.         return (1);
  125.     }
  126.  
  127.     /*
  128.      * Compute the discriminator.
  129.      */
  130.  
  131.     ka = -nbq / aq;
  132.     kb = cq / aq;
  133.  
  134.     disc = (ka * ka) - (kb);
  135.     if (disc < MIN_T)
  136.         return (0);
  137.  
  138.     t = ka - sqrt(disc);
  139.  
  140.     if (t < MIN_T)
  141.     {
  142.         t = ka + sqrt(disc);
  143.         if (t < MIN_T)
  144.             return (0);
  145.         inter->inside = 1;
  146.     }
  147.     else
  148.         inter->inside = 0;
  149.  
  150.     inter->t = t;
  151.     inter->obj = obj;
  152.  
  153.     return (1);
  154. }
  155.  
  156.  
  157. /*
  158.  * Quadric_normal()
  159.  * 
  160.  * Return the normal to a quadric at a given point along the surface.
  161.  */
  162.  
  163. Quadric_normal(q, ray, ip, normal)
  164. QUADRIC        *q;
  165. RAY            *ray;
  166. VECTOR         *ip;
  167. VECTOR         *normal;
  168. {
  169.  
  170.     normal->x = q->a * ip->x + q->b * ip->y + q->c * ip->z + q->d;
  171.     normal->y = q->b * ip->x + q->e * ip->y + q->f * ip->z + q->g;
  172.     normal->z = q->c * ip->x + q->f * ip->y + q->h * ip->z + q->i;
  173.  
  174.     VecNormalize(normal);
  175.  
  176.     if (VecDot(ray->dir, *normal) >= 0)
  177.         VecNegate(*normal);
  178. }
  179.