home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / radiance / simplerd.lha / simplerad / FinalFTP / Light / rad.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-21  |  32.4 KB  |  1,030 lines
  1. /**********************************************************************/
  2. /* rad.c                                                              */
  3. /*                                                                    */
  4. /* Radiosity solution step routines. Algorithms used:                 */
  5. /*                                                                    */
  6. /* General:                                                           */
  7. /*    Back face culling of source and receiver planes.                */
  8. /*    Receiver subdivision by: visibility, intensity gradients        */
  9. /* Chen89:                                                            */
  10. /*   progressive refinement                                           */
  11. /* Wallace89, etc:                                                    */
  12. /*    Hierarchical bounding volumes with object bv culling            */
  13. /*    Has ray casting for visiblity testing                           */
  14. /* Haines91: shaft culling                                            */
  15. /* (see documentation)                                                */
  16. /*                                                                    */
  17. /* Initial program shell taken from E. Chen in Graphics Gems II, 1991 */
  18. /*                                                                    */
  19. /* Copyright (C) 1992, Bernard Kwok                                   */
  20. /* All rights reserved.                                               */
  21. /* Revision 1.0                                                       */
  22. /* May, 1992                                                          */
  23. /**********************************************************************/
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <stdlib.h>
  27. #include <sys/time.h>
  28. #include <sys/resource.h>
  29. #include "geo.h"
  30. #include "struct.h"
  31. #include "rad.h"
  32. #include "display.h"
  33. #include "io.h"
  34. #include "misc.h"
  35. #include "ff.h"
  36. #include "bvol.h"
  37. #include "ray.h"
  38. #include "scull.h"
  39. #include "rtime.h"
  40.  
  41. #include "radmisc.h"
  42.  
  43. Polygon *shootPatch;          /* Current shooting patch               */
  44. Polygon *recPatch;            /* Current receiver patch               */
  45. Polygon *prev_Shooter = 0;    /* Previous shooting patch              */
  46. float avg_rec = 0.0;          /* Average percent receivers per 
  47.                  iteration */
  48. int num_src_samples;          /* Number of source samples             */
  49.  
  50. /**********************************************************************/
  51. /* Initialize list of possible receivers                              */
  52. /**********************************************************************/
  53. void Init_ReceiverList(p)
  54.      RadParams *p;
  55. {
  56.   int i;
  57.   Elist *el, *tmp;
  58.  
  59.   if (Option.debug) printf("\t> Initializing initial receiver list\n");
  60.   el = p->elements;
  61.   for (i=0;i<p->num_elements;i++) {
  62.     tmp = el->next;
  63.     free(el);
  64.     el = tmp;
  65.   }
  66.   p->num_elements = p->num_receivers = 0;
  67.   p->eltail = p->elements = 0;
  68. }
  69.  
  70. /**********************************************************************/
  71. /* Initialize structures for ray traced form factors */
  72. /**********************************************************************/
  73. void Init_RayTrace()
  74. {
  75.   if (Option.debug) printf("\t> Initializing for ray traced ff\n");
  76.   StartIntersect();
  77. }
  78.  
  79. /**********************************************************************/
  80. /* Find initial set of elements to display                            */
  81. /**********************************************************************/
  82. void Create_Initial_ReceiverList(theScene)
  83.      Scene *theScene;
  84. {
  85.   int i,j,k;
  86.   Scene *sptr = theScene;
  87.   Objectt *optr;
  88.   Mesh *mptr;
  89.   Polygon *pptr;
  90.   
  91.   Init_ReceiverList(&ReadLog);
  92.   for(i=0, optr=sptr->objects;i<sptr->num_objects;i++, optr++) 
  93.     for(j=0, mptr=optr->meshes;j<optr->num_meshes;j++, mptr++) 
  94.       for(k=0, pptr=mptr->polys;k<mptr->num_polys;k++, pptr++) 
  95.     Find_Receivers(pptr,(Polygon *)INITIAL_DISPLAY,1);
  96. }
  97.  
  98. /**********************************************************************/
  99. /* Initialize radioisity parameters, not set when reading scene       */
  100. /**********************************************************************/
  101. void Init_Rad(theScene)
  102.      Scene *theScene;
  103. {
  104.   if (Option.debug) printf("\n\t*** Initializing radiosity structures ***\n");
  105.   
  106.   /* Prepare the original scene */
  107.   Create_Initial_ReceiverList(theScene);
  108.  
  109.   /* Initialize for ray traced or hemicube form-factor calcs */
  110.   if (Option.ff_raytrace) {
  111.     Init_RayTrace();
  112.     if (FF_Options.shaft_cull) 
  113.       Init_ShaftStats();
  114.   } else 
  115.     InitHemicube();
  116.   
  117.   /* Prepare for display of intermediate PR steps if option set */
  118.   if (Option.show_pr_steps) {
  119.     if (Option.debug) {
  120.       print_View(stderr,&(ReadLog.displayView));
  121.       printf("\t> Displaying initial results\n");
  122.     }
  123.     Create_Display_Buffer(ReadLog.displayView.xRes,ReadLog.displayView.yRes);
  124.     if (Option.ambient)              /* Find initial ambient term */
  125.       ReadLog.ambient_term =  Pr_Ambient_Term();
  126.     DisplayResults(ReadLog.displayView);
  127.   } else {
  128.     if (Option.debug) 
  129.       print_View(stderr,&(ReadLog.displayView));
  130.   }
  131. }
  132.  
  133. int iteration = 1;
  134. /**********************************************************************/
  135. /* Print current iteration                                            */
  136. /**********************************************************************/
  137. int Print_Iter()
  138. {
  139.   int i;
  140.  
  141.   /* Print out intermediate results if specified */
  142.   if (Option.write_result) {
  143.     for(i=0;i<num_inter;i++) 
  144.       if ((iteration-1) == intermed[i]) {
  145.     if (test_optik)
  146.       Write_RadOptik(ReadLog, Option.OutSceneFilename, iteration);
  147.     else
  148.       Write_Rad(ReadLog, Option.OutSceneFilename, iteration);
  149.     i = num_inter+1;
  150.       }
  151.   }
  152.  
  153.   /* LogConv_Pt(iteration, ReadLog.totalEnergyLeft); */
  154.   LogConv_Pt(iteration, ReadLog.totalEnergyLeft / ReadLog.totalEnergy);
  155.   printf("\n\t*** Iteration %d. Energy Left = %g (%g %%) ***\n", 
  156.      iteration, ReadLog.totalEnergyLeft,
  157.      ReadLog.totalEnergyLeft / ReadLog.totalEnergy * 100.0);
  158.  
  159.   if (Option.device == FILES)
  160.     fprintf(Option.StatFile, 
  161.         "\n\t*** Iteration %d. Energy Left = %g (%g %%) ***\n", 
  162.         iteration, ReadLog.totalEnergyLeft,
  163.         ReadLog.totalEnergyLeft / ReadLog.totalEnergy * 100.0);
  164.   iteration++;
  165.   return iteration;
  166. }
  167.  
  168. /**********************************************************************/
  169. /* Form linkage between this patch and all others [Hanrahan91] (not   */
  170. /* implemented)                                                       */
  171. /**********************************************************************/
  172. void Form_Linkages(pptr)
  173.      Polygon *pptr;
  174. {
  175.   if (Option.debug) 
  176.     printf("\t> Linkages for %s%d\n", pptr->name, pptr->id);
  177. }
  178.  
  179. /**********************************************************************/
  180. /* Forming of linkages between patches using BF_Refinement [Hanra91]  */
  181. /* (not used)                                                         */
  182. /**********************************************************************/
  183. void Do_Linkage()
  184. {
  185.   int i;
  186.   Elist *Int_List;
  187.   int Int_ListSize;
  188.   Elist *elptr;
  189.   Polygon *pptr;
  190.  
  191.   if (Option.debug) {
  192.     printf("\n\t*** Forming patch linkages ***\n");
  193.   }
  194.   /* Setup list of patches to compare against others */
  195.   Create_Initial_ReceiverList(&RadScene);
  196.   Int_List = ReadLog.elements; /* This is wrong ! Need to fix ! */
  197.   Int_ListSize = ReadLog.num_elements;
  198.  
  199.   elptr = Int_List;
  200.   for(i=0;i<Int_ListSize;i++) {
  201.     pptr = elptr->element;
  202.     Init_ObjectList();
  203.     Make_ObjectList(pptr, &Scene_BVH, 1);
  204.     Make_ReceiverList(pptr);
  205.     Form_Linkages(pptr);
  206.   }
  207. }
  208.  
  209. /**********************************************************************/
  210. /* Cleanup stuff, and log statistics                                  */
  211. /**********************************************************************/
  212. void Do_Cleanup()
  213. {
  214.   /* Collect statistics for rays and shaft */
  215.   if (Option.tablelog) {
  216.     fprintf(rlogfile,"%d \\\\\n%g \\\\\n%g \\\\\n", 
  217.         iteration-1, ReadLog.totalEnergyLeft,
  218.         ReadLog.totalEnergyLeft / ReadLog.totalEnergy * 100.0);
  219.   }
  220.   /* Log CPU times */
  221.   Print_Times(Option.StatFile, FF_Options.shaft_cull);
  222.  
  223.   if (Option.device == PRINT) 
  224.     EndIntersect(stdout);
  225.   else 
  226.     EndIntersect(Option.StatFile);
  227.   if (FF_Options.shaft_cull) 
  228.     EndShaft(TRUE, Option.StatFile);
  229.  
  230.   /* Log percent receivers per iteration */
  231.   if (Option.tablelog) 
  232.     fprintf(rlogfile,"AR:%g \\\\\n", avg_rec);
  233.  
  234.   /* Log convergence results */
  235.   LogConv(1,"");
  236.  
  237.   /* Save final image, if any (only on IRIS) */
  238.   if (Option.show_pr_steps) {
  239.     if (Option.write_result) {
  240.       /* getsize(&xsize, &ysize);
  241.      Save_Screen_Image(xsize, ysize, 3, "rad.img"); */
  242.     }
  243.     /* Cleanup display buffers */
  244.     CleanUpBuffers();                    
  245.   }
  246. }
  247.  
  248. /**********************************************************************/
  249. /* Radiosity main algorithm                                           */
  250. /**********************************************************************/
  251. void Do_Rad()
  252. {
  253.   int xsize, ysize;
  254.   float iter_start, iter_end;
  255.  
  256.   if (Option.debug) printf("\n\t*** Doing radiosity computations ***\n");
  257.   
  258.   LogConv(1,Option.meshfilename); /* Start logging convergence data */
  259.  
  260.   /* Start recording iteration */
  261.   iter_start = Cpu_Time(&tstats.utime, &tstats.stime);  
  262.  
  263.   while(((shootPatch = FindShootPatch()) != 0) &&
  264.     (Print_Iter() <= ReadLog.max_iterations+1)) {
  265.     prev_Shooter = shootPatch;
  266.     Init_ObjectList();
  267.     if (Option.grid) {
  268.       Make_ObjectList(shootPatch, &Scene_BVH, 1);
  269.       if (Option.debug)
  270.     Print_ObjectList();
  271.     }
  272.     Make_ReceiverList(shootPatch);
  273.     ComputeFormfactors(shootPatch);
  274.     Distribute_Rad(shootPatch);
  275.     
  276.     if (ReadLog.num_receivers > 0) { /* Need at least 1 receiver */
  277.       if (Option.show_pr_steps) {
  278.     if (Option.ambient) {
  279.       ReadLog.ambient_term =  Pr_Ambient_Term();
  280.       print_Spectra(stdout, "Ambient term", ReadLog.ambient_term);
  281.     }
  282.     DisplayResults(ReadLog.displayView);
  283.       }
  284.     } /* At least 1 receiver */
  285.  
  286.     Subdivide_Elements();
  287.     
  288.     if (Option.ff_raytrace) {
  289.       if (Option.device == PRINT)
  290.     EndIntersect(stdout);
  291.     } else 
  292.       FreeSumFactors();
  293.   }
  294.  
  295.   /* Stop recording iteration */
  296.   iter_end = Cpu_Time(&tstats.utime, &tstats.stime); 
  297.   tstats.tot_time += (iter_end - iter_start);
  298.   tstats.avg_iter += (iter_end - iter_start) / (float) (iteration-1);
  299.   tstats.avg_ff = tstats.avg_ff / (float) (iteration-1);
  300.   tstats.avg_shaft /= (float) (iteration-1);
  301.   tstats.avg_ray /= (float) (iteration-1);
  302.   avg_rec /= (float) (iteration-1);
  303.   
  304.   if (Option.UseVCR) { /* Stop the video */
  305.     vcr_end_video();
  306.   }
  307.  
  308.   /* Cleanup program */
  309.   Do_Cleanup();
  310. }
  311.  
  312. /**********************************************************************/
  313. /* Check if at bottom of poly tri-quad tree, i.e. poly is leaf node   */
  314. /**********************************************************************/
  315. int IsLeaf(pptr, children)
  316.      Polygon *pptr;
  317.      int children[MAX_PATCH_CHILDREN];
  318. {
  319.   int i;
  320.   int no_children = 1;          /* Assume no children */
  321.   
  322.   for(i=0;i<MAX_PATCH_CHILDREN;i++)
  323.     if (pptr->child[i] != 0) {
  324.       no_children = 0;
  325.       children[i] = 1;
  326.     } else {
  327.       children[i] = 0;
  328.     }
  329.   return(no_children);
  330. }
  331.  
  332. /**********************************************************************/
  333. /* Add polygon to element list as possible receiver */
  334. /**********************************************************************/
  335. void AddTo_ReceiverList(pptr,patches_seen)
  336.      Polygon *pptr;
  337.      int patches_seen;
  338. {
  339.   Elist *elptr;
  340.   
  341.   /* printf("\t> Adding reciever[%d] %s%d = %d\n", ReadLog.num_receivers,
  342.      pptr->name, pptr->id, patches_seen); */
  343.   elptr = (Elist *) malloc(sizeof(Elist));
  344.   elptr->element = pptr;
  345.   elptr->next = 0;
  346.   elptr->ff = 0.0;                    /* Reset form-factor */
  347.   elptr->is_subdivided = 0;           /* Reset subdivision flag */
  348.   elptr->is_receiver = patches_seen;  /* Mark if element is a receiver */
  349.   
  350.   if (ReadLog.num_elements == 0) {    /* First possible receiver */
  351.     ReadLog.elements = elptr;
  352.     ReadLog.eltail = elptr;
  353.   } else {
  354.     ReadLog.eltail->next = elptr;     
  355.     ReadLog.eltail = elptr;
  356.   }
  357.   ReadLog.num_elements++;             /* Count total # of elements */
  358.   if (patches_seen)
  359.     ReadLog.num_receivers++;          /* Count # of receivers */
  360. }
  361.  
  362. /**********************************************************************/
  363. /* Check if any of sources unshot energy will be reflected by the     */
  364. /* receiver                                                           */
  365. /**********************************************************************/
  366. int Receiver_ReflectsEnergy(rec, src)
  367.      Polygon *rec, *src;
  368. {
  369.   int i;
  370.   Spectra rec_p;
  371.   Spectra src_uB;
  372.   double pB_sum = 0.0;
  373.   double p_sum = 0.0;
  374.  
  375.   rec_p = poly_reflect(rec); 
  376.   src_uB = src->unshot_B;
  377.   pB_sum = 0.0;
  378.   for (i=0;i<MAX_SPECTRA_SAMPLES;i++) {
  379.     pB_sum += rec_p.samples[i] * src_uB.samples[i];
  380.     p_sum += rec_p.samples[i];
  381.   }
  382.  
  383.   /* If reflects no unshot or has no reflectance, then it is 
  384.      not a receiver. */
  385.   if ((pB_sum == 0.0) || (p_sum == 0.0))
  386.     return 0;
  387.   else
  388.     return 1;
  389. }
  390.  
  391. /**********************************************************************/
  392. /* Find all elements for a poly = all leaves in tri/quad polygon tree */
  393. /**********************************************************************/
  394. void Find_Receivers(pptr,src,patches_seen)
  395.      Polygon *pptr;
  396.      Polygon *src;
  397.      int patches_seen;
  398. {
  399.   int i;
  400.   int children[MAX_PATCH_CHILDREN];
  401.   BoundingBoxType pbox;
  402.  
  403.   if (IsLeaf(pptr,children)) {
  404.  
  405.     if (src == INITIAL_DISPLAY) {    /* Display all, flag is meaningless */
  406.       AddTo_ReceiverList(pptr,1);
  407.     }
  408.  
  409.     /* Do checks on geometry to see if patch is a potential receiver 
  410.        for source (only if still a possible receiver) */
  411.     else if (patches_seen == 0) { 
  412.       AddTo_ReceiverList(pptr,0);
  413.     }
  414.      
  415.     else {
  416.       pbox = BoxPoly(pptr);
  417.       
  418.       if (src == pptr) {                /* Can't see yourself ! */
  419.     AddTo_ReceiverList(pptr,0); 
  420.       }
  421.       
  422.       else if ((pptr->Pfather != 0) &&
  423.            (src->Pfather == pptr->Pfather || pptr->Pfather == src)) {
  424.     AddTo_ReceiverList(pptr,0);  /* Can't see your own patch */
  425.        } 
  426.       
  427.       /* If receiver is behind, or on the plane of the source then 
  428.      the source cannot "see" it */
  429.       else if (Bounds_Behind_Plane(&pbox, src->normal[0], src->d)) {
  430.     AddTo_ReceiverList(pptr,0); 
  431.       }
  432.       
  433.       else 
  434.     AddTo_ReceiverList(pptr,patches_seen); 
  435.     }
  436.  
  437.   } else {
  438.     for(i=0;i<MAX_PATCH_CHILDREN;i++) {
  439.       if (children[i]) Find_Receivers(pptr->child[i],src,
  440.                       patches_seen);
  441.     }
  442.   }
  443. }
  444.  
  445. /**********************************************************************/
  446. /* Check if any of sources unshot energy will be reflected by the     */
  447. /* receiver                                                           */
  448. /**********************************************************************/
  449. int Object_ReflectsEnergy(rec_p, src)
  450.      Spectra rec_p;     
  451.      Polygon *src;
  452. {
  453.   int i;
  454.   Spectra src_uB;
  455.   double pB_sum = 0.0;
  456.  
  457.   src_uB = src->unshot_B;
  458.   pB_sum = 0.0;
  459.   for (i=0;i<MAX_SPECTRA_SAMPLES;i++)
  460.     pB_sum += rec_p.samples[i] * src_uB.samples[i];
  461.  
  462.   /* If reflects no unshot or has no reflectance, then is 
  463.      not a receiver ! */
  464.   if (pB_sum == 0.0)
  465.     return 0;
  466.   else
  467.     return 1;
  468. }
  469.  
  470. /**********************************************************************/
  471. /* Find potential receivers. Use some culling to speed up stuff       */
  472. /**********************************************************************/
  473. void Make_ReceiverList(shootPatch)
  474.      Polygon *shootPatch;
  475. {
  476.   int i,j,k;
  477.   int patches_seen;
  478.   Scene *sptr = &RadScene;
  479.   Objectt *optr;
  480.   Object_List *olptr;
  481.   Mesh *mptr;
  482.   Polygon *pptr;
  483.   Spectra rec_p;
  484.  
  485.   if (Option.debug) printf("\n\t> Finding receivers\n");
  486.  
  487.   /* Empty list of potential receivers, and then scan all 
  488.      polygons for potential receivers */
  489.  
  490.   Init_ReceiverList(&ReadLog);  
  491.  
  492.   /* Use bounding volumes to try and reduce number of receivers */
  493.   if (Option.grid) {                    
  494.     for(i=0, olptr=objlist;i<objlist_size;i++, olptr = olptr->next) {
  495.       optr = olptr->hbox->object;
  496.       
  497.       /* If object's bounding volume can be seen, and object will 
  498.      reflect some energy then assume polygon is a receiver */
  499.       if (olptr->is_receiver) {
  500.     rec_p = obj_reflect(optr);
  501.     olptr->is_receiver = Object_ReflectsEnergy(rec_p, shootPatch);
  502.       }
  503.       patches_seen = olptr->is_receiver;      
  504.       for(j=0, mptr=optr->meshes; j<optr->num_meshes;j++,mptr++) 
  505.     for(k=0, pptr=mptr->polys;k<mptr->num_polys;k++, pptr++) {
  506.       
  507.       /* Find potential receivers = elements of patches */
  508.       Find_Receivers(pptr,shootPatch,patches_seen);
  509.     }
  510.     }
  511.   }
  512.   
  513.   /* Don't use bounding volumes. Must scan all objects */
  514.   else {
  515.     for(i=0, optr=sptr->objects;i<sptr->num_objects;i++, optr++) { 
  516.  
  517.       /* if object can reflect some energy then assume polygon 
  518.      is a receiver */
  519.       rec_p = obj_reflect(optr);
  520.       patches_seen  = Object_ReflectsEnergy(rec_p, shootPatch);
  521.  
  522.       for(j=0, mptr=optr->meshes;j<optr->num_meshes;j++, mptr++)
  523.     for(k=0, pptr=mptr->polys;k<mptr->num_polys;k++, pptr++) {
  524.  
  525.       /* Find potential receivers = elements of patches */
  526.       Find_Receivers(pptr,shootPatch,patches_seen);
  527.     }
  528.     }
  529.   }
  530.   printf("\tNumber of receivers / total = %d / %d\n",
  531.      ReadLog.num_receivers, ReadLog.num_elements);
  532.   if (Option.statistics) {
  533.     if (Option.device == FILES)
  534.       fprintf(Option.StatFile,"\tNumber of receivers / total = %d / %d\n",
  535.           ReadLog.num_receivers, ReadLog.num_elements);
  536.     avg_rec += (float) ReadLog.num_receivers /  (float)ReadLog.num_elements;
  537.   }
  538. }
  539.  
  540. /**********************************************************************/
  541. /* Print current shooting patch                                       */
  542. /**********************************************************************/
  543. void Print_ShootPatch(pptr, fptr, energySum)
  544.      Polygon *pptr;
  545.      FILE *fptr;
  546.      double energySum;
  547. {
  548.   FILE *fp;
  549.   
  550.   if (Option.device == PRINT)
  551.     fp = stdout;
  552.   else fp = fptr;
  553.   fprintf(fp,"\tShooter %s%d unshot: %g %g %g -> %g\n",
  554.       pptr->name, pptr->id, pptr->unshot_B.samples[0],
  555.       pptr->unshot_B.samples[1], pptr->unshot_B.samples[2], energySum);
  556. }
  557.  
  558. /**********************************************************************/
  559. /* Find next shooting patch based on unshot energy for each patch */
  560. /* Return 0 if convergence is reached, otherwise return 1 */
  561. /**********************************************************************/
  562. Polygon *FindShootPatch()
  563. {
  564.   Polygon *shootPatch;
  565.   int i,j,k,l;
  566.   double energySum, error, maxEnergySum;
  567.   Scene *sptr = &RadScene;
  568.   Objectt *optr;
  569.   Mesh *mptr;
  570.   Polygon *pptr;
  571.  
  572.   if (Option.debug) printf("\n\t> Finding shooting patch\n");
  573.  
  574.   maxEnergySum = 0.0;
  575.   shootPatch = 0;
  576.   
  577.   /* Scan all polygons for energy levels, and potential receivers */
  578.   ReadLog.totalEnergyLeft = 0.0;
  579.   for(i=0, optr=sptr->objects;i<sptr->num_objects;i++, optr++) 
  580.     for(j=0, mptr=optr->meshes;j<optr->num_meshes;j++, mptr++) 
  581.       for(k=0, pptr=mptr->polys;k<mptr->num_polys;k++, pptr= pptr->next) {
  582.     
  583.     /* Find patch with maximum energy to shoot 
  584.        and compute total energy left in scene */
  585.     energySum = 0.0;
  586.     for (l=0;l<MAX_SPECTRA_SAMPLES; l++) {
  587.       energySum = energySum + (pptr->unshot_B.samples[l] * pptr->area);
  588.     }
  589.     ReadLog.totalEnergyLeft += energySum;
  590.     
  591.     if (energySum > maxEnergySum) {
  592.       if (pptr != prev_Shooter) {
  593.         shootPatch = pptr;
  594.         maxEnergySum = energySum;
  595.       }
  596.     }
  597.       }
  598.  
  599.   if (shootPatch == 0) /* No shooting patch found */
  600.     return(0);
  601.   
  602.   /* Check for convergence */
  603.   error = maxEnergySum / ReadLog.totalEnergy;
  604.   if (error < ReadLog.threshold) 
  605.     return(0);
  606.   else {
  607.     Print_ShootPatch(shootPatch, Option.StatFile, maxEnergySum); 
  608.     printf("\t\tShooting patch has energy %g\n", maxEnergySum);
  609.     if (Option.debug) { 
  610.       fprintf(stdout,"\t> Shooting Patch is:\n");
  611.       print_Polygon(stdout,"",shootPatch);
  612.     }
  613.     return(shootPatch);
  614.   }
  615. }
  616.  
  617. /**********************************************************************/
  618. /* Find approximate centre of a polygon. Average of vertices used.    */
  619. /**********************************************************************/
  620. Vector Poly_Center(pptr)
  621.      Polygon *pptr;
  622. {
  623.   int i;
  624.   Vector center;
  625.   double ratio;
  626.  
  627.   center.x = center.y = center.z = 0.0;
  628.   ratio = 1.0 / (double) pptr->numVert;
  629.   for(i=0;i<pptr->numVert;i++) {
  630.     center.x += pptr->vert[i]->pos.x;
  631.     center.y += pptr->vert[i]->pos.y;
  632.     center.z += pptr->vert[i]->pos.z;
  633.   }
  634.   center.x = center.x * ratio;
  635.   center.y = center.y * ratio;
  636.   center.z = center.z * ratio;
  637.   return(center);
  638. }
  639.  
  640. /**********************************************************************/
  641. /* Check for receiver subdivision due to:                             */
  642. /*  1) Form factor > unity.                                           */
  643. /*  2) Form factor gradient > gradient threshold.                     */
  644. /* Return whether did subdivision or not.                             */
  645. /**********************************************************************/
  646. int Receiver_Subdivision(ff, eptr, vtx_ff)
  647.      double *ff;
  648.      Polygon *eptr;
  649.      double vtx_ff[MAX_PATCH_VTX];
  650. {
  651.   int i,ii;
  652.   int do_subdiv = 0;
  653.   double ff_diff;
  654.   int level_area = 0;
  655.  
  656.   if (Option.debug) printf("\tDoing subdivision tests...\n");
  657.  
  658.   /* TEST 1 : Is ff > one */
  659.   if (*ff > 1.0) { 
  660.     if (Option.debug)
  661.       printf("\t++ Form factor %g > 1. Subdividing poly\n", *ff);
  662.     *ff = 1.0;
  663.     do_subdiv = 1;                          /* Do subdivision        */
  664.   }
  665.  
  666.   if (eptr->level >= FF_Options.max_levels) /* Only to max level     */
  667.     return 0;
  668.   if (eptr->area <= FF_Options.min_element_area) /* Only to min area */
  669.     return 0;
  670.  
  671.   /* Only for ray cast ff computations. */
  672.   if ((Option.ff_raytrace) && (do_subdiv==0) && 
  673.       (eptr->level < FF_Options.max_levels)) {
  674.     i = 0;
  675.     while ((do_subdiv==0) && (i<eptr->numVert)) {
  676.       ii = (i + 1) % eptr->numVert;
  677.  
  678.       /* TEST 2 : ff gradient is too large between vertices of
  679.      edges on the polygon, stop if too small */
  680.       ff_diff = fabs(vtx_ff[i] - vtx_ff[ii]);
  681.       if (ff_diff > FF_Options.F_diff_edge) {
  682.     if (Option.debug)
  683.       printf("\t++ F diff = %g - %g = %g. Subdividing poly\n",
  684.          vtx_ff[i], vtx_ff[ii], ff_diff);
  685.     do_subdiv = 1;               /* Do subdivision */
  686.       }
  687.       
  688.       else {
  689.     i++;
  690.       }
  691.     }
  692.   }
  693.   return (do_subdiv);
  694. }
  695.  
  696. /**********************************************************************/
  697. /* Use hemicube to calculate form-factor from shooting patch to       */
  698. /* receiver elements.                                                 */
  699. /* (Modified from Eric Chen (Graphics Gems II), 1991)                 */
  700. /**********************************************************************/
  701. void ff_Hemicube(shootPatch)
  702.      Polygon *shootPatch;
  703. {
  704.   long i;
  705.   Vector up[5];
  706.   Vector lookat[5];
  707.   Vector center;
  708.   Vector normal, tangentU, tangentV, vec;
  709.   int face;
  710.   double norm_length;
  711.   Polygon *sp;                    /* Shooting patch */
  712.   Elist *elptr;                    /* Receiving elements */
  713.   double *fp;
  714.   long receiver_num = 0;
  715.   
  716.   /* Get center of shootPatch */
  717.   sp = shootPatch;
  718.   center = Poly_Center(sp);
  719.   if (Option.debug)
  720.     printf("\t ** Poly center is %g,%g,%g\n",center.x, center.y, center.z);
  721.   normal = sp->normal[0];
  722.  
  723.   /* rotate hemicube along normal axis of patch randomly */
  724.   /* will reduce hemicube aliasing artifacts */
  725.   do {
  726.     vec.x = RandomFloat;
  727.     vec.y = RandomFloat;
  728.     vec.z = RandomFloat;
  729.     tangentU = cross(&normal, &vec);
  730.     norm_length = norm(&tangentU);
  731.   } while(norm_length == 0);
  732.  
  733.   /* compute tangentV */
  734.   tangentV = cross(&normal, &tangentU);
  735.  
  736.   /* assign lookats, and ups for each hemicube face */
  737.   lookat[0] = *vadd(¢er, &normal);
  738.   up[0] = tangentU;
  739.   lookat[1] = *vadd(¢er, &tangentU);
  740.   up[1] = normal;
  741.   lookat[2] = *vadd(¢er, &tangentV);
  742.   up[2] = normal;
  743.   lookat[3] = *vsub(¢er, &tangentU);
  744.   up[3] = normal;
  745.   lookat[4] = *vsub(¢er, &tangentV);
  746.   up[4] = normal;
  747.   
  748.   /* place cube at center of shooting patch */
  749.   hemicube.view.lookfrom = center;
  750.   
  751.   /* Allocate temporary ff space of size of (number of receivers in scene) */
  752.   /* Note: this # changes with subdivision */
  753.   InitSumFactors(ReadLog.num_receivers);
  754.  
  755.   for (face=0;face<5;face++) {
  756.     hemicube.view.lookat = lookat[face];
  757.     hemicube.view.lookup = up[face];
  758.  
  759.     /* Only draw receiver elements */
  760.     receiver_num = 0;
  761.     if (Option.show_pr_steps)
  762.       Begin_DrawHC(hemicube.view, kBackgroundItem, face); 
  763.     for (i=0, elptr=ReadLog.elements;i<ReadLog.num_elements;i++) {
  764.       if (elptr->is_receiver) { 
  765.     Draw_Element(elptr->element, (unsigned long) receiver_num++); 
  766.       } 
  767.       elptr = elptr->next;
  768.     }
  769.     if (Option.show_pr_steps)
  770.       End_DrawHC();
  771.  
  772.     /* if (Option.debug) printf("+ summing delta ff\n"); */
  773.     if (face==0)
  774.       SumFactors(hemicube.view.xRes, hemicube.view.yRes,
  775.          hemicube.view.buffer, hemicube.topFactors);
  776.     else
  777.       SumFactors(hemicube.view.xRes, hemicube.view.yRes/2,
  778.          hemicube.view.buffer, hemicube.topFactors);
  779.   }
  780.  
  781.   /* Compute reciprocal form-factors, save in receiver list. */
  782.   elptr = ReadLog.elements;
  783.   fp = formfactors;
  784.  
  785.   for (i=0;i<ReadLog.num_elements; i++) {
  786.     if (elptr->is_receiver) {
  787.       elptr->ff = *fp * sp->area / elptr->element->area;
  788.       
  789.       /* Check for receiver subdivision */
  790.       elptr->is_subdivided = 
  791.     Receiver_Subdivision(&elptr->ff, elptr->element, (double *)0);
  792.       fp++;
  793.     }
  794.     elptr = elptr->next;
  795.   }
  796. }      
  797.  
  798. /**********************************************************************/
  799. /* Use ray tracing to compute form factors */
  800. /**********************************************************************/
  801. void ff_RayTrace(shootPatch)
  802.      Polygon*shootPatch;
  803. {
  804.   int i;
  805.   Elist *elptr;                       /* Receiving elements          */
  806.   Polygon *recPatch;
  807.   BoundingBoxType sbox, rbox;         /* Source / Receiver bboxes */
  808.   Sampletype src_samples[MAX_SAMPLES]; /* Samples of source */
  809.   Vector src_corners[MAX_PATCH_VTX];
  810.   int patches_close = 0;
  811.   int tmp;
  812.   float shaft_start, shaft_end;
  813.  
  814.   /* Find corners of source area, and generate sample points on 
  815.      source patch */
  816.   for(i=0;i<shootPatch->numVert;i++) 
  817.     src_corners[i] = shootPatch->vert[i]->pos;
  818.   num_src_samples = 
  819.     ff_Generate_Patch_Samples(shootPatch,src_samples, 
  820.                   src_corners, shootPatch->area, 
  821.                   shootPatch->normal[0]);
  822.   
  823.   /* Create source bbox if shaft culling */
  824.   if (FF_Options.shaft_cull) {
  825.     StartShaft();
  826.     sbox = BoxPoly(shootPatch);
  827.   }
  828.  
  829.   /* Calculate element to shooting patch form factors for 
  830.      all potential receivers */
  831.   elptr = ReadLog.elements;
  832.   for (i=0;i<ReadLog.num_elements;i++) {
  833.     if (elptr->is_receiver) {
  834.       recPatch = elptr->element;
  835.  
  836.       /* Shaft cull bounding volume hierarchy of scene */
  837.       if (FF_Options.shaft_cull) {
  838.     Init_ObjectList();
  839.     rbox = BoxPoly(recPatch);
  840.  
  841.     shaft_start = Cpu_Time(&tstats.utime, &tstats.stime);  
  842.     ShaftCull(shootPatch, recPatch, sbox, rbox, &Scene_BVH);
  843.     shaft_end = Cpu_Time(&tstats.utime, &tstats.stime);  
  844.     tstats.avg_shaft += (shaft_end - shaft_start);
  845.  
  846.     if (Option.debug) { 
  847.       Print_Shaft(shootPatch, recPatch, stdout);
  848.       Print_CandidateList();
  849.     }
  850.       }
  851.       
  852.       /* Only do source, receiver plane culling */
  853.       else if (FF_Options.src_rec_cull) 
  854.     Cull_ObjectList(recPatch);
  855.  
  856.       if (FF_Options.use_analytic)
  857.     if (patches_close = Patches_Abut(shootPatch, recPatch)) {
  858.       tmp = FF_Options.fftype; FF_Options.fftype = ANALYTIC_FF;
  859.     }
  860.       
  861.       /* Use equal weight per vtx ff to get patch ff */
  862.       elptr->ff = ff_patches(shootPatch, recPatch, 
  863.                  elptr->vtx_ff, num_src_samples,
  864.                  src_samples, src_corners) /
  865.                    (double) elptr->element->numVert;
  866.       if (FF_Options.use_analytic)
  867.     if (patches_close) FF_Options.fftype = tmp;
  868.       
  869.       /* Check for receiver subdivision */
  870.       if (elptr->ff < 0.0) elptr->ff = 0.0;
  871.       elptr->is_subdivided = 
  872.     Receiver_Subdivision(&elptr->ff, recPatch, elptr->vtx_ff);
  873.       if (elptr->ff > 1.0) elptr->ff = 1.0;
  874.     }
  875.     elptr = elptr->next;
  876.   }
  877.  
  878.   if (FF_Options.shaft_cull) EndShaft(FALSE, Option.StatFile);
  879. }
  880.  
  881. /**********************************************************************/
  882. /* Use ray tracing to compute form factors */
  883. /**********************************************************************/
  884. void Subdivide_Elements()
  885. {
  886.   int i;
  887.   Elist *elptr;                       /* Receiving elements          */
  888.  
  889.   elptr = ReadLog.elements;
  890.   for (i=0;i<ReadLog.num_elements;i++, elptr=elptr->next) {
  891.     if (elptr->is_receiver)
  892.       if (elptr->is_subdivided)
  893.     Subdivide_Polygon(elptr->element);
  894.   }
  895. }
  896.  
  897. /**********************************************************************/
  898. /* Calc form-factors from shooting patch to every element */
  899. /**********************************************************************/
  900. void ComputeFormfactors(shootPatch)
  901.      Polygon *shootPatch;
  902. {
  903.   if (Option.debug) printf("\n\t> Computing form factors\n");
  904.  
  905.   /* Compute form-factors */
  906.   if (Option.ff_raytrace) 
  907.     ff_RayTrace(shootPatch);     /* Element to shooting patch factors */
  908.   else 
  909.     ff_Hemicube(shootPatch);     /* Shooting patch to element factors */
  910. }
  911.  
  912. /**********************************************************************/
  913. /* Distribute energy to elements                                      */
  914. /**********************************************************************/
  915. void Distribute_Rad(shootPatch)
  916.      Polygon *shootPatch;
  917. {
  918.   long i,j,k;
  919.   Polygon *ep;                          /* Receiver to update */
  920.   Elist *elptr;                         /* List of receivers */
  921.   Spectra delta_B;                      /* Change in patch radiosity */
  922.   Spectra vtx_delta_B[MAX_PATCH_VTX];   /* Change in vertex radiosity */
  923.   double weight;                        /* Weight of element to patch */
  924.   Spectra element_p;                    /* Elements reflectance */
  925.   double element_ff;                    /* Elements form factor */
  926.  
  927.   if (Option.debug) printf("\n\t> Distributing energy\n");
  928.     
  929.   /* Distribute unshot radiosity to every element */
  930.   elptr = ReadLog.elements;
  931.   for (i=0;i<ReadLog.num_elements;i++, elptr=elptr->next) {
  932.     if (elptr->is_receiver) {
  933.       ep = elptr->element;
  934.       element_ff = elptr->ff;
  935.       
  936.       if (element_ff != 0.0) {
  937.     element_p = poly_reflect(ep);
  938.  
  939.     /* Store vertex radiosity changes, and compute patch changes */
  940.     if (Option.ff_raytrace) {
  941.       weight = 1.0; /* (1.0 / ep->numVert); */
  942.       for(k=0;k<ep->numVert;k++) {
  943.         for (j=0;j<MAX_SPECTRA_SAMPLES;j++) {
  944.           delta_B.samples[j] = 0.0;
  945.           vtx_delta_B[k].samples[j] = shootPatch->unshot_B.samples[j] * 
  946.         elptr->vtx_ff[k] *  element_p.samples[j];
  947.         }
  948.         for (j=0;j<MAX_SPECTRA_SAMPLES;j++) {
  949.           delta_B.samples[j] += (weight * vtx_delta_B[k].samples[j]);
  950.         }
  951.       }
  952.     }
  953.     
  954.     /* Used hemicube: store patch changes only */
  955.     else {
  956.       for (j=0;j<MAX_SPECTRA_SAMPLES;j++) {
  957.         delta_B.samples[j] = shootPatch->unshot_B.samples[j] * 
  958.           element_ff * element_p.samples[j];
  959.         for(k=0;k<ep->numVert;k++) 
  960.           vtx_delta_B[k].samples[j] = shootPatch->unshot_B.samples[j] * 
  961.         elptr->vtx_ff[k] * element_p.samples[j];
  962.       }
  963.     }
  964.       
  965.     /* Update element's radiosity and patch's unshot radiosity */
  966.     if (ep->Pfather != 0)
  967.       weight = ep->area / ep->Pfather->area;
  968.     else
  969.       weight = 1.0;
  970.     for (j=0;j<MAX_SPECTRA_SAMPLES; j++) {
  971.       ep->B.samples[j] += delta_B.samples[j];
  972.  
  973.       for(k=0;k<ep->numVert;k++) {
  974.         ep->vtx_B[k].samples[j] += vtx_delta_B[k].samples[j];
  975.       }
  976.  
  977.       if (ep->Pfather != 0) {
  978.         ep->Pfather->unshot_B.samples[j] += 
  979.           (delta_B.samples[j] * weight);
  980.       } else {
  981.         ep->unshot_B.samples[j] += delta_B.samples[j];
  982.       }
  983.     }
  984.  
  985.     /* Some debugging statistics */
  986.     if (Option.debug) { 
  987.       printf("\t+ Rec %s%d: p = %g,%g,%g; ff = %g\n", ep->name, ep->id,
  988.          element_p.samples[0], element_p.samples[1],  
  989.          element_p.samples[2], element_ff);
  990.  
  991.       if (Option.ff_raytrace == 0) {
  992.         printf("\t++ delta_B = %g,%g,%g; Shoot unshot = %g,%g,%g\n", 
  993.            delta_B.samples[0], delta_B.samples[1], delta_B.samples[2],
  994.            shootPatch->unshot_B.samples[0], 
  995.            shootPatch->unshot_B.samples[1],
  996.            shootPatch->unshot_B.samples[2]);
  997.       } else {
  998.         for(k=0;k<ep->numVert;k++) {
  999.           printf("\t++ v[%d]: delta_B = %g,%g,%g, B = %g,%g,%g\n", 
  1000.              k, vtx_delta_B[k].samples[0], vtx_delta_B[k].samples[1], 
  1001.              vtx_delta_B[k].samples[2], ep->vtx_B[k].samples[0], 
  1002.              ep->vtx_B[k].samples[1], ep->vtx_B[k].samples[2]);
  1003.         }
  1004.       }
  1005.     }
  1006.     
  1007.       } /* ff != 0 */
  1008.     } /* is receiver */
  1009.   }
  1010.  
  1011.   /* reset shooting patch's unshot rad */
  1012.   for(i=0;i<MAX_SPECTRA_SAMPLES;i++)
  1013.     shootPatch->unshot_B.samples[i] = 0.0;
  1014. }
  1015.  
  1016. /**********************************************************************/
  1017. /* Clean up                                                           */
  1018. /**********************************************************************/
  1019. void CleanUp_Rad()
  1020. {
  1021.   if (Option.debug) printf("\t*** Cleaning up radiosity structures ***\n");
  1022.   
  1023.   if (Option.ff_raytrace == 0) {
  1024.     free(hemicube.topFactors);
  1025.     free(hemicube.sideFactors);
  1026.     free(hemicube.view.buffer);
  1027.     free(formfactors);
  1028.   }
  1029. }
  1030.