 Hi everyone. I’m trying to write a sky shader based on Nishita’s model that allows you to specify altitude as one of the sky parameters. I’m following this article here: http://www.scratchapixel.com/lessons/3d-advanced-lessons/simulating-the-colors-of-the-sky/atmospheric-scattering/

So far I’ve had some luck, and my shader compiles, but I can’t seem to set the sun direction properly - whatever vector I specify the sun always seems to end up below or at the horizon. I’m clearly doing something wrong with some vector maths somewhere, but I can’t track it down. Can anyone see what I’m doing wrong?

``````
#include &lt;stdosl.h&gt;

float Height = 1,
float Intensity = 20,

vector InSunDirection = vector(0,0,1),
vector InDirection = P,
output color Rayleigh = 0.0,
output color Mie = 0.0,
output color Sky = 0.0)

{
/* Set out my important variables */
float Re = 6360e3; // Radius of the earth
float Ra = 6420e3; // Radius of the atmosphere
float Hr = 7994; // Rayleigh scale height
float Hm = 1200; // Mie scale height
float g = 0.76; // Anisotropy term for Mie scattering.

vector BetaR = vector(5.8e-6,13.0e-6,22.4e-6);
vector BetaM = vector(21e-5);

vector SumR = vector(0);
vector SumM = vector(0);

/* Compute intersection point of camera ray with atmosphere*/

// Five possible situations:
// 1. No intersection, return black.
// 2. Camera inside planet. Return black.
// 3. Camera inside atmosphere. Return correct d value.
// 4. Camera inside atmoshere, facing earth, return d for collision with earth.
// 5. Camera outside atmosphere, can see through to other side. Return two d values one for each entry and exit.
// 6. Camera outside atmosphere, can't see through to other side. Return two d values, one for entry, one for earth.

// For this we need d, the distance along the camera ray we need to travel to reach the edge of the atmosphere.

vector Direction = normalize (InDirection);
vector SunDirection = normalize (InSunDirection);

// Our start and end points for our integral, these should be overwritten below.
vector Pu = vector(0,0,0);
vector Pv = vector(0,0,0);

if (Height &gt; 0) { // Checks height is above ground.
vector Pc = vector (0, 0, Re + Height);
Pu = Pc;
float b = 2 * dot (Direction, Pc);
float bsquared = pow(b,2);
float fourac = 4 * dot(Direction, Direction) * (1 - pow(Ra,2));

if (bsquared - fourac &gt; 0) { //Checks that our view direction intersects with the atmosphere.
float d1 = (-b - pow(bsquared - fourac, 0.5)) / 2;
float d2 = (-b + pow(bsquared - fourac, 0.5)) / 2;
// Check solutions for our outcomes:
if (Pc &gt; Ra) { //Checks if camera is outside atmosphere.
//We've already ruled out no solutions, so we expect two, they should both be positive.
if (d1 &gt; 0 && d2 &gt; 0) {
float d_entry = min(d1,d2);
float d_exit = max(d1,d2);
vector Pe = Pc + d_entry * Direction; // Pe is the entry point to the atmosphere.
Pu = Pe; // Overwrite Pu with Pe, its now our start point for the integral.
// Does the view pass through earth? If so find Pr.
float b_earth = 2 * dot (Direction, Pc);
float bsquared_earth = 2 * b_earth;
float fourac_earth = 4 * dot (Direction,Direction) * (1 - pow(Re, 2));
if (bsquared_earth - fourac_earth &gt; 0) { // Check if ray hits earth.
float d1_earth = (-b_earth - pow(bsquared_earth - fourac_earth, 0.5)) / 2;
float d2_earth = (-b_earth + pow(bsquared_earth - fourac_earth, 0.5)) / 2;
if (d1_earth &gt; 0 && d2_earth &gt; 0) { // Solutions aren't negative.
float d_entry = min(d1_earth,d2_earth);
vector Pr = Pc + d_entry * Direction;
Pv = Pr;
}
} else { // View doesn't hit earth, exits atmosphere.
vector Pa = Pc + d_exit * Direction;
Pv = Pa;
}
}
} else {
// Our camera is inside the atmosphere.
// We expect one negative solution, one positive, we only want the positive.
float d_exit = max(d1,d2);
vector Pa = Pc + d_exit * Direction;
Pv = Pa;
}
}
}

/* Now compute Rayleigh and Mie Phases, these are constant terms, not part of the integral.*/

float mu = dot(Direction, SunDirection);
float phaseR = (3/(16 * M_PI)) * (1 + mu*mu);
float phaseM = (3/(8 * M_PI)) * ((1 - g*g) * (1 + mu*mu) / ((2 + g*g) * pow(1 + g*g - 2 * g * mu, 1.5)));

/* Create samples along camera ray. */

int numSamples = 16;
int numSamplesL = 8;

float segmentLength = distance (Pu, Pv) / numSamples;

float opticalDepthR = 0;
float opticalDepthM = 0;

// Now begin taking samples for the integtral.
for (int i = 0; i &lt; numSamples; i++) {
vector samplePosition = Pu + segmentLength * Direction * (i + 0.5);
float sampleHeight = length(samplePosition) - Re;

/* Get optical depth for light at this sample: */
float Hr_sample = exp(-sampleHeight/Hr) * segmentLength;
float Hm_sample = exp(-sampleHeight/Hm) * segmentLength;

opticalDepthR += Hr_sample;
opticalDepthM += Hm_sample;

/* Find light sample termination point Ps for this sample: */

// There are 2 possibilities here (we now know we are in the atmosphere because only lines within the atmosphere are sampled):
// 1. The sun is visible from this position.
// 2. It isn't, because its behind the earth.

float opticalDepthLR = 0;
float opticalDepthLM = 0;

// Check that the vector in the direction of the sun from the sample position doesnt intersect the earth.
float b_s = 2 * dot (SunDirection,samplePosition);
float bsquared_s = pow(b_s, 2);
float fourac_s_earth = 4 * dot(SunDirection,SunDirection) * (1 - pow(Re,2));
int canseesun = 0;
if (bsquared_s - fourac_s_earth &gt; 0) { // Check if our vector towards the sun goes through the earth. If so dont bother sampling.
// There are Solutions! But if they're both negative we can carry on.
float d1_s_earth = (-b_s - pow(bsquared_s - fourac_s_earth, 0.5)) / 2;
float d2_s_earth = (-b_s + pow(bsquared_s - fourac_s_earth, 0.5)) / 2;

if (d1_s_earth &lt; 0 && d2_s_earth &lt; 0) {
//Both of our earth intersection vectors are behind the camera. So we can see the sun.
canseesun = 1;
}
} else { //Doesn't hit earth. Can see sun.
canseesun = 1;
}
if (canseesun) {
// Calculate distance to Ps and do samples.
float fourac_s = 4 * dot(SunDirection, SunDirection) * (1 - pow(Ra,2)); // This really shouldn't be negative at any time, that would mean we we inside a sphere but somehow the sun wasn't outside it...
float d1_s = (-b_s - pow(bsquared_s - fourac_s, 0.5)) / 2;
float d2_s = (-b_s + pow(bsquared_s - fourac_s, 0.5)) / 2;
// One of the two above should be negative, pick the positive(negative?) one.
float d_s = 0;
if (d1_s &gt; 0) {
d_s = d1_s;
} else {
d_s = d2_s;
}

vector Ps = samplePosition + d_s * SunDirection;

/* Create samples along sun ray*/
float segmentLengthL = d_s / numSamplesL;

for (int j = 0; j &lt; numSamplesL; j++) {
vector samplePositionL = samplePosition + segmentLengthL * SunDirection * (j + 0.5);

/* Test if light ray is below ground, and discard if it is. */
float sampleHeightL = length(samplePositionL) - Re;
if (sampleHeightL &gt; 0) {
// ignore sampleheights &lt; 0 they're inside the planet...
/* Get optical depth for light at this sample */
float Hlr_sample = exp(-sampleHeightL/Hr)*segmentLengthL;
float Hlm_sample = exp(-sampleHeightL/Hm)*segmentLengthL;

opticalDepthLR += Hlr_sample;
opticalDepthLM += Hlm_sample;
}
}
}

/* once we've done all our samples we can calculate our attenuation for the light */

vector tauR = BetaR * (opticalDepthR + opticalDepthLR);
vector tauM = BetaM * (opticalDepthM + opticalDepthLM);
vector attnR = vector (exp(-tauR), exp(-tauR), exp(-tauR));
vector attnM = vector (exp(-tauM));

SumR += Hr_sample * attnR;
SumM += Hm_sample * attnM;
}

//printf("%.2f
", Direction);
Rayleigh = color (SumR * phaseR * BetaR) * Intensity;
Mie = color (SumM * phaseM * BetaM) * Intensity;
Sky =  Rayleigh + Mie;

}

