home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / GRAPHICS / rayshade.lzh / cone.c < prev    next >
Text File  |  1990-05-08  |  7KB  |  296 lines
  1. /*
  2.  * cone.c
  3.  *
  4.  * Copyright (C) 1989, Craig E. Kolb
  5.  *
  6.  * This software may be freely copied, modified, and redistributed,
  7.  * provided that this copyright notice is preserved on all copies.
  8.  *
  9.  * There is no warranty or other guarantee of fitness for this software,
  10.  * it is provided solely .  Bug reports or fixes may be sent
  11.  * to the author, who may or may not act on them as he desires.
  12.  *
  13.  * You may not include this software in a program or other software product
  14.  * without supplying the source, or without informing the end-user that the
  15.  * source is available for no extra charge.
  16.  *
  17.  * If you modify this software, you should include a notice giving the
  18.  * name of the person performing the modification, the date of modification,
  19.  * and the reason for such modification.
  20.  *
  21.  * $Id: cone.c,v 3.0.1.5 90/04/09 14:30:08 craig Exp $
  22.  *
  23.  * $Log:    cone.c,v $
  24.  * Revision 3.0.1.5  90/04/09  14:30:08  craig
  25.  * patch5: Transformation information now stored locally.
  26.  * patch5: Canonical cone now truly canonical.
  27.  * 
  28.  * Revision 3.0.1.4  90/04/04  14:51:25  craig
  29.  * patch5: Fixed divide by zero problem in intcone().
  30.  * 
  31.  * Revision 3.0.1.3  90/02/12  13:20:58  craig
  32.  * patch4: Changes to avoid rotation about null axis.
  33.  * 
  34.  * Revision 3.0.1.2  89/12/06  16:33:44  craig
  35.  * patch2: Added calls to new error/warning routines.
  36.  * 
  37.  * Revision 3.0.1.1  89/11/18  14:08:09  craig
  38.  * patch1: Changes to reflect new names of transformation routines.
  39.  * 
  40.  * Revision 3.0  89/10/27  02:05:47  craig
  41.  * Baseline for first official release.
  42.  * 
  43.  */
  44. #include <stdio.h>
  45. #include <math.h>
  46. #include "typedefs.h"
  47. #include "funcdefs.h"
  48. #include "constants.h"
  49.  
  50. Object *
  51. makcone(surf, cent, ax, br, ar)
  52. char *surf;
  53. Vector *cent, *ax;
  54. double br, ar;
  55. {
  56.     Cone *cone;
  57.     Primitive *prim;
  58.     Object *newobj;
  59.     double tantheta, lprime, tlen, len, dtmp;
  60.     Vector axis, base, tmp;
  61.  
  62.     prim = mallocprim();
  63.     prim->surf = find_surface(surf);
  64.     prim->type = CONE;
  65.     newobj = new_object(NULL, CONE, (char *)prim, (Trans *)NULL);
  66.     cone = (Cone *)Malloc(sizeof(Cone));
  67.     prim->objpnt.p_cone = cone;
  68.  
  69.     /*
  70.      * To render a cone, we transform the desired cone into
  71.      * a canonical, Z-axis aligned, unit length, unit radius
  72.      * at the apex cone.
  73.      *
  74.      * Here, we construct the transformation matrix to map
  75.      * from cone<-->world space.
  76.      */
  77.  
  78.     /*
  79.      * The passed basepoint must be closer to the origin of the
  80.      * cone than the apex point, implying that the base radius
  81.      * must be smaller than the apex radius.  If the values passed
  82.      * reflect the opposite, we switch everything.
  83.      */
  84.     if(ar < br) {
  85.         tmp = *cent;
  86.         *cent = *ax;
  87.         *ax = tmp;
  88.         dtmp = br;
  89.         br = ar;
  90.         ar = dtmp;
  91.     } else if (equal(ar, br)) {
  92.         /*
  93.          * If the base and apex radii are equal, then we
  94.          * can treat the cone as a cylinder.
  95.          */
  96.         return makcyl(surf, cent, ax, br);
  97.     }
  98.     /*
  99.      * Find the axis and axis length.
  100.      */
  101.     vecsub(*ax, *cent, &axis);
  102.     len = normalize(&axis);
  103.     if (len < EPSILON) {
  104.         yywarning("Degenerate cone.\n");
  105.         free((char *)cone);
  106.         free((char *)prim);
  107.         free((char *)newobj);
  108.         return (Object *)0;
  109.     }
  110.     /*
  111.      * "tantheta" is the change in radius per unit length along
  112.      * the cone axis.
  113.      */
  114.     tantheta = (ar - br) / len;
  115.     /*
  116.      * lprime defines the distance along the axis where the first
  117.      * endcap should be placed.
  118.      */
  119.     lprime = br / tantheta;
  120.     /*
  121.      * Find the true base (origin) of the cone.
  122.      */
  123.     scalar_prod(-lprime, axis, &base);
  124.     vecadd(base, *cent, &base);
  125.     /*
  126.      * tlen is the total length of the cone.
  127.      */
  128.     tlen = lprime + len;
  129.     /*
  130.      * start_dist is the distance from the origin of the canonical
  131.      * cone at which the first endcap appears.
  132.      */
  133.     cone->start_dist = lprime / tlen;
  134.     /*
  135.      * Calculate transformation to map from cone space to world space.
  136.      */
  137.     if (abs(axis.z) == 1.) {
  138.         tmp.x = 1;
  139.         tmp.y = tmp.z = 0.;
  140.     } else {
  141.         tmp.x = axis.y;
  142.         tmp.y = -axis.x;
  143.         tmp.z = 0.;
  144.     }
  145.     init_trans(&cone->trans.obj2world);
  146.     RS_scale(&cone->trans.obj2world, ar, ar, tlen);
  147.     RS_rotate(&cone->trans.obj2world, &tmp, acos(axis.z));
  148.     RS_translate(&cone->trans.obj2world, &base);
  149.     invert_trans(&cone->trans.world2obj, &cone->trans.obj2world);
  150.  
  151.     return newobj;
  152. }
  153.  
  154. /*
  155.  * Ray-cone intersection test.  This routine is far from optimal, but
  156.  * it's straight-forward and it works...
  157.  */
  158. double
  159. intcone(pos, ray, obj)
  160. Vector           *pos, *ray;
  161. Primitive       *obj;
  162. {
  163.     double t1, t2, a, b, c, disc, zpos, et1, et2;
  164.     double x, y, distfact;
  165.     Ray newray;
  166.     Vector nray, npos;
  167.     extern double TransformRay();
  168.     extern unsigned long primtests[];
  169.     Cone *cone;
  170.  
  171.     primtests[CONE]++;
  172.     cone = obj->objpnt.p_cone;
  173.  
  174.     /*
  175.      * Transform ray from world to cone space.
  176.      */
  177.     newray.pos = *pos;
  178.     newray.dir = *ray;
  179.     distfact = TransformRay(&newray, &cone->trans.world2obj);
  180.     nray = newray.dir;
  181.     npos = newray.pos;
  182.  
  183.     a = nray.x * nray.x + nray.y * nray.y - nray.z*nray.z;
  184.     b = nray.x * npos.x + nray.y * npos.y - nray.z*npos.z;
  185.     c = npos.x*npos.x + npos.y*npos.y - npos.z*npos.z;
  186.  
  187.     if (equal(a, 0.)) {
  188.         /*
  189.          * Only one intersection point...
  190.          */
  191.         t1 = -c / b;
  192.         zpos = npos.z + t1 * nray.z;
  193.         if (t1 < EPSILON || zpos < cone->start_dist || zpos > 1.)
  194.             t1 = FAR_AWAY;
  195.         t2 = FAR_AWAY;
  196.     } else {
  197.         disc = b*b - a*c;
  198.         if(disc < 0.)
  199.             return 0.;        /* No possible intersection */
  200.         disc = sqrt(disc);
  201.         t1 = (-b + disc) / a;
  202.         t2 = (-b - disc) / a;
  203.         /*
  204.          * Clip intersection points.
  205.          */
  206.         zpos = npos.z + t1 * nray.z;
  207.         if (t1 < EPSILON || zpos < cone->start_dist || zpos > 1.)
  208.             t1 = FAR_AWAY;
  209.         zpos = npos.z + t2 * nray.z;
  210.         if (t2 < EPSILON || zpos < cone->start_dist || zpos > 1.)
  211.             t2 = FAR_AWAY;
  212.     }
  213.  
  214.     if (equal(nray.z, 0.)) {
  215.         t1 = min(t1, t2);
  216.         return (t1 == FAR_AWAY ? 0. : t1 / distfact);
  217.     }
  218.  
  219.     /*
  220.      * Find t for both endcaps.
  221.      */
  222.     et1 = (cone->start_dist - npos.z) / nray.z;
  223.     x = npos.x + et1 * nray.x;
  224.     y = npos.y + et1 * nray.y;
  225.     if (x*x + y*y > cone->start_dist*cone->start_dist)
  226.         et1 = FAR_AWAY;
  227.     et2 = (1. - npos.z) / nray.z;
  228.     x = npos.x + et2 * nray.x;
  229.     y = npos.y + et2 * nray.y;
  230.     if (x*x + y*y > 1.)
  231.         et2 = FAR_AWAY;
  232.  
  233.     t1 = min(t1, min(t2, min(et1, et2)));
  234.     return (t1 == FAR_AWAY ? 0. : t1 / distfact);
  235. }
  236.  
  237. /*
  238.  * Compute the normal to a cone at a given location on its surface.
  239.  */
  240. nrmcone(pos, obj, nrm)
  241. Vector *pos, *nrm;
  242. Primitive *obj;
  243. {
  244.     Cone *cone;
  245.     Vector npos;
  246.  
  247.     cone = obj->objpnt.p_cone;
  248.  
  249.     /*
  250.      * Transform intersection point to cone space.
  251.      */
  252.     npos = *pos;
  253.     transform_point(&npos, &cone->trans.world2obj);
  254.     
  255.     if (equal(npos.z, cone->start_dist)) {
  256.         nrm->x = nrm->y = 0.;
  257.         nrm->z = -1.;
  258.     } else if (equal(npos.z, 1)) {
  259.         nrm->x = nrm->y = 0.;
  260.         nrm->z = 1.;
  261.     } else {
  262.         /*
  263.          * The following is equal to
  264.          * (pos X (0, 0, 1)) X pos
  265.          */
  266.         nrm->x = npos.x * npos.z;
  267.         nrm->y = npos.y * npos.z;
  268.         nrm->z = -npos.x * npos.x - npos.y * npos.y;
  269.     }
  270.     /*
  271.      * Transform normal back to world space.
  272.      */
  273.     TransformNormal(nrm, &cone->trans.world2obj);
  274. }
  275.  
  276. /*
  277.  * Return the extent of a cone.
  278.  */
  279. coneextent(o, bounds)
  280. Primitive *o;
  281. double bounds[2][3];
  282. {
  283.     Cone *cone;
  284.  
  285.     cone = o->objpnt.p_cone;
  286.  
  287.     bounds[LOW][X] = bounds[LOW][Y] = -1;
  288.     bounds[HIGH][X] = bounds[HIGH][Y] = 1;
  289.     bounds[LOW][Z] = cone->start_dist;
  290.     bounds[HIGH][Z] = 1;
  291.     /*
  292.      * Transform bounding box to world space.
  293.      */
  294.     transform_bounds(&cone->trans.obj2world, bounds);
  295. }
  296.