home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2214 / cone.c next >
Encoding:
C/C++ Source or Header  |  1990-12-28  |  5.5 KB  |  331 lines
  1.  
  2. /*
  3.  * cone.c
  4.  * 
  5.  * This module conatins all of the code which relates to cones and cylinders.
  6.  * 
  7.  * I must admit, MTV's code helped to "inspire" this module.
  8.  *
  9.  * Copyright (C) 1990, Kory Hamzeh.
  10.  */
  11.  
  12.  
  13. #include <stdio.h>
  14. #include <math.h>
  15. #include "rt.h"
  16. #include "externs.h"
  17.  
  18.  
  19. int             Cone_intersect(), Cone_normal();
  20.  
  21.  
  22.  
  23. Build_cone(cd)
  24. CONE           *cd;
  25. {
  26.     OBJECT         *obj;
  27.     double          dmin, dmax, d, ftmp;
  28.     VECTOR          tmp;
  29.     int             i;
  30.  
  31.     if (nobjects == MAX_PRIMS)
  32.     {
  33.         fprintf(stderr, "%s: too many objects specified\n", my_name);
  34.         exit(1);
  35.     }
  36.  
  37.     if ((obj = (OBJECT *) malloc(sizeof(OBJECT))) == NULL)
  38.     {
  39.         fprintf(stderr, "%s: malloc failed\n", my_name);
  40.         exit(1);
  41.     }
  42.  
  43.     objects[nobjects++] = obj;
  44.  
  45.     obj->type = T_CONE;
  46.     obj->inter = Cone_intersect;
  47.     obj->normal = Cone_normal;
  48.     obj->surf = cur_surface;
  49.  
  50.     VecSub(cd->apex, cd->base, cd->w);
  51.     cd->height = VecNormalize(&cd->w);
  52.     cd->slope = (cd->apex_radius - cd->base_radius) / (cd->height);
  53.     cd->base_d = -VecDot(cd->base, cd->w);
  54.  
  55.     tmp.x = 0;
  56.     tmp.y = 0;
  57.     tmp.z = 1;
  58.  
  59.     if (1.0 - fabs(VecDot(tmp, cd->w)) < MIN_T)
  60.     {
  61.         tmp.y = 1;
  62.         tmp.x = 0;
  63.     }
  64.  
  65.     /*
  66.      * find two axes which are at right angles to cone_w
  67.      */
  68.  
  69.     VecCross(cd->w, tmp, cd->u);
  70.     VecCross(cd->u, cd->w, cd->v);
  71.  
  72.     VecNormalize(&cd->u);
  73.     VecNormalize(&cd->v);
  74.  
  75.     cd->min_d = VecDot(cd->w, cd->base);
  76.     cd->max_d = VecDot(cd->w, cd->apex);
  77.  
  78.     if (cd->max_d < cd->min_d)
  79.     {
  80.         ftmp = cd->max_d;
  81.         cd->max_d = cd->min_d;
  82.         cd->min_d = ftmp;
  83.     }
  84.  
  85.     obj->obj = cd;
  86.  
  87.     /*
  88.      * Create the bounding box for this puppy.
  89.      */
  90.  
  91.     dmin = HUGE;
  92.     dmax = -HUGE;
  93.  
  94.     /* first the X plane */
  95.     d = cd->base.x - cd->base_radius;
  96.  
  97.     if (d < dmin)
  98.         dmin = d;
  99.     if (d > dmax)
  100.         dmax = d;
  101.  
  102.     d = cd->base.x + cd->base_radius;
  103.  
  104.     if (d < dmin)
  105.         dmin = d;
  106.     if (d > dmax)
  107.         dmax = d;
  108.  
  109.     d = cd->apex.x - cd->apex_radius;
  110.     if (d < dmin)
  111.         dmin = d;
  112.     if (d > dmax)
  113.         dmax = d;
  114.  
  115.     d = cd->apex.x + cd->apex_radius;
  116.     if (d < dmin)
  117.         dmin = d;
  118.     if (d > dmax)
  119.         dmax = d;
  120.  
  121.     obj->b_min.x = dmin;
  122.     obj->b_max.x = dmax;
  123.  
  124.     /* now the Y plane */
  125.     dmin = HUGE;
  126.     dmax = -HUGE;
  127.  
  128.     d = cd->base.y - cd->base_radius;
  129.  
  130.     if (d < dmin)
  131.         dmin = d;
  132.     if (d > dmax)
  133.         dmax = d;
  134.  
  135.     d = cd->base.y + cd->base_radius;
  136.  
  137.     if (d < dmin)
  138.         dmin = d;
  139.     if (d > dmax)
  140.         dmax = d;
  141.  
  142.     d = cd->apex.y - cd->apex_radius;
  143.     if (d < dmin)
  144.         dmin = d;
  145.     if (d > dmax)
  146.         dmax = d;
  147.  
  148.     d = cd->apex.y + cd->apex_radius;
  149.     if (d < dmin)
  150.         dmin = d;
  151.     if (d > dmax)
  152.         dmax = d;
  153.  
  154.     obj->b_min.y = dmin;
  155.     obj->b_max.y = dmax;
  156.  
  157.     /* and finally the Z plane */
  158.     dmin = HUGE;
  159.     dmax = -HUGE;
  160.  
  161.     d = cd->base.z - cd->base_radius;
  162.  
  163.     if (d < dmin)
  164.         dmin = d;
  165.     if (d > dmax)
  166.         dmax = d;
  167.  
  168.     d = cd->base.z + cd->base_radius;
  169.  
  170.     if (d < dmin)
  171.         dmin = d;
  172.     if (d > dmax)
  173.         dmax = d;
  174.  
  175.     d = cd->apex.z - cd->apex_radius;
  176.     if (d < dmin)
  177.         dmin = d;
  178.     if (d > dmax)
  179.         dmax = d;
  180.  
  181.     d = cd->apex.z + cd->apex_radius;
  182.     if (d < dmin)
  183.         dmin = d;
  184.     if (d > dmax)
  185.         dmax = d;
  186.  
  187.     obj->b_min.z = dmin;
  188.     obj->b_max.z = dmax;
  189. }
  190.  
  191.  
  192. Cone_intersect(obj, ray, inter)
  193. OBJECT         *obj;
  194. RAY            *ray;
  195. INTERSECT      *inter;
  196. {
  197.     RAY             tray;
  198.     CONE           *cd;
  199.     VECTOR          v, p;
  200.     double          a, b, c, d, disc;
  201.     double          t1, t2;
  202.     int             nroots;
  203.  
  204.     cd = (CONE *) (obj->obj);
  205.  
  206.     /*
  207.      * First, we get the coordinates of the ray origin in the objects
  208.      * space....
  209.      */
  210.  
  211.     VecSub(ray->pos, cd->base, v);
  212.  
  213.     tray.pos.x = VecDot(v, cd->u);
  214.     tray.pos.y = VecDot(v, cd->v);
  215.     tray.pos.z = VecDot(v, cd->w);
  216.  
  217.     tray.dir.x = VecDot(ray->dir, cd->u);
  218.     tray.dir.y = VecDot(ray->dir, cd->v);
  219.     tray.dir.z = VecDot(ray->dir, cd->w);
  220.  
  221.  
  222.     a = tray.dir.x * tray.dir.x
  223.         + tray.dir.y * tray.dir.y
  224.         - cd->slope * cd->slope * tray.dir.z * tray.dir.z;
  225.  
  226.     b = 2.0 * (tray.pos.x * tray.dir.x + tray.pos.y * tray.dir.y -
  227.            cd->slope * cd->slope * tray.pos.z * tray.dir.z
  228.            - cd->base_radius * cd->slope * tray.dir.z);
  229.  
  230.     c = cd->slope * tray.pos.z + cd->base_radius;
  231.     c = tray.pos.x * tray.pos.x + tray.pos.y * tray.pos.y - (c * c);
  232.  
  233.     disc = b * b - 4.0 * a * c;
  234.  
  235.     if (disc < 0.0)
  236.         return (0);
  237.  
  238.     disc = sqrt(disc);
  239.     t1 = (-b - disc) / (2.0 * a);
  240.     t2 = (-b + disc) / (2.0 * a);
  241.  
  242.     if (t2 < MIN_T)
  243.         return (0);
  244.  
  245.     if (t1 < MIN_T)
  246.     {
  247.         nroots = 1;
  248.         t1 = t2;
  249.     }
  250.     else
  251.     {
  252.         nroots = 2;
  253.     }
  254.  
  255.     /*
  256.      * ensure that the points are between the two bounding planes...
  257.      */
  258.  
  259.     switch (nroots)
  260.     {
  261.     case 1:
  262.         VecAddS(t1, ray->dir, ray->pos, p);
  263.         d = VecDot(cd->w, p);
  264.  
  265.         if (d >= cd->min_d && d <= cd->max_d)
  266.         {
  267.             inter->t = t1;
  268.             inter->obj = obj;
  269.             inter->inside = 1;
  270.             return (1);
  271.         }
  272.         else
  273.         {
  274.             return (0);
  275.         }
  276.         break;
  277.  
  278.     case 2:
  279.         VecAddS(t1, ray->dir, ray->pos, p);
  280.         d = VecDot(cd->w, p);
  281.  
  282.         if (d >= cd->min_d && d <= cd->max_d)
  283.         {
  284.             inter->t = t1;
  285.             inter->obj = obj;
  286.             inter->inside = 0;
  287.             return (1);
  288.         }
  289.         else
  290.         {
  291.             VecAddS(t2, ray->dir, ray->pos, p);
  292.             d = VecDot(cd->w, p);
  293.             if (d >= cd->min_d && d <= cd->max_d)
  294.             {
  295.                 inter->t = t2;
  296.                 inter->obj = obj;
  297.                 inter->inside = 1;
  298.                 return (1);
  299.             }
  300.         }
  301.         return (0);
  302.     }
  303.     return (0);
  304. }
  305.  
  306.  
  307. Cone_normal(cone, ray, ip, normal)
  308. CONE           *cone;
  309. RAY            *ray;
  310. VECTOR         *ip;
  311. VECTOR         *normal;
  312. {
  313.     VECTOR          v;
  314.     double          t;
  315.  
  316.     /*
  317.      * fill in the real normal... Project the point onto the base plane.
  318.      * The normal is a vector from the basepoint through this point, plus
  319.      * the slope times the cone_w vector...
  320.      */
  321.  
  322.     t = -(VecDot(*ip, cone->w) + cone->base_d);
  323.     VecAddS(t, cone->w, *ip, v);
  324.     VecSub(v, cone->base, *normal);
  325.     VecNormalize(normal);
  326.     VecAddS(-cone->slope, cone->w, *normal, *normal);
  327.     VecNormalize(normal);
  328.     if (VecDot(ray->dir, *normal) >= 0)
  329.         VecNegate(*normal);
  330. }
  331.