Back to index

radiance  4R0+20100331
o_cone.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: o_cone.c,v 2.6 2004/06/28 10:07:17 greg Exp $";
00003 #endif
00004 /*
00005  *  o_cone.c - routine to determine ray intersection with cones.
00006  */
00007 
00008 #include "copyright.h"
00009 
00010 #include  "ray.h"
00011 #include  "otypes.h"
00012 #include  "rtotypes.h"
00013 #include  "cone.h"
00014 
00015 
00016 extern int
00017 o_cone(                     /* intersect ray with cone */
00018        OBJREC  *o,
00019        register RAY  *r
00020 )
00021 {
00022        FVECT  rox, rdx;
00023        double  a, b, c;
00024        double  root[2];
00025        int  nroots, rn;
00026        register CONE  *co;
00027        register int  i;
00028 
00029                                           /* get cone structure */
00030        co = getcone(o, 1);
00031 
00032        /*
00033         *     To intersect a ray with a cone, we transform the
00034         *  ray into the cone's normalized space.  This greatly
00035         *  simplifies the computation.
00036         *     For a cone or cup, normalization results in the
00037         *  equation:
00038         *
00039         *            x*x + y*y - z*z == 0
00040         *
00041         *     For a cylinder or tube, the normalized equation is:
00042         *
00043         *            x*x + y*y - r*r == 0
00044         *
00045         *     A normalized ring obeys the following set of equations:
00046         *
00047         *            z == 0               &&
00048         *            x*x + y*y >= r0*r0   &&
00049         *            x*x + y*y <= r1*r1
00050         */
00051 
00052                                    /* transform ray */
00053        multp3(rox, r->rorg, co->tm);
00054        multv3(rdx, r->rdir, co->tm);
00055                                    /* compute intersection */
00056 
00057        if (o->otype == OBJ_CONE || o->otype == OBJ_CUP) {
00058 
00059               a = rdx[0]*rdx[0] + rdx[1]*rdx[1] - rdx[2]*rdx[2];
00060               b = 2.0*(rdx[0]*rox[0] + rdx[1]*rox[1] - rdx[2]*rox[2]);
00061               c = rox[0]*rox[0] + rox[1]*rox[1] - rox[2]*rox[2];
00062 
00063        } else if (o->otype == OBJ_CYLINDER || o->otype == OBJ_TUBE) {
00064 
00065               a = rdx[0]*rdx[0] + rdx[1]*rdx[1];
00066               b = 2.0*(rdx[0]*rox[0] + rdx[1]*rox[1]);
00067               c = rox[0]*rox[0] + rox[1]*rox[1] - CO_R0(co)*CO_R0(co);
00068 
00069        } else { /* OBJ_RING */
00070 
00071               if (rdx[2] <= FTINY && rdx[2] >= -FTINY)
00072                      return(0);                  /* parallel */
00073               root[0] = -rox[2]/rdx[2];
00074               if (root[0] <= FTINY || root[0] >= r->rot)
00075                      return(0);                  /* distance check */
00076               b = root[0]*rdx[0] + rox[0];
00077               c = root[0]*rdx[1] + rox[1];
00078               a = b*b + c*c;
00079               if (a < CO_R0(co)*CO_R0(co) || a > CO_R1(co)*CO_R1(co))
00080                      return(0);                  /* outside radii */
00081               r->ro = o;
00082               r->rot = root[0];
00083               for (i = 0; i < 3; i++)
00084                      r->rop[i] = r->rorg[i] + r->rdir[i]*r->rot;
00085               VCOPY(r->ron, co->ad);
00086               r->rod = -rdx[2];
00087               r->rox = NULL;
00088               return(1);                         /* good */
00089        }
00090                                    /* roots for cone, cup, cyl., tube */
00091        nroots = quadratic(root, a, b, c);
00092 
00093        for (rn = 0; rn < nroots; rn++) {  /* check real roots */
00094               if (root[rn] <= FTINY)
00095                      continue;            /* too small */
00096               if (root[rn] >= r->rot)
00097                      break;               /* too big */
00098                                           /* check endpoints */
00099               for (i = 0; i < 3; i++) {
00100                      rox[i] = r->rorg[i] + root[rn]*r->rdir[i];
00101                      rdx[i] = rox[i] - CO_P0(co)[i];
00102               }
00103               b = DOT(rdx, co->ad); 
00104               if (b < 0.0)
00105                      continue;            /* before p0 */
00106               if (b > co->al)
00107                      continue;            /* after p1 */
00108               r->ro = o;
00109               r->rot = root[rn];
00110               VCOPY(r->rop, rox);
00111                                           /* get normal */
00112               if (o->otype == OBJ_CYLINDER)
00113                      a = CO_R0(co);
00114               else if (o->otype == OBJ_TUBE)
00115                      a = -CO_R0(co);
00116               else { /* OBJ_CONE || OBJ_CUP */
00117                      c = CO_R1(co) - CO_R0(co);
00118                      a = CO_R0(co) + b*c/co->al;
00119                      if (o->otype == OBJ_CUP) {
00120                             c = -c;
00121                             a = -a;
00122                      }
00123               }
00124               for (i = 0; i < 3; i++)
00125                      r->ron[i] = (rdx[i] - b*co->ad[i])/a;
00126               if (o->otype == OBJ_CONE || o->otype == OBJ_CUP)
00127                      for (i = 0; i < 3; i++)
00128                             r->ron[i] = (co->al*r->ron[i] - c*co->ad[i])
00129                                           /co->sl;
00130               a = DOT(r->ron, r->ron);
00131               if (a > 1.+FTINY || a < 1.-FTINY) {
00132                      c = 1./(.5 + .5*a);     /* avoid numerical error */
00133                      r->ron[0] *= c; r->ron[1] *= c; r->ron[2] *= c;
00134               }
00135               r->rod = -DOT(r->rdir, r->ron);
00136               r->pert[0] = r->pert[1] = r->pert[2] = 0.0;
00137               r->uv[0] = r->uv[1] = 0.0;
00138               r->rox = NULL;
00139               return(1);                  /* good */
00140        }
00141        return(0);
00142 }