home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / media.c < prev    next >
C/C++ Source or Header  |  2001-01-23  |  69KB  |  2,802 lines
  1. /****************************************************************************
  2. *                   media.c
  3. *
  4. *  This module contains all functions for participating media.
  5. *
  6. *  from Persistence of Vision(tm) Ray Tracer
  7. *  Copyright 1996,1999 Persistence of Vision Team
  8. *---------------------------------------------------------------------------
  9. *  NOTICE: This source code file is provided so that users may experiment
  10. *  with enhancements to POV-Ray and to port the software to platforms other
  11. *  than those supported by the POV-Ray Team.  There are strict rules under
  12. *  which you are permitted to use this file.  The rules are in the file
  13. *  named POVLEGAL.DOC which should be distributed with this file.
  14. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  15. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  16. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  17. *
  18. * This program is based on the popular DKB raytracer version 2.12.
  19. * DKBTrace was originally written by David K. Buck.
  20. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  21. *
  22. *****************************************************************************/
  23.  
  24. #include "frame.h"
  25. #include "vector.h"
  26. #include "povproto.h"
  27. #include "chi2.h"
  28. #include "colour.h"
  29. #include "povray.h"
  30. #include "texture.h"
  31. #include "pigment.h"
  32. #include "objects.h"
  33. #include "lighting.h"
  34. #include "matrices.h"
  35. #include "media.h"
  36. #include "pattern.h"
  37. #include "povray.h"
  38. #include "point.h"
  39. #include "texture.h"
  40. #include "ray.h"
  41. #include "lightgrp.h"
  42.  
  43. #include "photons.h"
  44.  
  45. /*****************************************************************************
  46. * Local preprocessor defines
  47. ******************************************************************************/
  48.  
  49. /*****************************************************************************
  50. * Local typedefs
  51. ******************************************************************************/
  52. #ifdef UseMediaAndLightCache
  53. extern long MaxMediaCacheDepth;
  54. extern LIT_INTERVAL **MediaLitIntervalCache;
  55. extern LIGHT_LIST **MediaLightListCache;
  56. extern MEDIA_INTERVAL **MediaIntervalCache;
  57. extern long *MediaIntervalCacheSize;
  58. extern DBL *s0, *s1;
  59. extern long MediaCacheIndex;
  60. #ifdef AccumulateCacheStatistics
  61. extern long MaxSimMediatRecCntr;
  62. #endif
  63. #endif
  64.  
  65. /*****************************************************************************
  66. * Local variables
  67. ******************************************************************************/
  68. static int sampCount_s;
  69. static LIGHT_SOURCE *currLight;
  70.  
  71.  
  72. /*****************************************************************************
  73. * Static functions
  74. ******************************************************************************/
  75.  
  76. /* NK samples */
  77. static void sample_media_rec(LIGHT_LIST *Light_List, RAY *Ray, IMEDIA **Media_List, MEDIA_INTERVAL *Interval, int light_ray,
  78.                              DBL d1, DBL d3,
  79.                              COLOUR Result, COLOUR C1, COLOUR C3,
  80.                              COLOUR ODResult, COLOUR od1, COLOUR od3,
  81.                              int depth, DBL Jitter, DBL aa_threshold
  82. #ifdef EmissionMethodPatch
  83.                              , COLOUR Em1, COLOUR Em3, COLOUR EmR
  84. #endif
  85.                              );
  86. /* NK ---- */
  87.  
  88. static void sample_media (LIGHT_LIST *, RAY *, IMEDIA **, MEDIA_INTERVAL *, int, DBL, COLOUR, COLOUR, int
  89. #ifdef EmissionMethodPatch
  90.                           , COLOUR
  91. #endif
  92.                           );
  93.  
  94. static void get_light_list (LIGHT_LIST *, RAY *, INTERSECTION *);
  95. static void get_lit_interval (int *, LIT_INTERVAL *, int, LIGHT_LIST *, INTERSECTION *);
  96. static int set_up_sampling_intervals (MEDIA_INTERVAL *,
  97.   int, LIT_INTERVAL *, IMEDIA *);
  98.  
  99. static int intersect_spotlight (RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL *d2);
  100. static int intersect_cylinderlight (RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL *d2);
  101.  
  102. static int CDECL compdoubles (CONST void *in_a, CONST void *in_b);
  103.  
  104. static void evaluate_density_pattern (IMEDIA *, VECTOR , COLOUR);
  105.  
  106.  
  107. void Backtrace_Simulate_Media(IMEDIA **Media_List, RAY *Ray, INTERSECTION *Inter, COLOUR Colour)
  108. {
  109.   int i, j, intervals, use_extinction, use_scattering, ignore_photons;
  110.   int lit_interval_entries;
  111.   COLOUR Od;
  112.   LIGHT_LIST *Light_List = NULL;
  113.   LIT_INTERVAL *Lit_Interval;
  114.   IMEDIA *IMedia, **Tmp, *Local;
  115.   MEDIA_INTERVAL *Media_Interval, *curr;
  116.   int minSamples;
  117.   /* NK samples */
  118.   DBL d0;
  119.   COLOUR C0;
  120.   COLOUR od0;
  121.   COLOUR PhotonColour;
  122.   /* NK ---- */
  123.   VECTOR TempPoint;
  124.   DBL mediaSpacingFactor;
  125. #ifdef EmissionMethodPatch
  126.   COLOUR Unused;
  127. #endif
  128.  
  129.   /* Why are we here? */
  130.  
  131. #ifdef UseMediaAndLightCache
  132.     int thesize;    
  133.  
  134.   MediaCacheIndex++;
  135. #ifdef AccumulateCacheStatistics
  136.   MaxSimMediatRecCntr=max(MaxSimMediatRecCntr,MediaCacheIndex);
  137. #endif
  138.  
  139.   if ((Media_List == NULL) || (Media_List[0] == NULL))  
  140.   {
  141.     MediaCacheIndex--;
  142.     return;  
  143.   }
  144. #else
  145.   if ((Media_List == NULL) || (Media_List[0] == NULL))
  146.   {
  147.     return;
  148.   }
  149. #endif
  150.  
  151.   /* Find media with the largest number of intervals. */
  152.  
  153.   intervals = 0;
  154.  
  155.   use_extinction = FALSE;
  156.   use_scattering = FALSE;
  157.   ignore_photons = TRUE;
  158.  
  159.   IMedia = Media_List[0];
  160.  
  161.   for (Tmp = Media_List; (*Tmp) != NULL; Tmp++)
  162.   {
  163.     for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
  164.     {
  165.       if (Local->Intervals > intervals)
  166.       {
  167.         intervals = Local->Intervals;
  168.  
  169.         IMedia = Local;
  170.       }
  171.  
  172.       use_extinction |= Local->use_extinction;
  173.       use_scattering |= Local->use_scattering;
  174.       ignore_photons &= Local->ignore_photons;
  175.       if (!Check_Media_Light_Group(Local,photonOptions.Light))
  176.         ignore_photons = TRUE;
  177.     }
  178.   }
  179.  
  180.   /* If this is a light ray and no extinction is used we can return. */
  181.  
  182.   if (!use_extinction)
  183.   {
  184. #ifdef UseMediaAndLightCache
  185.     MediaCacheIndex--;
  186. #endif
  187.     return;
  188.   }
  189.  
  190.   /*
  191.    * Prepare the Monte Carlo integration along the ray from P0 to P1.
  192.    */
  193.  
  194.   /* light_ray always true */
  195. #ifdef UseMediaAndLightCache
  196.   if ( MediaCacheIndex >= MaxMediaCacheDepth )
  197.   {
  198.     ResizeMediaMallocCaches(MaxMediaCacheDepth*2);
  199.   }
  200.  
  201.   Lit_Interval=MediaLitIntervalCache[MediaCacheIndex];
  202. #else
  203.   Lit_Interval = (LIT_INTERVAL *)POV_MALLOC(sizeof(LIT_INTERVAL), "lit interval");
  204. #endif
  205.  
  206.   lit_interval_entries = 1;
  207.  
  208.   Lit_Interval[0].lit = FALSE;
  209.  
  210.   Lit_Interval[0].s0 = 0.0;
  211.   Lit_Interval[0].s1 =
  212.   Lit_Interval[0].ds = Inter->Depth;
  213.  
  214.   /* Set up sampling intervals. */
  215.  
  216.   /* NK samples - make sure we will always have enough intervals 
  217.      getting half-way through a really long render and then having it stop with a
  218.      "too few intervals" error message is just not acceptable.
  219.      And for adaptive subdivision, the fewer intervals, the better
  220.   
  221.      By the way, I don't know everything about this media code, so this may not
  222.      work perfectly, but I haven't had any problems yet, and it doesn't look like
  223.      it should cause any in the future. */
  224. #ifdef UseMediaAndLightCache
  225.   if (lit_interval_entries>IMedia->Intervals)
  226.     thesize=lit_interval_entries*sizeof(MEDIA_INTERVAL);
  227.   else        
  228.     thesize=IMedia->Intervals*sizeof(MEDIA_INTERVAL);
  229.  
  230.   if ( MediaCacheIndex >= MaxMediaCacheDepth)    
  231.   {
  232.     ResizeMediaMallocCaches(MaxMediaCacheDepth*2);
  233.   }
  234.  
  235.   if ( thesize <= MediaIntervalCacheSize[MediaCacheIndex])
  236.     Media_Interval=MediaIntervalCache[MediaCacheIndex];
  237.   else        
  238.   {
  239.     POV_FREE(MediaIntervalCache[MediaCacheIndex]);
  240.     Media_Interval=MediaIntervalCache[MediaCacheIndex] = (MEDIA_INTERVAL *)POV_MALLOC(thesize, "media intervals");
  241.     MediaIntervalCacheSize[MediaCacheIndex]=thesize;
  242.   }
  243. #else
  244.   if (lit_interval_entries>IMedia->Intervals)
  245.   {
  246.     Media_Interval = (MEDIA_INTERVAL *)POV_MALLOC(lit_interval_entries*sizeof(MEDIA_INTERVAL), "media intervals");
  247.   }
  248.   else
  249.   {
  250.     Media_Interval = (MEDIA_INTERVAL *)POV_MALLOC(IMedia->Intervals*sizeof(MEDIA_INTERVAL), "media intervals");
  251.   }
  252. #endif
  253.   intervals = set_up_sampling_intervals(Media_Interval, lit_interval_entries, Lit_Interval, IMedia);
  254.  
  255.   minSamples = IMedia->Min_Samples;
  256.   if (IMedia->Sample_Method == 3)
  257.     minSamples*=2;
  258.  
  259.   /* -------- should only use one interval! ----------- */
  260.  
  261.   /* Sample all intervals. */
  262.   Od[0] = Colour[0];
  263.   Od[1] = Colour[1];
  264.   Od[2] = Colour[2];
  265.  
  266.   for (i = 0; i < intervals; i++)
  267.   {
  268.     minSamples=(int)
  269.       (Media_Interval[i].ds /
  270.       (photonOptions.photonSpread *
  271.       photonOptions.photonDepth *
  272.       photonOptions.mediaSpacingFactor));
  273.  
  274.     if (minSamples<=photonOptions.maxMediaSteps)
  275.     {
  276.       /* all's well */
  277.       mediaSpacingFactor = photonOptions.mediaSpacingFactor;
  278.     }
  279.     else
  280.     {
  281.       /* too many steps - use fewer steps and a bigger spacing factor */
  282.       minSamples = photonOptions.maxMediaSteps;
  283.       mediaSpacingFactor = 
  284.         (Media_Interval[i].ds /
  285.         (photonOptions.photonSpread *
  286.         photonOptions.photonDepth *
  287.         minSamples));
  288.     }
  289.     /* Sample current interval. */
  290.  
  291.     Increase_Counter(stats[Media_Intervals]);
  292.  
  293.     /* set static sample count */
  294.     sampCount_s = minSamples;
  295.  
  296.     for (j = 0; j < minSamples; j++)
  297.     {
  298.       d0 = (j + 0.5 + FRAND()*photonOptions.jitter - 0.5*photonOptions.jitter) / minSamples;
  299.       sample_media(Light_List, Ray, Media_List, &Media_Interval[i], TRUE, d0, C0, od0, 2 /* use method 2 */
  300. #ifdef EmissionMethodPatch
  301.         , Unused
  302. #endif
  303.         );
  304.  
  305.       if (use_scattering && !ignore_photons)
  306.       {
  307.         VScale(PhotonColour,Od, mediaSpacingFactor * 
  308.                                 photonOptions.photonSpread *
  309.                                 (photonOptions.photonDepth+d0*Media_Interval[i].ds+Media_Interval[i].s0));
  310.  
  311.         Od[0] = Colour[0]*exp(-Media_Interval[i].od[0]/(minSamples*2));
  312.         Od[1] = Colour[1]*exp(-Media_Interval[i].od[1]/(minSamples*2));
  313.         Od[2] = Colour[2]*exp(-Media_Interval[i].od[2]/(minSamples*2));
  314.  
  315.         VEvaluateRay(TempPoint, Ray->Initial, d0*Media_Interval[i].ds+Media_Interval[i].s0, Ray->Direction);
  316.         
  317.         addMediaPhoton(TempPoint, Ray->Initial, PhotonColour, d0*Media_Interval[i].ds+Media_Interval[i].s0);
  318.       }
  319.     }
  320.   }
  321.  
  322.   /* Sum the influences of all intervals. */
  323.   Make_Colour(Od, 0.0, 0.0, 0.0);
  324.  
  325.   curr = &Media_Interval[0];
  326.  
  327.   for (i = 0; i < intervals; i++)
  328.   {
  329.     /* Add optical depth of current interval. */
  330.  
  331.     Od[0] += curr->od[0] / (DBL)curr->samples;
  332.     Od[1] += curr->od[1] / (DBL)curr->samples;
  333.     Od[2] += curr->od[2] / (DBL)curr->samples;
  334.  
  335.     curr++;
  336.   }
  337.  
  338.   /* Add contribution estimated for the participating media. */
  339.   Colour[0] = Colour[0] * exp(-Od[0]);
  340.   Colour[1] = Colour[1] * exp(-Od[1]);
  341.   Colour[2] = Colour[2] * exp(-Od[2]);
  342.  
  343. #ifdef UseMediaAndLightCache
  344.     MediaCacheIndex--;
  345. #else
  346.   POV_FREE(Lit_Interval);
  347.   POV_FREE(Media_Interval);
  348. #endif
  349. }
  350. /*****************************************************************************
  351. *
  352. * FUNCTION
  353. *
  354. *   Simulate_Media
  355. *
  356. * INPUT
  357. *
  358. *   Ray       - Current ray, start point P0
  359. *   Inter     - Current intersection, end point P1
  360. *   Colour    - Color emitted at P1 towards P0
  361. *   light_ray - TRUE if we are looking at a light source ray
  362. *
  363. * OUTPUT
  364. *
  365. *   Colour    - Color arriving at the end point
  366. *
  367. * RETURNS
  368. *
  369. * AUTHOR
  370. *
  371. *   Dieter Bayer
  372. *
  373. * DESCRIPTION
  374. *
  375. *   Simulate participating media using volume sampling.
  376. *
  377. *   The effects of participating media on the light emitted at P1
  378. *   towards P0 are determined using Monte Carlo integration.
  379. *
  380. *   The effects include: emission, absoprtion and scattering.
  381. *
  382. *   Currently one global medium with constant coefficients is implemented.
  383. *
  384. *   Ideas for the atmospheric scattering were taken from:
  385. *
  386. *     - M. Inakage, "An Illumination Model for Atmospheric Environments", ..
  387. *
  388. *     - Nishita, T., Miyawaki, Y. and Nakamae, E., "A Shading Model for
  389. *       Atmospheric Scattering Considering Luminous Intensity Distribution
  390. *       of Light Sources", Computer Graphics, 21, 4, 1987, 303-310
  391. *
  392. * CHANGES
  393. *
  394. *   Nov 1994 : Creation.
  395. *
  396. *   Jan 1995 : Added support of cylindrical light sources. [DB]
  397. *
  398. *   Jun 1995 : Added code for alpha channel support. [DB]
  399. *
  400. ******************************************************************************/
  401.  
  402. void Simulate_Media(IMEDIA **Media_List, RAY *Ray, INTERSECTION *Inter, COLOUR Colour, int light_ray)
  403. {
  404.   int i, j, intervals, use_extinction;
  405.   int lit_interval_entries;
  406.   DBL n;
  407.   COLOUR Od, Te, Va;
  408.   LIGHT_LIST *Light_List = NULL;
  409.   LIT_INTERVAL *Lit_Interval;
  410.   IMEDIA *IMedia, **Tmp, *Local;
  411.   MEDIA_INTERVAL *Media_Interval, *curr;
  412.   /* NK samples */
  413.   DBL d0, dd;
  414.   COLOUR C0, C1, Result;
  415.   COLOUR ODResult;
  416.   COLOUR od0,od1;
  417.   /* NK ---- */
  418.   /* NK fast light_ray media calculation for constant media */
  419.   int all_constant_and_light_ray;  /* is all the media constant? */
  420.   /* NK ---- */
  421. #ifdef SampleSpacingPatch
  422.   DBL length;            /* JB */
  423. #endif
  424. #ifdef EmissionMethodPatch
  425.   COLOUR Em1, Em2, EmR;        /* JB */
  426. #endif
  427.   int minSamples;        /* JB */
  428.  
  429.   /* Why are we here? */
  430.  
  431. #ifdef UseMediaAndLightCache
  432.     int thesize;    
  433.  
  434.   MediaCacheIndex++;
  435. #ifdef AccumulateCacheStatistics
  436.   MaxSimMediatRecCntr=max(MaxSimMediatRecCntr,MediaCacheIndex);
  437. #endif
  438.  
  439.   if ((Media_List == NULL) || (Media_List[0] == NULL))  
  440.   {
  441.     MediaCacheIndex--;
  442.     return;  
  443.   }
  444. #else
  445.   if ((Media_List == NULL) || (Media_List[0] == NULL))
  446.   {
  447.     return;
  448.   }
  449. #endif
  450.   /* Find media with the largest number of intervals. */
  451.  
  452.   intervals = 0;
  453.  
  454.   use_extinction = FALSE;
  455.  
  456.   /* NK fast light_ray media calculation for constant media */
  457.   all_constant_and_light_ray = light_ray;
  458.   /* NK ---- */
  459.  
  460.   IMedia = Media_List[0];
  461.  
  462.   for (Tmp = Media_List; (*Tmp) != NULL; Tmp++)
  463.   {
  464.     for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
  465.     {
  466.       if (Local->Intervals > intervals)
  467.       {
  468.         intervals = Local->Intervals;
  469.  
  470.         IMedia = Local;
  471.       }
  472.  
  473.       use_extinction |= Local->use_extinction;
  474.  
  475.       /* NK fast light_ray media calculation for constant media */
  476.       if(Local->Density)
  477.       {
  478.         all_constant_and_light_ray &= (Local->Density->Type == PLAIN_PATTERN);
  479.       }
  480.       /* NK ---- */
  481.     }
  482.   }
  483.  
  484.   /* If this is a light ray and no extinction is used we can return. */
  485.  
  486.   if ((light_ray) && (!use_extinction))
  487.   {
  488. #ifdef UseMediaAndLightCache
  489.     MediaCacheIndex--;
  490. #endif
  491.     return;
  492.   }
  493.  
  494.   /* NK fast light_ray media calculation for constant media */
  495.   if (all_constant_and_light_ray)
  496.   {
  497.     intervals = 1;
  498.   }
  499.   /* NK ---- */
  500.  
  501.   /*
  502.    * Prepare the Monte Carlo integration along the ray from P0 to P1.
  503.    */
  504.  
  505.   if (light_ray || (Frame.Number_Of_Light_Sources==0))
  506.   {
  507. #ifdef UseMediaAndLightCache
  508.     if ( MediaCacheIndex >= MaxMediaCacheDepth )
  509.     {
  510.       ResizeMediaMallocCaches(MaxMediaCacheDepth*2);
  511.     }
  512.  
  513.     Lit_Interval=MediaLitIntervalCache[MediaCacheIndex];
  514. #else
  515.     Lit_Interval = (LIT_INTERVAL *)POV_MALLOC(sizeof(LIT_INTERVAL), "lit interval");
  516. #endif
  517.  
  518.     lit_interval_entries = 1;
  519.  
  520.     Lit_Interval[0].lit = FALSE;
  521.  
  522.     Lit_Interval[0].s0 = 0.0;
  523.     Lit_Interval[0].s1 =
  524.     Lit_Interval[0].ds = Inter->Depth;
  525.   }
  526.   else
  527.   {
  528.     /* Get light list. */
  529.  
  530. #ifdef UseMediaAndLightCache
  531.       if ( MediaCacheIndex >= MaxMediaCacheDepth)    
  532.     {
  533.       ResizeMediaMallocCaches(MaxMediaCacheDepth*2);
  534.     }
  535.       Light_List = MediaLightListCache[MediaCacheIndex];
  536. #else    
  537.       Light_List = (LIGHT_LIST *)POV_MALLOC(Frame.Number_Of_Light_Sources*sizeof(LIGHT_LIST), "light list");
  538. #endif
  539.     get_light_list(Light_List, Ray, Inter);
  540.  
  541.     /* Get lit intervals. */
  542. #ifdef UseMediaAndLightCache
  543.     /* already resized */
  544.         Lit_Interval=MediaLitIntervalCache[MediaCacheIndex];
  545. #else  
  546.     Lit_Interval = (LIT_INTERVAL *)POV_MALLOC((2*Frame.Number_Of_Light_Sources+1)*sizeof(LIT_INTERVAL), "lit interval");
  547. #endif
  548.     get_lit_interval(&lit_interval_entries, Lit_Interval, Frame.Number_Of_Light_Sources, Light_List, Inter);
  549.   }
  550.  
  551.   /* Set up sampling intervals. */
  552.  
  553.   /* NK samples - make sure we will always have enough intervals 
  554.      getting half-way through a really long render and then having it stop with a
  555.      "too few intervals" error message is just not acceptable.
  556.      And for adaptive subdivision, the fewer intervals, the better
  557.   
  558.      By the way, I don't know everything about this media code, so this may not
  559.      work perfectly, but I haven't had any problems yet, and it doesn't look like
  560.      it should cause any in the future. */
  561. #ifdef UseMediaAndLightCache
  562.   if (lit_interval_entries>IMedia->Intervals)
  563.     thesize=lit_interval_entries*sizeof(MEDIA_INTERVAL);
  564.   else        
  565.     thesize=IMedia->Intervals*sizeof(MEDIA_INTERVAL);
  566.  
  567.   if ( MediaCacheIndex >= MaxMediaCacheDepth)    
  568.   {
  569.     ResizeMediaMallocCaches(MaxMediaCacheDepth*2);
  570.   }
  571.  
  572.   if ( thesize <= MediaIntervalCacheSize[MediaCacheIndex])
  573.     Media_Interval=MediaIntervalCache[MediaCacheIndex];
  574.   else        
  575.   {
  576.     POV_FREE(MediaIntervalCache[MediaCacheIndex]);
  577.     Media_Interval=MediaIntervalCache[MediaCacheIndex] = (MEDIA_INTERVAL *)POV_MALLOC(thesize, "media intervals");
  578.     MediaIntervalCacheSize[MediaCacheIndex]=thesize;
  579.   }
  580. #else
  581.   if (lit_interval_entries>IMedia->Intervals)
  582.   {
  583.     Media_Interval = (MEDIA_INTERVAL *)POV_MALLOC(lit_interval_entries*sizeof(MEDIA_INTERVAL), "media intervals");
  584.   }
  585.   else
  586.   {
  587.     Media_Interval = (MEDIA_INTERVAL *)POV_MALLOC(IMedia->Intervals*sizeof(MEDIA_INTERVAL), "media intervals");
  588.   }
  589. #endif
  590.   intervals = set_up_sampling_intervals(Media_Interval, lit_interval_entries, Lit_Interval, IMedia);
  591.  
  592. #ifdef SampleSpacingPatch
  593.   /* JB +++ */
  594.   if (IMedia->Sample_Spacing > 0) {
  595.     VDist (length, Ray->Initial, Inter->IPoint);
  596.     /*    fprintf (stderr, "P0: <%f, %f, %f>\n P1: <%f, %f, %f>\n",
  597.       Ray->Initial[0], Ray->Initial[1], Ray->Initial[2],
  598.       Inter->IPoint[0], Inter->IPoint[1], Inter->IPoint[2]);
  599.       fprintf (stderr, "%f\n", Ray->Direction[2]);*/
  600.     minSamples = (int)(length/IMedia->Sample_Spacing+(FRAND()-0.5)*IMedia->Sample_Jitter)+1;
  601.     /*    fprintf (stderr, "length: %f, minSamples: %d, Therefore: ", length, minSamples);*/
  602.     if (minSamples < IMedia->Min_Samples) {
  603.       /*if (!light_ray)
  604.     fprintf (stderr, "Ideal min. samples: %d, actual min. samples: %d.\n", minSamples, IMedia->Min_Samples);*/
  605.       minSamples = IMedia->Min_Samples;
  606.     }
  607.     if (minSamples > IMedia->Max_Samples) {
  608.       /*fprintf (stderr, "Ideal min. samples: %d, actual min. samples: %d.\n", minSamples, IMedia->Max_Samples);*/
  609.       minSamples = IMedia->Max_Samples;
  610.     }
  611.     /*    fprintf(stderr, "%6f, %d\n", length, minSamples);*/
  612.   } 
  613.   else
  614.                  /* JB --- */
  615. #endif
  616.   minSamples = IMedia->Min_Samples;
  617.  
  618.   /* Sample all intervals. */
  619.  
  620.   if ((IMedia->Sample_Method == 3) && !all_constant_and_light_ray)
  621.   {
  622.     /* adaptive sampling */
  623.     int sampleCount; /* internal count for samples to take */
  624.     IMEDIA *AA_Search_Media=*Media_List;
  625.     DBL aa_threshold=1000;
  626.  
  627.     /* find smallest AA_Threshold */
  628.     while (AA_Search_Media) {
  629.     if (AA_Search_Media->AA_Threshold < aa_threshold)
  630.       aa_threshold=AA_Search_Media->AA_Threshold;
  631.       AA_Search_Media=AA_Search_Media->Next_Media;
  632.     }
  633.  
  634.     for (i = 0; i < intervals; i++)
  635.     {
  636.       /* Sample current interval. */
  637.  
  638.       Increase_Counter(stats[Media_Intervals]);
  639.  
  640.       sampleCount = (int)((minSamples) / 2.0);
  641.  
  642.       if (sampleCount < 2)
  643.         sampleCount = 1;
  644.  
  645.       /* set static sample count */
  646.       sampCount_s = sampleCount;
  647.  
  648.       if (sampleCount < 2)
  649.       {
  650.         /* always do at least three samples - one on each end and one in the middle */
  651.  
  652.         /* don't re-sample this if we've already done it in the previous interval */
  653.         sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray, 0.0+IMedia->Jitter*(FRAND()-0.5), C0, od0, 3
  654. #ifdef EmissionMethodPatch
  655.           , Em1
  656. #endif
  657.           );
  658.         sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray, 1.0+IMedia->Jitter*(FRAND()-0.5), C1, od1, 3
  659. #ifdef EmissionMethodPatch
  660.           , Em2
  661. #endif
  662.           );
  663.         sample_media_rec(Light_List, Ray, Media_List, &Media_Interval[i], light_ray,
  664.                          0.0, 1.0, Media_Interval[i].te, C0, C1,
  665.                          Media_Interval[i].od, od0, od1,
  666.                          IMedia->AA_Level-1, IMedia->Jitter,aa_threshold
  667. #ifdef EmissionMethodPatch
  668.                          , Em1, Em2, Media_Interval[i].Em2
  669. #endif
  670.                          );
  671.         Media_Interval[i].samples = 1;
  672.  
  673.         /* move c1 to c0 to go on to next sample/interval */
  674.         C0[0]=C1[0];
  675.         C0[1]=C1[1];
  676.         C0[2]=C1[2];
  677.         /* move od1 to od0 to go on to the next sample/interval */
  678.         od0[0]=od1[0];
  679.         od0[1]=od1[1];
  680.         od0[2]=od1[2];
  681.  
  682. #ifdef EmissionMethodPatch
  683.         Em1[0]=Em2[0];
  684.         Em1[1]=Em2[1];
  685.         Em1[2]=Em2[2];
  686. #endif
  687.       }
  688.       else
  689.       {
  690.         dd = 1.0 / (sampleCount+1);
  691.         d0 = 0.0;
  692.  
  693.           sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray,
  694.                      d0+dd*IMedia->Jitter*(FRAND()-0.5), C0, od0, 3
  695. #ifdef EmissionMethodPatch
  696.                      , Em1
  697. #endif
  698.           );
  699.  
  700.         /* clear out od & te */
  701.         Media_Interval[i].te[0] =
  702.         Media_Interval[i].te[1] =
  703.         Media_Interval[i].te[2] = 0.;
  704.  
  705.         Media_Interval[i].od[0] =
  706.         Media_Interval[i].od[1] =
  707.         Media_Interval[i].od[2] = 0.;
  708.  
  709. #ifdef EmissionMethodPatch
  710.         Media_Interval[i].Em2[0] =
  711.         Media_Interval[i].Em2[1] =
  712.         Media_Interval[i].Em2[2] = 0.;
  713. #endif
  714.  
  715.         for (j = 1, d0+=dd; j <= sampleCount; j++, d0+=dd)
  716.         {
  717.           sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray,
  718.                        d0+dd*IMedia->Jitter*(FRAND()-0.5), C1, od1, 3
  719. #ifdef EmissionMethodPatch
  720.                        , Em2
  721. #endif
  722.                        );
  723.           sample_media_rec(Light_List, Ray, Media_List, &Media_Interval[i], light_ray,
  724.                            d0-dd, d0, Result, C0, C1,
  725.                            ODResult, od0, od1,
  726.                            IMedia->AA_Level-1, IMedia->Jitter,aa_threshold
  727. #ifdef EmissionMethodPatch
  728.                            , Em1, Em2, EmR
  729. #endif
  730.                            );
  731.  
  732.           /* keep a sum of the results */
  733.           /* do some attenuation, too, since we are doing samples in order */
  734.           Media_Interval[i].te[0]+=Result[0] * exp(-Media_Interval[i].od[0]/sampleCount);  /* should this be sampleCount+1? */
  735.           Media_Interval[i].te[1]+=Result[1] * exp(-Media_Interval[i].od[1]/sampleCount);
  736.           Media_Interval[i].te[2]+=Result[2] * exp(-Media_Interval[i].od[2]/sampleCount);
  737.           /* move c1 to c0 to go on to next sample/interval */
  738.           C0[0]=C1[0];
  739.           C0[1]=C1[1];
  740.           C0[2]=C1[2];
  741.  
  742.           /* now do the same for optical depth */
  743.           Media_Interval[i].od[0]+=ODResult[0];
  744.           Media_Interval[i].od[1]+=ODResult[1];
  745.           Media_Interval[i].od[2]+=ODResult[2];
  746.  
  747. #ifdef EmissionMethodPatch
  748.           Media_Interval[i].Em2[0] += EmR[0];
  749.           Media_Interval[i].Em2[1] += EmR[1];
  750.           Media_Interval[i].Em2[2] += EmR[2];
  751. #endif
  752.           /* move od1 to od0 to go on to the next sample/interval */
  753.           od0[0]=od1[0];
  754.           od0[1]=od1[1];
  755.           od0[2]=od1[2];
  756.  
  757. #ifdef EmissionMethodPatch
  758.           Em1[0]=Em2[0];
  759.           Em1[1]=Em2[1];
  760.           Em1[2]=Em2[2];
  761. #endif
  762.         }
  763.         Media_Interval[i].samples = sampleCount;
  764.       }
  765.     }
  766.   }
  767.   else
  768.   {
  769.  
  770.   for (i = 0; i < intervals; i++)
  771.   {
  772.     /* Sample current interval. */
  773.  
  774.     Increase_Counter(stats[Media_Intervals]);
  775.  
  776.     /* set static sample count */
  777.     sampCount_s = minSamples;
  778.  
  779.     for (j = 0; j < minSamples; j++)
  780.     {
  781.       if(IMedia->Sample_Method==2)
  782.       {
  783.         d0 = (j+0.5) / minSamples +
  784.             (FRAND() * IMedia->Jitter / minSamples);
  785. #ifdef EmissionMethodPatch
  786.         sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray, d0, C0, od0, 2, Em2);
  787. #else
  788.         sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray, d0, C0, od0, 2);
  789. #endif
  790.       }
  791.       else
  792.       {
  793.         /* we may get here with media method 3 */
  794.         d0 = FRAND();
  795. #ifdef EmissionMethodPatch
  796.         sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray, d0, C0, od0, 1, Em2);
  797. #else
  798.         sample_media(Light_List, Ray, Media_List, &Media_Interval[i], light_ray, d0, C0, od0, 1);
  799. #endif
  800.       }
  801.  
  802.       if (all_constant_and_light_ray)
  803.         j = minSamples;
  804.     }
  805.   }
  806.  
  807.   /* Cast additional samples if necessary. */
  808.  
  809.   if ((!light_ray) && (IMedia->Max_Samples > minSamples))
  810.   {
  811.     curr = &Media_Interval[0];
  812.  
  813.     for (i = 0; i < intervals; i++)
  814.     {
  815.       if (curr->samples < IMedia->Max_Samples)
  816.       {
  817.         /* Get variance of samples. */
  818.  
  819.         n = (DBL)curr->samples;
  820.  
  821.         Va[0] = (curr->te2[0] / n - Sqr(curr->te[0] / n)) / n;
  822.         Va[1] = (curr->te2[1] / n - Sqr(curr->te[1] / n)) / n;
  823.         Va[2] = (curr->te2[2] / n - Sqr(curr->te[2] / n)) / n;
  824.  
  825.         /* Take additional samples until variance is small enough. */
  826.  
  827.         while ((Va[0] >= IMedia->Sample_Threshold[curr->samples-1]) ||
  828.                (Va[1] >= IMedia->Sample_Threshold[curr->samples-1]) ||
  829.                (Va[2] >= IMedia->Sample_Threshold[curr->samples-1]))
  830.         {
  831.           /* Sample current interval again. */
  832.  
  833.           sample_media(Light_List, Ray, Media_List, curr, light_ray, FRAND(), C0, od0, 1 /* treat all like samp meth 1 */
  834. #ifdef EmissionMethodPatch
  835.             , Em2
  836. #endif
  837.             );
  838.  
  839.           /* Have we reached maximum number of samples. */
  840.  
  841.           if (curr->samples > IMedia->Max_Samples)
  842.           {
  843.             break;
  844.           }
  845.  
  846.           /* Get variance of samples. */
  847.  
  848.           n = (DBL)curr->samples;
  849.  
  850.           Va[0] = (curr->te2[0] / n - Sqr(curr->te[0] / n)) / n;
  851.           Va[1] = (curr->te2[1] / n - Sqr(curr->te[1] / n)) / n;
  852.           Va[2] = (curr->te2[2] / n - Sqr(curr->te[2] / n)) / n;
  853.         }
  854.       }
  855.  
  856.       curr++;
  857.     }
  858.   }
  859.  
  860.   } /* end of if for adaptive sampling else */
  861.  
  862.   /* Sum the influences of all intervals. */
  863.  
  864.   Make_Colour(Od, 0.0, 0.0, 0.0);
  865.   Make_Colour(Te, 0.0, 0.0, 0.0);
  866. #ifdef EmissionMethodPatch
  867.   Make_Colour(Em2, 0.0, 0.0, 0.0); /* JB */
  868. #endif
  869.  
  870.   curr = &Media_Interval[0];
  871.  
  872.   for (i = 0; i < intervals; i++)
  873.   {
  874.     /* Add total emission. */
  875.  
  876.     Te[0] += curr->te[0] / (DBL)curr->samples * exp(-Od[0]);
  877.     Te[1] += curr->te[1] / (DBL)curr->samples * exp(-Od[1]);
  878.     Te[2] += curr->te[2] / (DBL)curr->samples * exp(-Od[2]);
  879.  
  880.     /* Add optical depth of current interval. */
  881.  
  882.     Od[0] += curr->od[0] / (DBL)curr->samples;
  883.     Od[1] += curr->od[1] / (DBL)curr->samples;
  884.     Od[2] += curr->od[2] / (DBL)curr->samples;
  885.  
  886. #ifdef EmissionMethodPatch
  887.                 /* JB +++ */
  888.     Em2[0] += curr->Em2[0] / (DBL)curr->samples;
  889.     Em2[1] += curr->Em2[1] / (DBL)curr->samples;
  890.     Em2[2] += curr->Em2[2] / (DBL)curr->samples;
  891.     /*if ((curr->Em2[0]>0.5)||(curr->Em2[1]>0.5)||(curr->Em2[2]>0.5))
  892.      fprintf(stderr, "<%6f, %6f, %6f>\n", curr->Em2[0], curr->Em2[1], curr->Em2[2]); */
  893.                 /* JB --- */
  894. #endif
  895.  
  896.     curr++;
  897.   }
  898.  
  899.   /* Add contribution estimated for the participating media. */
  900.  
  901. #ifdef EmissionMethodPatch
  902.   Colour[0] = Colour[0] * exp(-Od[0]) + Te[0] + Em2[0];
  903.   Colour[1] = Colour[1] * exp(-Od[1]) + Te[1] + Em2[1];
  904.   Colour[2] = Colour[2] * exp(-Od[2]) + Te[2] + Em2[2];
  905. #else
  906.   Colour[0] = Colour[0] * exp(-Od[0]) + Te[0];
  907.   Colour[1] = Colour[1] * exp(-Od[1]) + Te[1];
  908.   Colour[2] = Colour[2] * exp(-Od[2]) + Te[2];
  909. #endif
  910.  
  911. #ifdef UseMediaAndLightCache
  912.     MediaCacheIndex--;
  913. #else
  914.   if (!(light_ray || (Frame.Number_Of_Light_Sources==0)))
  915.   {
  916.     POV_FREE(Light_List);
  917.   }
  918.  
  919.   POV_FREE(Lit_Interval);
  920.  
  921.   POV_FREE(Media_Interval);
  922. #endif
  923. }
  924.  
  925. static void sample_media_rec(LIGHT_LIST *Light_List, RAY *Ray, IMEDIA **Media_List, 
  926.                              MEDIA_INTERVAL *Interval, int light_ray, 
  927.                              DBL d1, DBL d3, 
  928.                              COLOUR Result, COLOUR C1, COLOUR C3, 
  929.                              COLOUR ODResult, COLOUR od1, COLOUR od3, 
  930.                              int depth, DBL Jitter, DBL aa_threshold
  931. #ifdef EmissionMethodPatch
  932.                              , COLOUR Em1, COLOUR Em3, COLOUR EmR
  933. #endif
  934.                              )
  935. {
  936.   COLOUR C2, Result2;
  937. #ifdef EmissionMethodPatch
  938.   COLOUR Em2, EmR2;
  939. #endif
  940.   COLOUR od2, ODResult2;
  941.   DBL d2, jdist;
  942.  
  943.   /* d2 is between d1 and d3 (all in range of 0..1 */
  944.   d2 = 0.5 * (d1 + d3);
  945.   jdist = d2 + Jitter * (d3 - d1) * (FRAND() - 0.5);
  946.  
  947.   sample_media(Light_List, Ray, Media_List, Interval, light_ray, 
  948.                jdist, C2, od2, 3/* sample method*/
  949. #ifdef EmissionMethodPatch
  950.                , Em2 
  951. #endif
  952.                );
  953.  
  954.   /* if we're at max depth, then let's just use this last sample and average
  955.      it with the two end points */
  956.   if (depth<=0)
  957.   {
  958.     /* average colors */
  959.     Result[0]=(C1[0]+C3[0]+C2[0])/3.;
  960.     Result[1]=(C1[1]+C3[1]+C2[1])/3.;
  961.     Result[2]=(C1[2]+C3[2]+C2[2])/3.;
  962.  
  963.     /* average the optical depth */
  964.     ODResult[0]=(od1[0]+od3[0]+od2[0])/3.;
  965.     ODResult[1]=(od1[1]+od3[1]+od2[1])/3.;
  966.     ODResult[2]=(od1[2]+od3[2]+od2[2])/3.;
  967.  
  968. #ifdef EmissionMethodPatch
  969.     EmR[0] = (Em1[0]+Em2[0]+Em3[0])/3;
  970.     EmR[1] = (Em1[1]+Em2[1]+Em3[1])/3;
  971.     EmR[2] = (Em1[2]+Em2[2]+Em3[2])/3;
  972. #endif
  973.     /* bail out - we're done now */
  974.     return;
  975.   }
  976.  
  977.   /* check if we should sample between points 1 and 2 */
  978. #ifdef EmissionMethodPatch
  979.   if (( Colour_Distance(C1, C2) > aa_threshold ) || (Colour_Distance(Em1, Em2) > aa_threshold))
  980. #else
  981.   if ( Colour_Distance(C1, C2) > aa_threshold )
  982. #endif
  983.   {
  984.     /* recurse again */
  985.     sample_media_rec(Light_List, Ray, Media_List, Interval, light_ray, 
  986.                      d1, d2,
  987.                      Result2, C1, C2,
  988.                      ODResult2, od1, od2,
  989.                      depth-1, Jitter, aa_threshold
  990. #ifdef EmissionMethodPatch
  991.                      , Em1, Em2, EmR2
  992. #endif
  993.                      );
  994.  
  995.     /* average colors */
  996.     Result[0]=Result2[0]*0.5;
  997.     Result[1]=Result2[1]*0.5;
  998.     Result[2]=Result2[2]*0.5;
  999.  
  1000.     /* average the optical depth */
  1001.     ODResult[0]=ODResult2[0]*0.5;
  1002.     ODResult[1]=ODResult2[1]*0.5;
  1003.     ODResult[2]=ODResult2[2]*0.5;
  1004.  
  1005. #ifdef EmissionMethodPatch
  1006.     EmR[0] = EmR2[0]*0.5;
  1007.     EmR[1] = EmR2[1]*0.5;
  1008.     EmR[2] = EmR2[2]*0.5;
  1009. #endif
  1010.   }
  1011.   else
  1012.   {
  1013.     /* no new points needed - just average what we've got.
  1014.        use c1 weight = 1/3, c2 weight = 1/6 (we use c2, the middle point, again)
  1015.     */
  1016.     /* average colors */
  1017.     Result[0]=C1[0]/3.0 + C2[0]/6.0;
  1018.     Result[1]=C1[1]/3.0 + C2[1]/6.0;
  1019.     Result[2]=C1[2]/3.0 + C2[2]/6.0;
  1020.  
  1021.     /* average the optical depth */
  1022.     ODResult[0]=od1[0]/3.0 + od2[0]/6.0;
  1023.     ODResult[1]=od1[1]/3.0 + od2[1]/6.0;
  1024.     ODResult[2]=od1[2]/3.0 + od2[2]/6.0;
  1025.  
  1026. #ifdef EmissionMethodPatch
  1027.     EmR[0] = Em1[0]/3.0 + Em2[0]/6.0;
  1028.     EmR[1] = Em1[1]/3.0 + Em2[1]/6.0;
  1029.     EmR[2] = Em1[2]/3.0 + Em2[2]/6.0;
  1030. #endif
  1031.   }
  1032.  
  1033.   /* check if we should sample between points 2 and 3 */
  1034. #ifdef EmissionMethodPatch
  1035.   if (( Colour_Distance(C2, C3) > aa_threshold ) || (Colour_Distance(Em2, Em3) > aa_threshold))
  1036. #else
  1037.   if ( Colour_Distance(C2, C3) > aa_threshold )
  1038. #endif
  1039.   {
  1040.     /* recurse again */
  1041.     sample_media_rec(Light_List, Ray, Media_List, Interval, light_ray,
  1042.                      d2, d3,
  1043.                      Result2, C2, C3,
  1044.                      ODResult2, od2, od3,
  1045.                      depth-1, Jitter, aa_threshold
  1046. #ifdef EmissionMethodPatch
  1047.                      , Em2, Em3, EmR2
  1048. #endif
  1049.                      );
  1050.  
  1051.     /* average colors */
  1052.     Result[0]+=Result2[0]*0.5;
  1053.     Result[1]+=Result2[1]*0.5;
  1054.     Result[2]+=Result2[2]*0.5;
  1055.  
  1056.     /* average the optical depth */
  1057.     ODResult[0]+=ODResult2[0]*0.5;
  1058.     ODResult[1]+=ODResult2[1]*0.5;
  1059.     ODResult[2]+=ODResult2[2]*0.5;
  1060.  
  1061. #ifdef EmissionMethodPatch
  1062.     EmR[0] += EmR2[0]*0.5;
  1063.     EmR[1] += EmR2[1]*0.5;
  1064.     EmR[2] += EmR2[2]*0.5;
  1065. #endif
  1066.   }
  1067.   else
  1068.   {
  1069.     /* no new points needed - just average what we've got.
  1070.        use c1 weight = 1/3, c2 weight = 1/6 (we already used c2 once)
  1071.     */
  1072.     /* average colors */
  1073.     Result[0]+=C3[0]/3.0 + C2[0]/6.0;
  1074.     Result[1]+=C3[1]/3.0 + C2[1]/6.0;
  1075.     Result[2]+=C3[2]/3.0 + C2[2]/6.0;
  1076.  
  1077.     /* average the optical depth */
  1078.     ODResult[0]+=od3[0]/3.0 + od2[0]/6.0;
  1079.     ODResult[1]+=od3[1]/3.0 + od2[1]/6.0;
  1080.     ODResult[2]+=od3[2]/3.0 + od2[2]/6.0;
  1081.  
  1082. #ifdef EmissionMethodPatch
  1083.     EmR[0] += Em3[0]/3.0 + Em2[0]/6.0;
  1084.     EmR[1] += Em3[1]/3.0 + Em2[1]/6.0;
  1085.     EmR[2] += Em3[2]/3.0 + Em2[2]/6.0;
  1086. #endif
  1087.   }
  1088.   
  1089.  
  1090.  
  1091. }
  1092.  
  1093.  
  1094. /*****************************************************************************
  1095. *
  1096. * FUNCTION
  1097. *
  1098. *   sample_media
  1099. *
  1100. * INPUT
  1101. *
  1102. *   dist  - distance of current sample
  1103. *   Ray   - pointer to ray
  1104. *   IMedia - pointer to media to use
  1105. *
  1106. * OUTPUT
  1107. *
  1108. *   Col          - color of current sample
  1109. *
  1110. * RETURNS
  1111. *
  1112. * AUTHOR
  1113. *
  1114. *   Dieter Bayer
  1115. *
  1116. * DESCRIPTION
  1117. *
  1118. *   Calculate the color of the current media sample.
  1119. *
  1120. * CHANGES
  1121. *
  1122. *   Nov 1994 : Creation.
  1123. *
  1124. ******************************************************************************/
  1125.  
  1126. static void sample_media(LIGHT_LIST *Light_List, RAY *Ray, IMEDIA **Media_List, MEDIA_INTERVAL *Interval, int light_ray, DBL d0,
  1127.                          COLOUR SampCol, COLOUR SampOptDepth, int sample_method
  1128. #ifdef EmissionMethodPatch
  1129.                          , COLOUR SampEm2
  1130. #endif
  1131.                          )
  1132. {
  1133.   int i, n, use_scattering;
  1134.   /* NK samples - moved d0 to parameter list */
  1135.   DBL alpha, /*d0,*/ d1, len, k, g, g2;
  1136.   VECTOR P, H;
  1137.   COLOUR C0, Light_Colour, Te;
  1138.   COLOUR Em, Ex, Sc;
  1139.   RAY Light_Ray;
  1140.   IMEDIA **Tmp, *Local;
  1141.   /* NK phmap */
  1142.   int ignore_photons;
  1143.   /* NK ---- */
  1144. #ifdef EmissionMethodPatch
  1145.   COLOUR Em2;            /* JB */
  1146. #endif
  1147.  
  1148.   Increase_Counter(stats[Media_Samples]);
  1149.  
  1150.   /* Set up sampling location. */
  1151.  
  1152.   /*
  1153.   d0 = Interval->ds * FRAND();
  1154.   */
  1155.  
  1156.   /* NK intervals */
  1157.   d0 *= Interval->ds;
  1158.   /* NK ---- */
  1159.  
  1160.   d1 = Interval->s0 + d0;
  1161.  
  1162.   VEvaluateRay(H, Ray->Initial, d1, Ray->Direction);
  1163.  
  1164.   /* Get coefficients in current sample location. */
  1165.  
  1166.   Make_Colour(Em, 0.0, 0.0, 0.0);
  1167.   Make_Colour(Ex, 0.0, 0.0, 0.0);
  1168.   Make_Colour(Sc, 0.0, 0.0, 0.0);
  1169. #ifdef EmissionMethodPatch
  1170.   Make_Colour(Em2, 0.0, 0.0, 0.0); /* JB */
  1171. #endif
  1172.  
  1173.   use_scattering = FALSE;
  1174.  
  1175.   for (Tmp = Media_List; (*Tmp) != NULL; Tmp++)
  1176.   {
  1177.     for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
  1178.     {
  1179.       /* NK !!!!!!! check this for photon problems */
  1180.         if (light_ray && currLight && !Check_Media_Light_Group(Local,currLight))
  1181.             continue;
  1182.  
  1183.       Assign_Vector(P, H);
  1184.  
  1185.       evaluate_density_pattern(Local, P, C0);
  1186.  
  1187.       Ex[0] += C0[0] * Local->Extinction[0];
  1188.       Ex[1] += C0[1] * Local->Extinction[1];
  1189.       Ex[2] += C0[2] * Local->Extinction[2];
  1190.  
  1191. #ifdef EmissionMethodPatch
  1192.       if (Local->Emission_Method == 1) {
  1193.         DBL bw = (0.30*C0[0]*Local->Emission[0]+
  1194.         0.59*C0[1]*Local->Emission[1]+
  1195.         0.11*C0[2]*Local->Emission[2])*Local->Emission_Extinction;
  1196.         Ex[0] += bw;
  1197.         Ex[1] += bw;
  1198.         Ex[2] += bw;
  1199.       }
  1200. #endif
  1201.  
  1202.       if (!light_ray)
  1203.       {
  1204. #ifdef EmissionMethodPatch
  1205.         if (Local->Emission_Method == 1) 
  1206.         {
  1207.           Em2[0] += C0[0] * Local->Emission[0];
  1208.           Em2[1] += C0[1] * Local->Emission[1];
  1209.           Em2[2] += C0[2] * Local->Emission[2];
  1210.         } 
  1211.         else 
  1212.         {
  1213. #endif
  1214.           Em[0] += C0[0] * Local->Emission[0];
  1215.           Em[1] += C0[1] * Local->Emission[1];
  1216.           Em[2] += C0[2] * Local->Emission[2];
  1217. #ifdef EmissionMethodPatch
  1218.         }
  1219. #endif
  1220.  
  1221.         Sc[0] += C0[0] * Local->Scattering[0];
  1222.         Sc[1] += C0[1] * Local->Scattering[1];
  1223.         Sc[2] += C0[2] * Local->Scattering[2];
  1224.       }
  1225.  
  1226.       use_scattering |= Local->use_scattering;
  1227.     }
  1228.   }
  1229.  
  1230.   /* Get estimate for the total optical depth of the current interval. */
  1231.  
  1232.   SampOptDepth[0] = Ex[0] * Interval->ds;
  1233.   SampOptDepth[1] = Ex[1] * Interval->ds;
  1234.   SampOptDepth[2] = Ex[2] * Interval->ds;
  1235.  
  1236.   if(sample_method!=3)
  1237.   {
  1238.     Interval->od[0] += Ex[0] * Interval->ds;
  1239.     Interval->od[1] += Ex[1] * Interval->ds;
  1240.     Interval->od[2] += Ex[2] * Interval->ds;
  1241.   }
  1242.  
  1243.  
  1244.   /* Get estimate for the total emission of the current interval. */
  1245.  
  1246.   Te[0] = Em[0];
  1247.   Te[1] = Em[1];
  1248.   Te[2] = Em[2];
  1249.  
  1250.   if ((!light_ray) && (use_scattering) && (Interval->lit))
  1251.   {
  1252.     /* NK phmap */
  1253.     /* set up flags for photons */
  1254.     if(!backtraceFlag)
  1255.     {
  1256.     photonOptions.objectFlags = PH_FLAG_RFR_ON | PH_FLAG_RFL_ON;
  1257.  
  1258.     ignore_photons = TRUE;
  1259.     for (n = 0, Tmp = Media_List; (*Tmp) != NULL; n++, Tmp++)
  1260.       for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
  1261.         if (!Local->ignore_photons)
  1262.         {
  1263.           ignore_photons = FALSE;
  1264.           break;
  1265.         }
  1266.     if (ignore_photons)
  1267.       photonOptions.objectFlags |= PH_FLAG_IGNORE_PHOTONS;
  1268.  
  1269.     photonOptions.photonObject = NULL;
  1270.     }
  1271.     /* NK ---- */
  1272.  
  1273.     /* Process all light sources. */
  1274.  
  1275.     for (i = 0; i < Frame.Number_Of_Light_Sources; i++)
  1276.     {
  1277.       /* Use light only if active and within it's boundaries. */
  1278.  
  1279.       if ((Light_List[i].active) && (d1 >= Light_List[i].s0) && (d1 <= Light_List[i].s1))
  1280.       {
  1281.         photonOptions.lightFlags = Light_List[i].Light->Ph_Flags;
  1282.         if (!(Test_Shadow(Light_List[i].Light, &len, &Light_Ray, Ray, P, Light_Colour)))
  1283.         {
  1284.           /* Get attenuation due to scattering. */
  1285.  
  1286.           k = 0.0;
  1287.  
  1288.           for (n = 0, Tmp = Media_List; (*Tmp) != NULL; n++, Tmp++)
  1289.           {
  1290.             for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
  1291.             {
  1292.               if (!Check_Media_Light_Group(Local, Light_List[i].Light))
  1293.                   continue;
  1294.  
  1295.               switch (Local->Type)
  1296.               {
  1297.                 case RAYLEIGH_SCATTERING:
  1298.  
  1299.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1300.  
  1301.                   k += 0.799372013 * (1.0 + Sqr(alpha));
  1302.  
  1303.                   break;
  1304.  
  1305.                 case MIE_HAZY_SCATTERING:
  1306.  
  1307.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1308.  
  1309.                   k += 0.576655375 * (1.0 + 9.0 * pow(0.5 * (1.0 + alpha), 8.0));
  1310.  
  1311.                   break;
  1312.  
  1313.                 case MIE_MURKY_SCATTERING:
  1314.  
  1315.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1316.  
  1317.                   k += 0.495714547 * (1.0 + 50.0 * pow(0.5 * (1.0 + alpha), 32.0));
  1318.  
  1319.                   break;
  1320.  
  1321.                 case HENYEY_GREENSTEIN_SCATTERING:
  1322.  
  1323.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1324.  
  1325.                   g = Local->Eccentricity;
  1326.  
  1327.                   g2 = Sqr(g);
  1328.  
  1329.                   k += (1.0 - g2) / pow(1.0 + g2 - 2.0 * g * alpha, 1.5);
  1330.  
  1331.                   break;
  1332.  
  1333.                 case ISOTROPIC_SCATTERING:
  1334.                 default:
  1335.  
  1336.                   k += 1.0;
  1337.  
  1338.                 break;
  1339.               }
  1340.             }
  1341.           }
  1342.  
  1343.           k /= (DBL)n;
  1344.  
  1345.           Te[0] += k * Sc[0] * Light_Colour[0];
  1346.           Te[1] += k * Sc[1] * Light_Colour[1];
  1347.           Te[2] += k * Sc[2] * Light_Colour[2];
  1348.         }
  1349.       }
  1350.     }
  1351.   }
  1352.  
  1353. #if MEDIA_INTERACTION
  1354.   /* use all intervals for photons */
  1355.   if ((!light_ray) && (use_scattering))
  1356.   {
  1357.     if (photonOptions.mediaPhotonMap.numPhotons>0)
  1358.     {
  1359.       /* Process photons. */
  1360.       int n,tempn,j,step=0;
  1361.       int minGatherCount;
  1362.       DBL Size,r=0,tempr=0.0;
  1363.       DBL thisDensity=0;
  1364.       DBL prevDensity=0.0000000000000001; /* avoid div-by-zero error */
  1365.       int expanded=FALSE;
  1366.       COLOUR Colour2, TempCol;
  1367.  
  1368.       Size = photonOptions.mediaPhotonMap.minGatherRad;
  1369.       /* size /= 2; */
  1370.       Make_Colour(Colour2,0,0,0);
  1371.       n=-1;
  1372.       step=0;
  1373.       minGatherCount=photonOptions.minGatherCount;
  1374.       while(n<minGatherCount && step<photonOptions.mediaPhotonMap.gatherNumSteps)
  1375.       {
  1376.         Make_Colour(TempCol,0,0,0);
  1377.         tempr = 0;
  1378.  
  1379.         /* gather the photons */
  1380.         tempn=gatherPhotons(H, Size, &tempr,NULL,FALSE,&photonOptions.mediaPhotonMap);
  1381.  
  1382.         /* now go through these photons and add up their contribution */
  1383.         for(j=0; j<tempn; j++)
  1384.         {
  1385.           SMALL_COLOUR col;
  1386.           /*DBL theta,phi;*/
  1387.           int theta,phi;
  1388.  
  1389.           /* convert small color to normal color */
  1390.           col = photonOptions.photonGatherList[j]->Colour;
  1391.           UNPACK_COLOUR(Light_Colour, col, photonOptions.mediaPhotonMap);
  1392.  
  1393.           /* convert theta/phi to vector direction 
  1394.              Use a pre-computed array of sin/cos to avoid many calls to the
  1395.              sin() and cos() functions.  These arrays were initialized in
  1396.              InitBacktraceEverything.
  1397.           */
  1398.           theta = photonOptions.photonGatherList[j]->theta+127;
  1399.           phi = photonOptions.photonGatherList[j]->phi+127;
  1400.       
  1401.           Light_Ray.Direction[Y] = photonOptions.sinTheta[theta];
  1402.           Light_Ray.Direction[X] = photonOptions.cosTheta[theta];
  1403.  
  1404.           Light_Ray.Direction[Z] = Light_Ray.Direction[X]*photonOptions.sinTheta[phi];
  1405.           Light_Ray.Direction[X] = Light_Ray.Direction[X]*photonOptions.cosTheta[phi];
  1406.  
  1407.           VSub(Light_Ray.Initial, photonOptions.photonGatherList[j]->Loc, Light_Ray.Direction);
  1408.  
  1409.           /* Get attenuation due to scattering. */
  1410.  
  1411.           k = 0.0;
  1412.  
  1413.           for (n = 0, Tmp = Media_List; (*Tmp) != NULL; n++, Tmp++)
  1414.           {
  1415.             for (Local = *Tmp; Local != NULL; Local = Local->Next_Media)
  1416.             {
  1417.               switch (Local->Type)
  1418.               {
  1419.                 case RAYLEIGH_SCATTERING:
  1420.  
  1421.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1422.  
  1423.                   k += 0.799372013 * (1.0 + Sqr(alpha));
  1424.  
  1425.                   break;
  1426.  
  1427.                 case MIE_HAZY_SCATTERING:
  1428.  
  1429.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1430.  
  1431.                   k += 0.576655375 * (1.0 + 9.0 * pow(0.5 * (1.0 + alpha), 8.0));
  1432.  
  1433.                   break;
  1434.  
  1435.                 case MIE_MURKY_SCATTERING:
  1436.  
  1437.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1438.  
  1439.                   k += 0.495714547 * (1.0 + 50.0 * pow(0.5 * (1.0 + alpha), 32.0));
  1440.  
  1441.                   break;
  1442.  
  1443.                 case HENYEY_GREENSTEIN_SCATTERING:
  1444.  
  1445.                   VDot(alpha, Light_Ray.Direction, Ray->Direction);
  1446.  
  1447.                   g = Local->Eccentricity;
  1448.  
  1449.                   g2 = Sqr(g);
  1450.  
  1451.                   k += (1.0 - g2) / pow(1.0 + g2 - 2.0 * g * alpha, 1.5);
  1452.  
  1453.                   break;
  1454.  
  1455.                 case ISOTROPIC_SCATTERING:
  1456.                 default:
  1457.  
  1458.                   k += 1.0;
  1459.  
  1460.                 break;
  1461.               }
  1462.             }
  1463.           }
  1464.  
  1465.           k /= (DBL)n;
  1466.  
  1467.           TempCol[0] += k * Sc[0] * Light_Colour[0];
  1468.           TempCol[1] += k * Sc[1] * Light_Colour[1];
  1469.           TempCol[2] += k * Sc[2] * Light_Colour[2];
  1470.         }
  1471.         /* density of this search */
  1472.         thisDensity = tempn / (tempr*tempr*tempr*4.0/3.0);
  1473.  
  1474.         if(((thisDensity-prevDensity)/prevDensity < photonOptions.expandTolerance) 
  1475.            || (step==0) 
  1476.            || (tempn<photonOptions.minExpandCount && tempn>0))
  1477.         {
  1478.           /* it passes the tests, so use the new color */
  1479.       
  1480.           if (step>0)
  1481.             expanded = TRUE;
  1482.  
  1483.           prevDensity = thisDensity;
  1484.           if (prevDensity==0)
  1485.             prevDensity = 0.0000000000000001;  /* avoid div-by-zero error */
  1486.  
  1487.           Assign_Colour(Colour2, TempCol);
  1488.  
  1489.           r = tempr;
  1490.           n = tempn;
  1491.         }
  1492.  
  1493.         Size+=photonOptions.mediaPhotonMap.gatherRadStep;
  1494.         step++;
  1495.         /* force one step only for now */
  1496.         step = photonOptions.mediaPhotonMap.gatherNumSteps;
  1497.       }
  1498.       /* finish the photons equation */
  1499.       VScaleEq(Colour2, (DBL)(3.0)/(M_PI*r*r*r*4.0));
  1500.       
  1501.       Te[0] += Colour2[0];
  1502.       Te[1] += Colour2[1];
  1503.       Te[2] += Colour2[2];
  1504.     }
  1505. #endif
  1506.   }
  1507.  
  1508.   if (sample_method==3)
  1509.   {
  1510.     /* we're doing the samples in order, so we can attenuate correctly
  1511.        instead of assuming a constant absorption/extinction */
  1512.     DBL unused;
  1513.     DBL coef = modf(d0*sampCount_s/Interval->ds, &unused);
  1514.     if (coef < EPSILON) coef = 1.0;
  1515.     coef *= Interval->ds/sampCount_s;
  1516.     Te[0] *= Interval->ds * exp(-Ex[0] * coef/*Interval->ds/sampCount_s*/);
  1517.     Te[1] *= Interval->ds * exp(-Ex[1] * coef/*Interval->ds/sampCount_s*/);
  1518.     Te[2] *= Interval->ds * exp(-Ex[2] * coef/*Interval->ds/sampCount_s*/);
  1519.   }
  1520.   else
  1521.   {
  1522.     /* NOTE: this assumes constant absorption+extinction over the length
  1523.        of the interval */
  1524.     Te[0] *= Interval->ds * exp(-Ex[0] * d0);
  1525.     Te[1] *= Interval->ds * exp(-Ex[1] * d0);
  1526.     Te[2] *= Interval->ds * exp(-Ex[2] * d0);
  1527.   }
  1528.  
  1529. #ifdef EmissionMethodPatch
  1530.   Em2[0] *= Interval->ds;
  1531.   Em2[1] *= Interval->ds;
  1532.   Em2[2] *= Interval->ds;
  1533. #endif
  1534.  
  1535.   SampCol[0] = Te[0];
  1536.   SampCol[1] = Te[1];
  1537.   SampCol[2] = Te[2];
  1538.  
  1539. #ifdef EmissionMethodPatch
  1540.   SampEm2[0] = Em2[0];
  1541.   SampEm2[1] = Em2[1];
  1542.   SampEm2[2] = Em2[2];
  1543. #endif
  1544.  
  1545.   if(sample_method!=3)
  1546.   {
  1547.     /* Add emission. */
  1548.  
  1549.     Interval->te[0] += Te[0];
  1550.     Interval->te[1] += Te[1];
  1551.     Interval->te[2] += Te[2];
  1552.  
  1553.     Interval->te2[0] += Sqr(Te[0]);
  1554.     Interval->te2[1] += Sqr(Te[1]);
  1555.     Interval->te2[2] += Sqr(Te[2]);
  1556.  
  1557. #ifdef EmissionMethodPatch
  1558.                 /* JB +++ */
  1559.     Interval->Em2[0] += Em2[0];
  1560.     Interval->Em2[1] += Em2[1];
  1561.     Interval->Em2[2] += Em2[2];
  1562.                 /* JB --- */
  1563. #endif
  1564.   }
  1565.  
  1566.   Interval->samples++;
  1567. }
  1568.  
  1569.  
  1570.  
  1571. /*****************************************************************************
  1572. *
  1573. * FUNCTION
  1574. *
  1575. *   evaluate_density_pattern
  1576. *
  1577. * INPUT
  1578. *
  1579. * OUTPUT
  1580. *
  1581. * RETURNS
  1582. *
  1583. * AUTHOR
  1584. *
  1585. *   Dieter Bayer
  1586. *
  1587. * DESCRIPTION
  1588. *
  1589. * CHANGES
  1590. *
  1591. *   Dec 1996 : Creation.
  1592. *
  1593. ******************************************************************************/
  1594.  
  1595. static void evaluate_density_pattern(IMEDIA *IMedia, VECTOR P, COLOUR C)
  1596. {
  1597.   COLOUR Local_Color;
  1598.   PIGMENT *Temp=IMedia->Density;
  1599.   
  1600.   Make_Colour (C, 1.0, 1.0, 1.0);
  1601.  
  1602.   while (Temp != NULL)
  1603.   {
  1604.     Make_Colour (Local_Color, 0.0, 0.0, 0.0);
  1605.     
  1606.     Compute_Pigment (Local_Color, Temp, P, NULL);
  1607.  
  1608.     C[0] *= Local_Color[0];
  1609.     C[1] *= Local_Color[1];
  1610.     C[2] *= Local_Color[2];
  1611.  
  1612.     Temp=(PIGMENT *)Temp->Next;
  1613.   }
  1614. }
  1615.  
  1616.  
  1617.  
  1618. /*****************************************************************************
  1619. *
  1620. * FUNCTION
  1621. *
  1622. *   get_light_list
  1623. *
  1624. * INPUT
  1625. *
  1626. *   Light_List - array containing light source information
  1627. *   dist       - distance of current sample
  1628. *   Ray        - pointer to ray
  1629. *   IMedia      - pointer to media to use
  1630. *
  1631. * OUTPUT
  1632. *
  1633. *   Col        - color of current sample
  1634. *
  1635. * RETURNS
  1636. *
  1637. * AUTHOR
  1638. *
  1639. *   Dieter Bayer
  1640. *
  1641. * DESCRIPTION
  1642. *
  1643. * CHANGES
  1644. *
  1645. *   Nov 1994 : Creation.
  1646. *
  1647. ******************************************************************************/
  1648.  
  1649. static void get_light_list(LIGHT_LIST *Light_List, RAY *Ray, INTERSECTION *Inter)
  1650. {
  1651.   int i, insert;
  1652.   DBL t1, t2;
  1653.   LIGHT_SOURCE *Light;
  1654.  
  1655.   /* Get depths for all light sources and disconnected sampling intervals. */
  1656.  
  1657.   t1 = t2 = 0.0;
  1658.  
  1659.   for (i = 0, Light = Frame.Light_Sources; Light != NULL; Light = Light->Next_Light_Source, i++)
  1660.   {
  1661.     /* Init interval. */
  1662.  
  1663.     Light_List[i].active = FALSE;
  1664.     Light_List[i].s0     = 0.0;
  1665.     Light_List[i].s1     = Max_Distance;
  1666.     Light_List[i].Light  = NULL;
  1667.  
  1668.     insert = FALSE;
  1669.  
  1670.     Light_List[i].Light = Light;
  1671.  
  1672.     if (!Light->Media_Interaction)
  1673.     {
  1674.       continue;
  1675.     }
  1676.  
  1677.     switch (Light->Light_Type)
  1678.     {
  1679.       case CYLINDER_SOURCE:
  1680.  
  1681.         if (intersect_cylinderlight(Ray, Light, &t1, &t2))
  1682.         {
  1683.           if ((t1 < Inter->Depth) && (t2 > Small_Tolerance))
  1684.           {
  1685.             insert = TRUE;
  1686.           }
  1687.         }
  1688.  
  1689.         break;
  1690.  
  1691.       case POINT_SOURCE:
  1692.  
  1693.         t1 = 0.0;
  1694.         t2 = Inter->Depth;
  1695.  
  1696.         insert = TRUE;
  1697.  
  1698.         break;
  1699.  
  1700.       case SPOT_SOURCE:
  1701.  
  1702.         if (intersect_spotlight(Ray, Light, &t1, &t2))
  1703.         {
  1704.           if ((t1 < Inter->Depth) && (t2 > Small_Tolerance))
  1705.           {
  1706.             insert = TRUE;
  1707.           }
  1708.         }
  1709.  
  1710.         break;
  1711.     }
  1712.  
  1713.     /* Insert distances into sampling interval list. */
  1714.  
  1715.     if (insert)
  1716.     {
  1717.       /* Insert light source intersections into light list. */
  1718.  
  1719.       t1 = max(t1, 0.0);
  1720.       t2 = min(t2, Inter->Depth);
  1721.  
  1722.       Light_List[i].active = TRUE;
  1723.       Light_List[i].s0 = t1;
  1724.       Light_List[i].s1 = t2;
  1725.     }
  1726.   }
  1727. }
  1728.  
  1729.  
  1730.  
  1731. /*****************************************************************************
  1732. *
  1733. * FUNCTION
  1734. *
  1735. *   get_lit_interval
  1736. *
  1737. * INPUT
  1738. *
  1739. *   Light_List - array containing light source information
  1740. *   dist       - distance of current sample
  1741. *   Ray        - pointer to ray
  1742. *   IMedia      - pointer to media to use
  1743. *
  1744. * OUTPUT
  1745. *
  1746. *   Col        - color of current sample
  1747. *
  1748. * RETURNS
  1749. *
  1750. * AUTHOR
  1751. *
  1752. *   Dieter Bayer
  1753. *
  1754. * DESCRIPTION
  1755. *
  1756. * CHANGES
  1757. *
  1758. *   Nov 1994 : Creation.
  1759. *
  1760. ******************************************************************************/
  1761.  
  1762. static void get_lit_interval(int *number, LIT_INTERVAL *Lit_Interval, int entries, LIGHT_LIST *Light_List, INTERSECTION *Inter)
  1763. {
  1764.   int a, i, n;
  1765.   LIT_INTERVAL *curr, *prev;
  1766. #ifndef UseMediaAndLightCache
  1767.   DBL *s0, *s1;
  1768.  
  1769.   s0 = (DBL *)POV_MALLOC(Frame.Number_Of_Light_Sources*sizeof(DBL), "temp data");
  1770.   s1 = (DBL *)POV_MALLOC(Frame.Number_Of_Light_Sources*sizeof(DBL), "temp data");
  1771. #endif
  1772.  
  1773.   for (i = a = 0; i < entries; i++)
  1774.   {
  1775.     if (Light_List[i].active)
  1776.     {
  1777.       s0[a] = Light_List[i].s0;
  1778.       s1[a] = Light_List[i].s1;
  1779.  
  1780.       a++;
  1781.     }
  1782.   }
  1783.  
  1784.   n = 0;
  1785.  
  1786.   curr = &Lit_Interval[0];
  1787.  
  1788.   if (a)
  1789.   {
  1790.     QSORT((void *)(&s0[0]), (unsigned long)a, sizeof(DBL), compdoubles);
  1791.     QSORT((void *)(&s1[0]), (unsigned long)a, sizeof(DBL), compdoubles);
  1792.  
  1793.     if (s0[0] > 0.0)
  1794.     {
  1795.       curr->lit = FALSE;
  1796.  
  1797.       curr->s0 = 0.0;
  1798.       curr->s1 = s0[0];
  1799.  
  1800.       curr++;
  1801.  
  1802.       n++;
  1803.     }
  1804.  
  1805.     curr->lit = TRUE;
  1806.  
  1807.     curr->s0 = s0[0];
  1808.     curr->s1 = s1[0];
  1809.  
  1810.     prev = curr;
  1811.  
  1812.     curr++;
  1813.  
  1814.     n++;
  1815.  
  1816.     for (i = 1; i < a; i++)
  1817.     {
  1818.       if (s0[i] > prev->s1)
  1819.       {
  1820.         curr->lit = FALSE;
  1821.  
  1822.         curr->s0 = prev->s1;
  1823.         curr->s1 = s0[i];
  1824.  
  1825.         prev++;
  1826.         curr++;
  1827.  
  1828.         n++;
  1829.  
  1830.         curr->lit = TRUE;
  1831.  
  1832.         curr->s0 = s0[i];
  1833.         curr->s1 = s1[i];
  1834.  
  1835.         prev++;
  1836.         curr++;
  1837.  
  1838.         n++;
  1839.       }
  1840.       else
  1841.       {
  1842.         if (s1[i] > prev->s1)
  1843.         {
  1844.           prev->s1 = s1[i];
  1845.         }
  1846.       }
  1847.     }
  1848.  
  1849.     if (prev->s1 < Inter->Depth)
  1850.     {
  1851.       curr->lit = FALSE;
  1852.  
  1853.       curr->s0 = prev->s1;
  1854.       curr->s1 = Inter->Depth;
  1855.  
  1856.       curr++;
  1857.  
  1858.       n++;
  1859.     }
  1860.   }
  1861.   else
  1862.   {
  1863.     curr->lit = FALSE;
  1864.  
  1865.     curr->s0 = 0.0;
  1866.     curr->s1 = Inter->Depth;
  1867.  
  1868.     curr++;
  1869.  
  1870.     n++;
  1871.   }
  1872.  
  1873.   curr = &Lit_Interval[0];
  1874.  
  1875.   for (i = 0; i < n; i++)
  1876.   {
  1877.     curr->ds = curr->s1 - curr->s0;
  1878.  
  1879.     curr++;
  1880.   }
  1881.  
  1882. #ifndef UseMediaAndLightCache
  1883.   POV_FREE(s0);
  1884.   POV_FREE(s1);
  1885. #endif
  1886.  
  1887.   *number = n;
  1888. }
  1889.  
  1890.  
  1891.  
  1892. /*****************************************************************************
  1893. *
  1894. * FUNCTION
  1895. *
  1896. *   set_up_sampling_intervals
  1897. *
  1898. * INPUT
  1899. *
  1900. *   interval - array containing media intervals
  1901. *   number   - number of lit intervals
  1902. *   list     - array of lit intervals
  1903. *   media    - media to use
  1904. *
  1905. * OUTPUT
  1906. *
  1907. *   interval - array containing media intervals
  1908. *
  1909. * RETURNS
  1910. *
  1911. *   int      - number of media intervals created
  1912. *
  1913. * AUTHOR
  1914. *
  1915. *   Dieter Bayer
  1916. *
  1917. * DESCRIPTION
  1918. *
  1919. *   Distribute samples along an interval according to the maximum
  1920. *   number of samples and the ratio of samples in lit and unlit
  1921. *   areas as given by the participating media.
  1922. *
  1923. * CHANGES
  1924. *
  1925. *   Nov 1994 : Creation.
  1926. *
  1927. ******************************************************************************/
  1928.  
  1929. static int set_up_sampling_intervals(MEDIA_INTERVAL *interval, int number, LIT_INTERVAL *list, IMEDIA *media)
  1930. {
  1931.   int i, j, n, r, remaining, intervals;
  1932.   DBL delta, sum, weight;
  1933.   MEDIA_INTERVAL *curr;
  1934.   LIT_INTERVAL *entry;
  1935.  
  1936.   /* Set up sampling intervals. */
  1937.  
  1938.   /* NK samples - we will always have enough intervals 
  1939.      we always use the larger of the two numbers */
  1940.   if (media->Intervals < number)
  1941.     intervals = number;
  1942.   else
  1943.     intervals = media->Intervals;
  1944.  
  1945.   /* Use one interval if no lit intervals and constant media. */
  1946.  
  1947.   if ((number == 0) && (media->is_constant))
  1948.   {
  1949.     intervals = 1;
  1950.  
  1951.     delta = list[0].ds;
  1952.  
  1953.     curr = interval;
  1954.  
  1955.     curr->lit = TRUE;
  1956.  
  1957.     curr->samples = 0;
  1958.  
  1959.     curr->s0 = list[0].s0;
  1960.     curr->s1 = list[0].s0 + delta;
  1961.     curr->ds = delta;
  1962.  
  1963.     Make_Colour(curr->od,  0.0, 0.0, 0.0);
  1964.     Make_Colour(curr->te,  0.0, 0.0, 0.0);
  1965.     Make_Colour(curr->te2, 0.0, 0.0, 0.0);
  1966. #ifdef EmissionMethodPatch
  1967.     Make_Colour(curr->Em2, 0.0, 0.0, 0.0); /* JB */
  1968. #endif
  1969.  
  1970.     return(intervals);
  1971.   }
  1972.  
  1973.   /* Choose intervals. */
  1974.  
  1975.   if (number == 1)
  1976.   {
  1977.     /* Use uniform intervals. */
  1978.  
  1979.     delta = list[0].ds / (DBL)intervals;
  1980.  
  1981.     curr = interval;
  1982.  
  1983.     for (i = 0; i < intervals; i++)
  1984.     {
  1985.       curr->lit = TRUE;
  1986.  
  1987.       curr->samples = 0;
  1988.  
  1989.       curr->s0 = list[0].s0 + delta * (DBL)i;
  1990.       curr->s1 = list[0].s0 + delta * (DBL)(i + 1);
  1991.       curr->ds = delta;
  1992.  
  1993.       Make_Colour(curr->od,  0.0, 0.0, 0.0);
  1994.       Make_Colour(curr->te,  0.0, 0.0, 0.0);
  1995.       Make_Colour(curr->te2, 0.0, 0.0, 0.0);
  1996. #ifdef EmissionMethodPatch
  1997.       Make_Colour(curr->Em2, 0.0, 0.0, 0.0);
  1998. #endif
  1999.  
  2000.       curr++;
  2001.     }
  2002.   }
  2003.   else
  2004.   {
  2005.     /* Choose intervals according to the specified ratio. */
  2006.  
  2007.     if (number > intervals)
  2008.     {
  2009.       Error("Too few sampling intervals.\n");
  2010.     }
  2011.  
  2012.     sum = 0.0;
  2013.  
  2014.     entry = list;
  2015.  
  2016.     for (i = 0; i < number; i++)
  2017.     {
  2018. /*
  2019.       sum += ((entry->lit) ? (0.9) : (0.1)) * entry->ds;
  2020. */
  2021.       sum += ((entry->lit) ? (media->Ratio) : (1.0 - media->Ratio));
  2022.  
  2023.       entry++;
  2024.     }
  2025.  
  2026.     remaining = intervals;
  2027.  
  2028.     curr = interval;
  2029.  
  2030.     entry = list;
  2031.  
  2032.     for (i = 0; i < number; i++)
  2033.     {
  2034. /*
  2035.       weight = ((entry->lit) ? (0.9) : (0.1)) * entry->ds;
  2036. */
  2037.       weight = ((entry->lit) ? (media->Ratio) : (1.0 - media->Ratio));
  2038.  
  2039.       n = (int)(weight / sum * (DBL)intervals) + 1;
  2040.  
  2041.       r = remaining - number + i + 1;
  2042.  
  2043.       if (n > r)
  2044.       {
  2045.         n = r;
  2046.       }
  2047.  
  2048.       delta = entry->ds / (DBL)n;
  2049.  
  2050.       for (j = 0; j < n; j++)
  2051.       {
  2052.         curr->lit = entry->lit;
  2053.  
  2054.         curr->samples = 0;
  2055.  
  2056.         curr->s0 = entry->s0 + delta * (DBL)j;
  2057.         curr->s1 = entry->s0 + delta * (DBL)(j + 1);
  2058.         curr->ds = delta;
  2059.  
  2060.         Make_Colour(curr->od,  0.0, 0.0, 0.0);
  2061.         Make_Colour(curr->te,  0.0, 0.0, 0.0);
  2062.         Make_Colour(curr->te2, 0.0, 0.0, 0.0);
  2063. #ifdef EmissionMethodPatch
  2064.         Make_Colour(curr->Em2, 0.0, 0.0, 0.0);
  2065. #endif
  2066.  
  2067.         curr++;
  2068.       }
  2069.  
  2070.       remaining -= n;
  2071.  
  2072.       entry++;
  2073.     }
  2074.   }
  2075.  
  2076.   return(intervals);
  2077. }
  2078.  
  2079.  
  2080.  
  2081. /*****************************************************************************
  2082. *
  2083. * FUNCTION
  2084. *
  2085. *   intersect_spotlight
  2086. *
  2087. * INPUT
  2088. *
  2089. *   Ray    - current ray
  2090. *   Light  - current light source
  2091. *
  2092. * OUTPUT
  2093. *
  2094. *   d1, d2 - intersection depths
  2095. *
  2096. * RETURNS
  2097. *
  2098. *   int - TRUE, if hit
  2099. *
  2100. * AUTHOR
  2101. *
  2102. *   Dieter Bayer
  2103. *
  2104. * DESCRIPTION
  2105. *
  2106. *   Intersect a ray with the light cone of a spotlight.
  2107. *
  2108. * CHANGES
  2109. *
  2110. *   Nov 1994 : Creation.
  2111. *
  2112. ******************************************************************************/
  2113.  
  2114. static int intersect_spotlight(RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL  *d2)
  2115. {
  2116.   int viewpoint_is_in_cone;
  2117.   DBL a, b, c, d, m, l, l1, l2, t, t1, t2, k1, k2, k3, k4;
  2118.   VECTOR V1;
  2119.  
  2120.   /* Get cone's slope. Note that cos(falloff) is stored in Falloff! */
  2121.  
  2122.   a = acos(Light->Falloff);
  2123.  
  2124.   /* This only works for a < 180 degrees! */
  2125.  
  2126.   m = tan(a);
  2127.  
  2128.   m = 1.0 + Sqr(m);
  2129.  
  2130.   VSub(V1, Ray->Initial, Light->Center);
  2131.  
  2132.   VDot(k1, Ray->Direction, Light->Direction);
  2133.  
  2134.   VDot(k2, V1, Light->Direction);
  2135.  
  2136.   VLength(l, V1);
  2137.  
  2138.   if (l > EPSILON)
  2139.   {
  2140.     viewpoint_is_in_cone = (k2 / l >= Light->Falloff);
  2141.   }
  2142.   else
  2143.   {
  2144.     viewpoint_is_in_cone = FALSE;
  2145.   }
  2146.  
  2147.   if ((k1 <= 0.0) && (k2 < 0.0))
  2148.   {
  2149.     return (FALSE);
  2150.   }
  2151.  
  2152.   VDot(k3, V1, Ray->Direction);
  2153.  
  2154.   VDot(k4, V1, V1);
  2155.  
  2156.   a = 1.0 - Sqr(k1) * m;
  2157.  
  2158.   b = k3 - k1 * k2 * m;
  2159.  
  2160.   c = k4 - Sqr(k2) * m;
  2161.  
  2162.   if (a != 0.0)
  2163.   {
  2164.     d = Sqr(b) - a * c;
  2165.  
  2166.     if (d > EPSILON)
  2167.     {
  2168.       d = sqrt(d);
  2169.  
  2170.       t1 = (-b + d) / a;
  2171.       t2 = (-b - d) / a;
  2172.  
  2173.       if (t1 > t2)
  2174.       {
  2175.         t = t1; t1 = t2; t2 = t;
  2176.       }
  2177.  
  2178.       l1 = k2 + t1 * k1;
  2179.       l2 = k2 + t2 * k1;
  2180.  
  2181.       if ((l1 <= 0.0) && (l2 <= 0.0))
  2182.       {
  2183.         return (FALSE);
  2184.       }
  2185.  
  2186.       if ((l1 <= 0.0) || (l2 <= 0.0))
  2187.       {
  2188.         if (l1 <= 0.0)
  2189.         {
  2190.           if (viewpoint_is_in_cone)
  2191.           {
  2192.             t1 = 0.0;
  2193.             t2 = (t2 > 0.0) ? (t2) : (Max_Distance);
  2194.           }
  2195.           else
  2196.           {
  2197.             t1 = t2;
  2198.             t2 = Max_Distance;
  2199.           }
  2200.         }
  2201.         else
  2202.         {
  2203.           if (viewpoint_is_in_cone)
  2204.           {
  2205.             t2 = t1;
  2206.             t1 = 0.0;
  2207.           }
  2208.           else
  2209.           {
  2210.             t2 = Max_Distance;
  2211.           }
  2212.         }
  2213.       }
  2214.  
  2215.       *d1 = t1;
  2216.       *d2 = t2;
  2217.  
  2218.       return (TRUE);
  2219.     }
  2220.     else
  2221.     {
  2222.       if (d > -EPSILON)
  2223.       {
  2224.         if (viewpoint_is_in_cone)
  2225.         {
  2226.           *d1 = 0.0;
  2227.           *d2 = -b / a;
  2228.         }
  2229.         else
  2230.         {
  2231.           *d1 = -b / a;
  2232.           *d2 = Max_Distance;
  2233.         }
  2234.  
  2235.         return(TRUE);
  2236.       }
  2237.     }
  2238.   }
  2239.   else
  2240.   {
  2241.     if (viewpoint_is_in_cone)
  2242.     {
  2243.       *d1 = 0.0;
  2244.       *d2 = -c/b;
  2245.  
  2246.       return(TRUE);
  2247.     }
  2248.   }
  2249.  
  2250.   return (FALSE);
  2251. }
  2252.  
  2253.  
  2254.  
  2255. /*****************************************************************************
  2256. *
  2257. * FUNCTION
  2258. *
  2259. *   intersect_cylinderlight
  2260. *
  2261. * INPUT
  2262. *
  2263. *   Ray    - current ray
  2264. *   Light  - current light source
  2265. *
  2266. * OUTPUT
  2267. *
  2268. *   d1, d2 - intersection depths
  2269. *
  2270. * RETURNS
  2271. *
  2272. *   int - TRUE, if hit
  2273. *
  2274. * AUTHOR
  2275. *
  2276. *   Dieter Bayer
  2277. *
  2278. * DESCRIPTION
  2279. *
  2280. *   Intersect a ray with the light cylinder of a cylinderlight.
  2281. *
  2282. * CHANGES
  2283. *
  2284. *   Jan 1995 : Creation.
  2285. *
  2286. ******************************************************************************/
  2287.  
  2288. static int intersect_cylinderlight(RAY *Ray, LIGHT_SOURCE *Light, DBL *d1, DBL  *d2)
  2289. {
  2290.   DBL a, b, c, d, l1, l2, t, t1, t2, k1, k2, k3, k4;
  2291.   VECTOR V1;
  2292.  
  2293.   VSub(V1, Ray->Initial, Light->Center);
  2294.  
  2295.   VDot(k1, Ray->Direction, Light->Direction);
  2296.  
  2297.   VDot(k2, V1, Light->Direction);
  2298.  
  2299.   if ((k1 <= 0.0) && (k2 < 0.0))
  2300.   {
  2301.     return (FALSE);
  2302.   }
  2303.  
  2304.   a = 1.0 - Sqr(k1);
  2305.  
  2306.   if (a != 0.0)
  2307.   {
  2308.     VDot(k3, V1, Ray->Direction);
  2309.  
  2310.     VDot(k4, V1, V1);
  2311.  
  2312.     b = k3 - k1 * k2;
  2313.  
  2314.     c = k4 - Sqr(k2) - Sqr(Light->Falloff);
  2315.  
  2316.     d = Sqr(b) - a * c;
  2317.  
  2318.     if (d > EPSILON)
  2319.     {
  2320.       d = sqrt(d);
  2321.  
  2322.       t1 = (-b + d) / a;
  2323.       t2 = (-b - d) / a;
  2324.  
  2325.       if (t1 > t2)
  2326.       {
  2327.         t = t1; t1 = t2; t2 = t;
  2328.       }
  2329.  
  2330.       l1 = k2 + t1 * k1;
  2331.       l2 = k2 + t2 * k1;
  2332.  
  2333.       if ((l1 <= 0.0) && (l2 <= 0.0))
  2334.       {
  2335.         return (FALSE);
  2336.       }
  2337.  
  2338.       if ((l1 <= 0.0) || (l2 <= 0.0))
  2339.       {
  2340.         if (l1 <= 0.0)
  2341.         {
  2342.           t1 = 0.0;
  2343.         }
  2344.         else
  2345.         {
  2346.           t2 = (Max_Distance - k2) / k1;
  2347.         }
  2348.       }
  2349.  
  2350.       *d1 = t1;
  2351.       *d2 = t2;
  2352.  
  2353.       return (TRUE);
  2354.     }
  2355.   }
  2356.  
  2357.   return (FALSE);
  2358. }
  2359.  
  2360.  
  2361.  
  2362. /*****************************************************************************
  2363. *
  2364. * FUNCTION
  2365. *
  2366. *   Post_Media
  2367. *
  2368. * INPUT
  2369. *
  2370. * OUTPUT
  2371. *
  2372. * RETURNS
  2373. *
  2374. * AUTHOR
  2375. *
  2376. *   Dieter Bayer
  2377. *
  2378. * DESCRIPTION
  2379. *
  2380. *   Initialize media.
  2381. *
  2382. * CHANGES
  2383. *
  2384. *   Dec 1996 : Creation.
  2385. *
  2386. ******************************************************************************/
  2387.  
  2388. void Post_Media(IMEDIA *IMedia)
  2389. {
  2390.   int i;
  2391.   DBL t;
  2392.  
  2393.   if (IMedia == NULL)
  2394.   {
  2395.     return;
  2396.   }
  2397.   
  2398.   /* Get extinction coefficient. */
  2399.  
  2400.   CLinComb2(IMedia->Extinction, 1.0, IMedia->Absorption, IMedia->sc_ext, IMedia->Scattering);
  2401.  
  2402.   /* Determine used effects. */
  2403.  
  2404.   IMedia->use_absorption = (IMedia->Absorption[0] != 0.0) ||
  2405.                           (IMedia->Absorption[1] != 0.0) ||
  2406.                           (IMedia->Absorption[2] != 0.0);
  2407.  
  2408.   IMedia->use_emission = (IMedia->Emission[0] != 0.0) ||
  2409.                         (IMedia->Emission[1] != 0.0) ||
  2410.                         (IMedia->Emission[2] != 0.0);
  2411.  
  2412.   IMedia->use_scattering = (IMedia->Scattering[0] != 0.0) ||
  2413.                           (IMedia->Scattering[1] != 0.0) ||
  2414.                           (IMedia->Scattering[2] != 0.0);
  2415.  
  2416.   IMedia->use_extinction = IMedia->use_absorption || IMedia->use_scattering;
  2417.  
  2418.   IMedia->is_constant = (IMedia->Density == NULL);
  2419.  
  2420.   /* Init sample threshold array. */
  2421.  
  2422.   if (IMedia->Sample_Threshold != NULL)
  2423.   {
  2424.     POV_FREE(IMedia->Sample_Threshold);
  2425.   }
  2426.  
  2427.   /* Create list of thresholds for confidence test. */
  2428.  
  2429.   IMedia->Sample_Threshold = (DBL *)POV_MALLOC(IMedia->Max_Samples*sizeof(DBL), "sample threshold list");
  2430.  
  2431.   if (IMedia->Max_Samples > 1)
  2432.   {
  2433.     t = chdtri((DBL)(IMedia->Max_Samples-1), IMedia->Confidence);
  2434.  
  2435.     if (t > 0.0)
  2436.     {
  2437.       t = IMedia->Variance / t;
  2438.     }
  2439.     else
  2440.     {
  2441.       t = IMedia->Variance * EPSILON;
  2442.     }
  2443.  
  2444.     for (i = 0; i < IMedia->Max_Samples; i++)
  2445.     {
  2446.       IMedia->Sample_Threshold[i] = t * chdtri((DBL)(i+1), IMedia->Confidence);
  2447.  
  2448. /*
  2449.       fprintf(stderr, "threshold for n = %3d: %f\n", i+1, IMedia->Sample_Threshold[i]);
  2450. */
  2451.     }
  2452.   }
  2453.   else
  2454.   {
  2455.     IMedia->Sample_Threshold[0] = 0.0;
  2456.   }
  2457.  
  2458.   if (IMedia->Density != NULL)
  2459.   {
  2460.     Post_Pigment(IMedia->Density);
  2461.   }
  2462.   
  2463.   if (IMedia->Light_Group == NONE_GROUP) 
  2464.       IMedia->Light_Group = ALL_GROUP;
  2465.  
  2466.   Post_Media(IMedia->Next_Media);  
  2467. }
  2468.  
  2469.  
  2470.  
  2471. /*****************************************************************************
  2472. *
  2473. * FUNCTION
  2474. *
  2475. *   Create_Media
  2476. *
  2477. * INPUT
  2478. *
  2479. * OUTPUT
  2480. *
  2481. * RETURNS
  2482. *
  2483. *   IMEDIA * - created media
  2484. *
  2485. * AUTHOR
  2486. *
  2487. *   Dieter Bayer
  2488. *
  2489. * DESCRIPTION
  2490. *
  2491. *   Create a media.
  2492. *
  2493. * CHANGES
  2494. *
  2495. *   Dec 1994 : Creation.
  2496. *
  2497. ******************************************************************************/
  2498.  
  2499. IMEDIA *Create_Media()
  2500. {
  2501.   IMEDIA *New;
  2502.  
  2503.   New = (IMEDIA *)POV_MALLOC(sizeof(IMEDIA), "media");
  2504.  
  2505.   New->Type = ISOTROPIC_SCATTERING;
  2506.  
  2507.   New->Intervals      = 10;
  2508.   New->Min_Samples    = 1;
  2509.   New->Max_Samples    = 1;
  2510. #ifdef SampSpacingPatch
  2511.   New->Sample_Spacing = 0.0;
  2512.   New->Sample_Jitter  = 0.0;
  2513. #endif
  2514.   New->Eccentricity   = 0.0;
  2515.  
  2516.   Make_Colour(New->Absorption, 0.0, 0.0, 0.0);
  2517.   Make_Colour(New->Emission,   0.0, 0.0, 0.0);
  2518.   Make_Colour(New->Extinction, 0.0, 0.0, 0.0);
  2519.   Make_Colour(New->Scattering, 0.0, 0.0, 0.0);
  2520.  
  2521.   New->use_absorption = FALSE;
  2522.   New->use_emission   = FALSE;
  2523.   New->use_extinction = FALSE;
  2524.   New->use_scattering = FALSE;
  2525.  
  2526.   New->is_constant = FALSE;
  2527.  
  2528.   New->Light_Group            = NONE_GROUP;
  2529.   New->Invert_Light_Group    = FALSE;
  2530.  
  2531. #ifdef EmissionMethodPatch
  2532.   New->Emission_Method = 0;    /* JB */
  2533.   New->Emission_Extinction = 1.0; /* JB */
  2534. #endif
  2535.  
  2536.   New->sc_ext     = 1.0;
  2537.   New->Ratio      = 0.9;
  2538.   New->Confidence = 0.9;
  2539.   New->Variance   = 1.0 / 128.0;
  2540.  
  2541.   New->Sample_Threshold = NULL;
  2542.  
  2543.   New->Density = NULL;
  2544.  
  2545.   New->Next_Media = NULL;
  2546.  
  2547.   New->Sample_Method = 1;
  2548.   New->AA_Threshold = 0.1;
  2549.   New->AA_Level = 3;
  2550.   New->Jitter = 0.0;
  2551.  
  2552.   New->ignore_photons = FALSE;
  2553.  
  2554.   return (New);
  2555. }
  2556.  
  2557.  
  2558.  
  2559. /*****************************************************************************
  2560. *
  2561. * FUNCTION
  2562. *
  2563. *   Copy_Media
  2564. *
  2565. * INPUT
  2566. *
  2567. *   Old - media to copy
  2568. *
  2569. * OUTPUT
  2570. *
  2571. * RETURNS
  2572. *
  2573. *   IMEDIA * - new media
  2574. *
  2575. * AUTHOR
  2576. *
  2577. *   Dieter Bayer
  2578. *
  2579. * DESCRIPTION
  2580. *
  2581. *   Copy an media.
  2582. *
  2583. * CHANGES
  2584. *
  2585. *   Dec 1994 : Creation.
  2586. *
  2587. ******************************************************************************/
  2588.  
  2589. IMEDIA *Copy_Media(IMEDIA *Old)
  2590. {
  2591.   int i;
  2592.   IMEDIA *New, *First, *Previous, *Local_Media;
  2593.  
  2594.   Previous = First = NULL;
  2595.  
  2596.   if (Old != NULL)
  2597.   {
  2598.     for (Local_Media = Old; Local_Media != NULL; Local_Media = Local_Media->Next_Media)
  2599.     {
  2600.       New = Create_Media();
  2601.  
  2602.       *New = *Local_Media;
  2603.  
  2604.       if (Local_Media->Sample_Threshold != NULL)
  2605.       {
  2606.         if (New->Intervals > 0)
  2607.         {
  2608.           New->Sample_Threshold = (DBL *)POV_MALLOC(New->Intervals*sizeof(DBL), "sample threshold list");
  2609.  
  2610.           for (i = 0; i < New->Intervals; i++)
  2611.           {
  2612.             New->Sample_Threshold[i] =  Local_Media->Sample_Threshold[i];
  2613.           }
  2614.         }
  2615.       }
  2616.  
  2617.       New->Density = Copy_Pigment(Local_Media->Density);
  2618.  
  2619.       if (First == NULL)
  2620.       {
  2621.         First = New;
  2622.       }
  2623.  
  2624.       if (Previous != NULL)
  2625.       {
  2626.         Previous->Next_Media = New;
  2627.       }
  2628.  
  2629.       Previous = New;
  2630.     }
  2631.   }
  2632.  
  2633.   return(First);
  2634. }
  2635.  
  2636.  
  2637.  
  2638. /*****************************************************************************
  2639. *
  2640. * FUNCTION
  2641. *
  2642. *   Destroy_Media
  2643. *
  2644. * INPUT
  2645. *
  2646. *   IMedia - media to destroy
  2647. *
  2648. * OUTPUT
  2649. *
  2650. * RETURNS
  2651. *
  2652. * AUTHOR
  2653. *
  2654. *   Dieter Bayer
  2655. *
  2656. * DESCRIPTION
  2657. *
  2658. *   Destroy an media.
  2659. *
  2660. * CHANGES
  2661. *
  2662. *   Dec 1994 : Creation.
  2663. *
  2664. ******************************************************************************/
  2665.  
  2666. void Destroy_Media(IMEDIA *IMedia)
  2667. {
  2668.   IMEDIA *Local_Media, *Temp;
  2669.  
  2670.   if (IMedia != NULL)
  2671.   {
  2672.     Local_Media = IMedia;
  2673.  
  2674.     while (Local_Media != NULL)
  2675.     {
  2676.       if (Local_Media->Sample_Threshold != NULL)
  2677.       {
  2678.         POV_FREE(Local_Media->Sample_Threshold);
  2679.       }
  2680.  
  2681.       /* Note Destroy_Pigment also handles Density->Next */
  2682.       Destroy_Pigment(Local_Media->Density);
  2683.  
  2684.       Temp = Local_Media->Next_Media;
  2685.  
  2686.       POV_FREE(Local_Media);
  2687.  
  2688.       Local_Media = Temp;
  2689.     }
  2690.   }
  2691. }
  2692.  
  2693.  
  2694.  
  2695. /*****************************************************************************
  2696. *
  2697. * FUNCTION
  2698. *
  2699. *   compdoubles
  2700. *
  2701. * INPUT
  2702. *
  2703. *   in_a, in_b - Elements to compare
  2704. *
  2705. * OUTPUT
  2706. *
  2707. * RETURNS
  2708. *
  2709. *   int - result of comparison
  2710. *
  2711. * AUTHOR
  2712. *
  2713. *   Dieter Bayer
  2714. *
  2715. * DESCRIPTION
  2716. *
  2717. *   -
  2718. *
  2719. * CHANGES
  2720. *
  2721. *   Dec 1996 : Creation.
  2722. *
  2723. ******************************************************************************/
  2724.  
  2725. static int CDECL compdoubles(CONST void *in_a, CONST void *in_b)
  2726. {
  2727.   DBL *a, *b;
  2728.  
  2729.   a = (DBL *)in_a;
  2730.   b = (DBL *)in_b;
  2731.  
  2732.   if (*a < *b)
  2733.   {
  2734.     return (-1);
  2735.   }
  2736.   else
  2737.   {
  2738.     if (*a == *b)
  2739.     {
  2740.       return (0);
  2741.     }
  2742.     else
  2743.     {
  2744.       return (1);
  2745.     }
  2746.   }
  2747. }
  2748.  
  2749.  
  2750.  
  2751. /*****************************************************************************
  2752. *
  2753. * FUNCTION
  2754. *
  2755. *   Transform_Media
  2756. *
  2757. * INPUT
  2758. *
  2759. * OUTPUT
  2760. *
  2761. * RETURNS
  2762. *
  2763. * AUTHOR
  2764. *
  2765. *   Dieter Bayer
  2766. *
  2767. * DESCRIPTION
  2768. *
  2769. *   Transform media.
  2770. *
  2771. * CHANGES
  2772. *
  2773. *   Dec 1996 : Creation.
  2774. *
  2775. ******************************************************************************/
  2776.  
  2777. void Transform_Media(IMEDIA *IMedia, TRANSFORM *Trans)
  2778. {
  2779.   IMEDIA *Temp;
  2780.  
  2781.   if (IMedia != NULL)
  2782.   {
  2783.     for (Temp = IMedia; Temp != NULL; Temp = Temp->Next_Media)
  2784.     {
  2785.       Transform_Density(Temp->Density, Trans);
  2786.     }
  2787.   }
  2788. }
  2789.  
  2790. void Transform_Density(PIGMENT *Density, TRANSFORM *Trans)
  2791. {
  2792.   TPATTERN *Temp = (TPATTERN *)Density;
  2793.   
  2794.   while (Temp != NULL)
  2795.   {
  2796.     Transform_Tpattern(Temp, Trans);
  2797.     Temp = Temp->Next;
  2798.   }
  2799. }
  2800.  
  2801. void Set_Media_Light(LIGHT_SOURCE *light) {currLight=light;}
  2802.