Sign in to follow this  
vince35

atmospheric scattering theory

Recommended Posts

Hi, I'm trying to implement an atmospheric scattering model. I've read a few papers on the topic, but something always seems to be wrong. The full double nested integral as in equation 3.4 of this document http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/2554/pdf/imm2554.pdf makes sens. However, this complex equation is always simplified for shaders implementation. A common form of simplification I found is described in http://ati.amd.com/developer/dx9/ATI-LightScattering.pdf in the implementation section at the end of the document. The inscattering goes as follow: Lin = ((betaRtheta + betaMtheta) / (betaR + betaM)) * Esun * (1-exp((betaR+betaM)*S))) The problem I have with this equation is with S. They say S is the distance to the camera. Fine. I can get a decent day scene with it. But what about dusk and night? At dusk, the sky goes red at the horizon because the light travels more though the atmosphere. With S being the distance to the camera, it's constant regardless of the sun position, so the sun position as no effect other than concentrating the light with the phase function. This is not going to get me a red sunset! I could decide S is the distance from the viewpoint to the inscattering point, then to the edge of the atmosphere in the sun direction. I could get the desired sunset effet this way. But that wouldn't work for the terrain because I would always get inscattering light, even for a vertex at the viewpoint's position. A mix of both solution wouldn't work either because the transition between the sky and the groud would be obvious. So is there a trick with that system or am I just a fool not getting it? Another relation question is about the Rayleigh/Mie total scattering equations: betaR = 8*pow(PI,3)*pow(n*n-1,2)/(3*N*pow(lambda,4)) betaM = 0.434*PI*c*4*PI*PI*K/pow(lambda,2) I didn't find nowhere good numerical values for the constants with units. It seems with my choices betaM is too large compare to betaR. Does anybody know the constant values and units? (for example, is lambda in m, um, mm, other?) Thanks for taking the time to read that long post to the end :).

Share this post


Link to post
Share on other sites
ATI paper is an extension of http://www.cs.utah.edu/vissim/papers/sunsky/index.html. Check the paper and code for further details. i don't know if you convert color space, check http://www.eisscholle.de/articles/daysky.pdf.

void CalculateScatteringConstants()
{
const float n = 1.003f; // refractive index
const float N = 2.545e25f;
const float pn = 0.035f;

float fLambda[3],fLambda2[3], fLambda4[3];
fLambda[0] = 1/650e-9f; // red // note 650e-9 m = 650nm.
fLambda[1] = 1/570e-9f; // green
fLambda[2] = 1/475e-9f; // blue
for (int i=0; i < 3; i++)
{
fLambda2[i] = fLambda[i]*fLambda[i];
fLambda4[i] = fLambda2[i]*fLambda2[i];
}
D3DXVECTOR3 vLambda2(fLambda2[0], fLambda2[1], fLambda2[2]);
D3DXVECTOR3 vLambda4(fLambda4[0], fLambda4[1], fLambda4[2]);

// Rayleigh scattering constants.
float fTemp = D3DX_PI*D3DX_PI*(n*n-1)*(n*n-1)*(6+3*pn)/(6-7*pn)/N;
float fBeta = 8*fTemp*D3DX_PI/3;
vBetaR = fBeta * vLambda4;

float fBetaDash = fTemp/2;
vBetaDashR = fBetaDash * vLambda4;

// Mie scattering constants.
const float T = 2.0f;
float c = (6.544f*T - 6.51f)*1e-17f; // from page 57 of my thesis.
float fTemp2 = 0.434f*c*(2*D3DX_PI)*(2*D3DX_PI)*0.5f;
vBetaDashM = fTemp2 * vLambda2;

float K[3] = {0.685f, 0.679f, 0.670f}; // from pg 64 of my thesis.
float fTemp3 = 0.434f*c*D3DX_PI*(2*D3DX_PI)*(2*D3DX_PI);
D3DXVECTOR3 vBetaMieTemp(K[0]*fLambda2[0],K[1]*fLambda2[1], K[2]*fLambda2[2]);
vBetaM = fTemp3*vBetaMieTemp;

vBetaRM = vBetaR + vBetaM;
vOneOverBetaRM.x = 1.0f / vBetaRM.x;
vOneOverBetaRM.y = 1.0f / vBetaRM.y;
vOneOverBetaRM.z = 1.0f / vBetaRM.z;
}

