home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / unix / volume18 / mtvraytrace / part02 / cone.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-03-26  |  6.2 KB  |  275 lines
  1. /***********************************************************************
  2.  * $Author: markv $
  3.  * $Revision: 1.3 $
  4.  * $Date: 88/10/04 14:30:02 $
  5.  * $Log:    cone.c,v $
  6.  * Revision 1.3  88/10/04  14:30:02  markv
  7.  * Fixed bug reported by koblas@mips.
  8.  * 
  9.  * Revision 1.2  88/09/15  09:33:02  markv
  10.  * Fixed bug reported by koblas@mips.  Cones which were specified with
  11.  * axis coincident with -y axis had problems.  Caused by incorrect 
  12.  * handling when finding orthogonal axes for the local coordinate system.
  13.  * 
  14.  * Revision 1.1  88/09/11  11:00:38  markv
  15.  * Initial revision
  16.  * 
  17.  ***********************************************************************/
  18. #include <stdio.h>
  19. #include <math.h>
  20. #include "defs.h"
  21. #include "extern.h"
  22.  
  23. typedef struct t_conedata {
  24.     Vec         cone_base ;
  25.     Flt        cone_base_radius ;
  26.     Flt        cone_base_d ;
  27.     Vec        cone_apex ;
  28.     Flt        cone_apex_radius ;
  29.     Vec        cone_u ;
  30.     Vec        cone_v ;
  31.     Vec        cone_w ;
  32.     Flt        cone_height ;
  33.     Flt        cone_slope ;
  34.     Flt        cone_min_d ;
  35.     Flt        cone_max_d ;
  36. } ConeData ;
  37.  
  38. int ConePrint ();
  39. int ConeIntersect ();
  40. int ConeNormal ();
  41.  
  42. ObjectProcs ConeProcs = {
  43.     ConePrint,
  44.     ConeIntersect,
  45.     ConeNormal,
  46. } ;
  47.  
  48. int 
  49. ConePrint(obj)
  50.  Object *obj ;
  51. {
  52.     ConeData * cp ;
  53.  
  54.     cp = (ConeData *) obj -> o_data ;
  55.  
  56.     printf("c %g %g %g %g %g %g %g %g\n", cp -> cone_base[0], 
  57.                   cp -> cone_base[1],
  58.                   cp -> cone_base[2],
  59.                   cp -> cone_base_radius,
  60.                   cp -> cone_apex[0],
  61.                   cp -> cone_apex[1],
  62.                   cp -> cone_apex[2],
  63.                   cp -> cone_apex_radius) ;
  64. }
  65.  
  66. ConeIntersect(obj, ray, hit)
  67.  Object * obj ;
  68.  Ray * ray ;
  69.  Isect * hit ;
  70. {
  71.     Ray tray ;
  72.     ConeData * cd ;
  73.     Vec V, P ;
  74.     Flt a, b, c, d, disc ;
  75.     Flt t1, t2 ;
  76.     int nroots ;
  77.  
  78.     cd = (ConeData *) (obj -> o_data) ;
  79.  
  80.     /*
  81.      * First, we get the coordinates of the ray origin in 
  82.      * the objects space....
  83.      */
  84.     
  85.     VecSub(ray -> P, cd -> cone_base, V) ;
  86.  
  87.     tray.P[0] = VecDot(V, cd -> cone_u) ;
  88.     tray.P[1] = VecDot(V, cd -> cone_v) ;
  89.     tray.P[2] = VecDot(V, cd -> cone_w) ;
  90.  
  91.     /*
  92.     VecAdd(ray -> P, ray -> D, V) ;
  93.     VecSub(V, cd -> cone_base, V) ;
  94.     */
  95.  
  96.     tray.D[0] = VecDot(ray -> D, cd -> cone_u) ;
  97.     tray.D[1] = VecDot(ray -> D, cd -> cone_v) ;
  98.     tray.D[2] = VecDot(ray -> D, cd -> cone_w) ;
  99.  
  100.     /*
  101.     VecSub(tray.D, tray.P, tray.D) ;
  102.     */
  103.  
  104.     a = tray.D[0] * tray.D[0] 
  105.         + tray.D[1] * tray.D[1] 
  106.         - cd -> cone_slope * cd -> cone_slope * tray.D[2] * tray.D[2] ;
  107.  
  108.     b = 2.0 * (tray.P[0] * tray.D[0] + tray.P[1] * tray.D[1] - 
  109.         cd -> cone_slope * cd -> cone_slope * tray.P[2] * tray.D[2]
  110.         - cd -> cone_base_radius * cd -> cone_slope * tray.D[2]) ;
  111.  
  112.     c = cd -> cone_slope * tray.P[2] + cd -> cone_base_radius ;
  113.     c = tray.P[0] * tray.P[0] + tray.P[1] * tray.P[1] - (c * c) ;
  114.  
  115.     disc = b * b - 4.0 * a * c ;
  116.  
  117.     if (disc < 0.0)
  118.         return (0) ;
  119.     
  120.     disc = (Flt) sqrt(disc) ;
  121.     t1 = (-b - disc) / (2.0 * a) ;
  122.     t2 = (-b + disc) / (2.0 * a) ;
  123.  
  124.     if (t2 < rayeps)
  125.         return (0) ;
  126.     if (t1 < rayeps) {
  127.         nroots = 1 ;
  128.         t1 = t2 ;
  129.     } else {
  130.         nroots = 2 ;
  131.     }
  132.         
  133.     /*
  134.      * ensure that the points are between the two bounding planes...
  135.      */
  136.     
  137.     switch(nroots) {
  138.     case 1:
  139.         RayPoint(ray, t1, P) ;
  140.         d = VecDot(cd -> cone_w, P) ;
  141.         if (d >= cd -> cone_min_d && d <= cd -> cone_max_d) {
  142.             hit -> isect_t = t1 ;    
  143.             hit -> isect_prim = obj ;
  144.             hit -> isect_surf = obj -> o_surf ;
  145.             hit -> isect_enter = 0 ;
  146.             return(1) ;
  147.         } else {
  148.             return(0) ;
  149.         }
  150.         break ;
  151.     case 2:
  152.         RayPoint(ray, t1, P) ;
  153.         d = VecDot(cd -> cone_w, P) ;
  154.         if (d >= cd -> cone_min_d && d <= cd -> cone_max_d) {
  155.             hit -> isect_t = t1 ;    
  156.             hit -> isect_prim = obj ;
  157.             hit -> isect_surf = obj -> o_surf ;
  158.             hit -> isect_enter = 1 ;
  159.             return(1) ;
  160.         } else {
  161.             RayPoint(ray, t2, P) ;
  162.             d = VecDot(cd -> cone_w, P) ;
  163.             if (d >= cd -> cone_min_d && d <= cd -> cone_max_d) {
  164.                 hit -> isect_t = t2 ;    
  165.                 hit -> isect_prim = obj ;
  166.                 hit -> isect_surf = obj -> o_surf ;
  167.                 hit -> isect_enter = 0 ;
  168.                 return(1) ;
  169.             }
  170.         }
  171.         return(0) ;
  172.     }
  173.     return(0) ;
  174. }
  175.  
  176. int
  177. ConeNormal(obj, hit, P, N)
  178.  Object * obj ;
  179.  Isect * hit ;
  180.  Point P, N ;
  181. {
  182.     Flt t ;
  183.     Vec V ;
  184.     ConeData * cd ;
  185.  
  186.     cd = (ConeData *) obj -> o_data ;
  187.  
  188.     /*
  189.      * fill in the real normal...
  190.      * Project the point onto the base plane.  The normal is
  191.      * a vector from the basepoint through this point, plus the slope
  192.      * times the cone_w vector...
  193.      */
  194.  
  195.     t = - (VecDot(P, cd -> cone_w) + cd -> cone_base_d) ;
  196.     VecAddS(t, cd -> cone_w, P, V) ;
  197.     VecSub(V, cd -> cone_base, N) ;
  198.     VecNormalize(N) ;
  199.     VecAddS(- cd -> cone_slope, cd -> cone_w, N, N) ;
  200.     VecNormalize(N) ;
  201. }
  202.  
  203. Object *
  204. MakeCone(basepoint, baseradius, apexpoint, apexradius) 
  205.  Vec basepoint, apexpoint ;
  206.  Flt baseradius, apexradius ;
  207. {
  208.     Object * obj ;
  209.     ConeData * cd ;
  210.     Flt dmin, dmax, d , ftmp;
  211.     Vec tmp ;
  212.     int i ;
  213.  
  214.     obj = (Object *) malloc (sizeof (Object)) ;
  215.     obj -> o_type = T_CONE ;
  216.     obj -> o_procs = & ConeProcs ;
  217.     obj -> o_surf = CurrentSurface ;
  218.  
  219.     cd = (ConeData *) malloc (sizeof(ConeData)) ;
  220.  
  221.     VecCopy(basepoint, cd -> cone_base) ;
  222.     VecCopy(apexpoint, cd -> cone_apex) ;
  223.  
  224.     cd -> cone_base_radius = baseradius ;
  225.     cd -> cone_apex_radius = apexradius ;
  226.  
  227.  
  228.     VecSub(apexpoint, basepoint, cd -> cone_w) ;
  229.     cd -> cone_height = VecNormalize(cd -> cone_w) ;
  230.     cd -> cone_slope =  (cd -> cone_apex_radius - cd -> cone_base_radius) /
  231.                 (cd -> cone_height) ;
  232.     cd -> cone_base_d = - VecDot(basepoint, cd -> cone_w) ;
  233.  
  234.     MakeVector(0.0, 0.0, 1.0, tmp) ;
  235.  
  236.     if (1.0 - fabs(VecDot(tmp, cd ->  cone_w)) < rayeps) {
  237.         MakeVector(0.0, 1.0, 0.0, tmp) ;
  238.     }
  239.  
  240.     /* find two axes which are at right angles to cone_w
  241.      */
  242.  
  243.     VecCross(cd -> cone_w, tmp, cd -> cone_u) ;
  244.     VecCross(cd -> cone_u, cd -> cone_w, cd -> cone_v) ;
  245.  
  246.     cd -> cone_min_d = VecDot(cd -> cone_w, cd -> cone_base) ;
  247.     cd -> cone_max_d = VecDot(cd -> cone_w, cd -> cone_apex) ;
  248.  
  249.     if (cd -> cone_max_d < cd -> cone_min_d) {
  250.         ftmp = cd -> cone_max_d ;
  251.         cd -> cone_max_d = cd -> cone_min_d ;
  252.         cd -> cone_min_d = ftmp ;
  253.     }
  254.  
  255.     obj -> o_data = (void *) cd ;
  256.  
  257.     for (i = 0 ; i < NSLABS ; i ++) {
  258.         dmin = HUGE ;
  259.         dmax = -HUGE ;
  260.         d = VecDot(basepoint, Slab[i]) - baseradius ;
  261.         if (d < dmin) dmin = d ; if (d > dmax) dmax = d ;
  262.         d = VecDot(basepoint, Slab[i]) + baseradius ;
  263.         if (d < dmin) dmin = d ; if (d > dmax) dmax = d ;
  264.         d = VecDot(apexpoint, Slab[i]) - apexradius ;
  265.         if (d < dmin) dmin = d ; if (d > dmax) dmax = d ;
  266.         d = VecDot(apexpoint, Slab[i]) + apexradius ;
  267.         if (d < dmin) dmin = d ; if (d > dmax) dmax = d ;
  268.  
  269.         obj -> o_dmin[i] = dmin ;
  270.         obj -> o_dmax[i] = dmax ;
  271.     }
  272.  
  273.     return(obj) ;
  274. }
  275.