home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / photons.c < prev    next >
C/C++ Source or Header  |  2000-10-16  |  72KB  |  2,381 lines
  1. /*******************************************************
  2.  
  3. File:  photons.c
  4. Author: Nathan Kopp
  5.  
  6. This is part of the Photon Mapping custom version of
  7. POV-Ray.
  8.  
  9.  
  10. ********************************************************/
  11.  
  12.  
  13. #include "frame.h"
  14. #include "povray.h"
  15. #include "vector.h"
  16. #include "povproto.h"
  17. #include "texture.h"  /* for FRAND() */
  18. #include "parse.h"    /* for Warn() */
  19. #include "matrices.h"
  20. #include "objects.h"
  21. #include "csg.h"
  22. #include "photons.h"
  23. #include "ray.h"
  24.  
  25. #define BLACK_LEVEL 0.003
  26.  
  27. /* ------------------------------------------------------ */
  28. /* global variables */
  29. /* ------------------------------------------------------ */
  30. int backtraceFlag;
  31.  
  32. PHOTON_OPTIONS photonOptions;
  33.  
  34. int InitBacktraceWasCalled;
  35.  
  36. /* ------------------------------------------------------ */
  37. /* external variables */
  38. /* ------------------------------------------------------ */
  39. extern int Trace_Level;
  40.  
  41. #ifdef MotionBlurPatch
  42. extern int TimeStamp;
  43. #endif
  44.  
  45. /* ------------------------------------------------------ */
  46. /* static functions */
  47. /* ------------------------------------------------------ */
  48. static PHOTON* AllocatePhoton(PHOTON_MAP *map);
  49. static void FreePhotonMemory();
  50. static void InitPhotonMemory();
  51. static void sortAndSubdivide(int start, int end, int sorted);
  52. static void buildTree(PHOTON_MAP *map);
  53. static void ShootPhotonsAtObject(OBJECT *Object, LIGHT_SOURCE *Light, int count);
  54. static void SearchThroughObjects(OBJECT *Object, LIGHT_SOURCE *Light, int count);
  55. static int savePhotonMap(void);
  56. static int loadPhotonMap(void);
  57. static void swapPhotons(int a, int b);
  58. static void PQInsert(PHOTON *photon, DBL d);
  59. static void PQDelMax(void);
  60. static void gatherPhotonsRec(int start, int end);
  61. static void setGatherOptions(PHOTON_MAP *map, int mediaMap);
  62.  
  63. /* ------------------------------------------------------ */
  64. /* static variables */
  65. /* ------------------------------------------------------ */
  66. /* 
  67.   These static variables are used to conserve stack space during
  68.   extensive recursion when gathering photons.  All of these static
  69.   variables end in "_s".
  70. */
  71. static PHOTON **map_s;  /* photon map */
  72. static DBL size_sq_s;   /* search radius squared */
  73. static DBL Size_s;      /* search radius (static) */
  74. static DBL dmax_s;      /* dynamic search radius... current maximum */
  75. static int TargetNum_s; /* how many to gather */
  76. static DBL *pt_s;       /* point around which we are gathering */
  77. static int numfound_s;  /* number of photons found */
  78. static DBL *norm_s;     /* surface normal */
  79. static DBL flattenFactor; /* amount to flatten the spher to make it */
  80.                           /* an ellipsoid when gathering photons */
  81.                           /* zero = no flatten, one = regular */
  82. static DBL photonCountEstimate;
  83.  
  84. /*****************************************************************************
  85.  
  86.  FUNCTION
  87.  
  88.    CheckPassThru()
  89.  
  90.    Checks to make sure that pass-through, high-density, and refraction
  91.    are not simultaneously selected.  If all three are turned on, we need
  92.    to choose an appropriate one to turn off.
  93.  
  94.   Preconditions:
  95.     'o' is an initialized object
  96.     'flag' is PH_FLAG_PASSTHRU, ph_flag_target, or PH_FLAG_RFR_ON
  97.          (this is which flag was set most recently)
  98.  
  99.   Postconditions:
  100.     One of these flags in 'o' is turned off, since they cannot all be turned on.
  101.  
  102. ******************************************************************************/
  103. void CheckPassThru(OBJECT *o, int flag)
  104. {
  105.   if( (o->Ph_Flags & PH_FLAG_PASSTHRU) &&
  106.       (o->Ph_Flags & PH_FLAG_TARGET) &&
  107.       !(o->Ph_Flags & PH_FLAG_RFR_OFF) )
  108.   {
  109.     switch (flag)
  110.     {
  111.       case PH_FLAG_PASSTHRU:
  112.         Warning(0, "Cannot use pass through with refraction & high density photons.\nTurning off refraction.\n");
  113.         SET_PH_FLAG(o, PH_FLAG_RFR_OFF, PH_FLAG_RFR_ON);
  114.         break;
  115.  
  116.       case PH_FLAG_TARGET:
  117.         if(o->Ph_Flags & PH_FLAG_RFR_ON)
  118.         {
  119.           Warning(0, "Cannot use pass through with refraction & high density photons.\nTurning off pass through.\n");
  120.           o->Ph_Flags &= ~PH_FLAG_PASSTHRU;
  121.         }
  122.         else
  123.         {
  124.           Warning(0, "Cannot use pass through with refraction & high density photons.\nTurning off pass refraction.\n");
  125.           SET_PH_FLAG(o, PH_FLAG_RFR_OFF, PH_FLAG_RFR_ON);
  126.         }
  127.         break;
  128.  
  129.       case PH_FLAG_RFR_ON:
  130.         Warning(0, "Cannot use pass through with refraction & high density photons.\nTurning off pass through.\n");
  131.         o->Ph_Flags &= ~PH_FLAG_PASSTHRU;
  132.         break;
  133.     }
  134.   }
  135. }
  136.  
  137. /*****************************************************************************
  138.  
  139.  FUNCTION
  140.  
  141.    InitBacktraceEverything()
  142.  
  143.    Allocates memory.
  144.    Initializes all photon mapping stuff.
  145.    Does not create the photon map.
  146.  
  147.    Preconditions: InitBacktraceEverything() not yet called
  148.                     or
  149.                   both InitBacktraceEverything() and FreeBacktraceEverything() called
  150.  
  151.    Postconditions:
  152.       If photonOptions.photonsEnabled is true, then
  153.         memory for photon mapping is allocated.
  154.       else
  155.         nothing is done
  156.  
  157. ******************************************************************************/
  158. void InitBacktraceEverything()
  159. {
  160.   int i;
  161.   double theta;
  162.  
  163.   if (photonOptions.photonsEnabled)
  164.   {
  165.     InitBacktraceWasCalled = TRUE;
  166.  
  167.     photonOptions.photonMap.head = NULL;
  168.     photonOptions.photonMap.numPhotons  = 0;
  169.     photonOptions.photonMap.numBlocks  = 0;
  170.     photonOptions.globalPhotonMap.head = NULL;
  171.     photonOptions.globalPhotonMap.numPhotons  = 0;
  172.     photonOptions.globalPhotonMap.numBlocks  = 0;
  173.     photonOptions.mediaPhotonMap.head = NULL;
  174.     photonOptions.mediaPhotonMap.numPhotons  = 0;
  175.     photonOptions.mediaPhotonMap.numBlocks  = 0;
  176.  
  177.     photonOptions.photonGatherList = (PHOTON**)POV_MALLOC(sizeof(PHOTON *)*photonOptions.maxGatherCount, "Photon Map Info");
  178.     photonOptions.photonDistances = (DBL *)POV_MALLOC(sizeof(DBL)*photonOptions.maxGatherCount, "Photon Map Info");
  179.  
  180.     InitPhotonMemory();
  181.  
  182.     /* create the sin/cos arrays for speed */
  183.     /* range is -127..+127  =>  0..254 */
  184.     photonOptions.sinTheta = (DBL *)POV_MALLOC(sizeof(DBL)*255, "Photon Map Info");
  185.     photonOptions.cosTheta = (DBL *)POV_MALLOC(sizeof(DBL)*255, "Photon Map Info");
  186.     for(i=0; i<255; i++)
  187.     {
  188.       theta = (double)(i-127)*M_PI/127.0;
  189.       photonOptions.sinTheta[i] = sin(theta);
  190.       photonOptions.cosTheta[i] = cos(theta);
  191.     }
  192.   }
  193. }
  194.  
  195. /* savePhotonMap()
  196.  
  197.   Saves the caustic photon map to a file.
  198.  
  199.   Preconditions:
  200.     InitBacktraceEverything was called
  201.     the photon map has been built and balanced
  202.     photonOptions.fileName contains the filename to save
  203.  
  204.   Postconditions:
  205.     Returns 1 if success, 0 if failure.
  206.     If success, the photon map has been written to the file.
  207. */
  208. static int savePhotonMap()
  209. {
  210.   PHOTON *ph;
  211.   FILE *f;
  212.   int i, err;
  213.  
  214.   f = fopen(photonOptions.fileName, WRITE_BINFILE_STRING);
  215.   if (!f) return 0;
  216.  
  217.   /* caustic photons */
  218.   fwrite(&photonOptions.photonMap.numPhotons, sizeof(photonOptions.photonMap.numPhotons),1,f);
  219.   fwrite(&photonOptions.photonMap.rangeSelector, sizeof(photonOptions.photonMap.rangeSelector),1,f);
  220.   if (photonOptions.photonMap.numPhotons>0 && photonOptions.photonMap.head)
  221.   {
  222.     for(i=0; i<photonOptions.photonMap.numPhotons; i++)
  223.     {
  224.       ph = &(PHOTON_AMF(photonOptions.photonMap.head, i));
  225.       err = fwrite(ph, sizeof(PHOTON), 1, f);
  226.  
  227.       if (err<=0)
  228.       {
  229.         /* fwrite returned an error! */
  230.         fclose(f);
  231.         return 0;
  232.       }
  233.     }
  234.   }
  235.  
  236.   /* global photons */
  237.   fwrite(&photonOptions.globalPhotonMap.numPhotons, sizeof(photonOptions.globalPhotonMap.numPhotons),1,f);
  238.   fwrite(&photonOptions.globalPhotonMap.rangeSelector, sizeof(photonOptions.globalPhotonMap.rangeSelector),1,f);
  239.   if (photonOptions.globalPhotonMap.numPhotons>0 && photonOptions.globalPhotonMap.head)
  240.   {
  241.     for(i=0; i<photonOptions.globalPhotonMap.numPhotons; i++)
  242.     {
  243.       ph = &(PHOTON_AMF(photonOptions.globalPhotonMap.head, i));
  244.       err = fwrite(ph, sizeof(PHOTON), 1, f);
  245.  
  246.       if (err<=0)
  247.       {
  248.         /* fwrite returned an error! */
  249.         fclose(f);
  250.         return 0;
  251.       }
  252.     }
  253.   }
  254.  
  255.   /* media photons */
  256.   fwrite(&photonOptions.mediaPhotonMap.numPhotons, sizeof(photonOptions.mediaPhotonMap.numPhotons),1,f);
  257.   fwrite(&photonOptions.mediaPhotonMap.rangeSelector, sizeof(photonOptions.mediaPhotonMap.rangeSelector),1,f);
  258.   if (photonOptions.mediaPhotonMap.numPhotons>0 && photonOptions.mediaPhotonMap.head)
  259.   {
  260.     for(i=0; i<photonOptions.mediaPhotonMap.numPhotons; i++)
  261.     {
  262.       ph = &(PHOTON_AMF(photonOptions.mediaPhotonMap.head, i));
  263.       err = fwrite(ph, sizeof(PHOTON), 1, f);
  264.  
  265.       if (err<=0)
  266.       {
  267.         /* fwrite returned an error! */
  268.         fclose(f);
  269.         return 0;
  270.       }
  271.     }
  272.   }
  273.  
  274.   fclose(f);
  275.   return TRUE;
  276. }
  277.  
  278. /* loadPhotonMap()
  279.  
  280.   Loads the caustic photon map from a file.
  281.  
  282.   Preconditions:
  283.     InitBacktraceEverything was called
  284.     the photon map is empty
  285.     photonOptions.fileName contains the filename to load
  286.  
  287.   Postconditions:
  288.     Returns 1 if success, 0 if failure.
  289.     If success, the photon map has been loaded from the file.
  290.     If failure then the render should stop with an error
  291. */
  292. static int loadPhotonMap()
  293. {
  294.   int i;
  295.   int err;
  296.   PHOTON *ph;
  297.   FILE *f;
  298.   int numph;
  299.  
  300.   if (!photonOptions.photonsEnabled) return 0;
  301.  
  302.   f = fopen(photonOptions.fileName, READ_BINFILE_STRING);
  303.   if (!f) return 0;
  304.   fread(&numph, sizeof(numph),1,f);
  305.   fread(&photonOptions.photonMap.rangeSelector, sizeof(photonOptions.photonMap.rangeSelector),1,f);
  306.  
  307.   for(i=0; i<numph; i++)
  308.   {
  309.     ph = AllocatePhoton(&photonOptions.photonMap);
  310.     err = fread(ph, sizeof(PHOTON), 1, f);
  311.  
  312.     if (err<=0)
  313.     {
  314.       /* fread returned an error! */
  315.       fclose(f);
  316.       return 0;
  317.     }
  318.   }
  319.  
  320.   if (!feof(f)) /* for backwards file format compatibility */
  321.   {
  322.     /* global photons */
  323.     if (fread(&numph, sizeof(numph),1,f) > 0)
  324.     {
  325.       fread(&photonOptions.globalPhotonMap.rangeSelector, sizeof(photonOptions.globalPhotonMap.rangeSelector),1,f);
  326.       for(i=0; i<numph; i++)
  327.       {
  328.         ph = AllocatePhoton(&photonOptions.globalPhotonMap);
  329.         err = fread(ph, sizeof(PHOTON), 1, f);
  330.  
  331.         if (err<=0)
  332.         {
  333.           /* fread returned an error! */
  334.           fclose(f);
  335.           return 0;
  336.         }
  337.       }
  338.  
  339.       /* media photons */
  340.       fread(&numph, sizeof(numph),1,f);
  341.       fread(&photonOptions.mediaPhotonMap.rangeSelector, sizeof(photonOptions.mediaPhotonMap.rangeSelector),1,f);
  342.       for(i=0; i<numph; i++)
  343.       {
  344.         ph = AllocatePhoton(&photonOptions.mediaPhotonMap);
  345.         err = fread(ph, sizeof(PHOTON), 1, f);
  346.  
  347.         if (err<=0)
  348.         {
  349.           /* fread returned an error! */
  350.           fclose(f);
  351.           return 0;
  352.         }
  353.       }
  354.     }
  355.   }
  356.  
  357.   fclose(f);
  358.   return TRUE;
  359. }
  360.  
  361. /*****************************************************************************
  362.  
  363.  FUNCTION
  364.  
  365.   FreeBacktraceEverything()
  366.  
  367.   Preconditions:
  368.     if photonOptions.photonsEnabled is true, then InitBacktraceEverything()
  369.       must have been called
  370.     
  371.   PostConditions:
  372.     if photonOptions.photonsEnabled is true, then
  373.       photon memory is freed
  374.       sets photonOptions.photonsEnabled to false
  375.     else
  376.       does nothing
  377.  
  378. ******************************************************************************/
  379. void FreeBacktraceEverything()
  380. {
  381.   if (!InitBacktraceWasCalled) return;
  382.  
  383.   if (photonOptions.photonsEnabled)
  384.   {
  385.     /* free everything that we allocated */
  386.  
  387.     if(photonOptions.photonGatherList)
  388.       POV_FREE(photonOptions.photonGatherList);
  389.     photonOptions.photonGatherList = NULL;
  390.  
  391.     if(photonOptions.photonDistances)
  392.       POV_FREE(photonOptions.photonDistances);
  393.     photonOptions.photonDistances = NULL;
  394.  
  395.     if (photonOptions.sinTheta)
  396.       POV_FREE(photonOptions.sinTheta);
  397.     photonOptions.sinTheta = NULL;
  398.  
  399.     if (photonOptions.cosTheta)
  400.       POV_FREE(photonOptions.cosTheta);
  401.     photonOptions.cosTheta = NULL;
  402.  
  403.     FreePhotonMemory();
  404.     photonOptions.photonsEnabled = FALSE;
  405.   }
  406. }
  407.  
  408. /*****************************************************************************
  409.  
  410.  FUNCTION
  411.  
  412.   AllocatePhoton(PHOTON_MAP *map)
  413.     allocates a photon
  414.  
  415.     Photons are allocated in blocks.  map->head is a
  416.     dynamically-created array of these blocks.  The blocks themselves
  417.     are allocated as they are needed.
  418.  
  419.   Preconditions:
  420.     InitBacktraceEverything was called
  421.  
  422.   Postconditions:
  423.     Marks another photon as allocated (and allocates another block of
  424.     photons if necessary).
  425.     Returns a pointer to the new photon.
  426.     This will be the next available photon in array.
  427.  
  428. ******************************************************************************/
  429. PHOTON* AllocatePhoton(PHOTON_MAP *map)
  430. {
  431.   int i,j,k;
  432.  
  433.   /* array mapping funciton */
  434.   /* !!!!!!!!!!! warning
  435.      This code does the same function as the macro PHOTON_AMF
  436.      It is done here separatly instead of using the macro for
  437.      speed reasons (to avoid duplicate operations).  If the
  438.      macro is changed, this MUST ALSO BE CHANGED!
  439.   */
  440.   i=(map->numPhotons & PHOTON_BLOCK_MASK);
  441.   j=(map->numPhotons >> (PHOTON_BLOCK_POWER));
  442.  
  443.   /* new photon */
  444.   map->numPhotons++;
  445.  
  446.   if(j == map->numBlocks)
  447.   {
  448.     /* the base array is too small, we need to reallocate it */
  449.     PHOTON **newMap;
  450.     newMap = (PHOTON **)POV_MALLOC(sizeof(PHOTON *)*map->numBlocks*2, "photons");
  451.     map->numBlocks*=2;
  452.  
  453.     /* copy entries */
  454.     for(k=0; k<j; k++)
  455.       newMap[k] = map->head[k];
  456.  
  457.     /* set new entries to zero */
  458.     for(k=j; k<map->numBlocks; k++)
  459.       newMap[k] = NULL;
  460.  
  461.     /* free old map and put the new map in place */
  462.     POV_FREE(map->head);
  463.     map->head = newMap;
  464.   }
  465.  
  466.   if(map->head[j] == NULL)
  467.     /* allocate a new block of photons */
  468.     map->head[j] = (PHOTON *)POV_MALLOC(sizeof(PHOTON)*PHOTON_BLOCK_SIZE, "photons");
  469.  
  470.   return &(map->head[j][i]);
  471. }
  472.  
  473.  
  474. /*****************************************************************************
  475.  
  476.  FUNCTION
  477.  
  478.    InitPhotonMemory()
  479.  
  480.   Initializes photon memory.
  481.   Must only be called by InitBacktraceEverything().
  482.  
  483. ******************************************************************************/
  484. static void InitPhotonMemory()
  485. {
  486.   int k;
  487.  
  488.   /* allocate the base array */
  489.   photonOptions.photonMap.numPhotons = 0;
  490.   photonOptions.photonMap.numBlocks = INITIAL_BASE_ARRAY_SIZE;
  491.   photonOptions.photonMap.head = (PHOTON_BLOCK *)POV_MALLOC(sizeof(PHOTON_BLOCK *)*INITIAL_BASE_ARRAY_SIZE, "photons");
  492.  
  493.   /* zero the array */
  494.   for(k=0; k<photonOptions.photonMap.numBlocks; k++)
  495.     photonOptions.photonMap.head[k] = NULL;
  496.  
  497.   /* ------------ global photons ----------------*/
  498.   /* allocate the base array */
  499.   photonOptions.globalPhotonMap.numPhotons = 0;
  500.   photonOptions.globalPhotonMap.numBlocks = INITIAL_BASE_ARRAY_SIZE;
  501.   photonOptions.globalPhotonMap.head = (PHOTON_BLOCK *)POV_MALLOC(sizeof(PHOTON_BLOCK *)*INITIAL_BASE_ARRAY_SIZE, "photons");
  502.  
  503.   /* zero the array */
  504.   for(k=0; k<photonOptions.globalPhotonMap.numBlocks; k++)
  505.     photonOptions.globalPhotonMap.head[k] = NULL;
  506.  
  507.   /* ------------ media photons ----------------*/
  508.   /* allocate the base array */
  509.   photonOptions.mediaPhotonMap.numPhotons = 0;
  510.   photonOptions.mediaPhotonMap.numBlocks = INITIAL_BASE_ARRAY_SIZE;
  511.   photonOptions.mediaPhotonMap.head = (PHOTON_BLOCK *)POV_MALLOC(sizeof(PHOTON_BLOCK *)*INITIAL_BASE_ARRAY_SIZE, "photons");
  512.  
  513.   /* zero the array */
  514.   for(k=0; k<photonOptions.mediaPhotonMap.numBlocks; k++)
  515.     photonOptions.mediaPhotonMap.head[k] = NULL;
  516. }
  517.  
  518. /*****************************************************************************
  519.  
  520.  FUNCTION
  521.  
  522.   FreePhotonMemory()
  523.  
  524.   Frees all allocated blocks and the base array.
  525.   Must be called only by FreeBacktraceEverything()
  526.  
  527. ******************************************************************************/
  528. static void FreePhotonMemory()
  529. {
  530.   int j;
  531.  
  532.   /* if already freed then stop now */
  533.   if (photonOptions.photonMap.head==NULL)
  534.     return;
  535.     /*YS sept 17 2000 Memory leak */
  536.   if ( photonOptions.fileName)          /* file name to load or save caustic photon map */
  537.   {
  538.           POV_FREE(photonOptions.fileName);
  539.           photonOptions.fileName=NULL;
  540.   }
  541.   
  542.   /* free all non-NULL arrays */
  543.   for(j=0; j<photonOptions.photonMap.numBlocks; j++)
  544.   {
  545.     if(photonOptions.photonMap.head[j] != NULL)
  546.     {
  547.       POV_FREE(photonOptions.photonMap.head[j]);
  548.     }
  549.   }
  550.  
  551.   /* free the base array */
  552.   POV_FREE(photonOptions.photonMap.head);
  553.   photonOptions.photonMap.head = NULL;
  554.  
  555.   /* ---------------- global photons -------------- */
  556.   /* if already freed then stop now */
  557.   if (photonOptions.globalPhotonMap.head==NULL)
  558.     return;
  559.  
  560.   /* free all non-NULL arrays */
  561.   for(j=0; j<photonOptions.globalPhotonMap.numBlocks; j++)
  562.   {
  563.     if(photonOptions.globalPhotonMap.head[j] != NULL)
  564.     {
  565.       POV_FREE(photonOptions.globalPhotonMap.head[j]);
  566.     }
  567.   }
  568.  
  569.   /* free the base array */
  570.   POV_FREE(photonOptions.globalPhotonMap.head);
  571.   photonOptions.globalPhotonMap.head = NULL;
  572.  
  573.   /* ---------------- media photons -------------- */
  574.   /* if already freed then stop now */
  575.   if (photonOptions.mediaPhotonMap.head==NULL)
  576.     return;
  577.  
  578.   /* free all non-NULL arrays */
  579.   for(j=0; j<photonOptions.mediaPhotonMap.numBlocks; j++)
  580.   {
  581.     if(photonOptions.mediaPhotonMap.head[j] != NULL)
  582.     {
  583.       POV_FREE(photonOptions.mediaPhotonMap.head[j]);
  584.     }
  585.   }
  586.  
  587.   /* free the base array */
  588.   POV_FREE(photonOptions.mediaPhotonMap.head);
  589.   photonOptions.mediaPhotonMap.head = NULL;
  590.  
  591. }
  592.  
  593.  
  594. /*****************************************************************************
  595. *
  596. * FUNCTION
  597. *
  598. *   cubic_spline  (copied from point.c for use with light attenuation )
  599. *
  600. * INPUT
  601. *   
  602. * OUTPUT
  603. *   
  604. * RETURNS
  605. *   
  606. * AUTHOR
  607. *
  608. *   POV-Ray Team
  609. *   
  610. * DESCRIPTION
  611. *
  612. *   Cubic spline that has tangents of slope 0 at x == low and at x == high.
  613. *   For a given value "pos" between low and high the spline value is returned.
  614. *
  615. * CHANGES
  616. *
  617. *   -
  618. *
  619. ******************************************************************************/
  620. static DBL cubic_spline(DBL low, DBL  high, DBL  pos)
  621. {
  622.   /* Check to see if the position is within the proper boundaries. */
  623.   if (pos < low)
  624.     return(0.0);
  625.   else
  626.   {
  627.     if (pos >= high)
  628.       return(1.0);
  629.   }
  630.  
  631.   /* Normalize to the interval [0...1]. */
  632.   pos = (pos - low) / (high - low);
  633.  
  634.   /* See where it is on the cubic curve. */
  635.   return(3 - 2 * pos) * pos * pos;
  636. }
  637.  
  638. /*****************************************************************************
  639.  
  640.  FUNCTION
  641.  
  642.   SearchThroughObjects()
  643.  
  644.   Searches through 'object' and all siblings  and children of 'object' to
  645.   locate objects with PH_FLAG_TARGET set.  This flag means that the object
  646.   receives photons.
  647.  
  648.   Preconditions:
  649.     Photon mapping initialized (InitBacktraceEverything() called)
  650.     'Object' is a object (with or without siblings)
  651.     'Light' is a light source in the scene
  652.  
  653.   Postconditions:
  654.     Photons may have been shot at the object, its siblings, or its children.
  655.  
  656. ******************************************************************************/
  657. static void SearchThroughObjects(OBJECT *Object, LIGHT_SOURCE *Light, int count)
  658. {
  659.   OBJECT *Sib;
  660.  
  661.   /* check this object and all siblings */
  662.   for (Sib = Object; Sib != NULL; Sib = Sib -> Sibling)
  663.   {
  664.     if ((Sib->Ph_Flags & PH_FLAG_TARGET) &&
  665.         !(Sib->Type & LIGHT_SOURCE_OBJECT))
  666.       
  667.     {
  668.       ShootPhotonsAtObject(Sib, Light, count);
  669.     }
  670.     else
  671.     {
  672.       /* if it has children, check them too */
  673.       #ifdef MotionBlurPatch
  674.       if ((Sib->Type & COMPOUND_OBJECT) && 
  675.          !(Sib->Type & MOTION_BLUR_OBJECT))
  676.       #else
  677.       if ((Sib->Type & COMPOUND_OBJECT) )
  678.       #endif
  679.       {
  680.         SearchThroughObjects(((CSG *)Sib)->Children, Light, count);
  681.       }
  682.     }
  683.   }
  684. }
  685.  
  686. /*****************************************************************************
  687.  
  688.  FUNCTION
  689.  
  690.   ShootPhotonsAtObject()
  691.  
  692.   Shoots photons from 'Light' to 'Object'
  693.  
  694.   Preconditions:
  695.     Photon mapping initialized (InitBacktraceEverything() called)
  696.     'Object' is a object with PH_FLAG_TARGET set
  697.     'Light' is a light source in the scene
  698.  
  699.     Possible future expansion:
  700.       Object==NULL means create a global photon map.
  701.  
  702.   Postconditions:
  703.     Photons may have been shot at the object, depending on a variety
  704.     of flags
  705.  
  706. ******************************************************************************/
  707. static void ShootPhotonsAtObject(OBJECT *Object, LIGHT_SOURCE *Light, int count)
  708. {
  709.   RAY Ray;                       /* ray that we shoot */
  710.   COLOUR Colour, PhotonColour;   /* light color and photon color */
  711.   int i;                         /* counter */
  712.   DBL theta, phi;                /* rotation angles */
  713.   DBL dtheta, dphi;              /* deltas for theta and phi */
  714.   DBL jittheta, jitphi;          /* jittered versions of theta and phi */
  715.   DBL mintheta,maxtheta,minphi,maxphi;
  716.                                  /* these are minimum and maximum for theta and
  717.                                      phi for the spiral shooting */
  718.   DBL dist;                      /* distance from light to bounding sphere */
  719.   DBL rad;                       /* radius of bounding sphere */
  720.   DBL st,ct;                     /* cos(theta) & sin(theta) for rotation */
  721.   DBL Attenuation;               /* light attenuation for spotlight */
  722.   DBL costheta_spot;
  723.   VECTOR up, left, ctr, toctr, v; /* vectors to determine direction of shot */
  724.   TRANSFORM Trans;               /* transformation for rotation */
  725.   int mergedFlags=0;             /* merged flags to see if we should shoot photons */
  726.   int notComputed=TRUE;          /* have the ray containers been computed for this point yet?*/
  727.   int hitAtLeastOnce = FALSE;    /* have we hit the object at least once - for autostop stuff */
  728.  
  729.   /* get the light source colour */
  730.   Assign_Colour(Colour, Light->Colour);
  731.  
  732.   /* set global variable stuff */
  733.   photonOptions.Light = Light;
  734.   photonOptions.photonObject = Object;
  735.  
  736.   /* first, check on various flags... make sure all is a go for this object */
  737.   if(Object)
  738.   {
  739.     mergedFlags = Light->Ph_Flags | photonOptions.photonObject->Ph_Flags;
  740.     photonOptions.lightFlags = Light->Ph_Flags;
  741.   }
  742.   else
  743.   {
  744.     mergedFlags = photonOptions.lightFlags = PH_FLAG_RFR_ON | PH_FLAG_RFL_ON; /*Light->Ph_Flags;*/
  745.   }
  746.  
  747.   if (!(((mergedFlags & PH_FLAG_RFR_ON) && !(mergedFlags & PH_FLAG_RFR_OFF)) ||
  748.       ((mergedFlags & PH_FLAG_RFL_ON) && !(mergedFlags & PH_FLAG_RFL_OFF))))
  749.     /* it is a no-go for this object... bail out now */
  750.     return;
  751.  
  752.  
  753.   if(Object)
  754.   {
  755.     /* find bounding sphere based on bounding box */
  756.     ctr[X] = photonOptions.photonObject->BBox.Lower_Left[X] + photonOptions.photonObject->BBox.Lengths[X] / 2.0;
  757.     ctr[Y] = photonOptions.photonObject->BBox.Lower_Left[Y] + photonOptions.photonObject->BBox.Lengths[Y] / 2.0;
  758.     ctr[Z] = photonOptions.photonObject->BBox.Lower_Left[Z] + photonOptions.photonObject->BBox.Lengths[Z] / 2.0;
  759.     VSub(v, ctr,photonOptions.photonObject->BBox.Lower_Left);
  760.     VLength(rad, v);
  761.  
  762.     /* find direction from object to bounding sphere */
  763.     VSub(toctr, ctr, Light->Center);
  764.     VLength(dist, toctr);
  765.  
  766.     VNormalizeEq(toctr);
  767.     if ( fabs(fabs(toctr[Z])- 1.) < .1 ) {
  768.       /* too close to vertical for comfort, so use cross product with horizon */
  769.       up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
  770.     }
  771.     else
  772.     {
  773.       up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
  774.     }
  775.     VCross(left, toctr, up);  VNormalizeEq(left);
  776.     /* VCross(up, toctr, left);  VNormalizeEq(up); not needed */
  777.  
  778.     /*
  779.    light   dist         ctr
  780.     * ------------------ +
  781.          ---___          |
  782.                ---___    | rad
  783.                      ---_|
  784.  
  785.     */
  786.  
  787.     mintheta = 0;
  788.     if (dist>=rad)
  789.       maxtheta = atan(rad/dist);
  790.     else
  791.     {
  792.       maxtheta = M_PI;
  793.       if (fabs(dist)<EPSILON)
  794.       {
  795.         Make_Vector(up, 1,0,0);
  796.         Make_Vector(left, 0,1,0);
  797.         Make_Vector(toctr, 0,0,1);
  798.       }
  799.       dist = rad;
  800.     }
  801.  
  802.     /* set the photon density - calculate angular density from spacial */
  803.     photonOptions.photonSpread = photonOptions.photonObject->Ph_Density*photonOptions.surfaceSeparation / dist;
  804.  
  805.     if (count)
  806.     {
  807.       /* try to guess the number of photons */
  808.       DBL x=rad / (photonOptions.photonObject->Ph_Density*photonOptions.surfaceSeparation);
  809.       x=x*x*M_PI;
  810.  
  811. #if(1)
  812.       if ( ((mergedFlags & PH_FLAG_RFR_ON) && !(mergedFlags & PH_FLAG_RFR_OFF)) &&
  813.            ((mergedFlags & PH_FLAG_RFL_ON) && !(mergedFlags & PH_FLAG_RFL_OFF)) )
  814.       {
  815.         x *= 1.5;  /* assume 2 times as many photons with both reflection & refraction */
  816.       }
  817.  
  818.       if ( !(photonOptions.photonObject->Ph_Flags & PH_FLAG_IGNORE_PHOTONS) )
  819.       {
  820.         if ( ((mergedFlags & PH_FLAG_RFR_ON) && !(mergedFlags & PH_FLAG_RFR_OFF)) )
  821.         {
  822.           if ( ((mergedFlags & PH_FLAG_RFL_ON) && !(mergedFlags & PH_FLAG_RFL_OFF)) )
  823.             x *= 3;  /* assume 3 times as many photons if ignore_photons not used */
  824.           else
  825.             x *= 2;  /* assume less for only refraction */
  826.         }
  827.       }
  828.  
  829.       x *= 0.5;  /* assume 1/2 of photons hit target object */
  830. #endif
  831.  
  832.       photonCountEstimate += x;
  833.       return;
  834.     }
  835.   }
  836.   else
  837.   {
  838.     /* set up for global photon map */
  839.     mintheta = 0;
  840.     maxtheta = M_PI;
  841.     /* determine photonSpread from a number? */
  842.     photonOptions.photonSpread = Light->Ph_Density*photonOptions.globalSeparation;
  843.     Make_Vector(up, 1,0,0);
  844.     Make_Vector(left, 0,1,0);
  845.     Make_Vector(toctr, 0,0,1);
  846.     dist = 1.0;
  847.  
  848.     if (count)
  849.     {
  850.       /* try to guess the number of photons */
  851.       photonCountEstimate += M_PI*4.0/(photonOptions.photonSpread*photonOptions.photonSpread);
  852.       return;
  853.     }
  854.   }
  855.  
  856.  
  857.   if(Light->Area_Light && Light->Photon_Area_Light)
  858.   {
  859.     /* multiply spread if we are using an area light */
  860.     photonOptions.photonSpread *= sqrt(Light->Area_Size1*Light->Area_Size2);
  861.   }
  862.  
  863.   /* calculate delta theta */
  864.   dtheta = atan(photonOptions.photonSpread);
  865.  
  866.   /* if photonSpread <= 0.0, we must return or we'll get stuck in an infinite loop! */
  867.   if (photonOptions.photonSpread <= 0.0)
  868.     return;
  869.  
  870.   /* ---------------------------------------------
  871.          main ray-shooting loop 
  872.      --------------------------------------------- */
  873.   i = 0;
  874.   notComputed = TRUE;
  875.   for(theta=mintheta; theta<maxtheta; theta+=dtheta)
  876.   {
  877.     photonOptions.hitObject = FALSE;
  878.     if (theta<EPSILON)
  879.       dphi=2*M_PI;
  880.     else
  881.       dphi=dtheta/sin(theta);
  882.     minphi = -M_PI + dphi*FRAND()*0.5;
  883.     maxphi = M_PI - dphi/2 + (minphi+M_PI);
  884.     for(phi=minphi; phi<maxphi; phi+=dphi)
  885.     {
  886.       int x_samples,y_samples;
  887.       int area_x, area_y;
  888.       /* ------------------- shoot one photon ------------------ */
  889.  
  890.       /* jitter theta & phi */
  891.       jitphi = phi + (dphi)*(FRAND() - 0.5)*1.0*photonOptions.jitter;
  892.       jittheta = theta + (dtheta)*(FRAND() - 0.5)*1.0*photonOptions.jitter;
  893.  
  894.       /* actually, shoot multiple samples for area light */
  895.       if(Light->Area_Light && Light->Photon_Area_Light)
  896.       {
  897.         x_samples = Light->Area_Size1;
  898.         y_samples = Light->Area_Size2;
  899.       }
  900.       else
  901.       {
  902.         x_samples = 1;
  903.         y_samples = 1;
  904.       }
  905.  
  906.       for(area_x=0; area_x<x_samples; area_x++)
  907.       for(area_y=0; area_y<y_samples; area_y++)
  908.       {
  909.  
  910.         Assign_Vector(Ray.Initial,Light->Center);
  911.  
  912.         if (Light->Area_Light)
  913.         {
  914.           /* ------------- area light ----------- */
  915.           /* we need to make new up, left, and toctr vectors so we can
  916.              do proper rotations of theta and phi about toctr.  The
  917.              ray's initial point and ending points are both jittered to
  918.              produce the area-light effect. */
  919.           DBL Jitter_u, Jitter_v, ScaleFactor;
  920.           VECTOR NewAxis1, NewAxis2;
  921.  
  922.           /* we must recompute the media containers (new start point) */
  923.           notComputed = TRUE;
  924.  
  925.           /*
  926.           Jitter_u = (int)(FRAND()*Light->Area_Size1);
  927.           Jitter_v = (int)(FRAND()*Light->Area_Size2);
  928.           */
  929.           Jitter_u = area_x; /*+(0.5*FRAND() - 0.25);*/
  930.           Jitter_v = area_y; /*+(0.5*FRAND() - 0.25);*/
  931.  
  932.           if (Light->Area_Size1 > 1)
  933.           {
  934.             ScaleFactor = Jitter_u/(DBL)(Light->Area_Size1 - 1) - 0.5;
  935.             VScale (NewAxis1, Light->Axis1, ScaleFactor)
  936.           }
  937.           else
  938.           {
  939.             Make_Vector(NewAxis1, 0.0, 0.0, 0.0);
  940.           }
  941.  
  942.           if (Light->Area_Size2 > 1)
  943.           {
  944.             ScaleFactor = Jitter_v/(DBL)(Light->Area_Size2 - 1) - 0.5;
  945.             VScale (NewAxis2, Light->Axis2, ScaleFactor)
  946.           }
  947.           else
  948.           {
  949.             Make_Vector(NewAxis2, 0.0, 0.0, 0.0);
  950.           }
  951.  
  952.           /* need a new toctr & left */
  953.           VAddEq(Ray.Initial, NewAxis1);
  954.           VAddEq(Ray.Initial, NewAxis2);
  955.  
  956.           VSub(toctr, ctr, Ray.Initial);
  957.           VLength(dist, toctr);
  958.  
  959.           VNormalizeEq(toctr);
  960.           if ( fabs(fabs(toctr[Z])- 1.) < .1 ) {
  961.             /* too close to vertical for comfort, so use cross product with horizon */
  962.             up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
  963.           }
  964.           else
  965.           {
  966.             up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
  967.           }
  968.           VCross(left, toctr, up);  VNormalizeEq(left);
  969.  
  970.           if (fabs(dist)<EPSILON)
  971.           {
  972.             Make_Vector(up, 1,0,0);
  973.             Make_Vector(left, 0,1,0);
  974.             Make_Vector(toctr, 0,0,1);
  975.           }
  976.         }
  977.  
  978.         /* rotate toctr by theta around up */
  979.         st = sin(jittheta);
  980.         ct = cos(jittheta);
  981.         /* use fast rotation */
  982.         v[X] = -st*left[X] + ct*toctr[X];
  983.         v[Y] = -st*left[Y] + ct*toctr[Y];
  984.         v[Z] = -st*left[Z] + ct*toctr[Z];
  985.  
  986.         /* then rotate by phi around toctr */
  987.         /* use POV funcitons... slower but easy */
  988.         Compute_Axis_Rotation_Transform(&Trans,toctr,jitphi);
  989.         MTransPoint(Ray.Direction, v, &Trans);
  990.  
  991.         /* ------ attenuation for spot/cylinder (copied from point.c) ---- */ 
  992.         Attenuation = 1.0;
  993.  
  994.         VDot(costheta_spot, Ray.Direction, Light->Direction);
  995.         /*costheta_spot *= -1.0; - not needed?*/
  996.  
  997.         /* ---------- spot light --------- */
  998.         if (Light->Light_Type == SPOT_SOURCE)
  999.         {
  1000.           if (costheta_spot > 0.0)
  1001.           {
  1002.             Attenuation = pow(costheta_spot, Light->Coeff);
  1003.  
  1004.             if (Light->Radius > 0.0)
  1005.               Attenuation *= cubic_spline(Light->Falloff, Light->Radius, costheta_spot);
  1006.  
  1007.           }
  1008.           else
  1009.             Attenuation = 0.0;
  1010.         }
  1011.         /* ---------- cylinder light ----------- */
  1012.         else if (Light->Light_Type == CYLINDER_SOURCE)
  1013.         {
  1014.           VECTOR P;
  1015.           DBL k, len;
  1016.  
  1017.           VDot(k, Ray.Direction, Light->Direction);
  1018.  
  1019.           if (k > 0.0)
  1020.           {
  1021.             VLinComb2(P, 1.0, Ray.Direction, -k, Light->Direction);
  1022.             VLength(len, P);
  1023.  
  1024.             if (len < Light->Falloff)
  1025.             {
  1026.               len = 1.0 - len / Light->Falloff;
  1027.               Attenuation = pow(len, Light->Coeff);
  1028.  
  1029.               if (Light->Radius > 0.0)
  1030.                 Attenuation *= cubic_spline(1.0 - Light->Radius / Light->Falloff, 1.0, len);
  1031.  
  1032.             }
  1033.             else
  1034.               Attenuation = 0.0;
  1035.           }
  1036.           else
  1037.             Attenuation = 0.0;
  1038.         }
  1039.  
  1040.         /* set up defaults for reflection, refraction */
  1041.         photonOptions.passThruPrev = TRUE;
  1042.         photonOptions.passThruThis = FALSE;
  1043.  
  1044.         photonOptions.photonDepth = 0.0;
  1045.         Trace_Level = 1;
  1046.         Total_Depth = 0.0;
  1047.         Increase_Counter(stats[Number_Of_Photons_Shot]);
  1048.  
  1049.         /* attenuate for area light extra samples */
  1050.         Attenuation/=(x_samples*y_samples);
  1051.  
  1052.         /* compute photon color from light source & attenuation */
  1053.         if (Attenuation<BLACK_LEVEL) continue;
  1054.  
  1055.         PhotonColour[0] = Colour[0]*Attenuation;
  1056.         PhotonColour[1] = Colour[1]*Attenuation;
  1057.         PhotonColour[2] = Colour[2]*Attenuation;
  1058.         PhotonColour[3] = 0.0;
  1059.         PhotonColour[4] = 0.0;
  1060.  
  1061.         /* As mike said, "fire photon torpedo!" */
  1062.         Initialize_Ray_Containers(&Ray);
  1063.         initialize_ray_container_state(&Ray, notComputed);
  1064.         notComputed = FALSE;
  1065. #ifdef MotionBlurPatch
  1066.         /*TimeStamp = opts.photonTimeStamp;*/
  1067.         TimeStamp = (opts.motionBlurCount+1)/2; /* must be non-zero */
  1068. #endif
  1069.         Trace(&Ray, PhotonColour, 1.0, NULL);
  1070.  
  1071.         /* display here */
  1072.         i++;
  1073.         if ((i%100) == 0)
  1074.         {
  1075. #ifndef ProgressOnOneLine        
  1076.           Status_Info ("\rBuilding Photon Maps... this=%d, total=%d, m=%d, area %d, %d : %d", i, photonOptions.photonMap.numPhotons, photonOptions.mediaPhotonMap.numPhotons,x_samples,y_samples,(int)photonOptions.photonMap.rangeSelector);
  1077. #else
  1078.           Status_Info ("\rthis=%5otal=%5d, m=%5d, area %5d, %5d : %5d            ", i, photonOptions.photonMap.numPhotons, photonOptions.mediaPhotonMap.numPhotons,x_samples,y_samples,(int)photonOptions.photonMap.rangeSelector);
  1079. #endif
  1080.         }
  1081.  
  1082.       } /* end of multiple samples */
  1083.     }
  1084.  
  1085.     /* if we didn't hit anything and we're past the autostop angle, then
  1086.        we should stop 
  1087.        
  1088.        as per suggestion from Smellenberg, changed autostop to a percentage
  1089.        of the object's bounding sphere. */
  1090.  
  1091.     /* suggested by Pabs, we only use autostop if we have it it once */
  1092.     if (photonOptions.hitObject) hitAtLeastOnce=TRUE;
  1093.  
  1094.     if (hitAtLeastOnce && !photonOptions.hitObject && photonOptions.photonObject)
  1095.       if (theta>photonOptions.autoStopPercent*maxtheta) break;
  1096.   } /* end of rays loop */
  1097.  
  1098. #ifndef ProgressOnOneLine        
  1099.   Status_Info ("\rBuilding Photon Maps... this=%d, total=%d, m=%d : %d", i, photonOptions.photonMap.numPhotons, photonOptions.mediaPhotonMap.numPhotons,(int)photonOptions.photonMap.rangeSelector);
  1100. #else
  1101.   Status_Info ("\rthis=%5otal=%5d, m=%5d : %5d            ", i, photonOptions.photonMap.numPhotons, photonOptions.mediaPhotonMap.numPhotons,(int)photonOptions.photonMap.rangeSelector);
  1102. #endif
  1103. }
  1104.  
  1105. /*****************************************************************************
  1106.  
  1107.   FUNCTION
  1108.  
  1109.   BuildPhotonMaps()
  1110.  
  1111.   This is the primary function for building photon maps.
  1112.  
  1113.   Preconditions:
  1114.     Photon memory is allocated (InitBacktraceEverything has been called).
  1115.     The entire scene has been parsed.
  1116.     The ray-tracer has been completely initialized.
  1117.  
  1118.   Postconditions:
  1119.     The photon map is built based on options specified in the scene file.
  1120.  
  1121. ******************************************************************************/
  1122. void BuildPhotonMaps(void)
  1123. {
  1124.   LIGHT_SOURCE *Light;  /* light source to use */
  1125.   int old_mtl; /* saved max_trace_level */
  1126.   DBL old_adc; /* saved adc_bailout */
  1127.  
  1128.   /* if not enabled, then return */
  1129.   if (!photonOptions.photonsEnabled)
  1130.     return;
  1131.  
  1132.   /* should we load the photon map instead of building? */
  1133.   if (photonOptions.fileName && photonOptions.loadFile)
  1134.   {
  1135.     /* status bar for user */
  1136.     Status_Info ("\nLoading Photon Maps...");
  1137.     if (!loadPhotonMap())
  1138.     {
  1139.       Error("Could not load photon map (%s)\n",photonOptions.fileName);
  1140.     }
  1141.  
  1142.     /* don't build */
  1143.     /* but do set photon options automatically */
  1144.  
  1145.     /* ----------- surface photons ------------- */
  1146.     if (photonOptions.photonMap.numPhotons>0)
  1147.     {
  1148.       setGatherOptions(&photonOptions.photonMap, FALSE);
  1149.     }
  1150.  
  1151.     /* ----------- global photons ------------- */
  1152.     if (photonOptions.globalPhotonMap.numPhotons>0)
  1153.     {
  1154.       setGatherOptions(&photonOptions.globalPhotonMap, FALSE);
  1155.     }
  1156.  
  1157.     /* ----------- media photons ------------- */
  1158.     if (photonOptions.mediaPhotonMap.numPhotons>0)
  1159.     {
  1160.       setGatherOptions(&photonOptions.mediaPhotonMap, TRUE);
  1161.     }
  1162.  
  1163.     return;
  1164.   }
  1165.  
  1166.   /* set flag so POV knows we're in backtrace step */
  1167.   backtraceFlag = 1;
  1168.  
  1169.   /* status bar for user */
  1170.   Status_Info ("\nBuilding Photon Maps...");
  1171.  
  1172.   /* save adc_bailout and max_trace_level */
  1173.   old_adc = ADC_Bailout;
  1174.   old_mtl = Max_Trace_Level;
  1175.  
  1176.   /* use the photon-specific ones if they were specified */
  1177.   if (photonOptions.Max_Trace_Level>=0)
  1178.     Max_Trace_Level = photonOptions.Max_Trace_Level;
  1179.  
  1180.   if (photonOptions.ADC_Bailout>=0)
  1181.     ADC_Bailout = photonOptions.ADC_Bailout;
  1182.  
  1183.   /*  COUNT THE PHOTONS  */
  1184.   if(photonOptions.surfaceCount>0)
  1185.   {
  1186.     DBL factor;
  1187.     photonCountEstimate = 0.0;
  1188.  
  1189.     for (Light = Frame.Light_Sources;
  1190.          Light != NULL;
  1191.          Light = Light->Next_Light_Source)
  1192.     if (Light->Light_Type != FILL_LIGHT_SOURCE)
  1193.     {
  1194.       SearchThroughObjects(Frame.Objects, Light, TRUE);
  1195.     }
  1196.     factor = (DBL)photonCountEstimate/photonOptions.surfaceCount;
  1197.     factor = sqrt(factor);
  1198.     photonOptions.surfaceSeparation *= factor;
  1199.   }
  1200.  
  1201.   /*  COUNT THE GLOBAL PHOTONS  */
  1202.   if(photonOptions.globalCount>0)
  1203.   {
  1204.     DBL factor;
  1205.     photonCountEstimate = 0.0;
  1206.  
  1207.     for (Light = Frame.Light_Sources;
  1208.          Light != NULL;
  1209.          Light = Light->Next_Light_Source)
  1210.     if (Light->Light_Type != FILL_LIGHT_SOURCE)
  1211.     {
  1212.       ShootPhotonsAtObject(NULL, Light, TRUE);
  1213.     }
  1214.     factor = (DBL)photonCountEstimate/photonOptions.globalCount;
  1215.     factor = sqrt(factor);
  1216.     photonOptions.globalSeparation *= factor;
  1217.   }
  1218.  
  1219.   /*  loop through light sources  */
  1220.   for (Light = Frame.Light_Sources;
  1221.        Light != NULL;
  1222.        Light = Light->Next_Light_Source)
  1223.   if (Light->Light_Type != FILL_LIGHT_SOURCE)
  1224.   {
  1225.     /* give a warning for any parallel light, since they do not
  1226.     currently function properly */
  1227.     if (Light->Parallel)
  1228.     {
  1229.       Warning(0,"Parallel lights not supported by photons. Using a point source instead.\n");
  1230.     }
  1231.     if (Light->Light_Type == CYLINDER_SOURCE)
  1232.     {
  1233.       Warning(0,"Cylinder lights do not work perfectly with photons.\n");
  1234.     }
  1235.  
  1236.     /* do global lighting here if it is ever implemented */
  1237.     if ((Light->Ph_Flags & PH_FLAG_TARGET) && (photonOptions.globalCount>0))
  1238.       ShootPhotonsAtObject(NULL, Light, FALSE);
  1239.  
  1240.     /* do object-specific lighting */
  1241.     SearchThroughObjects(Frame.Objects, Light, FALSE);
  1242.   }
  1243.  
  1244.   /* clear this flag */
  1245.   backtraceFlag = 0;
  1246.  
  1247.   /* restore saved variables */
  1248.   ADC_Bailout = old_adc;
  1249.   Max_Trace_Level = old_mtl;
  1250.  
  1251.   /* now actually build the kd-tree by sorting the array of photons */
  1252.   if (photonOptions.photonMap.numPhotons>0)
  1253.   {
  1254.     buildTree(&photonOptions.photonMap);
  1255.     setGatherOptions(&photonOptions.photonMap, FALSE);
  1256.   }
  1257.  
  1258.   /* ----------- global photons ------------- */
  1259.   if (photonOptions.globalPhotonMap.numPhotons>0)
  1260.   {
  1261.     buildTree(&photonOptions.globalPhotonMap);
  1262.     setGatherOptions(&photonOptions.globalPhotonMap, FALSE);
  1263.   }
  1264.  
  1265.   /* ----------- media photons ------------- */
  1266.   if (photonOptions.mediaPhotonMap.numPhotons>0)
  1267.   {
  1268.     buildTree(&photonOptions.mediaPhotonMap);
  1269.     setGatherOptions(&photonOptions.mediaPhotonMap, TRUE);
  1270.   }
  1271.  
  1272.   if (photonOptions.photonMap.numPhotons+
  1273.       photonOptions.globalPhotonMap.numPhotons+
  1274.       photonOptions.mediaPhotonMap.numPhotons > 0)
  1275.   {
  1276.     /* should we load the photon map now that it is built? */
  1277.     if (photonOptions.fileName && !photonOptions.loadFile)
  1278.     {
  1279.       /* status bar for user */
  1280.       Status_Info ("\nSaving Photon Maps...");
  1281.       if (!savePhotonMap())
  1282.       {
  1283.         Warn(0,"Could not save photon map.\n");
  1284.       }
  1285.     }
  1286.   }
  1287.   else
  1288.   {
  1289.     if (photonOptions.fileName && !photonOptions.loadFile)
  1290.     {
  1291.       Warn(0,"Could not save photon map - no photons!\n");
  1292.     }
  1293.   }
  1294.  
  1295.  
  1296. }
  1297.  
  1298. /*****************************************************************************
  1299.  
  1300.   FUNCTION
  1301.  
  1302.   addSurfacePhoton()
  1303.  
  1304.   Adds a photon to the array of photons.
  1305.  
  1306.   Preconditions:
  1307.     InitBacktraceEverything() was called
  1308.     'Point' is the intersection point to store the photon
  1309.     'Origin' is the origin of the light ray
  1310.     'LightCol' is the color of the light propogated through the scene
  1311.     'RawNorm' is the raw normal of the surface intersected
  1312.  
  1313.   Postconditions:
  1314.     Another photon is allocated (by AllocatePhoton())
  1315.     The information passed in (as well as photonOptions.photonDepth)
  1316.       is stored in the photon data structure.
  1317.  
  1318. ******************************************************************************/
  1319.  
  1320. void addSurfacePhoton(VECTOR Point, VECTOR Origin, COLOUR LightCol, VECTOR RawNorm)
  1321. {
  1322.   PHOTON *Photon;
  1323.   int intensity;
  1324.   COLOUR LightCol2;
  1325.   DBL Attenuation;
  1326.   VECTOR d;
  1327.   DBL d_len, phi, theta;
  1328.   PHOTON_MAP *map;
  1329.  
  1330.   /* first, compensate for POV's weird light attenuation */
  1331.   if ((photonOptions.Light->Fade_Power > 0.0) && (fabs(photonOptions.Light->Fade_Distance) > EPSILON))
  1332.   {
  1333.     Attenuation = 2.0 / (1.0 + pow(photonOptions.photonDepth / photonOptions.Light->Fade_Distance, photonOptions.Light->Fade_Power));
  1334.   }
  1335.   else
  1336.     Attenuation = 1;
  1337.  
  1338.   VScale(LightCol2, LightCol, Attenuation);
  1339.   VScaleEq(LightCol2, photonOptions.photonDepth*photonOptions.photonDepth);
  1340.   VScaleEq(LightCol2, photonOptions.photonSpread*photonOptions.photonSpread);
  1341.  
  1342.   /* if too dark, maybe we should stop here */
  1343.  
  1344.   if(photonOptions.photonObject==NULL)
  1345.   {
  1346.     map = &photonOptions.globalPhotonMap;
  1347.     Increase_Counter(stats[Number_Of_Global_Photons_Stored]);
  1348.   }
  1349.   else
  1350.   {
  1351.     map = &photonOptions.photonMap;
  1352.     Increase_Counter(stats[Number_Of_Photons_Stored]);
  1353.   }
  1354.  
  1355.  
  1356.  
  1357.   /* allocate the photon */
  1358.   Photon = AllocatePhoton(map);
  1359.  
  1360.   /* adaptively choose a range selector for RGBI format */
  1361.   /* only do this if we didn't use it yet */
  1362.   /* and if user didn't give a range selector */
  1363.   if (map->rangeSelector == 0)
  1364.   {
  1365.     map->rangeSelector = 2.0/max3(LightCol2[0],LightCol2[1],LightCol2[2]);
  1366.   }
  1367.  
  1368.   /* convert photon from three floats to two bytes
  1369.      this is currently stored in RGBI format (designed by Nathan Kopp)
  1370.      RGBE format (by Greg Ward) could be used instead
  1371.   */
  1372.   intensity = (int)(256.0/(max3(LightCol2[0],LightCol2[1],LightCol2[2])*map->rangeSelector))-1;
  1373.   if (intensity<0) intensity = 0;
  1374.   if (intensity>255) intensity = 255;
  1375.  
  1376.   Photon->Colour.i = intensity;
  1377.   Photon->Colour.r = (int)(LightCol2[0]*(map->rangeSelector*(intensity+1)));
  1378.   Photon->Colour.g = (int)(LightCol2[1]*(map->rangeSelector*(intensity+1)));
  1379.   Photon->Colour.b = (int)(LightCol2[2]*(map->rangeSelector*(intensity+1)));
  1380.  
  1381.   /* store the location */
  1382.   Assign_SNGL_Vect(Photon->Loc, Point);
  1383.  
  1384.   /* now determine rotation angles */
  1385.   VSub(d,Origin, Point);
  1386.   VNormalizeEq(d);
  1387.   d_len = sqrt(d[X]*d[X]+d[Z]*d[Z]);
  1388.  
  1389.   phi = acos(d[X]/d_len);
  1390.   if (d[Z]<0) phi = -phi;
  1391.  
  1392.   theta = acos(d_len);
  1393.   if (d[Y]<0) theta = -theta;
  1394.  
  1395.   /* cram these rotation angles into two signed bytes */
  1396.   Photon->theta=(char)(theta*127.0/M_PI);
  1397.   Photon->phi=(char)(phi*127.0/M_PI);
  1398.  
  1399. }
  1400.  
  1401. /*****************************************************************************
  1402.  
  1403.   FUNCTION
  1404.  
  1405.   addMediaPhoton()
  1406.  
  1407.   Adds a photon to the array of photons.
  1408.  
  1409.   Preconditions:
  1410.     InitBacktraceEverything() was called
  1411.     'Point' is the intersection point to store the photon
  1412.     'Origin' is the origin of the light ray
  1413.     'LightCol' is the color of the light propogated through the scene
  1414.  
  1415.   Postconditions:
  1416.     Another photon is allocated (by AllocatePhoton())
  1417.     The information passed in (as well as photonOptions.photonDepth)
  1418.       is stored in the photon data structure.
  1419.  
  1420. ******************************************************************************/
  1421.  
  1422. void addMediaPhoton(VECTOR Point, VECTOR Origin, COLOUR LightCol, DBL depthDiff)
  1423. {
  1424.   PHOTON *Photon;
  1425.   int intensity;
  1426.   COLOUR LightCol2;
  1427.   DBL Attenuation;
  1428.   VECTOR d;
  1429.   DBL d_len, phi, theta;
  1430.  
  1431.   /* first, compensate for POV's weird light attenuation */
  1432.   if ((photonOptions.Light->Fade_Power > 0.0) && (fabs(photonOptions.Light->Fade_Distance) > EPSILON))
  1433.   {
  1434.     Attenuation = 2.0 / (1.0 + pow((photonOptions.photonDepth+depthDiff) / photonOptions.Light->Fade_Distance, photonOptions.Light->Fade_Power));
  1435.   }
  1436.   else
  1437.     Attenuation = 1;
  1438.  
  1439. #if 0
  1440.   VScale(LightCol2, LightCol, photonOptions.photonSpread*photonOptions.photonSpread);
  1441.   VScaleEq(LightCol2, (photonOptions.photonDepth+depthDiff)*(photonOptions.photonDepth+depthDiff)*Attenuation);
  1442. #else
  1443.   VScale(LightCol2, LightCol, Attenuation);
  1444.   VScaleEq(LightCol2, (photonOptions.photonDepth+depthDiff) *
  1445.                       (photonOptions.photonDepth+depthDiff));
  1446.   VScaleEq(LightCol2, photonOptions.photonSpread*photonOptions.photonSpread);
  1447. #endif
  1448.  
  1449.   /* if too dark, maybe we should stop here */
  1450.  
  1451.   /* allocate the photon */
  1452.   if(photonOptions.photonObject==NULL) return;
  1453.  
  1454.   Increase_Counter(stats[Number_Of_Media_Photons_Stored]);
  1455.  
  1456.   Photon = AllocatePhoton(&photonOptions.mediaPhotonMap);
  1457.  
  1458.   /* adaptively choose a range selector for RGBI format */
  1459.   /* only do this if we didn't use it yet */
  1460.   /* and if user didn't give a range selector */
  1461.   if (photonOptions.mediaPhotonMap.rangeSelector == 0)
  1462.   {
  1463.     photonOptions.mediaPhotonMap.rangeSelector = 2.0/max3(LightCol2[0],LightCol2[1],LightCol2[2]);
  1464.   }
  1465.  
  1466.   /* convert photon from three floats to two bytes
  1467.      this is currently stored in RGBI format (designed by Nathan Kopp)
  1468.      RGBE format (by Greg Ward) could be used instead
  1469.   */
  1470.   intensity = (int)(256.0/(max3(LightCol2[0],LightCol2[1],LightCol2[2])*photonOptions.mediaPhotonMap.rangeSelector))-1;
  1471.   if (intensity<0) intensity = 0;
  1472.   if (intensity>255) intensity = 255;
  1473.  
  1474.   Photon->Colour.i = intensity;
  1475.   Photon->Colour.r = (int)(LightCol2[0]*(photonOptions.mediaPhotonMap.rangeSelector*(intensity+1)));
  1476.   Photon->Colour.g = (int)(LightCol2[1]*(photonOptions.mediaPhotonMap.rangeSelector*(intensity+1)));
  1477.   Photon->Colour.b = (int)(LightCol2[2]*(photonOptions.mediaPhotonMap.rangeSelector*(intensity+1)));
  1478.  
  1479.   /* store the location */
  1480.   Assign_SNGL_Vect(Photon->Loc, Point);
  1481.  
  1482.   /* now determine rotation angles */
  1483.   VSub(d,Origin, Point);
  1484.   VNormalizeEq(d);
  1485.   d_len = sqrt(d[X]*d[X]+d[Z]*d[Z]);
  1486.  
  1487.   phi = acos(d[X]/d_len);
  1488.   if (d[Z]<0) phi = -phi;
  1489.  
  1490.   theta = acos(d_len);
  1491.   if (d[Y]<0) theta = -theta;
  1492.  
  1493.   /* cram these rotation angles into two signed bytes */
  1494.   Photon->theta=(char)(theta*127.0/M_PI);
  1495.   Photon->phi=(char)(phi*127.0/M_PI);
  1496.  
  1497. }
  1498.  
  1499. /* ====================================================================== */
  1500. /* ====================================================================== */
  1501. /*                              KD - TREE                                 */
  1502. /* ====================================================================== */
  1503. /* ====================================================================== */
  1504.  
  1505. /*****************************************************************************
  1506.  
  1507.   FUNCTION
  1508.  
  1509.   swapPhotons
  1510.  
  1511.   swaps two photons
  1512.  
  1513.   Precondition:
  1514.     photon memory initialized
  1515.     static map_s points to the photon map
  1516.     'a' and 'b' are indexes within the range of photons in map_s
  1517.       (NO ERROR CHECKING IS DONE)
  1518.  
  1519.   Postconditions:
  1520.     the photons indexed by 'a' and 'b' are swapped
  1521.  
  1522. *****************************************************************************/
  1523.  
  1524. static void swapPhotons(int a, int b)
  1525. {
  1526.   int ai,aj,bi,bj;
  1527.   PHOTON tmp;
  1528.  
  1529.   /* !!!!!!!!!!! warning
  1530.      This code does the same function as the macro PHOTON_AMF
  1531.      It is done here separatly instead of using the macro for
  1532.      speed reasons (to avoid duplicate operations).  If the
  1533.      macro is changed, this MUST ALSO BE CHANGED!
  1534.   */
  1535.   ai = a & PHOTON_BLOCK_MASK;
  1536.   aj = a >> PHOTON_BLOCK_POWER;
  1537.   bi = b & PHOTON_BLOCK_MASK;
  1538.   bj = b >> PHOTON_BLOCK_POWER;
  1539.  
  1540.   tmp = map_s[aj][ai];
  1541.   map_s[aj][ai] = map_s[bj][bi];
  1542.   map_s[bj][bi] = tmp;
  1543. }
  1544.  
  1545. /*****************************************************************************
  1546.  
  1547.   FUNCTION
  1548.  
  1549.   insertSort
  1550.   (modified from Data Structures textbook)
  1551.  
  1552.   Preconditions:
  1553.     photon memory initialized
  1554.     static map_s points to the photon map
  1555.     'start' is the index of the first photon
  1556.     'end' is the index of the last photon
  1557.     'd' is the dimension to sort on (X, Y, or Z)
  1558.  
  1559.   Postconditions:
  1560.     photons from 'start' to 'end' in map_s are sorted in
  1561.     ascending order on dimension d
  1562. ******************************************************************************/
  1563. static void insertSort(int start, int end, int d)
  1564. {
  1565.   int j,k;
  1566.   PHOTON tmp;
  1567.  
  1568.   for(k=end-1; k>=start; k--)
  1569.   {
  1570.     j=k+1;
  1571.     tmp = PHOTON_AMF(map_s, k);
  1572.     while ( (tmp.Loc[d] > PHOTON_AMF(map_s,j).Loc[d]) )
  1573.     {
  1574.       PHOTON_AMF(map_s,j-1) = PHOTON_AMF(map_s,j);
  1575.       j++;
  1576.       if (j>end) break;
  1577.     }
  1578.     PHOTON_AMF(map_s,j-1) = tmp;
  1579.   }
  1580. }
  1581.  
  1582. /*****************************************************************************
  1583.  
  1584.   FUNCTION
  1585.  
  1586.   quickSortRec
  1587.   (modified from Data Structures textbook)
  1588.  
  1589.   Recursive part of the quicksort routine
  1590.   This does not sort all the way.  once this is done, insertSort
  1591.   should be called to finish the sorting process!
  1592.  
  1593.   Preconditions:
  1594.     photon memory initialized
  1595.     static map_s points to the photon map
  1596.     'left' is the index of the first photon
  1597.     'right' is the index of the last photon
  1598.     'd' is the dimension to sort on (X, Y, or Z)
  1599.  
  1600.   Postconditions:
  1601.     photons from 'left' to 'right' in map_s are MOSTLY sorted in
  1602.     ascending order on dimension d
  1603. ******************************************************************************/
  1604. static void quickSortRec(int left, int right, int d)
  1605. {
  1606.   int j,k;
  1607.   if(left<right)
  1608.   {
  1609.     swapPhotons(((left+right)>>1), left+1);
  1610.     if(PHOTON_AMF(map_s,left+1).Loc[d] > PHOTON_AMF(map_s,right).Loc[d])
  1611.       swapPhotons(left+1,right);
  1612.     if(PHOTON_AMF(map_s,left).Loc[d] > PHOTON_AMF(map_s,right).Loc[d])
  1613.       swapPhotons(left,right);
  1614.     if(PHOTON_AMF(map_s,left+1).Loc[d] > PHOTON_AMF(map_s,left).Loc[d])
  1615.       swapPhotons(left+1,left);
  1616.  
  1617.     j=left+1; k=right;
  1618.     while(j<=k)
  1619.     {
  1620.       for(j++; ((j<=right)&&(PHOTON_AMF(map_s,j).Loc[d]<PHOTON_AMF(map_s,left).Loc[d])); j++);
  1621.       for(k--; ((k>=left)&&(PHOTON_AMF(map_s,k).Loc[d]>PHOTON_AMF(map_s,left).Loc[d])); k--);
  1622.  
  1623.       if(j<k)
  1624.         swapPhotons(j,k);
  1625.     }
  1626.  
  1627.     swapPhotons(left,k);
  1628.     if(k-left > 10)
  1629.     {
  1630.       quickSortRec(left,k-1,d);
  1631.     }
  1632.     if(right-k > 10)
  1633.     {
  1634.       quickSortRec(k+1,right,d);
  1635.     }
  1636.     /* leave the rest for insertSort */
  1637.   }
  1638. }
  1639.  
  1640. /*****************************************************************************
  1641.  
  1642.   FUNCTION
  1643.  
  1644.   sortAndSubdivide
  1645.   
  1646.   Finds the dimension with the greates range, sorts the photons on that
  1647.   dimension.  Then it recurses on the left and right halves (keeping
  1648.   the median photon as a pivot).  This produces a balanced kd-tree.
  1649.  
  1650.   Preconditions:
  1651.     photon memory initialized
  1652.     static map_s points to the photon map
  1653.     'start' is the index of the first photon
  1654.     'end' is the index of the last photon
  1655.     'sorted' is the dimension that was last sorted (so we don't sort again)
  1656.  
  1657.   Postconditions:
  1658.     photons from 'start' to 'end' in map_s are in a valid kd-tree format
  1659. ******************************************************************************/
  1660. static void sortAndSubdivide(int start, int end, int sorted)
  1661. {
  1662.   int i,j;             /* counters */
  1663.   SNGL_VECT min,max;   /* min/max vectors for finding range */
  1664.   int DimToUse;        /* which dimesion has the greatest range */
  1665.   int mid;             /* index of median (middle) */
  1666.   int len;             /* length of the array we're sorting */
  1667.  
  1668.   if (end==start) 
  1669.   {
  1670.     PHOTON_AMF(map_s, start).info = 0;
  1671.     return;
  1672.   }
  1673.  
  1674.   if(end<start) return;
  1675.  
  1676.   /*
  1677.    * loop and find greatest range
  1678.    */
  1679.   Make_Vector(min, 1/EPSILON, 1/EPSILON, 1/EPSILON);
  1680.   Make_Vector(max, -1/EPSILON, -1/EPSILON, -1/EPSILON);
  1681.  
  1682.   for(i=start; i<=end; i++)
  1683.   {
  1684.     for(j=X; j<=Z; j++)
  1685.     {
  1686.       PHOTON *ph = &(PHOTON_AMF(map_s,i));
  1687.  
  1688.       if (ph->Loc[j] < min[j])
  1689.         min[j]=ph->Loc[j];
  1690.       if (ph->Loc[j] > max[j])
  1691.         max[j]=ph->Loc[j];
  1692.     }
  1693.   }
  1694.  
  1695.   /* choose which dimension to use */
  1696.   DimToUse = X;
  1697.   if((max[Y]-min[Y])>(max[DimToUse]-min[DimToUse]))
  1698.     DimToUse=Y;
  1699.   if((max[Z]-min[Z])>(max[DimToUse]-min[DimToUse]))
  1700.     DimToUse=Z;
  1701.  
  1702.   if(DimToUse != sorted)
  1703.   {
  1704.     /* start with quicksort */
  1705.     len = end-start;
  1706.     if (len>=8)
  1707.     {
  1708.       /* only display status every so often */
  1709.       if(len > 1000)
  1710. #ifndef ProgressOnOneLine        
  1711.         Status_Info ("\nSorting photons... %d",end);
  1712. #else
  1713.         Status_Info ("\rSorting photons... %8d",end);
  1714. #endif
  1715.  
  1716.       quickSortRec(start,end,DimToUse);
  1717.     }
  1718.  
  1719.     /* finish up with quicksort */
  1720.     insertSort(start,end,DimToUse);
  1721.   }
  1722.  
  1723.   /* find midpoint and set DimToUse for it */
  1724.   mid = (end+start)>>1;
  1725.   PHOTON_AMF(map_s, mid).info = DimToUse;
  1726.  
  1727.   /* now recurse to continue building the kd-tree */
  1728.   sortAndSubdivide(start, mid - 1, DimToUse);
  1729.   sortAndSubdivide(mid + 1, end, DimToUse);
  1730. }
  1731.  
  1732. /*****************************************************************************
  1733.  
  1734.   FUNCTION
  1735.  
  1736.   buildTree
  1737.   
  1738.   Builds the kd-tree by calling sortAndSubdivide().  Sets the static
  1739.   variable map_s.
  1740.  
  1741.   Preconditions:
  1742.     photon memory initialized
  1743.     'map' is a pointer to a photon map containing an array of unsorted
  1744.          photons
  1745.  
  1746.   Postconditions:
  1747.     photons are in a valid kd-tree format
  1748. ******************************************************************************/
  1749. static void buildTree(PHOTON_MAP *map)
  1750. {
  1751. #ifdef ProgressOnOneLine
  1752.   Status_Info ("\n");
  1753. #endif
  1754.   map_s = map->head;
  1755.   sortAndSubdivide(0,map->numPhotons-1,X+Y+Z/*this is not X, Y, or Z*/);
  1756. }
  1757.  
  1758. /*****************************************************************************
  1759.  
  1760.   FUNCTION
  1761.  
  1762.   setGatherOptions
  1763.   
  1764.   determines gather options
  1765.  
  1766.   Preconditions:
  1767.     photon memory initialized
  1768.     'map' points to an already-built (and sorted) photon map
  1769.     'mediaMap' is true if 'map' contians media photons, and false if
  1770.          'map' contains surface photons
  1771.  
  1772.   Postconditions:
  1773.     gather gather options are set for this map
  1774. ******************************************************************************/
  1775. static void setGatherOptions(PHOTON_MAP *map, int mediaMap)
  1776. {
  1777.   DBL r;
  1778.   DBL density;
  1779.   VECTOR Point;
  1780.   int numToSample;
  1781.   int n,i,j;
  1782.   DBL mind,maxd,avgd;
  1783.   DBL sum,sum2;
  1784.   DBL saveDensity;
  1785.   int greaterThan;
  1786.  
  1787.   /* if user did not set minimum gather radius, 
  1788.     then we should calculate it */
  1789.   if (map->minGatherRad <= 0.0)
  1790.   {
  1791.     mind=10000000.0;
  1792.     maxd=avgd=sum=sum2=0.0;
  1793.  
  1794.     /* use 5% of photons, min 100, max 10000 */
  1795.     numToSample = map->numPhotons/20;
  1796.     if (numToSample>1000) numToSample = 1000;
  1797.     if (numToSample<100) numToSample = 100;
  1798.  
  1799.     for(i=0; i<numToSample; i++)
  1800.     {
  1801.       j = rand() % map->numPhotons;
  1802.  
  1803.       Assign_SNGL_Vect(Point,(PHOTON_AMF(map->head, j)).Loc);
  1804.  
  1805.       n=gatherPhotons(Point, 10000000.0, &r, NULL, FALSE, map);
  1806.  
  1807.       if(mediaMap)
  1808.         density = 3.0 * n / (4.0*M_PI*r*r*r); /* should that be 4/3? */
  1809.       else
  1810.         density = n / (M_PI*r*r);
  1811.  
  1812.  
  1813.       if (density>maxd) maxd=density;
  1814.       if (density<mind) mind=density;
  1815.       sum+=density;
  1816.       sum2+=density*density;
  1817.     }
  1818.     avgd = sum/numToSample;
  1819.  
  1820.     /* try using some standard deviation stuff instead of this */
  1821.     saveDensity = avgd;
  1822.  
  1823.     greaterThan = 0;
  1824.     for(i=0; i<numToSample; i++)
  1825.     {
  1826.       j = rand() % map->numPhotons;
  1827.  
  1828.       Assign_SNGL_Vect(Point,(PHOTON_AMF(map->head, j)).Loc);
  1829.  
  1830.       n=gatherPhotons(Point, 10000000.0, &r, NULL, FALSE, map);
  1831.  
  1832.       if(mediaMap)
  1833.         density = 3.0 * n / (4.0*M_PI*r*r*r); /* should that be 4/3? */
  1834.       else
  1835.         density = n / (M_PI*r*r);
  1836.  
  1837.       if (density>saveDensity)
  1838.         greaterThan++;
  1839.     }
  1840.  
  1841.     density = saveDensity * (DBL)greaterThan / numToSample;
  1842.  
  1843.     if(mediaMap)
  1844.     {
  1845.       map->minGatherRad = pow(3.0 * photonOptions.maxGatherCount / (density*M_PI*4.0), 0.3333);
  1846.     }
  1847.     else
  1848.       map->minGatherRad = sqrt(photonOptions.maxGatherCount / (density*M_PI));
  1849.   }
  1850.  
  1851.   if(mediaMap)
  1852.   {
  1853.     /* double the radius if it is a media map */
  1854.     map->minGatherRad *= 2;
  1855.   }
  1856.  
  1857.   /* always do this! */
  1858.   map->gatherRadStep = map->minGatherRad*1.0;
  1859.  
  1860.   /* somehow we could automatically determine the number of steps */
  1861. }
  1862.  
  1863.  
  1864. /**************************************************************
  1865.  
  1866.   =========== PRIORITY QUEUES ===============
  1867.   Each priority stores its data in the static variables below (such as
  1868.   numfound_s) and in the global variables 
  1869.  
  1870.   Preconditions:
  1871.  
  1872.   static DBL size_sq_s; - maximum radius given squared
  1873.   static DBL Size_s;  - radius
  1874.   static DBL dmax_s;    - square of radius used so far
  1875.   static int TargetNum_s; - target number
  1876.   static DBL *pt_s;       - center point
  1877.   static numfound_s;        - number of photons in priority queue
  1878.  
  1879.   these must be allocated:
  1880.     photonOptions.photonGatherList - array of photons in priority queue
  1881.     photonOptions.photonDistances  - corresponding priorities(distances)
  1882.  
  1883.   *Each priority queue has the following functions:
  1884.  
  1885.   function PQInsert(PHOTON *photon, DBL d)
  1886.  
  1887.     Inserts 'photon' into the priority queue with priority (distance
  1888.     from target point) 'd'.
  1889.  
  1890.   void PQDelMax()
  1891.   
  1892.     Removes the photon with the greates distance (highest priority)
  1893.     from the queue.
  1894.  
  1895. ********************************************************************/
  1896.  
  1897. /* try different priority queues */
  1898.  
  1899. #define ORDERED   0
  1900. #define UNORDERED 1
  1901. #define HEAP      2
  1902.  
  1903. #define PRI_QUE HEAP
  1904.  
  1905. /* -------------- ordered list implementation ----------------- */
  1906. #if (PRI_QUE == ORDERED)
  1907. static void PQInsert(PHOTON *photon, DBL d)
  1908. {
  1909.   int i,j;
  1910.  
  1911.   Increase_Counter(stats[Priority_Queue_Insert]);
  1912.   /* save this in order, remove maximum, save new dmax_s */
  1913.  
  1914.   /* store in array and shift - assumption is that we don't have
  1915.      to shift often */
  1916.   for (i=0; photonOptions.photonDistances[i]<d && i<(numfound_s); i++);
  1917.   for (j=numfound_s; j>i; j--)
  1918.   {
  1919.     photonOptions.photonGatherList[j] = photonOptions.photonGatherList[j-1];
  1920.     photonOptions.photonDistances[j] = photonOptions.photonDistances[j-1];
  1921.   }
  1922.  
  1923.   numfound_s++;
  1924.   photonOptions.photonGatherList[i] = photon;
  1925.   photonOptions.photonDistances[i] = d;
  1926.   if (numfound_s==TargetNum_s)
  1927.     dmax_s=photonOptions.photonDistances[numfound_s-1];
  1928.  
  1929. }
  1930.  
  1931. static void PQDelMax()
  1932. {
  1933.   Increase_Counter(stats[Priority_Queue_Remove]);
  1934.   numfound_s--;
  1935. }
  1936. #endif
  1937.  
  1938. /* -------------- unordered list implementation ----------------- */
  1939. #if (PRI_QUE == UNORDERED)
  1940. static void PQInsert(PHOTON *photon, DBL d)
  1941. {
  1942.   Increase_Counter(stats[Priority_Queue_Insert]);
  1943.  
  1944.   photonOptions.photonGatherList[numfound_s] = photon;
  1945.   photonOptions.photonDistances[numfound_s] = d;
  1946.  
  1947.   if (d>dmax_s)
  1948.     dmax_s=d;
  1949.  
  1950.   numfound_s++;
  1951. }
  1952.  
  1953. static void PQDelMax()
  1954. {
  1955.   int i,max;
  1956.  
  1957.   Increase_Counter(stats[Priority_Queue_Remove]);
  1958.  
  1959.   max=0;
  1960.   /* find max */
  1961.   for(i=1; i<numfound_s; i++)
  1962.     if (photonOptions.photonDistances[i]>photonOptions.photonDistances[max]) max = i;
  1963.  
  1964.   /* remove it, shifting the photons */
  1965.   for(i=max+1; i<numfound_s; i++)
  1966.   {
  1967.     photonOptions.photonGatherList[i-1] = photonOptions.photonGatherList[i];
  1968.     photonOptions.photonDistances[i-1] = photonOptions.photonDistances[i];
  1969.   }
  1970.  
  1971.   numfound_s--;
  1972.  
  1973.   /* find a new dmax_s */
  1974.   dmax_s=photonOptions.photonDistances[0];
  1975.   for(i=1; i<numfound_s; i++)
  1976.     if (photonOptions.photonDistances[i]>dmax_s) dmax_s = photonOptions.photonDistances[i];
  1977. }
  1978. #endif
  1979.  
  1980. /* -------------- heap implementation ----------------- */
  1981. /* this is from Sejwick (spelling?) */
  1982. #if (PRI_QUE == HEAP)
  1983.  
  1984. static void fixUp()
  1985. {
  1986.   int k,k2,k3;
  1987.   DBL d;
  1988.   PHOTON *ph;
  1989.   k=numfound_s;
  1990.  
  1991.   while(k>1 && photonOptions.photonDistances[(k>>1)-1]<photonOptions.photonDistances[k-1])
  1992.   {
  1993.     /* exchange k & k/2  */
  1994.     k2=k-1;
  1995.     k3=(k>>1)-1;
  1996.     ph = photonOptions.photonGatherList[k2];
  1997.     d = photonOptions.photonDistances[k2];
  1998.     photonOptions.photonGatherList[k2] = photonOptions.photonGatherList[k3];
  1999.     photonOptions.photonDistances[k2] = photonOptions.photonDistances[k3];
  2000.     photonOptions.photonGatherList[k3] = ph;
  2001.     photonOptions.photonDistances[k3] = d;
  2002.  
  2003.     k=k>>1;
  2004.   }
  2005. }
  2006.  
  2007. static void fixDown()
  2008. {
  2009.   DBL d;
  2010.   PHOTON *ph;
  2011.   int j,j2;
  2012.   int k,k2;
  2013.   k=1;
  2014.  
  2015.   while (2*k <= numfound_s)
  2016.   {
  2017.     j=2*k;
  2018.     if (j<numfound_s && photonOptions.photonDistances[j-1]<photonOptions.photonDistances[j+1-1]) j++;
  2019.  
  2020.     k2=k-1;
  2021.     j2=j-1;
  2022.     if (!(photonOptions.photonDistances[k2]<photonOptions.photonDistances[j2])) break;
  2023.  
  2024.     /* exchange k & j  */
  2025.     ph = photonOptions.photonGatherList[k2];
  2026.     d = photonOptions.photonDistances[k2];
  2027.     photonOptions.photonGatherList[k2] = photonOptions.photonGatherList[j2];
  2028.     photonOptions.photonDistances[k2] = photonOptions.photonDistances[j2];
  2029.     photonOptions.photonGatherList[j2] = ph;
  2030.     photonOptions.photonDistances[j2] = d;
  2031.  
  2032.     k=j;
  2033.   }
  2034. }
  2035.  
  2036.  
  2037. static void PQInsert(PHOTON *photon, DBL d)
  2038. {
  2039.   Increase_Counter(stats[Priority_Queue_Insert]);
  2040.  
  2041.   numfound_s++;
  2042.  
  2043.   photonOptions.photonGatherList[numfound_s-1] = photon;
  2044.   photonOptions.photonDistances[numfound_s-1] = d;
  2045.  
  2046.   fixUp();
  2047.  
  2048.   if (d>dmax_s)
  2049.     dmax_s=d;
  2050. }
  2051.  
  2052. static void PQDelMax()
  2053. {
  2054.   DBL d;
  2055.   PHOTON *ph;
  2056.   int i;
  2057.  
  2058.   Increase_Counter(stats[Priority_Queue_Remove]);
  2059.  
  2060.   /* exchange 1 & numfound_s */
  2061.   i = numfound_s-1;
  2062.   ph = photonOptions.photonGatherList[1-1];
  2063.   d = photonOptions.photonDistances[1-1];
  2064.   photonOptions.photonGatherList[1-1] = photonOptions.photonGatherList[i];
  2065.   photonOptions.photonDistances[1-1] = photonOptions.photonDistances[i];
  2066.   photonOptions.photonGatherList[i] = ph;
  2067.   photonOptions.photonDistances[i] = d;
  2068.  
  2069.   numfound_s--;
  2070.   fixDown();
  2071.  
  2072.   /* find a new dmax_s */
  2073.   dmax_s=photonOptions.photonDistances[0];
  2074.   /*
  2075.   this is not needed. we know that the biggest is at the top of the heap
  2076.   for(i=1; i<numfound_s; i++)
  2077.     if (photonOptions.photonDistances[i]>dmax_s) dmax_s = photonOptions.photonDistances[i];
  2078.   */
  2079. }
  2080.  
  2081. #endif
  2082.  
  2083. /*****************************************************************************
  2084.  
  2085.   FUNCTION
  2086.  
  2087.   gatherPhotonsRec()
  2088.  
  2089.   Recursive part of gatherPhotons
  2090.   Searches the kd-tree with range start..end (midpoint is pivot)
  2091.   
  2092.   Preconditions:
  2093.     same preconditions as priority queue functions
  2094.     static variable map_s points to the map to use
  2095.     'start' is the first photon in this range
  2096.     'end' is the last photon in this range
  2097.  
  2098.     the range 'start..end' must have been used in building photon map!!!
  2099.  
  2100.   Postconditions:
  2101.     photons within the range of start..end are added to the priority
  2102.     queue (photons may be delted from the queue to make room for photons
  2103.     of lower priority)
  2104.  
  2105. ******************************************************************************/
  2106.  
  2107. static void gatherPhotonsRec(int start, int end)
  2108. {
  2109.   DBL sqrt_dmax_s, delta;
  2110.   int DimToUse;
  2111.   DBL d,dx,dy,dz;
  2112.   int mid;
  2113.   PHOTON *photon;
  2114.   VECTOR ptToPhoton;
  2115.   DBL discFix;   /* use disc(ellipsoid) for gathering instead of sphere */
  2116.  
  2117.   /* find midpoint */
  2118.   mid = (end+start)>>1;
  2119.   photon = &(PHOTON_AMF(map_s, mid));
  2120.  
  2121.   DimToUse = photon->info & PH_MASK_XYZ;
  2122.  
  2123.   /*
  2124.    * check this photon
  2125.    */
  2126.  
  2127.   /* find distance from pt */
  2128.   ptToPhoton[X] = - pt_s[X] + photon->Loc[X];
  2129.   ptToPhoton[Y] = - pt_s[Y] + photon->Loc[Y];
  2130.   ptToPhoton[Z] = - pt_s[Z] + photon->Loc[Z];
  2131.   /* all distances are squared */
  2132.   dx = ptToPhoton[X]*ptToPhoton[X];
  2133.   dy = ptToPhoton[Y]*ptToPhoton[Y];
  2134.   dz = ptToPhoton[Z]*ptToPhoton[Z];
  2135.  
  2136.   if (!(  ((dx>dmax_s) && ((DimToUse)==X)) ||
  2137.           ((dy>dmax_s) && ((DimToUse)==Y)) ||
  2138.           ((dz>dmax_s) && ((DimToUse)==Z)) ))
  2139.   {
  2140.     /* it fits manhatten distance - maybe we can use this photon */
  2141.  
  2142.     /* find euclidian distance (squared) */
  2143.     d = dx + dy + dz;
  2144.  
  2145.     /* now fix this distance so that we gather using an ellipsoid
  2146.        alligned with the surface normal instead of a sphere.  This
  2147.        minimizes false bleeding of photons at sharp corners
  2148.  
  2149.        dmax_s is square of radius of major axis
  2150.        dmax_s/16 is  "   "   "     " minor  "    (1/6 of major axis)
  2151.      */
  2152.     /*
  2153.     VDot(discFix,norm_s,ptToPhoton);
  2154.     discFix*=discFix*(dmax_s/1000.0-dmax_s);
  2155.     */
  2156.     if (flattenFactor!=0.0)
  2157.     {
  2158.       VDot(discFix,norm_s,ptToPhoton);
  2159.       d += flattenFactor*fabs(discFix);
  2160.     }
  2161.     /* this will add zero if on the plane, and will double distance from
  2162.     point to photon if it is ptToPhoton is perpendicular to the surface */
  2163.  
  2164.     if(d < dmax_s)
  2165.     {
  2166.       if (numfound_s+1>TargetNum_s)
  2167.       { PQDelMax(); }
  2168.  
  2169.       PQInsert(photon, d);
  2170.     }
  2171.   }
  2172.  
  2173.   /* now go left & right if appropriate - if going left or right goes out
  2174.       the current range, then don't go that way. */
  2175.   /*
  2176.   delta=pt_s[DimToUse]-photon->Loc[DimToUse];
  2177.   if(delta<0)
  2178.   {
  2179.     if (end>=mid+1) gatherPhotonsRec(start, mid - 1);
  2180.     if (delta*delta < dmax_s )
  2181.       if (mid-1>=start) gatherPhotonsRec(mid + 1, end);
  2182.   }
  2183.   else
  2184.   {
  2185.     if (mid-1>=start) gatherPhotonsRec(mid+1,end);
  2186.     if (delta*delta < dmax_s )
  2187.       if (end>=mid+1) gatherPhotonsRec(start, mid - 1);
  2188.   }
  2189.   */
  2190.   sqrt_dmax_s = sqrt(dmax_s);
  2191.   delta=pt_s[DimToUse]-photon->Loc[DimToUse];
  2192.   if(delta<0)
  2193.   {
  2194.     /* on left - go left first */
  2195.     if (pt_s[DimToUse]-sqrt_dmax_s < photon->Loc[DimToUse])
  2196.     {
  2197.       if (mid-1>=start)
  2198.         gatherPhotonsRec(start, mid - 1);
  2199.     }
  2200.     if (pt_s[DimToUse]+sqrt_dmax_s > photon->Loc[DimToUse])
  2201.     {
  2202.       if(end>=mid+1)
  2203.         gatherPhotonsRec(mid + 1, end);
  2204.     }
  2205.   }
  2206.   else
  2207.   {
  2208.     /* on right - go right first */
  2209.     if (pt_s[DimToUse]+sqrt_dmax_s > photon->Loc[DimToUse])
  2210.     {
  2211.       if(end>=mid+1)
  2212.         gatherPhotonsRec(mid + 1, end);
  2213.     }
  2214.     if (pt_s[DimToUse]-sqrt_dmax_s < photon->Loc[DimToUse])
  2215.     {
  2216.       if (mid-1>=start)
  2217.         gatherPhotonsRec(start, mid - 1);
  2218.     }
  2219.   }
  2220. }
  2221.  
  2222. /*****************************************************************************
  2223.  
  2224.   FUNCTION
  2225.  
  2226.   gatherPhotons()
  2227.  
  2228.   gathers photons from the global photon map
  2229.  
  2230.   Preconditons:
  2231.  
  2232.     photonOptions.photonGatherList and photonOptions.photonDistances
  2233.       are allocated and are each maxGatherCount in length
  2234.  
  2235.     'Size' - maximum search radius
  2236.     'r' points to a double
  2237.  
  2238.     BuildPhotonMaps() has been called for this scene.
  2239.  
  2240.   Postconditions:
  2241.  
  2242.     *r is radius actually used for gathereing (maximum value is 'Size')
  2243.     photonOptions.photonGatherList and photonOptions.photonDistances
  2244.       contain the gathered photons
  2245.     returns number of photons gathered
  2246.  
  2247. ******************************************************************************/
  2248. int gatherPhotons(VECTOR pt, DBL Size, DBL *r, VECTOR norm, int flatten, PHOTON_MAP *map)
  2249. {
  2250.   if (map->numPhotons<=0) return 0; /* no crashes, please... */
  2251.  
  2252.   /* set the static variables */
  2253.   numfound_s=0;
  2254.   size_sq_s = Size*Size;
  2255.   dmax_s = size_sq_s;
  2256.   norm_s = norm;
  2257.  
  2258.   if(flatten)
  2259.   {
  2260.     flattenFactor = 1.0;
  2261.   }
  2262.   else
  2263.   {
  2264.     flattenFactor = 0.0;
  2265.   }
  2266.  
  2267.   Size_s = Size;
  2268.   TargetNum_s = photonOptions.maxGatherCount;
  2269.   pt_s = pt;
  2270.  
  2271.   map_s = map->head;
  2272.  
  2273.   /* now search the kd-tree recursively */
  2274.   gatherPhotonsRec(0, map->numPhotons-1);
  2275.  
  2276.   /* set the radius variable */
  2277.   *r = sqrt(dmax_s);
  2278.  
  2279.   /* return the number of photons found */
  2280.   return(numfound_s);
  2281. }
  2282.  
  2283.  
  2284. /******************************************************************
  2285. stuff grabbed from radiosit.h & radiosit.c
  2286. ******************************************************************/
  2287. #ifndef MACOS
  2288. typedef struct byte_xyz BYTE_XYZ;
  2289.  
  2290. struct byte_xyz {
  2291.   unsigned char x, y, z;
  2292. };
  2293. #endif
  2294.  
  2295. extern BYTE_XYZ rad_samples[];
  2296.  
  2297. static void VUnpack(VECTOR dest_vec, BYTE_XYZ * pack_vec)
  2298. {
  2299.   dest_vec[X] = ((double)pack_vec->x * (1./ 255.))*2.-1.;
  2300.   dest_vec[Y] = ((double)pack_vec->y * (1./ 255.))*2.-1.;
  2301.   dest_vec[Z] = ((double)pack_vec->z * (1./ 255.));
  2302.  
  2303.   VNormalizeEq(dest_vec);   /* already good to about 1%, but we can do better */
  2304. }
  2305.  
  2306. /******************************************************************
  2307. ******************************************************************/
  2308. void ChooseRay(RAY *NewRay, VECTOR Normal, RAY *Ray, VECTOR Raw_Normal, int WhichRay)
  2309. {
  2310.   VECTOR random_vec, up, n2, n3;
  2311.   int i;
  2312.   DBL /*n,*/ NRay_Direction;
  2313.  
  2314. #define REFLECT_FOR_RADIANCE 0
  2315. #if (REFLECT_FOR_RADIANCE)
  2316.   /* Get direction of reflected ray. */
  2317.   n = -2.0 * (Ray->Direction[X] * Normal[X] + Ray->Direction[Y] * Normal[Y] + Ray->Direction[Z] * Normal[Z]);
  2318.  
  2319.   VLinComb2(NewRay->Direction, n, Normal, 1.0, Ray->Direction);
  2320.  
  2321.   VDot(NRay_Direction, NewRay->Direction, Raw_Normal);
  2322.   if (NRay_Direction < 0.0)
  2323.   {
  2324.     /* subtract 2*(projection of NRay.Direction onto Raw_Normal)
  2325.        from NRay.Direction */
  2326.     DBL Proj;
  2327.     Proj = NRay_Direction * -2;
  2328.     VAddScaledEq(NewRay->Direction, Proj, Raw_Normal);
  2329.   }
  2330.   return;
  2331. #else
  2332.   Assign_Vector(NewRay->Direction, Normal);
  2333. #endif
  2334.  
  2335.   if ( fabs(fabs(NewRay->Direction[Z])- 1.) < .1 ) {
  2336.     /* too close to vertical for comfort, so use cross product with horizon */
  2337.     up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
  2338.   }
  2339.   else
  2340.   {
  2341.     up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
  2342.   }
  2343.  
  2344.   VCross(n2, NewRay->Direction, up);  VNormalizeEq(n2);
  2345.   VCross(n3, NewRay->Direction, n2);  VNormalizeEq(n3);
  2346.  
  2347.   /*i = (int)(FRAND()*1600);*/
  2348.   i = WhichRay;
  2349.   WhichRay = (WhichRay + 1) % 1600;
  2350.  
  2351.   VUnpack(random_vec, &rad_samples[i]);
  2352.  
  2353.   if ( fabs(NewRay->Direction[Z] - 1.) < .001 )         /* pretty well straight Z, folks */
  2354.   {
  2355.     /* we are within 1/20 degree of pointing in the Z axis. */
  2356.     /* use all vectors as is--they're precomputed this way */
  2357.     Assign_Vector(NewRay->Direction, random_vec);
  2358.   }
  2359.   else
  2360.   {
  2361.     NewRay->Direction[X] = n2[X]*random_vec[X] + n3[X]*random_vec[Y] + NewRay->Direction[X]*random_vec[Z];
  2362.     NewRay->Direction[Y] = n2[Y]*random_vec[X] + n3[Y]*random_vec[Y] + NewRay->Direction[Y]*random_vec[Z];
  2363.     NewRay->Direction[Z] = n2[Z]*random_vec[X] + n3[Z]*random_vec[Y] + NewRay->Direction[Z]*random_vec[Z];
  2364.   }
  2365.  
  2366.   /* if our new ray goes through, flip it back across raw_normal */
  2367.  
  2368.   VDot(NRay_Direction, NewRay->Direction, Raw_Normal);
  2369.   if (NRay_Direction < 0.0)
  2370.   {
  2371.     /* subtract 2*(projection of NRay.Direction onto Raw_Normal)
  2372.        from NRay.Direction */
  2373.     DBL Proj;
  2374.     Proj = NRay_Direction * -2;
  2375.     VAddScaledEq(NewRay->Direction, Proj, Raw_Normal);
  2376.   }
  2377.  
  2378.   VNormalizeEq(NewRay->Direction);
  2379. }
  2380.  
  2381.