Share this post


Link to post
Share on other sites
Hi,

Thanks for your post filousnt. Your coefficient code snippet was very usefull. It's the same as what I was doing, except the refractive index I found on various source is n = 1.0002926. I assume you just forgot a zero, unless you have a reason to use a different value.

Quote:
Original post by filousnt
const float n = 1.003f; // refractive index


I'm still digging this scattering problem and I foun something interesting. Or at least it's interesting because it doesn't make sens to me yet.

I'm concentrating on the extinction factor now. I assume the extinction factor have to be between 0 and 1. But my result is much larger than 1. I'm doing an integral by part of the equation 3.3 and 3.4 of the following document http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/2554/pdf/imm2554.pdf. It isn't necessary to read the document, but I find the symbolic form of the equation easier to read.

The extinction factor is the integral from the vertex to the viewpoint of
exp(-tau(s,lambda)) ds
extinction = integral(exp(-tau(s,lambda))*ds)
tau = betaR*integral(rhoR(l)*dl) + betaR*integral(rhoR(l)*dl)

If I understand those equations properly, if I take a very small ds, than I get a very small tau (<<1). At the limit, tau is 0. That means:
exp(-tau(s,lambda)) = 1.
That means the extinction integral would give me s, the length of the path. The smaller my elements are, the more it is true. I must be wrong, but where?

For those of you willing to look at the code snippet:
Quote:

for (nSeg=1; nSeg<=nSegment; nSeg++)
{
float a = -1./8300.;
float integ = abs((1./(a*fCosPhiS))*(exp(a*altSegStart)-exp(a*altSegEnd)));
extinction = fRhoR*fBetaR*integ + fRhoM*fBetaM*integ;
terrainExtinction += exp(-extinction)*seg.len;
}


where fCosPhiS is the path angle w/r to the vertical and altSegStart and altSegEnd are respectively the altitude of the start and end point of this segment.

Share this post


Link to post
Share on other sites
" 5.3 Areal Perspective
5.3.1 ProposedMethod for Calculating Areal Perspective

[...] To calculate the optical depth between a point in the terrain and the viewpoint accurately, it is necessary to solve the integral from equation (3.3).
This is not feasible for a real time approach.

The two simplest ways of estimating the average density are to either calculate
the average altitude of the viewing path and then use the density at
the average altitude, or to calculate the density at both ends of the viewing

path and use the average of these two values. Figure 5.2 shows a comparison
of the two methods for a terrain point at 1500 meters and an observer
altitude between 1000 and 30000 meters with a vertical viewing path.

As can be seen from figure 5.2, using the average altitude provides the
most precise results. [...]"

Vector eye to vertex "v = vert.xyz − eye.xyz"
length of v "l = lenght(v)"

molecular and aerosol density at the vertex position.
"Pm = e(−vert.y / 15000)"
"Pa = e(−vert.y / 5000)"

Calculate average molecular and aerosol density.
Pm,a = (Pm,a + (Pm,a at Eye pos)) / 2

Extinction Term "x = e(BetaMie * ta + BetaRay * tm)"
Optical Depth "tm,a = Pm,a . l"

ReRead the paper, check the shaders code at the end and there is some explanation on computing terms page 85. The symbolic form of the equations are not the final ones used in the shaders. Don't mix Br, Bm, Optical Depth calculation from different papers.

For a more accurate OpticalDepth, check ONeil (GPU Gems II) and Nishita's papers (Using integral with N samples).

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this