home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / D-G / GemsII / radiosity / rad.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  13.0 KB  |  514 lines  |  [TEXT/MPS ]
  1. /*****************************************************************************
  2. *  rad.c
  3. *
  4. *  This program contains three functions that should be called in sequence to
  5. *  perform radiosity rendering:
  6. *  InitRad(): Initialize radiosity.
  7. *  DoRad(): Perform the main radiosity iteration loop.
  8. *  CleanUpRad(): Clean up.
  9. *
  10. *  The following routines are assumed to be provided by the user:
  11. *  BeginDraw()
  12. *  DrawPolygon()
  13. *  EndDraw()
  14. *  Refer to rad.h for details
  15. *
  16. *  Copyright (C) 1990-1991 Apple Computer, Inc.
  17. *  All rights reserved.
  18. *
  19. *  12/1990 S. Eric Chen    
  20. ******************************************************************************/
  21.  
  22. #include "rad.h"
  23. #include <math.h>
  24. #include <stdlib.h>
  25.  
  26. #define kMaxPolyPoints    255
  27. #define PI    3.1415926
  28. #define AddVector(c,a,b) (c).x=(a).x+(b).x, (c).y=(a).y+(b).y, (c).z=(a).z+(b).z
  29. #define SubVector(c,a,b) (c).x=(a).x-(b).x, (c).y=(a).y-(b).y, (c).z=(a).z-(b).z
  30. #define CrossVector(c,a,b)    (c).x = (a).y*(b).z - (a).z*(b).y, \
  31.                                 (c).y = (a).z*(b).x - (a).x*(b).z, \
  32.                                 (c).z = (a).x*(b).y - (a).y*(b).x
  33. #define DotVector(a,b) (a).x*(b).x + (a).y*(b).y + (a).z*(b).z
  34. #define ScaleVector(c,s) (c).x*=(s), (c).y*=(s), (c).z*=(s)
  35. #define NormalizeVector(n,a)     ((n)=sqrt(DotVector(a,a)), \
  36.                                 (n)?((a).x/=n, (a).y/=n, (a).z/=n):0)
  37.  
  38. typedef struct {
  39.     TView    view;        /* we only need to store one face of the hemi-cube */
  40.     double*    topFactors;    /* delta form-factors(weight for each pixel) of the top face */
  41.     double*    sideFactors; /* delta form-factors of the side faces */
  42. } THemicube;
  43.  
  44. static TRadParams *params;    /* input parameters */
  45. static THemicube hemicube;    /* one hemi-cube */
  46. static double *formfactors;    /* a form-factor array which has the same length as the number of elements */
  47. static double totalEnergy;    /* total emitted energy; used for convergence checking */
  48.  
  49. static const TSpectra black = { 0, 0, 0 };    /* for initialization */
  50. static int FindShootPatch(unsigned long *shootPatch);
  51. static void SumFactors(double* formfactors, int xRes, int yRes, 
  52.     unsigned long* buf, double* deltaFactors);
  53. static void MakeTopFactors(int hres, double* deltaFactors);
  54. static void MakeSideFactors(int hres, double* deltaFactors);
  55. static void ComputeFormfactors(unsigned long shootPatch);
  56. static void DistributeRad(unsigned long shootPatch);
  57. static void DisplayResults(TView* view);
  58. static void DrawElement(TElement* ep, unsigned long color);
  59. static TColor32b SpectraToRGB(TSpectra* spectra);
  60.  
  61.  
  62. /* Initialize radiosity based on the input parameters p */
  63. void InitRad(TRadParams *p)
  64. {
  65.     int n;
  66.     int hRes;
  67.     unsigned long i;
  68.     int j;
  69.     TPatch*    pp;
  70.     TElement* ep;
  71.     
  72.     params = p;
  73.     
  74.     /* initialize hemi-cube */
  75.     hemicube.view.fovx = 90;
  76.     hemicube.view.fovy = 90;
  77.     /* make sure hemicube resolution is an even number */
  78.     hRes = ((int)(params->hemicubeRes/2.0+0.5))*2;
  79.     hemicube.view.xRes = hemicube.view.yRes = hRes;
  80.     n = hRes*hRes;
  81.     hemicube.view.buffer = calloc(n, sizeof(unsigned long));
  82.     hemicube.view.wid=0;
  83.     hemicube.view.near = params->worldSize*0.001;
  84.     hemicube.view.far = params->worldSize;
  85.     
  86.     /* take advantage of the symmetry in the delta form-factors */
  87.     hemicube.topFactors= calloc(n/4, sizeof(double));
  88.     hemicube.sideFactors= calloc(n/4, sizeof(double));
  89.     MakeTopFactors(hRes/2, hemicube.topFactors);
  90.     MakeSideFactors(hRes/2, hemicube.sideFactors);
  91.     
  92.     formfactors = calloc(params->nElements, sizeof(double));
  93.     
  94.     /* initialize radiosity */
  95.     pp = params->patches;
  96.     for (i=params->nPatches; i--; pp++)
  97.         pp->unshotRad = *(pp->emission);
  98.     ep = params->elements;
  99.     for (i=params->nElements; i--; ep++)
  100.         ep->rad = *(ep->patch->emission);
  101.  
  102.     /* compute total energy */
  103.     totalEnergy = 0;
  104.     pp = params->patches;
  105.     for (i=params->nPatches; i--; pp++)
  106.         for (j=0; j<kNumberOfRadSamples; j++)
  107.             totalEnergy += pp->emission->samples[j] * pp->area;
  108.  
  109.     DisplayResults(¶ms->displayView); 
  110.  
  111. }
  112.  
  113. /* Main iterative loop */
  114. void DoRad()
  115. {
  116.     unsigned long shootPatch;
  117.     
  118.     while (FindShootPatch(&shootPatch)) 
  119.     {
  120.         ComputeFormfactors(shootPatch);
  121.         DistributeRad(shootPatch);
  122.         DisplayResults(¶ms->displayView);
  123.     }
  124.     
  125. }
  126.  
  127. /* Clean up */
  128. void CleanUpRad()
  129. {
  130.     free(hemicube.topFactors);
  131.     free(hemicube.sideFactors);
  132.     free(hemicube.view.buffer);
  133.     free(formfactors);
  134.  
  135. }
  136.  
  137. /* Find the next shooting patch based on the unshot energy of each patch */
  138. /* Return 0 if convergence is reached; otherwise, return 1 */
  139. static int FindShootPatch(unsigned long *shootPatch)
  140. {
  141.     int i, j;
  142.     double energySum, error, maxEnergySum=0;
  143.     TPatch* ep;
  144.  
  145.     ep = params->patches;
  146.     for (i=0; i< params->nPatches; i++, ep++)
  147.     {
  148.         energySum =0;
  149.         for (j=0; j<kNumberOfRadSamples; j++)
  150.             energySum += ep->unshotRad.samples[j] * ep->area;
  151.         
  152.         if (energySum > maxEnergySum) 
  153.         {
  154.             *shootPatch = i;
  155.             maxEnergySum = energySum;
  156.         }
  157.     }
  158.  
  159.     error = maxEnergySum / totalEnergy;
  160.     /* check convergence */
  161.     if (error < params->threshold)
  162.         return (0);        /* converged */
  163.     else
  164.         return (1);
  165.     
  166.  
  167. }
  168.  
  169. /* Find out the index to the delta form-factors arrary */
  170. #define Index(i)    ((i)<hres? i: (hres-1- ((i)%hres)))
  171.  
  172. /* Use the largest 32bit unsigned long for background */
  173. #define kBackgroundItem 0xffffffff
  174.  
  175. /* Convert a hemi-cube face to form-factors */
  176. static void SumFactors(
  177. double* formfactors, /* output */
  178. int xRes, int yRes, /* resolution of the hemi-cube face */
  179. unsigned long* buf, /* we only need the storage of the top hemi-cube face */
  180. double* deltaFactors /* delta form-factors for each hemi-cube pixel */
  181. )
  182. {
  183.     int i, j;
  184.     int ii, jj;
  185.     unsigned long *ip=buf;
  186.     int hres = xRes/2;
  187.     for (i=0; i<yRes; i++) 
  188.     {
  189.         ii= Index(i)*hres;
  190.           for (j=0; j<xRes; j++, ip++) 
  191.             if ((*ip) != kBackgroundItem) 
  192.             {
  193.                 jj = Index(j);
  194.                 formfactors[*ip] += deltaFactors[ii+jj];
  195.             }
  196.     }
  197. }
  198.  
  199. /* Create the delta form-factors for the top face of hemi-cube */
  200. /* Only need to compute 1/4 of the form-factors because of the 4-way symmetry */
  201. static void MakeTopFactors(
  202. int hres, /* half resolution of the face */
  203. double* deltaFactors /* output */
  204. )
  205. {
  206.     int j,k;
  207.     double xSq , ySq, xy1Sq;
  208.     double n= hres;
  209.     double* wp;
  210.     double dj, dk;
  211.     
  212.     wp = deltaFactors;
  213.     for (j=0; j<hres; j++)
  214.     {
  215.         dj = (double)j;
  216.         ySq = (n - (j+0.5)) / n;
  217.            ySq *= ySq;
  218.            for ( k=0 ; k<hres ; k++ )
  219.            {
  220.             dk = (double)k;
  221.                xSq = ( n - (k + 0.5) ) / n;
  222.             xSq *= xSq;
  223.             xy1Sq =  xSq + ySq + 1.0 ;
  224.             xy1Sq *= xy1Sq;
  225.             *wp++ = 1.0 / (xy1Sq * PI * n * n);
  226.            }
  227.     }
  228. }
  229.  
  230. /* Create the delta form-factors for the side face of hemi-cube */
  231. /* Only need to compute 1/4 of the form-factors because of the 4-way symmetry */
  232. static void MakeSideFactors(
  233. int hres, /* half resolution of the face */
  234. double* deltaFactors /* output */
  235. )
  236. {
  237.     int j,k;
  238.     double x, xSq , y, ySq, xy1, xy1Sq;
  239.     double n= hres;
  240.     double* wp;
  241.     double dj, dk;
  242.     
  243.     wp = deltaFactors;
  244.     for (j=0; j<hres; j++)
  245.     {
  246.         dj = (double)j;
  247.         y = (n - (dj+0.5)) / n;
  248.            ySq = y*y;
  249.            for ( k=0 ; k<hres ; k++ )
  250.            {
  251.             dk = (double)k;
  252.                x = ( n - (dk + 0.5) ) / n;
  253.             xSq = x*x;
  254.             xy1 =  xSq + ySq + 1.0 ;
  255.             xy1Sq = xy1*xy1;
  256.             *wp++ = y / (xy1Sq * PI * n * n);
  257.            }
  258.     }
  259. }
  260.  
  261. /* Use drand48 instead if it is supported */
  262. #define RandomFloat ((float)(rand())/(float)RAND_MAX)
  263.  
  264. /* Compute form-factors from the shooting patch to every elements */
  265. static void ComputeFormfactors(unsigned long shootPatch)
  266. {
  267.     unsigned long i;
  268.     TVector3f    up[5]; 
  269.     TPoint3f    lookat[5];
  270.     TPoint3f    center;
  271.     TVector3f    normal, tangentU, tangentV, vec;
  272.     int face;
  273.     double        norm;
  274.     TPatch*        sp;
  275.     double*        fp;
  276.     TElement*    ep;
  277.  
  278.     /* get the center of shootPatch */
  279.     sp = &(params->patches[shootPatch]);
  280.     center = sp->center;
  281.     normal = sp->normal;
  282.     
  283.     /* rotate the hemi-cube along the normal axis of the patch randomly */
  284.     /* this will reduce the hemi-cube aliasing artifacts */
  285.     do {
  286.         vec.x = RandomFloat;
  287.         vec.y = RandomFloat;
  288.         vec.z = RandomFloat;
  289.         /* get a tangent vector */
  290.         CrossVector(tangentU, normal, vec);
  291.         NormalizeVector(norm, tangentU);
  292.     } while (norm==0);    /* bad choice of the radom vector */
  293.     
  294.     /* compute tangentV */
  295.     CrossVector(tangentV, normal, tangentU);
  296.     
  297.     /* assign the lookats and ups for each hemicube face */
  298.     AddVector(lookat[0], center, normal);
  299.     up[0] = tangentU;
  300.     AddVector(lookat[1], center, tangentU);
  301.     up[1] = normal;
  302.     AddVector(lookat[2], center, tangentV);
  303.     up[2] = normal;
  304.     SubVector(lookat[3], center, tangentU);
  305.     up[3] = normal;
  306.     SubVector(lookat[4], center, tangentV);
  307.     up[4] = normal;
  308.     
  309.     /* position the hemicube slightly above the center of the shooting patch */
  310.     ScaleVector(normal, params->worldSize*0.0001);
  311.     AddVector(hemicube.view.camera, center, normal);
  312.     
  313.     /* clear the formfactors */
  314.     fp = formfactors;
  315.     for (i=params->nElements; i--; fp++)
  316.         *fp = 0.0;
  317.         
  318.     for (face=0; face < 5; face++)
  319.     {
  320.         hemicube.view.lookat = lookat[face];
  321.         hemicube.view.up = up[face];
  322.  
  323.         /* draw elements */
  324.         BeginDraw(&(hemicube.view), kBackgroundItem);
  325.         for (i=0; i< params->nElements; i++)
  326.             DrawElement(¶ms->elements[i], i);    
  327.             /* color element i with its index */
  328.         EndDraw();
  329.         
  330.         /* get formfactors */
  331.         if (face==0)
  332.             SumFactors(formfactors, hemicube.view.xRes, hemicube.view.yRes, 
  333.                 hemicube.view.buffer, hemicube.topFactors);
  334.         else
  335.             SumFactors(formfactors, hemicube.view.xRes, hemicube.view.yRes/2, 
  336.                 hemicube.view.buffer, hemicube.sideFactors);
  337.     }
  338.     
  339.     /* compute reciprocal form-factors */
  340.     ep = params->elements;
  341.     fp = formfactors;
  342.     for (i=params->nElements; i--; ep++, fp++)
  343.     {
  344.         *fp *= sp->area / ep->area;
  345.  
  346.         /* This is a potential source of hemi-cube aliasing */
  347.         /* To do this right, we need to subdivide the shooting patch
  348.         and reshoot. For now we just clip it to unity */
  349.         if ((*fp) > 1.0)     *fp ==1.0;    
  350.     }
  351.  
  352. }
  353.  
  354. /* Distribute radiosity form shootPatch to every element */
  355. /* Reset the shooter's unshot radiosity to 0 */
  356. static void DistributeRad(unsigned long shootPatch)
  357. {
  358.     unsigned long i;
  359.     int j;
  360.     TPatch* sp;
  361.     TElement* ep;
  362.     double* fp;
  363.     TSpectra deltaRad;
  364.     double w;
  365.  
  366.     sp = &(params->patches[shootPatch]);
  367.     
  368.     /* distribute unshotRad to every element */
  369.     ep = params->elements;
  370.     fp = formfactors;
  371.     for (i=params->nElements; i--; ep++, fp++)
  372.     {
  373.         if ((*fp) != 0.0) 
  374.         {
  375.             for (j=0; j<kNumberOfRadSamples; j++)
  376.                  deltaRad.samples[j] =     sp->unshotRad.samples[j] * (*fp) * 
  377.                                         ep->patch->reflectance->samples[j];
  378.  
  379.             /* incremental element's radiosity and patch's unshot radiosity */
  380.             w = ep->area/ep->patch->area;
  381.             for (j=0; j<kNumberOfRadSamples; j++) 
  382.             {
  383.                 ep->rad.samples[j] += deltaRad.samples[j];
  384.                 ep->patch->unshotRad.samples[j] += deltaRad.samples[j] * w;
  385.             }
  386.         }
  387.     }
  388.     
  389.     /* reset shooting patch's unshot radiosity */
  390.     sp->unshotRad = black;
  391. }
  392.  
  393. /* Convert a TSpectra (radiosity) to a TColor32b (rgb color) */
  394. /* Assume the first three samples of the spectra are the r, g, b colors */
  395. /* More elaborated color space transformation could be performed here */
  396. static TColor32b
  397. SpectraToRGB(TSpectra* spectra)
  398. {
  399.     TColor32b    c;
  400.     TSpectra    r;
  401.     double     max=1.0;
  402.     int k;
  403.  
  404.     for (k=kNumberOfRadSamples; k--;) {
  405.         if (spectra->samples[k] > max)
  406.             max = spectra->samples[k];
  407.     }
  408.     /* Clip the intensity*/
  409.     r = *spectra;
  410.     if (max>1.0) {
  411.         for (k=kNumberOfRadSamples; k--; )
  412.             r.samples[k] /= max;
  413.     }
  414.     
  415.     /* Convert to a 32-bit color; Assume the first 3 samples in TSpectra 
  416.     are the r, g, b colors we want. Otherwise, do color conversion here */
  417.     c.a= 0;
  418.     c.r= (unsigned char) (r.samples[0] * 255.0 + 0.5);
  419.     c.g= (unsigned char) (r.samples[1] * 255.0 + 0.5);
  420.     c.b= (unsigned char) (r.samples[2] * 255.0 + 0.5);
  421.     
  422.     return c;
  423. }
  424.  
  425. static void
  426. GetAmbient(TSpectra* ambient)
  427. {
  428.     TPatch* p;
  429.     unsigned long i;
  430.     int k;
  431.     static int first = 1;
  432.     static TSpectra baseSum; 
  433.     TSpectra uSum;
  434.  
  435.     uSum=black;
  436.     if (first) {
  437.         double areaSum;
  438.         TSpectra rSum;
  439.         areaSum=0;
  440.         rSum=black;
  441.         /* sum area and (area*reflectivity) */
  442.         p= params->patches;
  443.         for (i=params->nPatches; i--; p++) {
  444.             areaSum += p->area;
  445.             for (k=kNumberOfRadSamples; k--; )
  446.                 rSum.samples[k] += p->reflectance->samples[k]* p->area;
  447.         }
  448.         for (k=kNumberOfRadSamples; k--; )
  449.             baseSum.samples[k] = areaSum - rSum.samples[k];
  450.         first = 0;
  451.     }
  452.  
  453.     /* sum (unshot radiosity * area) */
  454.     p= params->patches;
  455.     for (i=params->nPatches; i--; p++) {
  456.         for (k=kNumberOfRadSamples; k--; )
  457.             uSum.samples[k] += p->unshotRad.samples[k] * p->area;
  458.     }
  459.     
  460.     /* compute ambient */
  461.     for (k=kNumberOfRadSamples; k--; )
  462.         ambient->samples[k] = uSum.samples[k] / baseSum.samples[k];
  463.  
  464. }
  465.  
  466. static void
  467. DisplayResults(TView* view)
  468. {
  469.     unsigned long i;
  470.     register TElement* ep;
  471.     TSpectra ambient;
  472.     GetAmbient(&ambient);
  473.     
  474.     BeginDraw(view, 0);
  475.     ep = params->elements;
  476.     for (i=0; i< params->nElements; i++, ep++) {
  477.         TColor32b    c;
  478.         TSpectra  s;
  479.         int k;
  480.         /* add ambient approximation */
  481.         if (params->addAmbient) {
  482.             for (k=kNumberOfRadSamples; k--;)
  483.                 s.samples[k] = (ep->rad.samples[k] + (ambient.samples[k]*
  484.                     ep->patch->reflectance->samples[k]))*params->intensityScale;
  485.         } else {
  486.             for (k=kNumberOfRadSamples; k--; )
  487.                 s.samples[k] = ep->rad.samples[k]*params->intensityScale;
  488.         }
  489.         /* quantize color */
  490.         c = SpectraToRGB(&s);
  491.         DrawElement(ep, *(unsigned long*)&c);
  492.     }
  493.             
  494.     EndDraw();
  495.  
  496. }
  497.  
  498. static void
  499. DrawElement(TElement* ep, unsigned long color)
  500. {
  501.     static TPoint3f pts[kMaxPolyPoints];
  502.     int nPts = ep->nVerts;
  503.     int j;
  504.     for (j=0; j<nPts; j++)
  505.         pts[j] = params->points[ep->verts[j]];
  506.     
  507.     DrawPolygon(nPts, pts, &ep->normal, color);
  508.  
  509. }
  510.  
  511.  
  512.  
  513.  
  514.