home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / rad386 / radiosit / src / objects.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-22  |  28.7 KB  |  1,149 lines
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <ctype.h>
  4. #include <math.h>
  5. #include <malloc.h>
  6. #include "GraphicsGems.h"
  7. #include "data_structure.h"
  8. #include "objects.h"
  9.  
  10. #include "raddecl.h"
  11.  
  12. /************** Begin object specific intersections **************/
  13. static double intersect_polygon(p,ray_start,ray_direction,pu,pv)
  14. Polygon *p;
  15. Point3 *ray_start;
  16. Vector3 *ray_direction;
  17. double *pu,*pv;
  18. {
  19.     double ray_polygon();
  20.     *pu = *pv = 0.5;
  21.     return(
  22.         ray_polygon(ray_start,ray_direction,
  23.             &(p->plane_normal), p->plane_constant,
  24.             p->dominant_axis_flag,
  25.             p->nvertices,
  26.             p->vertex_list)
  27.     );
  28. }
  29.  
  30. static double intersect_quadrilateral(q,ray_start,ray_direction,pu,pv)
  31. Quadrilateral *q;
  32. Point3 *ray_start;
  33. Vector3 *ray_direction;
  34. double *pu,*pv;
  35. {
  36.     double ray_quadrilateral();
  37.     return(
  38.         ray_quadrilateral(ray_start,ray_direction,
  39.             &(q->plane_normal), q->plane_constant,
  40.             q->Du0,q->Du1,q->Du2,q->Dux,q->Duy,
  41.             q->Dv0,q->Dv1,q->Dv2,q->Dvx,q->Dvy,
  42.             &(q->Qux),&(q->Quy),&(q->Qvx),&(q->Qvy),
  43.             &(q->Na),&(q->Nb),&(q->Nc),
  44.             pu,pv)
  45.     );
  46. }
  47.  
  48. static double intersect_sphere(s,ray_start,ray_direction,pu,pv)
  49. Sphere *s;
  50. Point3 *ray_start;
  51. Vector3 *ray_direction;
  52. double *pu,*pv;
  53. {
  54.     double ray_sphere();
  55.     return(
  56.         ray_sphere(ray_start,ray_direction,
  57.             s->radius, &(s->center), 0, pu,pv)
  58.     );
  59. }
  60.  
  61. static double intersect_bubble(b,ray_start,ray_direction,pu,pv)
  62. Sphere *b;
  63. Point3 *ray_start;
  64. Vector3 *ray_direction;
  65. double *pu,*pv;
  66. {
  67.     double ray_sphere();
  68.     return(
  69.         ray_sphere(ray_start,ray_direction,
  70.             b->radius, &(b->center), 1, pu,pv)
  71.     );
  72. }
  73.  
  74. static void ray_transform(s_in,d_in,m,s_out,d_out)
  75. Point3 *s_in,*s_out;
  76. Vector3 *d_in,*d_out;
  77. Matrix4 *m;
  78. {
  79.         *d_out = transform_vector(d_in,m);
  80.         *s_out = *s_in;
  81.         V3MulPointByMatrix(s_out,m);
  82. }
  83.  
  84. static double intersect_cylinder(c,ray_start,ray_direction,pu,pv)
  85. Cylinder *c;
  86. Point3 *ray_start;
  87. Vector3 *ray_direction;
  88. double *pu,*pv;
  89. {
  90.     double ray_cylinder();
  91.     Point3 s;
  92.     Vector3 d;
  93.     ray_transform(ray_start,ray_direction,&(c->world_to_object),&s,&d);
  94.     return(
  95.         ray_cylinder(&s,&d,c->radius,c->height,0,pu,pv)
  96.     );
  97. }
  98.  
  99. static double intersect_tube(tube,ray_start,ray_direction,pu,pv)
  100. Cylinder *tube;
  101. Point3 *ray_start;
  102. Vector3 *ray_direction;
  103. double *pu,*pv;
  104. {
  105.     double ray_cylinder();
  106.     Point3 s;
  107.     Vector3 d;
  108.     ray_transform(ray_start,ray_direction,&(tube->world_to_object),&s,&d);
  109.     return(
  110.         ray_cylinder(&s,&d,tube->radius,tube->height,1,pu,pv)
  111.     );
  112. }
  113.  
  114. static double intersect_cone(c,ray_start,ray_direction,pu,pv)
  115. Cone *c;
  116. Point3 *ray_start;
  117. Vector3 *ray_direction;
  118. double *pu,*pv;
  119. {
  120.     double ray_cone();
  121.     Point3 s;
  122.     Vector3 d;
  123.     ray_transform(ray_start,ray_direction,&(c->world_to_object),&s,&d);
  124.     return(
  125.         ray_cone(&s,&d,c->distance,0,pu,pv)
  126.     );
  127. }
  128.  
  129. static double intersect_cup(c,ray_start,ray_direction,pu,pv)
  130. Cone *c;
  131. Point3 *ray_start;
  132. Vector3 *ray_direction;
  133. double *pu,*pv;
  134. {
  135.     double ray_cone();
  136.     Point3 s;
  137.     Vector3 d;
  138.     ray_transform(ray_start,ray_direction,&(c->world_to_object),&s,&d);
  139.     return(
  140.         ray_cone(&s,&d,c->distance,1,pu,pv)
  141.     );
  142. }
  143.  
  144. static double intersect_ring(r,ray_start,ray_direction,pu,pv)
  145. Ring *r;
  146. Point3 *ray_start;
  147. Vector3 *ray_direction;
  148. double *pu,*pv;
  149. {
  150.     double ray_ring();
  151.     Point3 s;
  152.     Vector3 d;
  153.     ray_transform(ray_start,ray_direction,&(r->world_to_object),&s,&d);
  154.     return(
  155.         ray_ring(&s,&d,r->inner_radius,r->outer_radius,pu,pv)
  156.     );
  157. }
  158.  
  159. static double intersect_other_geometry(dummy,ray_start,ray_direction,pu,pv)
  160. void *dummy;
  161. Point3 *ray_start;
  162. Vector3 *ray_direction;
  163. double *pu,*pv;
  164. {
  165.     error("This Geometry is unimplemented.");
  166. }
  167.  
  168.  
  169. /************* End object specific intersections **************/
  170. /************* Begin Object Specific Preprocess  **************/
  171. /*
  172. get_view_trans_matrix(U,V,N,m)
  173. Vector3 *U,*V,*N;
  174. Matrix4 *m;
  175. {
  176.         m->element[0][0]=U->x;
  177.         m->element[1][0]=U->y;
  178.         m->element[2][0]=U->z;
  179.         m->element[0][1]=V->x;
  180.         m->element[1][1]=V->y;
  181.         m->element[2][1]=V->z;
  182.         m->element[0][2]=N->x;
  183.         m->element[1][2]=N->y;
  184.         m->element[2][2]=N->z;
  185.  
  186.     m->element[0][3]= m->element[1][3]= m->element[2][3]=
  187.         m->element[3][0]= m->element[3][1]= m->element[3][2]=0.;
  188.     m->element[3][3]=1.;
  189. }
  190. */
  191. #define update_scene_bound(bb,b) {\
  192.     b->min.x = MIN(bb->min.x,b->min.x); b->max.x = MAX(bb->max.x,b->max.x);\
  193.     b->min.y = MIN(bb->min.y,b->min.y); b->max.y = MAX(bb->max.y,b->max.y);\
  194.     b->min.z = MIN(bb->min.z,b->min.z); b->max.z = MAX(bb->max.z,b->max.z);\
  195. }
  196.  
  197. static void null_preprocess(grid,un,vn,spec_struc,bb,b)
  198. Surface_Grid_structure *grid;
  199. int un,vn;
  200. void *spec_struc;
  201. Box3 *bb,*b;
  202. {
  203.     error("No Preprocess Implemented Yet for this Object.");
  204. }
  205. /*
  206. static double area (p1,p2,p3,p4)
  207. Point3 *p1,*p2,*p3,*p4;
  208. {
  209.     Vector3 v1,v2,N;
  210.     double a;
  211.     V3Sub(p2,p1,&v1); V3Sub(p4,p1,&v2); V3Cross(&v1,&v2,&N);
  212.     a = V3Length(&N);
  213.     V3Sub(p4,p3,&v1);V3Sub(p2,p3,&v2); V3Cross(&v1,&v2,&N);
  214.     a += V3Length(&N);
  215.     return(a/2.0);
  216.     }
  217. */
  218. static int get_dominant_axis(N)
  219. Vector3 *N;
  220. {
  221.     double maxval;
  222.     Point3 anorm;
  223.         /* Computes Dominant axis.
  224.            0 - for X-axis, 1 - for Y-axis, 2 - for Z axis */
  225.     anorm.x = fabs(N->x);
  226.     anorm.y = fabs(N->y);
  227.     anorm.z = fabs(N->z);
  228.     maxval  = MAX(anorm.x, MAX(anorm.y,anorm.z));
  229.     return((maxval == anorm.x)?0:((maxval == anorm.y)?1:2));
  230. }
  231. static double trapez_area(n,vlist,N,dominant_axis)
  232. /* Polygon area computation by trapezium area summation method */
  233. int n;
  234. Point3 vlist[];
  235. Vector3 *N;
  236. int dominant_axis;
  237. {
  238.     int i;
  239.     double area = 0.;
  240.     Point2 cur,last;
  241.     
  242.     project_point(last,vlist[n-1],dominant_axis);
  243.     for(i=0;i<n;i++,last=cur){
  244.         project_point(cur,vlist[i],dominant_axis);
  245.         area += (cur.x-last.x)*(cur.y+last.y);
  246.     }
  247.     area /= ((dominant_axis==0)?N->x:((dominant_axis==1)?N->y:N->z));
  248.     return(0.5*fabs(area));
  249. }
  250.  
  251. static void quadrilateral_preprocess(grid,un,vn,q,bb,b)
  252. Surface_Grid_structure *grid;
  253. int un,vn;
  254. Quadrilateral *q;
  255. Box3 *bb,*b;
  256. {
  257.     int i,j,n;
  258.     Point3 Pa, Pb,Pc,Pd;
  259.     double u,v,du=1.0/un,dv=1.0/vn;
  260.     double area();
  261.     /* Compute Bounding Box */
  262.     bb->max.x=MAX(MAX(q->P00.x,q->P10.x),MAX(q->P11.x,q->P01.x));
  263.         bb->max.y=MAX(MAX(q->P00.y,q->P10.y),MAX(q->P11.y,q->P01.y));
  264.     bb->max.z=MAX(MAX(q->P00.z,q->P10.z),MAX(q->P11.z,q->P01.z));
  265.  
  266.     bb->min.x=MIN(MIN(q->P00.x,q->P10.x),MIN(q->P11.x,q->P01.x));
  267.         bb->min.y=MIN(MIN(q->P00.y,q->P10.y),MIN(q->P11.y,q->P01.y));
  268.     bb->min.z=MIN(MIN(q->P00.z,q->P10.z),MIN(q->P11.z,q->P01.z));
  269.     update_scene_bound(bb,b);
  270.     /* Compute Other parameters. */
  271.            Pa.x=q->P00.x-q->P10.x+q->P11.x-q->P01.x;
  272.            Pa.y=q->P00.y-q->P10.y+q->P11.y-q->P01.y;
  273.            Pa.z=q->P00.z-q->P10.z+q->P11.z-q->P01.z;
  274.            Pb.x=q->P10.x-q->P00.x;
  275.            Pb.y=q->P10.y-q->P00.y;
  276.            Pb.z=q->P10.z-q->P00.z;
  277.            Pc.x=q->P01.x-q->P00.x;
  278.            Pc.y=q->P01.y-q->P00.y;
  279.            Pc.z=q->P01.z-q->P00.z;
  280.            Pd = q->P00;
  281.            V3Cross(&Pb,&Pc,&(q->plane_normal));
  282.     V3Normalize(&(q->plane_normal));
  283.     /* Compute grid area. */
  284.     for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
  285.         Point3 f[4];
  286.         f[0] = get_bilinear(q,(double)0.,v);
  287.         f[3] = get_bilinear(q,(double)0.,v+dv);
  288.         for(u=0.0,j=0;j<un;j++,u+=du,n++){
  289.             f[1] = get_bilinear(q,u+du,v);
  290.             f[2] = get_bilinear(q,u+du,v+dv);
  291.             grid[n].area=trapez_area(4,f,&(q->plane_normal),
  292.                 get_dominant_axis(&(q->plane_normal)));
  293.             f[0] = f[1]; f[3] = f[2];
  294.         }
  295.     }
  296.         /* Compute the Coordinate Transformation Matrix */
  297.     /*
  298.     {
  299.         Vector3 U,V,W;
  300.         Matrix4 m;
  301.         V=Pc;
  302.         V3Normalize(&V);
  303.         V3Cross(&V,&(q->plane_normal),&U);
  304.         get_view_trans_matrix(&U,&V,&(q->plane_normal),&m);
  305.         U.x=V.y=W.z=1.0;
  306.         U.y=U.z=V.x=V.z=W.x=W.y=0.0;
  307.         U = transform_vector(&U,&m);
  308.         V = transform_vector(&V,&m);
  309.         W = transform_vector(&W,&m);
  310.         get_view_trans_matrix(&U,&V,&W,&(q->local_trans));
  311.     }
  312.     */
  313.     align_Z_to_axis(&(q->plane_normal),&(q->local_trans)); 
  314.  
  315.            q->plane_constant = -V3Dot(&(q->plane_normal),&(q->P00));
  316.  
  317.            V3Cross(&Pa,&(q->plane_normal),&(q->Na));
  318.            V3Cross(&Pb,&(q->plane_normal),&(q->Nb));
  319.            V3Cross(&Pc,&(q->plane_normal),&(q->Nc));
  320.  
  321.            q->Du0 = V3Dot(&(q->Nc),&Pd);
  322.            q->Du1 = V3Dot(&(q->Na),&Pd)+V3Dot(&(q->Nc),&Pb);
  323.            q->Du2 = V3Dot(&(q->Na),&Pb);
  324.            if (fabs(q->Du2) > EPSILON){
  325.                       double frac = 1.0/q->Du2;
  326.                       q->Quy.x = -q->Nc.x*frac; 
  327.         q->Quy.y = -q->Nc.y*frac; 
  328.         q->Quy.z = -q->Nc.z*frac;
  329.                 q->Duy=q->Du0*frac;
  330.                       frac /= 2.0;
  331.                 q->Qux.x=q->Na.x*frac; 
  332.         q->Qux.y=q->Na.y*frac; 
  333.         q->Qux.z=q->Na.z*frac;
  334.                 q->Dux = -q->Du1*frac;
  335.            }
  336.            q->Dv0 = V3Dot(&(q->Nb),&Pd);
  337.            q->Dv1 = V3Dot(&(q->Na),&Pd)+V3Dot(&(q->Nb),&Pc);
  338.            q->Dv2 = V3Dot(&(q->Na),&Pc);
  339.            if (fabs(q->Dv2) > EPSILON){
  340.                       double frac = 1.0/q->Dv2;
  341.                       q->Qvy.x = -q->Nb.x*frac; 
  342.         q->Qvy.y = -q->Nb.y*frac; 
  343.         q->Qvy.z = -q->Nb.z*frac;
  344.                       q->Dvy=q->Dv0*frac;
  345.                       frac /= 2.0;
  346.                       q->Qvx.x=q->Na.x*frac; 
  347.         q->Qvx.y=q->Na.y*frac; 
  348.         q->Qvx.z=q->Na.z*frac;
  349.                       q->Dvx = -q->Dv1*frac;
  350.            }
  351. }
  352.  
  353. static void sphere_preprocess(grid,un,vn,s,bb,b)
  354. Surface_Grid_structure *grid;
  355. int un,vn;
  356. Sphere *s;
  357. Box3 *bb,*b;
  358. {
  359.     int i,j,n;
  360.     double area,mult,prevcos=1.0,nextcos;
  361.     double v,du=1.0/un,dv=1.0/vn;
  362.     /* Compute Bounding Box */
  363.     bb->min.x = s->center.x - s->radius;
  364.     bb->min.y = s->center.y - s->radius;
  365.     bb->min.z = s->center.z - s->radius;
  366.  
  367.     bb->max.x = s->center.x + s->radius;
  368.     bb->max.y = s->center.y + s->radius;
  369.     bb->max.z = s->center.z + s->radius;
  370.     update_scene_bound(bb,b);
  371.     /* Compute Area of the Grid Elements. */
  372.     mult = 2.*PI*s->radius*s->radius*du;
  373.     for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
  374.         nextcos = cos(PI*(v+dv));
  375.         area=mult*fabs(prevcos-nextcos);
  376.         prevcos=nextcos;
  377.         for(j=0;j<un;j++,n++)
  378.             grid[n].area=area;
  379.     }
  380. }
  381.  
  382. void translate(dx,dy,dz,m)
  383. double dx,dy,dz;
  384. Matrix4 *m;
  385. {
  386.     m->element[0][0] =
  387.         m->element[1][1] =
  388.             m->element[2][2] =
  389.                 m->element[3][3] = 1.;
  390.     m->element[0][1] = m->element[0][2] = m->element[0][3] =
  391.         m->element[1][0] = m->element[1][2] = m->element[1][3] =
  392.             m->element[2][0] = m->element[2][1] = m->element[2][3] = 0.;
  393.     m->element[3][0] = dx;
  394.     m->element[3][1] = dy;
  395.     m->element[3][2] = dz;
  396. }
  397.  
  398. void scale(sx,sy,sz,m)
  399. double sx,sy,sz;
  400. Matrix4 *m;
  401. {
  402.     m->element[0][0] = sx;
  403.     m->element[1][1] = sy;
  404.     m->element[2][2] = sz;
  405.  
  406.     m->element[3][3] = 1.;
  407.     m->element[0][1] = m->element[0][2] = m->element[0][3] =
  408.         m->element[1][0] = m->element[1][2] = m->element[1][3] =
  409.             m->element[2][0] = m->element[2][1] = m->element[2][3] =
  410.                 m->element[3][0] = m->element[3][1] = m->element[3][2] = 0.;
  411. }
  412.  
  413. void bound_transform(b,m)
  414. Box3 *b;
  415. Matrix4 *m;
  416. {
  417.     Point3 p;
  418.     Box3 bb;
  419.     bb = *b; 
  420.     b->min.x = b->min.y = b->min.z =  LARGE;
  421.     b->max.x = b->max.y = b->max.z = -LARGE;
  422.     p = bb.min; V3MulPointByMatrix(&p,m);
  423.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  424.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  425.     p = bb.min; p.x = bb.max.x; V3MulPointByMatrix(&p,m);
  426.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  427.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  428.     p = bb.min; p.y = bb.max.y; V3MulPointByMatrix(&p,m);
  429.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  430.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  431.     p = bb.min; p.x = bb.max.x; p.y = bb.max.y; V3MulPointByMatrix(&p,m);
  432.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  433.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  434.  
  435.     p = bb.max; V3MulPointByMatrix(&p,m);
  436.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  437.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  438.     p = bb.max; p.x = bb.min.x; V3MulPointByMatrix(&p,m);
  439.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  440.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  441.     p = bb.max; p.y = bb.min.y; V3MulPointByMatrix(&p,m);
  442.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  443.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  444.     p = bb.max; p.x = bb.min.x; p.y = bb.min.y; V3MulPointByMatrix(&p,m);
  445.     b->min.x = MIN(b->min.x,p.x); b->min.y = MIN(b->min.y,p.y); b->min.z = MIN(b->min.z,p.z);
  446.     b->max.x = MAX(b->max.x,p.x); b->max.y = MAX(b->max.y,p.y); b->max.z = MAX(b->max.z,p.z);
  447. }
  448.  
  449. static void cylinder_preprocess(grid,un,vn,c,bb,b)
  450. Surface_Grid_structure *grid;
  451. int un,vn;
  452. Cylinder *c;
  453. Box3 *bb,*b;
  454. {
  455.     Vector3 dir;
  456.     Matrix4 m,m1;
  457.     double du=1.0/un,dv=1.0/vn,area;
  458.     int i,j,n;
  459.  
  460.     /* Compute derived data */
  461.     c->height = V3DistanceBetween2Points(&(c->base_point1),&(c->base_point2));
  462.     /* Compute Transformations */
  463.     V3Sub(&(c->base_point2),&(c->base_point1),&dir); V3Normalize(&dir);
  464.     translate(-c->base_point1.x,-c->base_point1.y,-c->base_point1.z,&m);
  465.     align_axis_to_Z(&dir,&m1);
  466.     V3MatMul(&m,&m1,&(c->world_to_object));
  467.     align_Z_to_axis(&dir,&m);
  468.     translate(c->base_point1.x,c->base_point1.y,c->base_point1.z,&m1);
  469.     V3MatMul(&m,&m1,&(c->object_to_world));
  470.  
  471.     /* Compute Bounding Box */
  472.     bb->min.x = bb->min.y = -c->radius; bb->min.z = 0.;
  473.     bb->max.x = bb->max.y = c->radius; bb->max.z = c->height;
  474.     bound_transform(bb,&(c->object_to_world));
  475.     update_scene_bound(bb,b);
  476.     /* Compute Area of the Grid Elements */
  477.     area = 2.*PI*du*c->radius*c->height*dv;
  478.     for(n=0,i=0;i<vn;i++)
  479.         for(j=0;j<un;j++,n++)
  480.             grid[n].area=area;
  481. }
  482.  
  483. static void cone_preprocess(grid,un,vn,c,bb,b)
  484. Surface_Grid_structure *grid;
  485. int un,vn;
  486. Cone *c;
  487. Box3 *bb,*b;
  488. {
  489.     Point3 origin;
  490.     Matrix4 m,m1,m2,m3;
  491.     double H1,H2,tan_theta;
  492.     double v,du=1.0/un,dv=1.0/vn,mult;
  493.     int i,j,n;
  494.  
  495.     /* Make sure 0 <= apex_radius <= base_radius. */
  496.     if (c->apex_radius > c->base_radius){
  497.         /* interchange */
  498.         Point3 p;
  499.         double r;
  500.         p=c->base_point; c->base_point = c->apex_point; c->apex_point = p;
  501.         r=c->base_radius;c->base_radius=c->apex_radius; c->apex_radius= r;
  502.     }
  503.  
  504.     /* Compute derived data */
  505.     H1 = V3DistanceBetween2Points(&(c->base_point),&(c->apex_point));
  506.     tan_theta = (c->base_radius-c->apex_radius)/H1;
  507.     H2 = c->apex_radius/tan_theta;
  508.     c->distance = H2/(H1+H2);
  509.     V3Sub(&(c->base_point),&(c->apex_point),&(c->axis)); V3Normalize(&(c->axis));
  510.  
  511.     /* Compute Transformations */
  512.     origin.x = c->apex_point.x - H2*c->axis.x;
  513.     origin.y = c->apex_point.y - H2*c->axis.y;
  514.     origin.z = c->apex_point.z - H2*c->axis.z;
  515.     translate(-origin.x,-origin.y,-origin.z,&m);
  516.     align_axis_to_Z(&(c->axis),&m1);
  517.     scale(1./c->base_radius,1./c->base_radius,1./(H1+H2),&m2);
  518.     V3MatMul(&m,&m1,&m3);V3MatMul(&m3,&m2,&(c->world_to_object));
  519.     
  520.     scale(c->base_radius,c->base_radius,(H1+H2),&m);
  521.     align_Z_to_axis(&(c->axis),&m1);
  522.     translate(origin.x,origin.y,origin.z,&m2);
  523.     V3MatMul(&m,&m1,&m3);V3MatMul(&m3,&m2,&(c->object_to_world));
  524.  
  525.     /* Compute Bounding Box */
  526.     bb->min.x = bb->min.y = -1.; bb->min.z = c->distance;
  527.     bb->max.x = bb->max.y = bb->max.z = 1.;
  528.     bound_transform(bb,&(c->object_to_world));
  529.     update_scene_bound(bb,b);
  530.  
  531.     /* Compute Area of the Grid Elements */
  532.     mult = PI*du*sqrt((c->base_radius-c->apex_radius)*(c->base_radius-c->apex_radius)+H1*H1);
  533.     for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
  534.         double area=mult*((c->base_radius-c->apex_radius)*((v+dv)*(v+dv)-v*v)+2.0*c->apex_radius*dv);
  535.         for(j=0;j<un;j++,n++)
  536.             grid[n].area=area;
  537.     }
  538. }
  539.  
  540. static void ring_preprocess(grid,un,vn,r,bb,b)
  541. Surface_Grid_structure *grid;
  542. int un,vn;
  543. Ring *r;
  544. Box3 *bb,*b;
  545. {
  546.     Matrix4 m,m1;
  547.     double extent;
  548.     double v,du=1.0/un,dv=1.0/vn,mult;
  549.     int i,j,n;
  550.  
  551.     if (r->inner_radius > r->outer_radius){
  552.         double rad = r->outer_radius;
  553.         r->outer_radius = r->inner_radius;
  554.         r->inner_radius = rad;
  555.     }
  556.     V3Normalize(&(r->normal));
  557.  
  558.     /* Compute Transformations */
  559.     translate(-r->center.x,-r->center.y,-r->center.z,&m);
  560.     align_axis_to_Z(&(r->normal),&m1);
  561.     V3MatMul(&m,&m1,&(r->world_to_object));
  562.  
  563.     align_Z_to_axis(&(r->normal),&m);
  564.     translate(r->center.x,r->center.y,r->center.z,&m1);
  565.     V3MatMul(&m,&m1,&(r->object_to_world));
  566.  
  567.     /* Compute Bounding Box */
  568.     extent = r->outer_radius * sqrt(1. - r->normal.x*r->normal.x);
  569.     bb->min.x = r->center.x - extent; 
  570.     bb->max.x = r->center.x + extent;
  571.  
  572.     extent = r->outer_radius * sqrt(1. - r->normal.y*r->normal.y);
  573.     bb->min.y = r->center.y - extent;
  574.     bb->max.y = r->center.y + extent;
  575.  
  576.     extent = r->outer_radius * sqrt(1. - r->normal.z*r->normal.z);
  577.     bb->min.z = r->center.z - extent;
  578.     bb->max.z = r->center.z + extent;
  579.     update_scene_bound(bb,b);
  580.     /* Compute Area of the Grid Elements */
  581.     mult = PI*du*(r->outer_radius*r->outer_radius - r->inner_radius*r->inner_radius);
  582.     for(n=0,v=0.0,i=0;i<vn;i++,v+=dv){
  583.         double area=mult*((v+dv)*(v+dv)-v*v);
  584.         for(j=0;j<un;j++,n++)
  585.             grid[n].area=area;
  586.     }
  587. }
  588.  
  589. static void polygon_preprocess(grid,un,vn,p,bb,b)
  590. Surface_Grid_structure *grid;
  591. int un,vn;
  592. Polygon *p;
  593. Box3 *bb,*b;
  594. {
  595.     int i,j;
  596.     Point3 *point=p->vertex_list;
  597.     /* Compute derived data */
  598.         /* Normal by Newels Method */
  599.     p->plane_normal.x = p->plane_normal.y = p->plane_normal.z = 0.;
  600.     for(i=0,j=1;i<p->nvertices;i++,j++){
  601.         if (j==p->nvertices)j=0;
  602.         p->plane_normal.x += (point[i].y - point[j].y)*(point[i].z+point[j].z);
  603.         p->plane_normal.y += (point[i].z - point[j].z)*(point[i].x+point[j].x);
  604.         p->plane_normal.z += (point[i].x - point[j].x)*(point[i].y+point[j].y);
  605.     }
  606.     V3Normalize(&(p->plane_normal));
  607.     p->plane_constant = -V3Dot(&(p->plane_normal),point);
  608.     p->dominant_axis_flag = get_dominant_axis(&(p->plane_normal));
  609.     align_Z_to_axis(&(p->plane_normal),&(p->local_trans));
  610.     /* Compute Bounding Box */
  611.     bb->min = bb->max = *point;
  612.     for(i=1;i<p->nvertices;i++){
  613.         bb->min.x = MIN(bb->min.x,point[i].x);
  614.         bb->max.x = MAX(bb->max.x,point[i].x);
  615.         bb->min.y = MIN(bb->min.y,point[i].y);
  616.         bb->max.y = MAX(bb->max.y,point[i].y);
  617.         bb->min.z = MIN(bb->min.z,point[i].z);
  618.         bb->max.z = MAX(bb->max.z,point[i].z);
  619.     }
  620.     update_scene_bound(bb,b);
  621.     /* Compute Area of the Grid Elements */
  622.     if (un > 1 || vn > 1)
  623.         error("Does not know any way to grid arbitrary polygon.");
  624.     grid[0].area=
  625.     trapez_area(p->nvertices,p->vertex_list,&(p->plane_normal),p->dominant_axis_flag);
  626. }
  627. #undef update_scene_bound
  628. /************** End Object Specific Preprocess   **************/
  629. /************** Begin Object Specific Input      **************/
  630. /*
  631.         for quadrilateral
  632.             P00 
  633.             P10 
  634.             P11 
  635.             P01 (Vertices Counter Clockwise) 
  636.         for sphere and bubble
  637.             centre
  638.             radius
  639.         for cylinder and tube
  640.             base point1
  641.             base point2
  642.             radius
  643.         for cone and cup
  644.             apex point
  645.             apex radius
  646.             base point
  647.             base radius
  648.         for ring
  649.             center
  650.             normal
  651.             inner radius
  652.             outer radius
  653.         for polygon
  654.             number of vertices
  655.             vertices
  656.             .
  657.             .
  658.         for others
  659.             nill
  660. */
  661. static void *quadrilateral_input()
  662. {
  663.     Quadrilateral *q=(Quadrilateral *)malloc(sizeof(Quadrilateral));
  664.  
  665.     if (q==NULL)error("Cannot allocate Quadrilateral.");
  666.     scanf("%lf%lf%lf",&(q->P00.x),&(q->P00.y),&(q->P00.z));
  667.     scanf("%lf%lf%lf",&(q->P10.x),&(q->P10.y),&(q->P10.z));
  668.     scanf("%lf%lf%lf",&(q->P11.x),&(q->P11.y),&(q->P11.z));
  669.     scanf("%lf%lf%lf",&(q->P01.x),&(q->P01.y),&(q->P01.z));
  670.  
  671.     return((void *)q);
  672. }
  673.  
  674. static void *polygon_input()
  675. {
  676.     int i;
  677.     Point3 *point;
  678.     Polygon *p=(Polygon *)malloc(sizeof(Polygon));
  679.  
  680.     if (p==NULL)error("Cannot allocate Polygon.");
  681.     scanf("%d",&(p->nvertices));
  682.     point = p->vertex_list = (Point3 *)malloc(sizeof(Point3) * p->nvertices);
  683.     if (point == NULL) error("cannot allocate Polygon Vertices.");
  684.     for (i=0; i < p->nvertices; i++,point++)
  685.         scanf("%lf%lf%lf",&(point->x),&(point->y),&(point->z));
  686.     return((void *)p);
  687. }
  688.  
  689. static void *sphere_input()
  690. {
  691.     Sphere *s=(Sphere *)malloc(sizeof(Sphere));
  692.  
  693.     if (s==NULL)error("Cannot allocate Sphere.");
  694.     scanf("%lf%lf%lf",&(s->center.x),&(s->center.y),&(s->center.z));
  695.     scanf("%lf",&(s->radius));
  696.  
  697.     return((void *)s);
  698. }
  699.  
  700. static void *cylinder_input()
  701. {
  702.     Cylinder *c=(Cylinder *)malloc(sizeof(Cylinder));
  703.  
  704.     if (c==NULL)error("Cannot allocate Cylinder/Tube.");
  705.     scanf("%lf%lf%lf",&(c->base_point1.x),&(c->base_point1.y),&(c->base_point1.z));
  706.     scanf("%lf%lf%lf",&(c->base_point2.x),&(c->base_point2.y),&(c->base_point2.z));
  707.     scanf("%lf",&(c->radius));
  708.  
  709.     return((void *)c);
  710. }
  711.  
  712.  
  713. static void *cone_input()
  714. {
  715.     Cone *c=(Cone *)malloc(sizeof(Cone));
  716.  
  717.     if (c==NULL)error("Cannot allocate Cone/Cup.");
  718.     scanf("%lf%lf%lf%lf",&(c->apex_point.x),&(c->apex_point.y),&(c->apex_point.z),
  719.             &(c->apex_radius));
  720.     scanf("%lf%lf%lf%lf",&(c->base_point.x),&(c->base_point.y),&(c->base_point.z),
  721.             &(c->base_radius));
  722.  
  723.     return((void *)c);
  724. }
  725.  
  726. static void *ring_input()
  727. {
  728.     Ring *r=(Ring *)malloc(sizeof(Ring));
  729.  
  730.     if (r==NULL)error("Cannot allocate Ring.");
  731.     scanf("%lf%lf%lf",&(r->center.x),&(r->center.y),&(r->center.z));
  732.     scanf("%lf%lf%lf",&(r->normal.x),&(r->normal.y),&(r->normal.z));
  733.     scanf("%lf%lf",&(r->inner_radius),&(r->outer_radius));
  734.  
  735.     return((void *)r);
  736. }
  737.  
  738. static void *no_input()
  739. {
  740.     error("This type Not implemented.");
  741. }
  742.  
  743. /************** End Object Specific Input        **************/
  744. /************** Begin Object Specific Normal     **************/
  745.  
  746. static Vector3 quadrilateral_normal(q,u,v)
  747. Quadrilateral *q;
  748. double u,v;
  749. {
  750.     return(q->plane_normal);
  751. }
  752.  
  753. static Vector3 polygon_normal(p,u,v)
  754. Polygon *p;
  755. double u,v;
  756. {
  757.     return(p->plane_normal);
  758. }
  759.  
  760. static Vector3 sphere_normal(s,u,v)
  761. Sphere *s;
  762. double u,v;
  763. {
  764.     Vector3 N;
  765.     double phi = u * 2. * PI;
  766.     double theta = v * PI;
  767.     double sin_theta = sin(theta);
  768.  
  769.     N.x = sin_theta*cos(phi);
  770.     N.y = sin_theta*sin(phi);
  771.     N.z = cos(theta);
  772.  
  773.     return(N);
  774. }
  775.  
  776. static Vector3 bubble_normal(s,u,v)
  777. Sphere *s;
  778. double u,v;
  779. {
  780.     Vector3 N;
  781.  
  782.     N=sphere_normal(s,u,v);
  783.     V3Negate(&N);
  784.  
  785.     return(N);
  786. }
  787.  
  788. static Vector3 cylinder_normal(c,u,v)
  789. Cylinder *c;
  790. double u,v;
  791. {
  792.     Vector3 N;
  793.     double phi = u * 2. * PI;
  794.  
  795.     N.x = cos(phi);
  796.     N.y = sin(phi);
  797.     N.z = 0.;
  798.     N = transform_vector(&N,&(c->object_to_world));
  799.  
  800.     return(N);
  801. }
  802.  
  803. static Vector3 tube_normal(c,u,v)
  804. Cylinder *c;
  805. double u,v;
  806. {
  807.         Vector3 N;
  808.  
  809.     N=cylinder_normal(c,u,v);
  810.     V3Negate(&N);
  811.  
  812.     return(N);
  813. }
  814.  
  815. static Vector3 cone_normal(c,u,v)
  816. Cone *c;
  817. double u,v;
  818. {
  819.         Vector3 N,v1,v2;
  820.     double phi = u * 2. * PI;
  821.  
  822.     v1.x = cos(phi); v1.y = sin(phi); v1.z = 1.;
  823.     v1 = transform_vector(&v1,&(c->object_to_world));
  824.     V3Cross(&v1,&(c->axis),&v2);
  825.     V3Cross(&v1,&v2,&N);
  826.     V3Normalize(&N);
  827.  
  828.         return(N);
  829. }
  830.  
  831. static Vector3 cup_normal(c,u,v)
  832. Cone *c;
  833. double u,v;
  834. {
  835.         Vector3 N;
  836.  
  837.     N=cone_normal(c,u,v);
  838.     V3Negate(&N);
  839.  
  840.         return(N);
  841. }
  842.  
  843. static Vector3 ring_normal(r,u,v)
  844. Ring *r;
  845. double u,v;
  846. {
  847.         return(r->normal);
  848. }
  849.  
  850. static Vector3 no_normal(s,u,v)
  851. void *s;
  852. double u,v;
  853. {
  854.     error("No Normal Yet.");
  855. }
  856.  
  857. /************** End Object Specific Normal       **************/
  858. /************** Begin Object Specific Local Matrix ************/
  859.  
  860. static Matrix4 quadrilateral_matrix(N,q)
  861. Vector3 *N;
  862. Quadrilateral *q;
  863. {
  864.     return(q->local_trans);
  865. }
  866.  
  867.  
  868. static Matrix4 polygon_matrix(N,p)
  869. Vector3 *N;
  870. Polygon *p;
  871. {
  872.     return(p->local_trans);
  873. }
  874.  
  875. static Matrix4 align_matrix(N,s)
  876. Vector3 *N;
  877. void *s;
  878. {
  879.     Matrix4 m;
  880.  
  881.     align_Z_to_axis(N,&m);
  882.  
  883.     return(m);
  884. }
  885.  
  886. static Matrix4 ring_matrix(N,r)
  887. /*
  888.     CAUTION: This matrix also has a nonzero translation vector.
  889. */ 
  890. Vector3 *N;
  891. Ring *r;
  892. {
  893.     return(r->object_to_world);
  894. }
  895.  
  896. static Matrix4 no_matrix(N,s)
  897. Vector3 *N;
  898. void *s;
  899. {
  900.     error("No Matrix Yet.");
  901. }
  902.  
  903. /************** End Object Specific Local Matrix   ************/
  904. /************** Begin Object Specific Area Sampling ***********/
  905.  
  906. static Point3 sphere_point(s,u,v)
  907. Sphere *s;
  908. double u,v;
  909. {
  910.     Point3 p;
  911.     double R =s->radius;
  912.     Point3 C;
  913.     C =s->center;
  914.     p.x = R*sin(PI*v)*cos(2.*PI*u);
  915.     p.y = R*sin(PI*v)*sin(2.*PI*u);
  916.     p.z = R*cos(PI*v);
  917.     V3Add(&p,&C,&p);
  918.     return(p);
  919. }
  920.  
  921. static Point3 cylinder_point(c,u,v)
  922. Cylinder *c;
  923. double u,v;
  924. {
  925.     Point3 p;
  926.     double R = c->radius;
  927.     double H = c->height;
  928.     double phi = 2.*PI*u;
  929.     p.x = R*cos(phi);
  930.     p.y = R*sin(phi);
  931.     p.z = v*H;
  932.     V3MulPointByMatrix(&p,&(c->object_to_world));
  933.     return(p);
  934. }
  935. static Point3 cone_point(c,u,v)
  936. Cone *c;
  937. double u,v;
  938. {
  939.     Point3 posn;
  940.     double phi=u * 2. * PI;
  941.     Point3 p;
  942.     p.z = (1.-c->distance)*v+c->distance;
  943.     p.x = p.z * cos(phi); p.y = p.z * sin(phi);
  944.     V3MulPointByMatrix(&p,&(c->object_to_world));
  945.     return(p);
  946. }
  947. static Point3 ring_point(r,u,v)
  948. Ring *r;
  949. double u,v;
  950. {
  951.     Point3 p;
  952.     double radius=r->inner_radius+
  953.         v*(r->outer_radius-r->inner_radius);
  954.     double phi=2. * u * PI;
  955.     p.x = radius * cos(phi); p.y = radius * sin(phi); p.z = 0.;
  956.     V3MulPointByMatrix(&p,&(r->object_to_world));
  957.     return(p);
  958. }
  959. static Point3 triangle_point(p,u,v)
  960. Point3 *p;
  961. double u,v;
  962. {
  963.     Point3 posn;
  964.     posn.x = p[0].x + (1.-v)*u*(p[1].x-p[0].x)+ v*(p[2].x-p[0].x);
  965.     posn.y = p[0].y + (1.-v)*u*(p[1].y-p[0].y)+ v*(p[2].y-p[0].y);
  966.     posn.z = p[0].z + (1.-v)*u*(p[1].z-p[0].z)+ v*(p[2].z-p[0].z);
  967.     return(posn);
  968. }
  969. static Point3 polygon_point(p,u,v)
  970. Polygon *p;
  971. double u,v;
  972. {
  973.     return(p->vertex_list[0]); /* Arbitrary */
  974. }
  975. static Point3 no_object_point(p,u,v)
  976. void *p;
  977. double u,v;
  978. {
  979.     error("Cannot compute point for this object.");
  980. }
  981. /************** End Object Specific Area Sampling   ***********/
  982.  
  983. /************** Begin Grid Centre Point             ***********/
  984. #define get_U ((uindex+0.5) * 1./object[n].grid_h_reso)
  985. #define get_V ((vindex+0.5) * 1./object[n].grid_v_reso)
  986. static Point3 quadrilateral_grid_centre_point(n,vindex,uindex)
  987. int n;
  988. int vindex,uindex;
  989. {
  990.     double u = get_U;
  991.     double v = get_V;
  992.     return(get_bilinear(object[n].object_specific_structure,u,v));
  993. }
  994. static Point3 sphere_grid_centre_point(n,vindex,uindex)
  995. int n;
  996. int vindex,uindex;
  997. {
  998.     double u = get_U;
  999.     double v = get_V;
  1000.     return(sphere_point(object[n].object_specific_structure,u,v));
  1001. }
  1002. static Point3 cylinder_grid_centre_point(n,vindex,uindex)
  1003. int n;
  1004. int vindex,uindex;
  1005. {
  1006.     double u = get_U;
  1007.     double v = get_V;
  1008.     return(cylinder_point(object[n].object_specific_structure,u,v));
  1009. }
  1010. static Point3 cone_grid_centre_point(n,vindex,uindex)
  1011. int n;
  1012. int vindex,uindex;
  1013. {
  1014.     double u = get_U;
  1015.     double v = get_V;
  1016.     return(cone_point(object[n].object_specific_structure,u,v));
  1017. }
  1018. static Point3 ring_grid_centre_point(n,vindex,uindex)
  1019. int n;
  1020. int vindex,uindex;
  1021. {
  1022.     double u = get_U;
  1023.     double v = get_V;
  1024.     return(ring_point(object[n].object_specific_structure,u,v));
  1025. }
  1026. static Point3 polygon_grid_centre_point(n,vindex,uindex)
  1027. int n;
  1028. int vindex,uindex;
  1029. {
  1030.     Polygon *p=(Polygon *)object[n].object_specific_structure;
  1031.     if (p->nvertices==3){/* Triangle */
  1032.         double u = get_U;
  1033.         double v = get_V;
  1034.         return(triangle_point(p->vertex_list,u,v));
  1035.     }
  1036.     else /* We donot have a parameterisation method for
  1037.         arbitrary polygon. So I can send any point 
  1038.         on the surface. */
  1039.         return(p->vertex_list[0]);
  1040. }
  1041. static Point3 no_grid_centre_point(n,vindex,uindex)
  1042. int n;
  1043. int vindex,uindex;
  1044. {
  1045.     error("Grid center Point Not Implemented.");
  1046. }
  1047. #undef get_U
  1048. #undef get_V
  1049. /************** End Grid Centre Point               ***********/
  1050.  
  1051. ObjectFunctions ofunc[MAX_GEOMETRY_TYPE]={
  1052. /*
  1053.     NOTE : Each of the Normal routines must return Normalized Vectors. 
  1054. */
  1055. {
  1056.     quadrilateral_input,
  1057.     quadrilateral_normal,
  1058.     intersect_quadrilateral,
  1059.     quadrilateral_preprocess,
  1060.     quadrilateral_matrix,
  1061.     quadrilateral_grid_centre_point,
  1062.     get_bilinear
  1063. },
  1064. {
  1065.     sphere_input,
  1066.     sphere_normal,
  1067.     intersect_sphere,
  1068.     sphere_preprocess,
  1069.     align_matrix,
  1070.     sphere_grid_centre_point,
  1071.     sphere_point
  1072. },
  1073. {
  1074.     sphere_input,
  1075.     bubble_normal,
  1076.     intersect_bubble,
  1077.     sphere_preprocess,
  1078.     align_matrix,
  1079.     sphere_grid_centre_point,
  1080.     sphere_point
  1081. },
  1082. {
  1083.     cylinder_input,
  1084.     cylinder_normal,
  1085.     intersect_cylinder,
  1086.     cylinder_preprocess,
  1087.     align_matrix,
  1088.     cylinder_grid_centre_point,
  1089.     cylinder_point
  1090. },
  1091. {
  1092.     cylinder_input,
  1093.     tube_normal,
  1094.     intersect_tube,
  1095.     cylinder_preprocess,
  1096.     align_matrix,
  1097.     cylinder_grid_centre_point,
  1098.     cylinder_point,
  1099. },
  1100. {
  1101.     cone_input,
  1102.     cone_normal,
  1103.     intersect_cone,
  1104.     cone_preprocess,
  1105.     align_matrix,
  1106.     cone_grid_centre_point,
  1107.     cone_point
  1108. },
  1109. {
  1110.     cone_input,
  1111.     cup_normal,
  1112.     intersect_cup,
  1113.     cone_preprocess,
  1114.     align_matrix,
  1115.     cone_grid_centre_point,
  1116.     cone_point
  1117. },
  1118. {
  1119.     ring_input,
  1120.     ring_normal,
  1121.     intersect_ring,
  1122.     ring_preprocess,
  1123.     ring_matrix,
  1124.     ring_grid_centre_point,
  1125.     ring_point
  1126. },
  1127. {
  1128.     polygon_input,
  1129.     polygon_normal,
  1130.     intersect_polygon,
  1131.     polygon_preprocess,
  1132.     polygon_matrix,
  1133.     polygon_grid_centre_point,
  1134.     polygon_point
  1135. },
  1136. {
  1137.     no_input,
  1138.     no_normal,
  1139.     intersect_other_geometry,
  1140.     null_preprocess,
  1141.     no_matrix,
  1142.     no_grid_centre_point,
  1143.     no_object_point
  1144. }
  1145. };
  1146. #undef get_V
  1147. #undef get_U
  1148. #undef update_scene_bound
  1149.