``````

Here’s what the output looks like currently (you’re looking at the sky thorough a fisheye lens here):

1 Like

Woo! Some progress. Re-wrote some bits of the shader for calculating the intersection of rays with the earth/atmosphere a bit more neatly and got this!

(It’s a equirectangular projection this time. Makes a bit more sense to look at.

Here’s the updated shader. Was trying to put those four functions at the start into a single function that returned a struct, but my osl skills aren’t up to that… was just getting syntax errors I didn’t know how to correct.

``````#include &lt;stdosl.h&gt;

int has_solutions (vector StartPoint, vector UnitVector, float Radius) {
float b = 2 * dot (UnitVector, StartPoint);
float bsquared = pow(b,2);
float c = dot (StartPoint,StartPoint) - pow(Radius,2);
float fourac = 4 * c;

if (bsquared &gt; fourac) {
return 1;
} else {
return 0;
}
}

int linetype (vector StartPoint, vector UnitVector, float Radius) {
float b = 2 * dot (UnitVector, StartPoint);
float bsquared = pow(b,2);
float c = dot (StartPoint,StartPoint) - pow(Radius,2);
float fourac = 4 * c;

if (bsquared &gt; fourac) {
float d1 = (-b + pow(bsquared - fourac, 0.5)) / 2;
float d2 = (-b - pow(bsquared - fourac, 0.5)) / 2;
if (d1 &lt; 0 && d2 &lt; 0) {
return 2;
} else if (d1 &lt; 0 || d2 &lt; 0) {
return 1;
} else {
return 0;
}
} else {
return 3;
}
}

vector lmin (vector StartPoint, vector UnitVector, float Radius) {
float b = 2 * dot (UnitVector, StartPoint);
float bsquared = pow(b,2);
float c = dot (StartPoint,StartPoint) - pow(Radius,2);
float fourac = 4 * c;

if (bsquared &gt; fourac) {
float d1 = (-b + pow(bsquared - fourac, 0.5)) / 2;
float d2 = (-b - pow(bsquared - fourac, 0.5)) / 2;
return StartPoint + UnitVector * min(d1, d2);
} else {
return vector(0);
}
}

vector lmax (vector StartPoint, vector UnitVector, float Radius) {
float b = 2 * dot (UnitVector, StartPoint);
float bsquared = pow(b,2);
float c = dot (StartPoint,StartPoint) - pow(Radius,2);
float fourac = 4 * c;

if (bsquared &gt; fourac) {
float d1 = (-b + pow(bsquared - fourac, 0.5)) / 2;
float d2 = (-b - pow(bsquared - fourac, 0.5)) / 2;
return StartPoint + UnitVector * max(d1, d2);
} else {
return vector(0);
}
}

struct linesample {
int has_solutions;
int type;
vector min;
vector max;
};

float Height = 1,
float Intensity = 20,

vector SunVector = vector(0,0,1),
vector Vector = P,
output color Rayleigh = 0.0,
output color Mie = 0.0,
output color Sky = 0.0)

