]> git.mxchange.org Git - simgear.git/blob - simgear/scene/sky/clouds3d/SkyCloud.cpp
Melchior FRANZ: Make the clouds3d code valgrind clean
[simgear.git] / simgear / scene / sky / clouds3d / SkyCloud.cpp
1 //------------------------------------------------------------------------------
2 // File : SkyCloud.cpp
3 //------------------------------------------------------------------------------
4 // SkyWorks : Adapted from skyworks program writen by Mark J. Harris and
5 //                                              The University of North Carolina at Chapel Hill
6 //                                      : by J. Wojnaroski Sep 2002
7 //------------------------------------------------------------------------------
8 // Permission to use, copy, modify, distribute and sell this software and its 
9 // documentation for any purpose is hereby granted without fee, provided that 
10 // the above copyright notice appear in all copies and that both that copyright 
11 // notice and this permission notice appear in supporting documentation. 
12 // Binaries may be compiled with this software without any royalties or 
13 // restrictions. 
14 //
15 // The author(s) and The University of North Carolina at Chapel Hill make no 
16 // representations about the suitability of this software for any purpose. 
17 // It is provided "as is" without express or 
18 // implied warranty.
19 /**
20  * @file SkyCloud.cpp
21  * 
22  * Implementation of class SkyCloud.
23  */
24
25 // warning for truncation of template name for browse info
26 #pragma warning( disable : 4786)
27
28 #include <plib/ul.h>
29
30 #include "SkyCloud.hpp"
31 #include "SkyRenderableInstance.hpp" 
32 #include "SkyContext.hpp"
33 #include "SkyMaterial.hpp"
34 #include "SkyLight.hpp"
35 #include "SkyTextureManager.hpp"
36 #include "SkySceneManager.hpp"
37 #include <algorithm>
38
39 //! The version used for cloud archive files.
40 #define CLOUD_ARCHIVE_VERSION 0.1f
41
42 //------------------------------------------------------------------------------
43 // Static initialization
44 //------------------------------------------------------------------------------
45 SkyMaterial*  SkyCloud::s_pMaterial         = NULL;
46 SkyMaterial*  SkyCloud::s_pShadeMaterial    = NULL;
47 unsigned int  SkyCloud::s_iShadeResolution  = 32;
48 float         SkyCloud::s_rAlbedo           = 0.9f;
49 float         SkyCloud::s_rExtinction       = 80.0f;
50 float         SkyCloud::s_rTransparency     = exp(-s_rExtinction);
51 float         SkyCloud::s_rScatterFactor    = s_rAlbedo * s_rExtinction * SKY_INV_4PI;
52 float         SkyCloud::s_rSortAngleErrorTolerance = 0.8f;
53 float         SkyCloud::s_rSortSquareDistanceTolerance = 100;
54
55 //------------------------------------------------------------------------------
56 // Function               : SkyCloud::SkyCloud
57 // Description      : 
58 //------------------------------------------------------------------------------
59 /**
60  * @fn SkyCloud::SkyCloud()
61  * @brief Constructor.
62  */ 
63 SkyCloud::SkyCloud()
64 : SkyRenderable(),
65   _bUsePhaseFunction(true),
66   _vecLastSortViewDir(Vec3f::ZERO),
67   _vecLastSortCamPos(Vec3f::ZERO)
68 {
69   if (!s_pShadeMaterial)
70   {
71     s_pShadeMaterial = new SkyMaterial;
72     s_pShadeMaterial->SetAmbient(Vec4f(0.1f, 0.1f, 0.1f, 1));
73     s_pShadeMaterial->EnableDepthTest(false);
74     s_pShadeMaterial->SetBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
75     s_pShadeMaterial->EnableBlending(true);
76     s_pShadeMaterial->SetAlphaFunc(GL_GREATER);
77     s_pShadeMaterial->SetAlphaRef(0);
78     s_pShadeMaterial->EnableAlphaTest(true);
79     s_pShadeMaterial->SetColorMaterialMode(GL_DIFFUSE);
80     s_pShadeMaterial->EnableColorMaterial(true);
81     s_pShadeMaterial->EnableLighting(false);
82     s_pShadeMaterial->SetTextureApplicationMode(GL_MODULATE);
83   }
84   if (!s_pMaterial)
85   {
86     s_pMaterial = new SkyMaterial;
87     s_pMaterial->SetAmbient(Vec4f(0.3f, 0.3f, 0.3f, 1));
88     s_pMaterial->SetDepthMask(false);
89     s_pMaterial->SetBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
90     s_pMaterial->EnableBlending(true);
91     s_pMaterial->SetAlphaFunc(GL_GREATER);
92     s_pMaterial->SetAlphaRef(0);
93     s_pMaterial->EnableAlphaTest(true);
94     s_pMaterial->SetColorMaterialMode(GL_DIFFUSE);
95     s_pMaterial->EnableColorMaterial(true);
96     s_pMaterial->EnableLighting(false);
97     s_pMaterial->SetTextureApplicationMode(GL_MODULATE);
98     _CreateSplatTexture(32); // will assign the texture to both static materials
99   }
100 }
101
102
103 //------------------------------------------------------------------------------
104 // Function               : SkyCloud::~SkyCloud
105 // Description      : 
106 //------------------------------------------------------------------------------
107 /**
108  * @fn SkyCloud::~SkyCloud()
109  * @brief Destructor.
110  */ 
111 SkyCloud::~SkyCloud()
112 {
113 }
114
115
116 //------------------------------------------------------------------------------
117 // Function               : SkyCloud::Update
118 // Description      : 
119 //------------------------------------------------------------------------------
120 /**
121 * @fn SkyCloud::Update(const Camera &cam, SkyRenderableInstance* pInstance)
122 * @brief Currently does nothing.
123 */ 
124 SKYRESULT SkyCloud::Update(const Camera &cam, SkyRenderableInstance* pInstance)
125 {
126   return SKYRESULT_OK;
127 }
128
129
130
131 //------------------------------------------------------------------------------
132 // Function               : DrawQuad
133 // Description      : 
134 //------------------------------------------------------------------------------
135 /**
136  * @fn DrawQuad(Vec3f pos, Vec3f x, Vec3f y, Vec4f color)
137  * @brief Draw a quad.
138  */ 
139 inline void DrawQuad(Vec3f pos, Vec3f x, Vec3f y, Vec4f color)
140 {
141   glColor4fv(&(color.x));
142   Vec3f left  = pos;  left   -= y; 
143   Vec3f right = left; right  += x; 
144   left  -= x;
145   glTexCoord2f(0, 0); glVertex3fv(&(left.x));
146   glTexCoord2f(1, 0); glVertex3fv(&(right.x));
147   left  += y;  left  += y;
148   right += y;  right += y;
149   glTexCoord2f(1, 1); glVertex3fv(&(right.x));
150   glTexCoord2f(0, 1); glVertex3fv(&(left.x));
151 }
152
153
154 //------------------------------------------------------------------------------
155 // Function               : SkyCloud::Display
156 // Description      : 
157 //------------------------------------------------------------------------------
158 /**
159  * @fn SkyCloud::Display(const Camera &camera, SkyRenderableInstance *pInstance)
160  * @brief Renders the cloud.
161  * 
162  * The cloud is rendered by splatting the particles from back to front with respect
163  * to @a camera.  Since instances of clouds each have their own particles, which 
164  * are pre-transformed into world space, @a pInstance is not used.
165  *
166  * An alternative method is to store the particles untransformed, and transform the 
167  * camera and light into cloud space for rendering.  This is more complicated, 
168  * and not as straightforward.  Since I have to store the particles with each instance
169  * anyway, I decided to pre-transform them instead.
170  */ 
171 SKYRESULT SkyCloud::Display(const Camera &camera, SkyRenderableInstance *pInstance)
172 {
173   // copy the current camera
174   Camera cam(camera);
175   
176   // This cosine computation, along with the if() below, are an optimization.  The goal
177   // is to avoid sorting when it will make no visual difference.  This will be true when the 
178   // cloud particles are almost sorted for the current viewpoint.  This is the case most of the 
179   // time, since the viewpoint does not move very far in a single frame.  Each time we sort,
180   // we cache the current view direction.  Then, each time the cloud is displayed, if the 
181   // current view direction is very close to the current view direction (dot product is nearly 1)
182   // then we do not resort the particles.
183   float rCosAngleSinceLastSort = 
184       _vecLastSortViewDir * cam.ViewDir(); // dot product
185
186   float rSquareDistanceSinceLastSort = 
187     (cam.Orig - _vecLastSortCamPos).LengthSqr();
188
189   if (rCosAngleSinceLastSort < s_rSortAngleErrorTolerance || 
190       rSquareDistanceSinceLastSort > s_rSortSquareDistanceTolerance)
191   {
192     // compute the sort position for particles.
193     // don't just use the camera position -- if it is too far away from the cloud, then
194     // precision limitations may cause the STL sort to hang.  Instead, put the sort position
195     // just outside the bounding sphere of the cloud in the direction of the camera.
196     _vecSortPos = -cam.ViewDir();
197     _vecSortPos *= (1.1 * _boundingBox.GetRadius());
198     _vecSortPos += _boundingBox.GetCenter();
199     
200     // sort the particles from back to front wrt the camera position.
201     _SortParticles(cam.ViewDir(), _vecSortPos, SKY_CLOUD_SORT_TOWARD);
202
203     //_vecLastSortViewDir = GLVU::GetCurrent()->GetCurrentCam()->ViewDir();
204     //_vecLastSortCamPos = GLVU::GetCurrent()->GetCurrentCam()->Orig;
205     _vecLastSortViewDir = cam.ViewDir();
206     _vecLastSortCamPos = cam.Orig;
207   } 
208
209   // set the material state / properties that clouds use for rendering:
210   // Enables blending, with blend func (ONE, ONE_MINUS_SRC_ALPHA).
211   // Enables alpha test to discard completely transparent fragments.
212   // Disables depth test.
213   // Enables texturing, with modulation, and the texture set to the shared splat texture.
214   s_pMaterial->Activate();
215   
216   Vec4f color;
217   Vec3f eyeDir;
218
219   // Draw the particles using immediate mode.
220   glBegin(GL_QUADS);
221
222   int i = 0;
223   for (ParticleIterator iter = _particles.begin(); iter != _particles.end(); iter++)
224   {
225     i++;
226     SkyCloudParticle *p = *iter;
227
228     // Start with ambient light
229     color = p->GetBaseColor();
230     
231     if (_bUsePhaseFunction) // use the phase function for anisotropic scattering.
232     { 
233       eyeDir = cam.Orig;
234       eyeDir -= p->GetPosition();     
235       eyeDir.Normalize();
236       float pf;
237           
238       // add the color contribution to this particle from each light source, modulated by
239       // the phase function.  See _PhaseFunction() documentation for details.
240       for (int i = 0; i < p->GetNumLitColors(); i++)
241       {
242         pf = _PhaseFunction(_lightDirections[i], eyeDir);
243         // expand this to avoid temporary vector creation in the inner loop
244         color.x += p->GetLitColor(i).x * pf;
245         color.y += p->GetLitColor(i).y * pf;
246         color.z += p->GetLitColor(i).z * pf;
247       }
248     }
249     else  // just use isotropic scattering instead.
250     {   
251       for (int i = 0; i < (*iter)->GetNumLitColors(); ++i)
252       {
253         color += p->GetLitColor(i);
254       }
255     }
256     
257     // Set the transparency independently of the colors
258     color.w = 1 - s_rTransparency;
259     
260     // draw the particle as a textured billboard.
261     DrawQuad((*iter)->GetPosition(), cam.X * p->GetRadius(), cam.Y * p->GetRadius(), color);
262   }
263   glEnd();
264  
265   return SKYRESULT_OK;
266 }
267
268
269 //------------------------------------------------------------------------------
270 // Function               : SkyCloud::DisplaySplit
271 // Description      : 
272 //------------------------------------------------------------------------------
273 /**
274  * @fn SkyCloud::DisplaySplit(const Camera &camera, const Vec3f &vecSplitPoint, bool bBackHalf, SkyRenderableInstance *pInstance)
275  * @brief The same as Display(), except it displays only the particles in front of or behind the split point.
276  * 
277  * This is used to render clouds into two impostor images for displaying clouds that contain objects.
278  *
279  * @see SkyRenderableInstanceCloud
280  */ 
281 SKYRESULT SkyCloud::DisplaySplit(const Camera           &camera, 
282                                  const Vec3f            &vecSplitPoint,
283                                  bool                   bBackHalf,
284                                  SkyRenderableInstance  *pInstance /* = NULL */)
285 {
286   // copy the current camera
287   Camera cam(camera);
288
289   Vec3f vecCloudSpaceSplit = vecSplitPoint;
290
291   if (bBackHalf) // only sort when rendering the back half.  Reuse sort for front half.
292   {
293     // compute the sort position for particles.
294     // don't just use the camera position -- if it is too far away from the cloud, then
295     // precision limitations may cause the STL sort to hang.  Instead, put the sort position
296     // just outside the bounding sphere of the cloud in the direction of the camera.
297     _vecSortPos = -cam.ViewDir();
298     _vecSortPos *= (1.1 * _boundingBox.GetRadius());
299     _vecSortPos += _boundingBox.GetCenter();
300     
301     // sort the particles from back to front wrt the camera position.
302     _SortParticles(cam.ViewDir(), _vecSortPos, SKY_CLOUD_SORT_TOWARD);
303
304     // we can't use the view direction optimization when the cloud is split, or we get a lot
305     // of popping of objects in and out of cloud cover.  For consistency, though, we need to update
306     // the cached sort direction, since we just sorted the particles.
307     // _vecLastSortViewDir = GLVU::GetCurrent()->GetCurrentCam()->ViewDir(); 
308     // _vecLastSortCamPos = GLVU::GetCurrent()->GetCurrentCam()->Orig;
309
310     // compute the split distance.
311     vecCloudSpaceSplit  -= _vecSortPos;
312     _rSplitDistance = vecCloudSpaceSplit * cam.ViewDir();
313   }
314
315   // set the material state / properties that clouds use for rendering:
316   // Enables blending, with blend func (ONE, ONE_MINUS_SRC_ALPHA).
317   // Enables alpha test to discard completely transparent fragments.
318   // Disables depth test.
319   // Enables texturing, with modulation, and the texture set to the shared splat texture.
320   s_pMaterial->Activate();
321   
322   Vec4f color;
323   Vec3f eyeDir;
324
325   // Draw the particles using immediate mode.
326   glBegin(GL_QUADS);
327
328   // if bBackHalf is false, then we just continue where we left off.  If it is true, we
329   // reset the iterator to the beginning of the sorted list.
330   static ParticleIterator iter;
331   if (bBackHalf) 
332     iter = _particles.begin();
333   
334   // iterate over the particles and render them.
335   for (; iter != _particles.end(); ++iter)
336   {
337     SkyCloudParticle *p = *iter;
338
339     if (bBackHalf && (p->GetSquareSortDistance() < _rSplitDistance))
340       break;
341
342     // Start with ambient light
343     color = p->GetBaseColor();
344     
345     if (_bUsePhaseFunction) // use the phase function for anisotropic scattering.
346     {
347       eyeDir = cam.Orig;
348       eyeDir -= p->GetPosition();
349       eyeDir.Normalize();
350       float pf;
351           
352       // add the color contribution to this particle from each light source, modulated by
353       // the phase function.  See _PhaseFunction() documentation for details.
354       for (int i = 0; i < p->GetNumLitColors(); i++)
355       {
356         pf = _PhaseFunction(_lightDirections[i], eyeDir);
357         // expand this to avoid temporary vector creation in the inner loop
358         color.x += p->GetLitColor(i).x * pf;
359         color.y += p->GetLitColor(i).y * pf;
360         color.z += p->GetLitColor(i).z * pf;
361       }
362     }
363     else  // just use isotropic scattering instead.
364     {      
365       for (int i = 0; i < p->GetNumLitColors(); ++i)
366       {
367         color += p->GetLitColor(i);
368       }
369     }
370
371     // set the transparency independently of the colors.
372     color.w = 1 - s_rTransparency;
373     
374     // draw the particle as a textured billboard.
375     DrawQuad((*iter)->GetPosition(), cam.X * p->GetRadius(), cam.Y * p->GetRadius(), color);
376   }
377   glEnd();
378   
379   return SKYRESULT_OK;
380 }
381
382 //------------------------------------------------------------------------------
383 // Function               : SkyCloud::Illuminate
384 // Description      : 
385 //------------------------------------------------------------------------------
386 /**
387  * @fn SkyCloud::Illuminate(SkyLight *pLight, SkyRenderableInstance* pInstance, bool bReset)
388  * @brief Compute the illumination of the cloud by the lightsource @a pLight
389  * 
390  * This method uses graphics hardware to compute multiple forward scattering at each cloud
391  * in the cloud of light from the directional light source @a pLight.  The algorithm works
392  * by successively subtracting "light" from an initially white (fully lit) frame buffer by
393  * using hardware blending and read back.  The method stores the illumination from each light
394  * source passed to it separately at each particle, unless @a bReset is true, in which case
395  * the lists of illumination in the particles are reset before the lighting is computed.
396  * 
397  */ 
398 SKYRESULT SkyCloud::Illuminate(SkyLight *pLight, SkyRenderableInstance* pInstance, bool bReset)
399 {
400   int iOldVP[4];        
401
402   glGetIntegerv(GL_VIEWPORT, iOldVP);
403   glViewport(0, 0, s_iShadeResolution, s_iShadeResolution);
404
405   Vec3f vecDir(pLight->GetDirection());
406
407   // if this is the first pass through the lights, reset will be true, and the cached light 
408   // directions should be updated.  Light directions are cached in cloud space to accelerate 
409   // computation of the phase function, which depends on light direction and view direction.
410   if (bReset)
411     _lightDirections.clear();
412   _lightDirections.push_back(vecDir); // cache the (unit-length) light direction
413
414   // compute the light/sort position for particles from the light direction.
415   // don't just use the camera position -- if it is too far away from the cloud, then
416   // precision limitations may cause the STL sort to hang.  Instead, put the sort position
417   // just outside the bounding sphere of the cloud in the direction of the camera.
418   Vec3f vecLightPos(vecDir);
419   vecLightPos *= (1.1*_boundingBox.GetRadius());
420   vecLightPos += _boundingBox.GetCenter();
421
422   // Set up a camera to look at the cloud from the light position.  Since the sun is an infinite
423   // light source, this camera will use an orthographic projection tightly fit to the bounding
424   // sphere of the cloud.  
425   Camera cam;
426   
427   // Avoid degenerate camera bases.
428   Vec3f vecUp(0, 1, 0);
429   if (fabs(vecDir * vecUp) - 1 < 1e-6) // check that the view and up directions are not parallel.
430     vecUp.Set(1, 0, 0);
431   
432   cam.LookAt(vecLightPos, _boundingBox.GetCenter(), vecUp);
433
434   // sort the particles away from the light source.
435   _SortParticles(cam.ViewDir(), vecLightPos, SKY_CLOUD_SORT_AWAY);
436
437   // projected dist to cntr along viewdir
438   float DistToCntr = (_boundingBox.GetCenter() - vecLightPos) * cam.ViewDir();
439   
440   // calc tight-fitting near and far distances for the orthographic frustum
441   float rNearDist = DistToCntr - _boundingBox.GetRadius();
442   float rFarDist  = DistToCntr + _boundingBox.GetRadius();
443
444   // set the modelview matrix from this camera.
445   glMatrixMode(GL_MODELVIEW);
446   glPushMatrix();
447   float M[16];
448  
449  
450   cam.GetModelviewMatrix(M);
451   glLoadMatrixf(M);
452   
453   // switch to parallel projection
454   glMatrixMode(GL_PROJECTION);
455   glPushMatrix();
456   glLoadIdentity();
457   glOrtho(-_boundingBox.GetRadius(), _boundingBox.GetRadius(),
458           -_boundingBox.GetRadius(), _boundingBox.GetRadius(),
459           rNearDist, rFarDist);
460   
461   // set the material state / properties that clouds use for shading:
462   // Enables blending, with blend func (ONE, ONE_MINUS_SRC_ALPHA).
463   // Enables alpha test to discard completely transparent fragments.
464   // Disables depth test.
465   // Enables texturing, with modulation, and the texture set to the shared splat texture.
466   s_pShadeMaterial->Activate();
467
468   // these are used for projecting the particle position to determine where to read pixels.
469   double MM[16], PM[16];
470   int VP[4] = { 0, 0, s_iShadeResolution, s_iShadeResolution };
471   glGetDoublev(GL_MODELVIEW_MATRIX, MM);
472   glGetDoublev(GL_PROJECTION_MATRIX, PM);
473     
474   // initialize back buffer to all white -- modulation darkens areas where cloud particles
475   // absorb light, and lightens it where they scatter light in the forward direction.
476   glClearColor(1, 1, 1, 1);
477   glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
478
479   float rPixelsPerLength = s_iShadeResolution / (2 * _boundingBox.GetRadius());
480
481   // the solid angle over which we will sample forward-scattered light.
482   float rSolidAngle = 0.09;
483   int i = 0;
484   int iNumFailed = 0;
485   for (ParticleIterator iter = _particles.begin(); iter != _particles.end(); ++iter, ++i)
486   {
487     Vec3f vecParticlePos = (*iter)->GetPosition();
488
489     Vec3f vecOffset(vecLightPos);
490     vecOffset -= vecParticlePos;
491    
492     // compute the pixel area to read back in order to integrate the illumination of the particle
493     // over a constant solid angle.
494     float rDistance  = fabs(cam.ViewDir() * vecOffset) - rNearDist;
495
496     float rArea      = rSolidAngle * rDistance * rDistance;
497     int   iPixelDim  = sqrt(rArea) * rPixelsPerLength;
498     int   iNumPixels = iPixelDim * iPixelDim;
499     if (iNumPixels < 1) 
500     {
501       iNumPixels = 1;
502       iPixelDim = 1;
503     }
504
505     // the scale factor to convert the read back pixel colors to an average illumination of the area.
506     float rColorScaleFactor = rSolidAngle / (iNumPixels * 255.0f);
507
508     unsigned char *c = new unsigned char[4 * iNumPixels];
509     
510     Vec3d vecWinPos;
511     
512     // find the position in the buffer to which the particle position projects.
513     if (!gluProject(vecParticlePos.x, vecParticlePos.y, vecParticlePos.z, 
514                     MM, PM, VP, 
515                     &(vecWinPos.x), &(vecWinPos.y), &(vecWinPos.z)))
516     {
517       FAIL_RETURN_MSG(SKYRESULT_FAIL, 
518                       "Error: SkyCloud::Illuminate(): failed to project particle position.");
519     }
520    
521     // offset the projected window position by half the size of the readback region.
522     vecWinPos.x -= 0.5 * iPixelDim;
523     if (vecWinPos.x < 0) vecWinPos.x = 0;
524     vecWinPos.y -= 0.5 * iPixelDim;
525     if (vecWinPos.y < 0) vecWinPos.y = 0;
526     
527     // read back illumination of this particle from the buffer.
528     glReadBuffer(GL_BACK);
529     glReadPixels(vecWinPos.x, vecWinPos.y, iPixelDim, iPixelDim, GL_RGBA, GL_UNSIGNED_BYTE, c);
530     
531     // scattering coefficient vector.
532     Vec4f vecScatter(s_rScatterFactor, s_rScatterFactor, s_rScatterFactor, 1);
533
534     // add up the read back pixels (only need one component -- its grayscale)
535     int iSum = 0;
536     for (int k = 0; k < 4 * iNumPixels; k+=4)
537       iSum += c[k];
538     delete [] c;
539     
540     // compute the amount of light scattered to this particle by particles closer to the light.
541     // this is the illumination over the solid angle that we measured (using glReadPixels) times
542     // the scattering coefficient (vecScatter);
543     Vec4f vecScatteredAmount(iSum * rColorScaleFactor, 
544                              iSum * rColorScaleFactor, 
545                              iSum * rColorScaleFactor, 
546                              1 - s_rTransparency);
547     vecScatteredAmount &= vecScatter;
548
549     // the color of th particle (iter) contributed by this light source (pLight) is the 
550     // scattered light from the part of the cloud closer to the light, times the diffuse color
551     // of the light source.  The alpha is 1 - the uniform transparency of all particles (modulated
552     // by the splat texture).
553     Vec4f vecColor = vecScatteredAmount;
554     vecColor      &= pLight->GetDiffuse();  
555     vecColor.w     = 1 - s_rTransparency;
556     
557     // add this color to the list of lit colors for the particle.  The contribution from each light
558     // is kept separate because the phase function we apply at runtime depends on the light vector
559     // for each light source separately.  This view-dependent effect is impossible without knowing 
560     // the amount of light contributed for each light.  This, of course, assumes the clouds will
561     // be lit by a reasonably small number of lights (The sun plus some simulation of light reflected
562     // from the sky and / or ground.) This technique works very well for simulating anisotropic 
563     // illumination by skylight.
564     if (bReset)
565     {
566       (*iter)->SetBaseColor(s_pMaterial->GetAmbient());  
567       (*iter)->ClearLitColors();
568       (*iter)->AddLitColor(vecColor);
569     }
570     else
571     {
572       (*iter)->AddLitColor(vecColor);
573     }
574     
575     // the following computation (scaling of the scattered amount by the phase function) is done
576     // after the lit color is stored so we don't add the scattering to this particle twice.
577     vecScatteredAmount *= 1.5; // rayleigh scattering phase function for angle of zero or 180 = 1.5!
578     
579     // clamp the color
580     if (vecScatteredAmount.x > 1) vecScatteredAmount.x = 1;
581     if (vecScatteredAmount.y > 1) vecScatteredAmount.y = 1;
582     if (vecScatteredAmount.z > 1) vecScatteredAmount.z = 1;
583     vecScatteredAmount.w = 1 - s_rTransparency;
584     
585     vecScatteredAmount.x = 0.50;  vecScatteredAmount.y = 0.60;  vecScatteredAmount.z = 0.70; 
586
587     // Draw the particle as a texture billboard.  Use the scattered light amount as the color to 
588     // simulate forward scattering of light by this particle.
589     glBegin(GL_QUADS);
590     DrawQuad(vecParticlePos, cam.X * (*iter)->GetRadius(), cam.Y * (*iter)->GetRadius(), vecScatteredAmount);
591     glEnd();
592
593     //glutSwapBuffers(); // Uncomment this swap buffers to visualize cloud illumination computation.
594   }
595   
596   // Note: here we could optionally store the current back buffer as a shadow image
597   // to be projected from the light position onto the scene.  This way we can have clouds shadow
598   // the environment.
599   
600   // restore matrix stack and viewport.
601   glMatrixMode(GL_PROJECTION);
602   glPopMatrix();
603   glMatrixMode(GL_MODELVIEW);
604   glPopMatrix();
605   glViewport(iOldVP[0], iOldVP[1], iOldVP[2], iOldVP[3]);
606
607   return SKYRESULT_OK; 
608 }
609
610
611 //------------------------------------------------------------------------------
612 // Function               : SkyCloud::CopyBoundingVolume
613 // Description      : 
614 //------------------------------------------------------------------------------
615 /**
616  * @fn SkyCloud::CopyBoundingVolume() const
617  * @brief Returns a new copy of the SkyMinMaxBox for this cloud.
618  */ 
619 SkyMinMaxBox* SkyCloud::CopyBoundingVolume() const
620 {
621   SkyMinMaxBox *pBox = new SkyMinMaxBox();
622   pBox->SetMax(_boundingBox.GetMax());
623   pBox->SetMin(_boundingBox.GetMin());
624   return pBox; 
625 }
626
627 SKYRESULT SkyCloud::Load(const SkyArchive &archive, 
628                          float rScale, /* = 1.0f */ 
629                          double latitude, double longitude )
630 {
631   unsigned int iNumParticles;
632   Vec3f vecCenter = Vec3f::ZERO;
633   //Vec3f vecCenter;
634   //float rRadius;
635   //archive.FindVec3f("CldCenter", &vecCenter);
636   //archive.FindFloat32("CldRadius", &rRadius);
637
638   //_boundingBox.SetMin(vecCenter - Vec3f(rRadius, rRadius, rRadius));
639   //_boundingBox.SetMax(vecCenter + Vec3f(rRadius, rRadius, rRadius));
640
641   archive.FindUInt32("CldNumParticles", &iNumParticles);
642   //if (!bLocal)
643     archive.FindVec3f("CldCenter", &vecCenter);
644
645   Vec3f *pParticlePositions = new Vec3f[iNumParticles];
646   float *pParticleRadii     = new float[iNumParticles];
647   Vec4f *pParticleColors    = new Vec4f[iNumParticles];
648
649   unsigned int iNumBytes;
650   archive.FindData("CldParticlePositions", ANY_TYPE, (void**const)&pParticlePositions, &iNumBytes);
651   archive.FindData("CldParticleRadii",     ANY_TYPE, (void**const)&pParticleRadii,     &iNumBytes);
652   archive.FindData("CldParticleColors",    ANY_TYPE, (void**const)&pParticleColors,    &iNumBytes);
653   
654   for (unsigned int i = 0; i < iNumParticles; ++i)
655   {
656     SkyCloudParticle *pParticle = new SkyCloudParticle((pParticlePositions[i] + vecCenter) * rScale,
657                                                        pParticleRadii[i] * rScale,
658                                                        pParticleColors[i]);
659     _boundingBox.AddPoint(pParticle->GetPosition());
660     
661     _particles.push_back(pParticle);
662   }
663   // this is just a bad hack to align cloud field from skyworks with local horizon at KSFO
664   // this "almost" works not quite the right solution okay to get some up and running
665   // we need to develop our own scheme for loading and positioning clouds
666   Mat33f R;
667   Vec3f  moveit;
668
669   R.Set( 0, 1, 0,
670                          1, 0, 0,
671                          0, 0, 1);
672   // clouds sit in the y-z plane and x-axis is the vertical cloud height
673   Rotate( R ); 
674   
675 // rotate the cloud field about the fgfs z-axis based on initial longitude
676 float ex = 0.0;
677 float ey = 0.0;
678 float ez = 1.0;
679 float phi = longitude / 57.29578;
680 float one_min_cos = 1 - cos(phi);
681   
682 R.Set(
683 cos(phi) + one_min_cos*ex*ex, one_min_cos*ex*ey - ez*sin(phi), one_min_cos*ex*ez + ey*sin(phi),
684 one_min_cos*ex*ey + ez*sin(phi), cos(phi) + one_min_cos*ey*ey, one_min_cos*ey*ez - ex*sin(phi),
685 one_min_cos*ex*ez - ey*sin(phi), one_min_cos*ey*ez + ex*sin(phi), cos(phi) + one_min_cos*ez*ez );
686                         
687 Rotate( R );
688
689 // okay now that let's rotate about a vector for latitude where longitude forms the 
690 // components of a unit vector in the x-y plane
691 ex = sin( longitude  / 57.29578 );
692 ey = -cos( longitude  / 57.29578 );
693 ez = 0.0;
694 phi = latitude / 57.29578;
695 one_min_cos = 1 - cos(phi);
696
697 R.Set(
698 cos(phi) + one_min_cos*ex*ex, one_min_cos*ex*ey - ez*sin(phi), one_min_cos*ex*ez + ey*sin(phi),
699 one_min_cos*ex*ey + ez*sin(phi), cos(phi) + one_min_cos*ey*ey, one_min_cos*ey*ez - ex*sin(phi),
700 one_min_cos*ex*ez - ey*sin(phi), one_min_cos*ey*ez + ex*sin(phi), cos(phi) + one_min_cos*ez*ez );
701                         
702 Rotate( R );
703 // need to calculate an offset to place the clouds at ~3000 feet MSL  ATM this is an approximation 
704 // to move the clouds to some altitude above sea level. At some locations this could be underground
705 // will need a better scheme to position clouds per user preferences
706 float cloud_level_msl = 3000.0f;
707
708 float x_offset = ex * cloud_level_msl;
709 float y_offset = ey * cloud_level_msl; 
710 float z_offset = cloud_level_msl * 0.5;
711 moveit.Set( x_offset, y_offset, z_offset  );
712   
713   Translate( moveit );
714   
715   return SKYRESULT_OK;
716 }
717
718
719 //------------------------------------------------------------------------------
720 // Function               : SkyCloud::Save
721 // Description      : 
722 //------------------------------------------------------------------------------
723 /**
724 * @fn SkyCloud::Save(SkyArchive &archive) const
725 * @brief Saves the cloud data to @a archive.
726
727 * @todo <WRITE EXTENDED SkyCloud::Save FUNCTION DOCUMENTATION>
728 */ 
729 SKYRESULT SkyCloud::Save(SkyArchive &archive) const
730 {
731   SkyArchive myArchive("Cloud");
732   //myArchive.AddVec3f("CldCenter", _center);
733   //myArchive.AddFloat32("CldRadius", _boundingBox.GetRadius());
734   myArchive.AddUInt32("CldNumParticles", _particles.size());
735   
736   // make temp arrays
737   Vec3f *pParticlePositions = new Vec3f[_particles.size()];
738   float *pParticleRadii     = new float[_particles.size()];
739   Vec4f *pParticleColors    = new Vec4f[_particles.size()];
740   
741   unsigned int i = 0;
742
743   for (ParticleConstIterator iter = _particles.begin(); iter != _particles.end(); ++iter, ++i)
744   {
745     pParticlePositions[i] = (*iter)->GetPosition(); // position around origin
746     pParticleRadii[i]     = (*iter)->GetRadius();
747     pParticleColors[i]    = (*iter)->GetBaseColor();
748   }
749   
750   myArchive.AddData("CldParticlePositions", 
751                     ANY_TYPE, 
752                     pParticlePositions, 
753                     sizeof(Vec3f), 
754                     _particles.size());
755   
756   myArchive.AddData("CldParticleRadii", 
757                     ANY_TYPE, 
758                     pParticleRadii, 
759                     sizeof(float), 
760                     _particles.size());
761   
762   myArchive.AddData("CldParticleColors", 
763                     ANY_TYPE, 
764                     pParticleColors, 
765                     sizeof(Vec3f), 
766                     _particles.size());
767   
768   archive.AddArchive(myArchive);
769
770   return SKYRESULT_OK;
771 }
772
773
774 //------------------------------------------------------------------------------
775 // Function               : SkyCloud::Rotate
776 // Description      : 
777 //------------------------------------------------------------------------------
778 /**
779  * @fn SkyCloud::Rotate(const Mat33f& rot)
780  * @brief @todo <WRITE BRIEF SkyCloud::Rotate DOCUMENTATION>
781  * 
782  * @todo <WRITE EXTENDED SkyCloud::Rotate FUNCTION DOCUMENTATION>
783  */ 
784 void SkyCloud::Rotate(const Mat33f& rot)
785 {
786   _boundingBox.Clear();
787   for (int i = 0; i < _particles.size(); ++i)
788   {
789     _particles[i]->SetPosition(rot * _particles[i]->GetPosition());
790     _boundingBox.AddPoint(_particles[i]->GetPosition());
791   }
792 }
793
794
795 //------------------------------------------------------------------------------
796 // Function               : SkyCloud::Translate
797 // Description      : 
798 //------------------------------------------------------------------------------
799 /**
800  * @fn SkyCloud::Translate(const Vec3f& trans)
801  * @brief @todo <WRITE BRIEF SkyCloud::Translate DOCUMENTATION>
802  * 
803  * @todo <WRITE EXTENDED SkyCloud::Translate FUNCTION DOCUMENTATION>
804  */ 
805 void SkyCloud::Translate(const Vec3f& trans)
806 {
807   for (int i = 0; i < _particles.size(); ++i)
808   {
809     _particles[i]->SetPosition(_particles[i]->GetPosition() + trans);
810   }
811   _boundingBox.SetMax(_boundingBox.GetMax() + trans);
812   _boundingBox.SetMin(_boundingBox.GetMin() + trans);
813 }
814
815
816 //------------------------------------------------------------------------------
817 // Function               : SkyCloud::Scale
818 // Description      : 
819 //------------------------------------------------------------------------------
820 /**
821  * @fn SkyCloud::Scale(const float scale)
822  * @brief @todo <WRITE BRIEF SkyCloud::Scale DOCUMENTATION>
823  * 
824  * @todo <WRITE EXTENDED SkyCloud::Scale FUNCTION DOCUMENTATION>
825  */ 
826 void SkyCloud::Scale(const float scale)
827 {
828   _boundingBox.Clear();
829   for (int i = 0; i < _particles.size(); ++i)
830   {
831     _particles[i]->SetPosition(_particles[i]->GetPosition() * scale);
832     _boundingBox.AddPoint(_particles[i]->GetPosition());
833   }
834 }
835
836
837 //------------------------------------------------------------------------------
838 // Function               : SkyCloud::_SortParticles
839 // Description      : 
840 //------------------------------------------------------------------------------
841 /**
842  * @fn SkyCloud::_SortParticles(const Vec3f& vecViewDir, const Vec3f& sortPoint, SortDirection dir)
843  * @brief Sorts the cloud particles in the direction specified by @a dir.
844  * 
845  * @vecSortPoint is assumed to already be transformed into the basis space of the cloud.
846  */ 
847 void SkyCloud::_SortParticles(const Vec3f& vecViewDir,
848                               const Vec3f& vecSortPoint, 
849                               SortDirection dir)
850 {
851   Vec3f partPos;
852         for (int i = 0; i < _particles.size(); ++i)
853         {
854                 partPos = _particles[i]->GetPosition();
855                 partPos -= vecSortPoint;
856                 _particles[i]->SetSquareSortDistance(partPos * vecViewDir);//partPos.LengthSqr());
857         }
858
859   switch (dir)
860   {
861   case SKY_CLOUD_SORT_TOWARD:
862     std::sort(_particles.begin(), _particles.end(), _towardComparator);
863     break;
864   case SKY_CLOUD_SORT_AWAY:
865     std::sort(_particles.begin(), _particles.end(), _awayComparator);
866     break;
867   default:
868     break;
869   }
870 }
871
872
873
874 //------------------------------------------------------------------------------
875 // Function               : EvalHermite
876 // Description      : 
877 //------------------------------------------------------------------------------
878 /**
879  * EvalHermite(float pA, float pB, float vA, float vB, float u)
880  * @brief Evaluates Hermite basis functions for the specified coefficients.
881  */ 
882 inline float EvalHermite(float pA, float pB, float vA, float vB, float u)
883 {
884   float u2=(u*u), u3=u2*u;
885   float B0 = 2*u3 - 3*u2 + 1;
886   float B1 = -2*u3 + 3*u2;
887   float B2 = u3 - 2*u2 + u;
888   float B3 = u3 - u;
889   return( B0*pA + B1*pB + B2*vA + B3*vB );
890 }
891
892 // NORMALIZED GAUSSIAN INTENSITY MAP (N must be a power of 2)
893
894 //------------------------------------------------------------------------------
895 // Function               : CreateGaussianMap
896 // Description      : 
897 //------------------------------------------------------------------------------
898 /**
899  * CreateGaussianMap(int N)
900  * 
901  * Creates a 2D gaussian image using a hermite surface.
902  */ 
903 unsigned char* CreateGaussianMap(int N)
904 {
905   float *M = new float[2*N*N];
906   unsigned char *B = new unsigned char[4*N*N];
907   float X,Y,Y2,Dist;
908   float Incr = 2.0f/N;
909   int i=0;  
910   int j = 0;
911   Y = -1.0f;
912   for (int y=0; y<N; y++, Y+=Incr)
913   {
914     Y2=Y*Y;
915     X = -1.0f;
916     for (int x=0; x<N; x++, X+=Incr, i+=2, j+=4)
917     {
918       Dist = (float)sqrt(X*X+Y2);
919       if (Dist>1) Dist=1;
920       M[i+1] = M[i] = EvalHermite(0.4f,0,0,0,Dist);// * (1 - noise);
921       B[j+3] = B[j+2] = B[j+1] = B[j] = (unsigned char)(M[i] * 255);
922     }
923   }
924   SAFE_DELETE_ARRAY(M);
925   return(B);
926 }              
927
928
929 //------------------------------------------------------------------------------
930 // Function               : SkyCloud::_CreateSplatTexture
931 // Description      : 
932 //------------------------------------------------------------------------------
933 /**
934  * @fn SkyCloud::_CreateSplatTexture(unsigned int iResolution)
935  * @brief Creates the texture map used for cloud particles.
936  */ 
937 void SkyCloud::_CreateSplatTexture(unsigned int iResolution)
938 {
939   unsigned char *splatTexture = CreateGaussianMap(iResolution);
940   SkyTexture texture;
941   TextureManager::InstancePtr()->Create2DTextureObject(texture, iResolution, iResolution, 
942                                                        GL_RGBA, splatTexture);
943
944   s_pMaterial->SetTexture(0, GL_TEXTURE_2D, texture);
945   s_pShadeMaterial->SetTexture(0, GL_TEXTURE_2D, texture);
946   s_pMaterial->SetTextureParameter(0, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
947   s_pShadeMaterial->SetTextureParameter(0, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
948   s_pMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
949   s_pShadeMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
950   s_pMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
951   s_pShadeMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
952   s_pMaterial->EnableTexture(0, true);
953   s_pShadeMaterial->EnableTexture(0, true);
954
955   SAFE_DELETE_ARRAY(splatTexture);
956 }     
957
958
959 //------------------------------------------------------------------------------
960 // Function               : SkyCloud::_PhaseFunction
961 // Description      : 
962 //------------------------------------------------------------------------------
963 /**
964  * @fn SkyCloud::_PhaseFunction(const Vec3f& vecLightDir, const Vec3f& vecViewDir)
965  * @brief Computes the phase (scattering) function of the given light and view directions.
966  * 
967  * A phase function is a transfer function that determines, for any angle between incident 
968  * and outgoing directions, how much of the incident light intensity will be 
969  * scattered in the outgoing direction.  For example, scattering by very small 
970  * particles such as those found in clear air, can be approximated using <i>Rayleigh 
971  * scattering</i>.  The phase function for Rayleigh scattering is 
972  * p(q) = 0.75*(1 + cos<sup>2</sup>(q)), where q  is the angle between incident 
973  * and scattered directions.  Scattering by larger particles is more complicated.  
974  * It is described by Mie scattering theory.  Cloud particles are more in the regime 
975  * of Mie scattering than Rayleigh scattering.  However, we obtain good visual 
976  * results by using the simpler Rayleigh scattering phase function as an approximation.
977  */ 
978 float SkyCloud::_PhaseFunction(const Vec3f& vecLightDir, const Vec3f& vecViewDir)
979 {
980   float rCosAlpha = vecLightDir * vecViewDir;
981   return .75f * (1 + rCosAlpha * rCosAlpha); // rayleigh scattering = (3/4) * (1+cos^2(alpha))
982 }