Please note, that in your code the bottom and top images seem to be interchanged! (a libpng effect? Kind of known to me)
Yes, it's a libpng issue. I'll fix that so the code makes a bit more sense.
Did you fix the LUT, meanwhile?
Not yet--I've actually been working on the inscattering LUT, which is what the GLSL shader would eventually use. I'm going see what happens when I switch from an extinction LUT to an optical depth LUT; linear interpolation artifacts are the only thing I can think of that could be causing the problems I'm seeing. Also, I notice that the Schafhitzel paper uses optical depth instead of extinction as well.
I am presently comparing your code with the paper by Nishita et al. and Cornette&Shanks (Appl. Optics).
Why didn't you use the improved
Henyey-Greenstein(HG) function including Cornette's corrections for Mie scattering? You apparently took the original
HG phase function.
double phaseHenyeyGreenstein(double cosTheta, double g)
return (1.0 - g * g) / pow(1.0 + g * g - 2.0 * g * cosTheta, 1.5);
Since the shader must execute quickly, I was anticipating using the Schlick approximation to the original HG phase function. The phase function could be encoded in a 1D texture, but that adds some complexity, since Celestia would have to generate a new texture for each value of g used. Regardless of all that, we should support Schlick, the original HG, and the CS-corrected version of HG in scattersim so that we can compare differences between them. If the visual differences ends up being minimal, we can take the simple approach (Schlick) in Celestia.
Unlike the original HG phase function, the good thing about the CS-corrected version is that for g=0 one recovers exactly the Rayleigh phase function. That is however necessary, when the size parameter, a = 2*PI*r/ lambda is significantly
less than one and where r is the particle radius. After all, Mie scattering is the more general framework that must reproduce Rayleigh scattering as a limiting case for a<<1.
I actually have that latter CornetteShamks paper in Applied Optics (it costs money). I modified your code to include the CS-improved phase function. It doesn't take much longer.
Good. Please check in this change, and I'll add a new cfg file parameter to select the phase function.
I think in the Mie-backscattering strength something is wrong. If I take exactly 180 degs between source and camera, there is very little backscatter independent of the MieAsymmetry parameter! This might be connected with your uncorrected HG-Phasefunction?
Yes, something is fishy there I agree . . . I'm looking into it.
I also wonder whether one might better introduce more physical parameters to describe atmospheric scattering, like (using Nishita's nomenclature now)
for molecular scattering (Rayleigh):
- (effective) index of refraction of air n,
- the molecular density parameter of (standard) atmosphere N_s,
of course one may directy take 4*Pi8k/lambda^4, but using the individual parameters may be not so bad in various case...
- the actual (molecular) atmosphere heights instead of the Rayleigh scale height
The latter can usually be looked up in or estimated from planetary databases.
- the wavelength lambda of incident light along with parametrizing explicitly the Rayleight 1/lambda^4 constraint into the corresponding RGB colors, etc.
I think we should stick to a 1/lambda^4 Rayleigh constraint in RGB as long as possible.
Yes, we should make it easy to create an atmopshere by setting these physical parameters. However, for experimentation, I still want it to be possible to set the derived quantities directly, i.e. the attenuation coefficients due to outscattering and absorption. Especially in Celestia, I don't believe that we should enforce the 1/lambda^4 constraint, in order to allow some flexibility in atmosphere modeling.
for aerosol scattering (Mie)
- The quantity u -> g(u) determining the asymmetry factor g, u=0.5..1.0
Is the sign convention of your asymmetry factor in scattersim the same as in Celestia??
- The actual aerosol atmosphere height
Sure . . . although, if we specify an aerosol atmosphere height, how does the density vary with h? Which brings to mind another possibility: rather than simply using exp(-h / H0) for density, we could use a lookup table from some standard atmosphere model. I'm not sure that it's worth the effort given the other simplifications in scattersim (e.g. no multiple scattering), but it might be interesting. For my part, I'm going to focus on getting this all working on the GPU.
A possibly tricky issue might be the accuracy of the numerical integration of the optical depth. The trapezoidal rule with varying intervals ~ log(rho) might well be satisfactory often but not always. I know from my own experience long ago that due to the exponential variation of the atmospheric density care is necessary. Did you test the accuracies explicitly in your code?? I noted though that you allow for using different numbers of integration steps in the command line.
The only accuracy tests I performed were visual ones: how do things look with 10 vs 20 integration steps, etc. Since the ultimate aim is producing an image, this seemed adequate, though it would be nice to quantify the the differences.