{
//Variables:
float Re = 6360e3; // Radius of the earth
float Ra = 6420e3; // Radius of the atmosphere
float Hr = 7994; // Rayleigh scale height
float Hm = 1200; // Mie scale height
float g = 0.76; // Anisotropy term for Mie scattering.

vector BetaR = vector(5.8e-6,13.0e-6,22.4e-6);
vector BetaM = vector(21e-6);

//Clean up inputs:
vector SunDirection = normalize(SunVector);
vector Direction = normalize(Vector);

//Calcualte Rayleigh and mie phases.
float mu = dot(Direction, SunDirection);
float phaseR = (3/(16 * M_PI)) * (1 + mu*mu);
float phaseM = (3/(8 * M_PI)) * ((1 - g*g) * (1 + mu*mu) / ((2 + g*g) * pow(1 + g*g - 2 * g * mu, 1.5)));

//Stuff that will eventually help form my output.
vector SumR = vector(0);
vector SumM = vector(0);

//Set up start and end points for my integrals.
vector Pu = vector(0,0,0);
vector Pv = vector(0,0,0);

//Get start and end point for View Line.
if (Height &gt; 0) {
vector Pc = vector (0,0,Re + Height);
Pu = Pc;

linesample groundline;
groundline.has_solutions = has_solutions(Pc, Direction, Re);
groundline.type = linetype(Pc, Direction, Re);
groundline.min = lmin(Pc, Direction, Re);
groundline.max = lmax(Pc, Direction, Re);

linesample atmosphereline;
atmosphereline.has_solutions = has_solutions(Pc, Direction, Ra);
atmosphereline.type = linetype(Pc, Direction, Ra);
atmosphereline.min = lmin(Pc, Direction, Ra);
atmosphereline.max = lmax(Pc, Direction, Ra);

Pv = atmosphereline.max;
if (atmosphereline.has_solutions && atmosphereline.type == 0) {
//Line starts outside atmoshere, so Pu becomes Line.min
Pu = atmosphereline.min;
}
if (groundline.has_solutions && groundline.type == 0) {
//Line hits ground, Pv becomes point at which line hits ground, which will be Line.min
Pv = groundline.min;
}
}

// Create samples along view ray.
int numSamples = 16;
int numSamplesL = 8;

float segmentLength = distance (Pu, Pv) / numSamples;
float opticalDepthR = 0;
float opticalDepthM = 0;

for (int i = 0; i &lt; numSamples; i++) {
vector Px = Pu + segmentLength * Direction * (i + 0.5);
float sampleHeight = length(Px) - Re;

/* Get optical depth for light at this sample: */
float Hr_sample = exp(-sampleHeight/Hr) * segmentLength;
float Hm_sample = exp(-sampleHeight/Hm) * segmentLength;

opticalDepthR += Hr_sample;
opticalDepthM += Hm_sample;

//Get light depth to sun at this sample.
float opticalDepthLR = 0;
float opticalDepthLM = 0;

linesample atmosphereline_sample;
atmosphereline_sample.has_solutions = has_solutions(Px, SunDirection, Ra);
atmosphereline_sample.type = linetype(Px, SunDirection, Ra);
atmosphereline_sample.min = lmin(Px, SunDirection, Ra);
atmosphereline_sample.max = lmax(Px, SunDirection, Ra);

linesample groundline_sample;
groundline_sample.has_solutions = has_solutions(Px, SunDirection, Re);
groundline_sample.type = linetype(Px, SunDirection, Re);
groundline_sample.min = lmin(Px, SunDirection, Re);
groundline_sample.max = lmax(Px, SunDirection, Re);

vector Ps = atmosphereline_sample.max;
if (groundline_sample.has_solutions && groundline_sample.type == 0) {
//Light is below horizon from this point. No sample needed.
float opticalDepthLR = 0;
float opticalDepthLM = 0;
} else {
float segmentLengthL = distance(Px, Ps) / numSamplesL;

for (int j = 0; j &lt; numSamplesL; j++) {
vector Pl = Px + segmentLengthL * SunDirection * (j + 0.5);
float sampleHeightL = length(Pl) - Re;
if (sampleHeightL &gt; 0) {
// ignore sampleheights &lt; 0 they're inside the planet...
float Hlr_sample = exp(-sampleHeightL/Hr)*segmentLengthL;
float Hlm_sample = exp(-sampleHeightL/Hm)*segmentLengthL;

opticalDepthLR += Hlr_sample;
opticalDepthLM += Hlm_sample;
}
}
}

// With our light samples done we can calculate the attenuation.

vector tauR = BetaR * (opticalDepthR + opticalDepthLR);
vector tauM = BetaM * (opticalDepthM + opticalDepthLM);
vector attnR = vector (exp(-tauR), exp(-tauR), exp(-tauR));
vector attnM = vector (exp(-tauM));

SumR += Hr_sample * attnR;
SumM += Hm_sample * attnM;
}

Rayleigh = color (SumR * phaseR * BetaR) * Intensity;
Mie = color (SumM * phaseM * BetaM) * Intensity;
Sky =  Rayleigh + Mie;
}
``````

Very Nice! Hey Ben,

This is what happens if I test it:

Isn’t this what’s it supposed to look like?

Although, if I put my sun below the horizon (negative Z), I still get a blue sky. Should the algorithm take that into account?

Regards Nick

PS: I don’t think that I ever said this but I really like your book a lot. Thanks for taking the time to write it.

Thanks Nick! Yes it’s working much better now, tough still some bugs to fix. The sun below the horizon causing a blue sky being the main one. You might have noticed the sunset looks a bit white rather than orange too, I’ve fixed that one already, will post some better images tomorrow.

Pushed the shader to github so I can keep on top of updating it better.

Some more renders…

Improved sunset:

Sunny day, low altitude:

High altitude shot (space above, earth below):

This really a nice script! great results.

Thanks secrop! I’ve now managed to fix the issue with the sun going below the horizon too.

Here’s a sunset: And some desktop background res shots!

Congratulations, what amazing result!

Here’s a simple test using the sky for some lighting (with a real sun lamp for shadows).

Great!
Could you explain what is the input parameter “Vector”?
Thanks.

Sorry for my bad English … is the fault of Google Translator !!! Is it possible to change the script, so it works as shader for meshes? F.e. use it as an “modeled” planet atmosphere shader?

Incredible Shader man. On Blender 2.72 the shader gives an error as soon you use it (Missing stdols.h), but if I remove the #include it worked. It’s useful to create some simple Environment maps to use.

Thanks!

Awesome script Thank you.

How can I use a Sun object or maybe better light direction of sun object as input for your “sun vector”?
I tried to ad Atribute node, but without success.
Thank you for help.

Hi Ben,

the shader is pretty darn good. However, I don’t understand how the SunVector works. Is it like xyz coordinates of the sun with respect to global (0,0,0)? or some sort of rotation around the global origin? Can you elaborate on this?
Same question goes for the other vector input.
I tried to understand what these inputs are actually doing but I have the impression that both behave a bit inconsistent/arbitrary.

I would be glad if you could shed some light on the vector inputs.

Best regards,
Meikel.

If we could add some clouds also - just a wish
thank you