# Having trouble implementing a sky colouring algorithm...

###
#1
Members - Reputation: **122**

Posted 19 June 2002 - 03:08 AM

###
#2
Members - Reputation: **130**

Posted 19 June 2002 - 04:32 AM

1. generate a very large data set of intensities from model or photograph

2. calculate the weights. Save the weights and not the original intensities

Then when you need to display the sky,

3. (re)generate the intensities from the weights

You will not get exactly the same sky as if you had lugged around the whole large data set, but if you used enough expansion terms then it should be a good approximation. The DNKY paper discusses steps 2 and 3, not step 1.

This, at least, is my understanding. I have not looked at the paper closely or implemented it.

If you are trying to generate the sky intensities from a model, I suggest considering "A practical analytic model for daylight" (Preetham et al). They very considerately include parameter values and have source code available at their website. The paper has a nice review of the literature too. I have not implemented this either, though.

[edited by - greeneggs on June 19, 2002 12:24:57 PM]

###
#3
Members - Reputation: **1802**

Posted 19 June 2002 - 05:26 AM

The best way to get the raw data is to use photographic sources (panoramic 360° sky views). You could also use equations (1) and (2) of the mentioned paper, but results won''t be too realistic, unless you use a very accurate approximation of the integrals (long computation time...).

I implemented "A practical analytic model for daylight" some time ago, and I was somewhat dissappointed by the results. The sky colours are not bad, but not close to photorealistic, and take a rather long time to compute. But if you don''t need a photorealistic sky, but a nice approximation, then this paper is worth implementing - it will give you a fully procedural sky model. If you combine that with Nishita''s weighted basis functions, it might be an interesting and fast way to get the gradients.

/ Yann

###
#4
Members - Reputation: **122**

Posted 19 June 2002 - 05:40 AM

I''m definately after photorealism though.. Does anyone know of any analytical models that produce good values? Or should i go for sample data?

###
#6
Members - Reputation: **1802**

Posted 19 June 2002 - 07:02 AM

[edited by - Yann L on June 19, 2002 2:16:50 PM]

###
#7
Members - Reputation: **548**

Posted 19 June 2002 - 07:08 AM

There was a quite popular sky thread a little while back.

http://www.gamedev.net/community/forums/topic.asp?topic_id=86024

It''s not the same technique but it deals with sky anyway and I just wanted to mention it. Yann posts a few screenshots of his implementation with Perlin Noise octaves.

###
#8
Anonymous Poster_Anonymous Poster_*
Guests - Reputation:

Posted 24 June 2002 - 06:33 AM

quote:Original post by bishop_pass

Yann L, would you be interested in posting some screenshots of your implemetation? I would be very interested in seeing what your implemetation has produced. There are certain phenomena that I would be curious to see if it reproduces.

Looking at the paper more closely, it is actually all very simple. For skylight, I mean -- I suppose Hoffman & Preetham''s ATI article "Rendering outdoor light scattering in real time" is probably better for the aerial perspective stuff.

I don''t understand one thing, though. There are very simple formulas for computing chromaticities x and y, and luminance Y. But what should be done with these numbers? I guess you can convert to X, Y, Z, and then invert the RGB basis matrix to get RGB basis coordinates. (And then do gamma correction?)

The authors say "the luminance Y and chromaticities x and y can be converted to spectral radiance on the fly using the CIR daylight curve method described in the Appendix." What is going on here? What is the point? What does this give you more than immediately converting to RGB? Sorry, I am not familiar enough with graphics techniques to follow this step.

###
#9
Members - Reputation: **122**

Posted 26 June 2002 - 04:31 AM

Yann, itd be cool to see some of your screenshots of the algorithm at work.

-----------------------------

OpenMountains... an open source snowboarding game... One day...

###
#10
Members - Reputation: **1802**

Posted 26 June 2002 - 06:19 AM

quote:

Well I''m gonna have a crack at implementing it without the accompanied source code to try and understand it completely. If i''m successful i''ll write a tutorial or something

Yann, itd be cool to see some of your screenshots of the algorithm at work.

I can''t find the source anymore, sorry for that It must be somewhere on a dusty archival DDS tape, between gigabytes of useless junk...

But as I mentioned, it was more or less the reference implementation. And it wasn''t too realistic, I''m not sure if you''ll be really happy with the results. I would suggest to test it on the authors source (perhaps even without understanding it), and if you are satisfied with the images, then dive deeper into it and implement your own version.

Good luck anyway

/ Yann

###
#11
Anonymous Poster_Anonymous Poster_*
Guests - Reputation:

Posted 26 June 2002 - 08:01 AM

#include <math.h>

#define PI 3.14159265359f

#define radians(x) ((x) / 180.0f * PI) // longitude and latitude are inputted in degrees

// All angles in radians, theta angles measured from normal

inline float RiAngleBetween(float thetav, float phiv, float theta, float phi) {

float cospsi = sin(thetav) * sin(theta) * cos(phi-phiv) + cos(thetav) * cos(theta);

if (cospsi > 1) return 0;

if (cospsi < -1) return PI;

return acos(cospsi);

}

/**********************************************************

South = x, East = y, up = z

All times in decimal form (6.25 = 6:15 AM)

All angles in Radians

From IES Lighting Handbook pg 361

***********************************************************/

static float[3] toSun;

static float latitude,

longitude,

julianDay,

timeOfDay,

standardMeridian;

static float turbidity;

static float zenith_Y,

zenith_x,

zenith_y;

static float[5] perez_Y,

perez_x,

perez_y;

static void initSunThetaPhi() {

float solarTime, solarDeclination, solarAltitude,

opp, adj,

solarAzimuth;

solarTime = timeOfDay +

(0.170f*sin(4*PI*(julianDay - 80)/373) - 0.129f*sin(2*PI*(julianDay - 8)/355)) +

(standardMeridian - longitude)/15.0;

solarDeclination = (0.4093f*sin(2*PI*(julianDay - 81)/368));

solarAltitude = asin(sin(radians(latitude)) * sin(solarDeclination) -

cos(radians(latitude)) * cos(solarDeclination) * cos(PI*solarTime/12));

opp = -cos(solarDeclination) * sin(PI*solarTime/12);

adj = -(cos(radians(latitude)) * sin(solarDeclination) +

sin(radians(latitude)) * cos(solarDeclination) * cos(PI*solarTime/12));

solarAzimuth = atan2(opp, adj);

phiS = -solarAzimuth;

thetaS = PI / 2.0f - solarAltitude;

}

void skycolorInitialize(

float lat,

float longi,

int sm, // standardMeridian

int jd, // julianDay

float tOfDay, // timeOfDay

float turb // turbidity

)

{

float Ttheta2, theta3, T2, chi;

latitude = lat;

longitude = longi;

julianDay = jd;

timeOfDay = tOfDay;

standardMeridian = sm * 15.0; // sm is actually timezone number (east to west, zero based...)

turbidity = turb;

initSunThetaPhi();

toSun[0] = cos(phiS) * sin(thetaS);

toSun[1] = sin(phiS) * sin(thetaS);

toSun[2] = cos(thetaS);

theta2 = thetaS * thetaS;

theta3 = theta2 * thetaS;

T = turb;

T2 = turb * turb;

chi = (4.0f / 9.0f - T / 120.0f) * (PI - 2 * thetaS);

zenith_Y = (4.0453f * T - 4.9710f) * tan(chi) - .2155f * T + 2.4192f;

zenith_Y *= 1000; // conversion from kcd/m^2 to cd/m^2

zenith_x =

( 0.00165f * theta3 - 0.00375f * theta2 + 0.00209f * thetaS + 0.0f) * T2 +

(-0.02903f * theta3 + 0.06377f * theta2 - 0.03202f * thetaS + 0.00394f) * T +

( 0.11693f * theta3 - 0.21196f * theta2 + 0.06052f * thetaS + 0.25886f);

zenith_y =

( 0.00275f * theta3 - 0.00610f * theta2 + 0.00317f * thetaS + 0.0f) * T2 +

(-0.04214f * theta3 + 0.08970f * theta2 - 0.04153f * thetaS + 0.00516f) * T +

( 0.15346f * theta3 - 0.26756f * theta2 + 0.06670f * thetaS + 0.26688f);

perez_Y[1] = 0.17872f * T - 1.46303f;

perez_Y[2] = -0.35540f * T + 0.42749f;

perez_Y[3] = -0.02266f * T + 5.32505f;

perez_Y[4] = 0.12064f * T - 2.57705f;

perez_Y[5] = -0.06696f * T + 0.37027f;

perez_x[1] = -0.01925f * T - 0.25922f;

perez_x[2] = -0.06651f * T + 0.00081f;

perez_x[3] = -0.00041f * T + 0.21247f;

perez_x[4] = -0.06409f * T - 0.89887f;

perez_x[5] = -0.00325f * T + 0.04517f;

perez_y[1] = -0.01669f * T - 0.26078f;

perez_y[2] = -0.09495f * T + 0.00921f;

perez_y[3] = -0.00792f * T + 0.21023f;

perez_y[4] = -0.04405f * T - 1.65369f;

perez_y[5] = -0.01092f * T + 0.05291f;

}

static inline float PerezFunction(const float *lam, float theta, float gamma, float lvz) const {

float den = ((1 + lam[1]*exp(lam[2])) *

(1 + lam[3]*exp(lam[4]*thetaS) + lam[5]*cos(thetaS)*cos(thetaS)));

float num = ((1 + lam[1]*exp(lam[2]/cos(theta))) *

(1 + lam[3]*exp(lam[4]*gamma) + lam[5]*cos(gamma)*cos(gamma)));

return lvz* num/den;

}

RiSpectrum getSkySpectralRadiance(float theta, float phi) const {

float gamma = RiAngleBetween(theta, phi, thetaS, phiS);

// Compute xyY values

float x = PerezFunction(perez_x, theta, gamma, zenith_x);

float y = PerezFunction(perez_y, theta, gamma, zenith_y);

float Y = PerezFunction(perez_Y, theta, gamma, zenith_Y);

RiSpectrum spect = ChromaticityToSpectrum(x, y);

// A simple luminance function would be more efficient.

return Y * spect / RiColorXYZV(spect).Y();

}

I believe this is all the code you need (assuming you aren''t doing aerial perspective using this model, which you probably shouldn''t). Then render a mesh calling getSkySpectralRadiance for each vertex. Then render the sun itself.

I still don''t understand how they are getting the final color -- what is ChromaticityToSpectrum doing. What does the comment "A simple luminance function would be more efficient" mean? (use another Perez function, I suppose) See my earlier post. Any help?

###
#12
Members - Reputation: **130**

Posted 26 June 2002 - 08:48 AM

Maybe Yann knows a better way? (or can explain this way?) If Yann was able to reproduce the shots of the "practical.." paper, then he must have figured out something.

[edited by - greeneggs on June 26, 2002 3:54:10 PM]

###
#13
Members - Reputation: **1802**

Posted 27 June 2002 - 04:45 PM

quote:

Maybe Yann knows a better way? (or can explain this way?) If Yann was able to reproduce the shots of the "practical.." paper, then he must have figured out something.

I didn't use the spectral radiance conversion (with M1 and M2) at all. I directly converted the luminance Y and the two chromaticities (x,y) to CIE XYZ tristimulus values, and then to RGB (standard rec 709 RGB primaries and white point) using a simple matrix multiply.

As far as I see it, the Yxy->spectral radiance conversion step doesn't really add something to the image. It splits up the chromatic components into a spectral distribution, which might be handy if you want to do colour adaption later on (on different colour temperatures or display media). But you could also do that by adjusting the primaries or white point of the XYZ->RGB conversion. And it's a lot faster that way.

If anyone is interested, I can post the equations to go from (Y,xy) to RGB.

/ Yann

[edited by - Yann L on June 27, 2002 11:57:23 PM]

###
#14
Members - Reputation: **122**

Posted 27 June 2002 - 10:08 PM

quote:Original post by Yann L

If anyone is interested, I can post the equations to go from (Y,xy) to RGB.

Please! :O)

-----------------------------

OpenMountains... an open source snowboarding game... One day...

###
#15
Members - Reputation: **130**

Posted 28 June 2002 - 12:15 PM

quote:Original post by Yann L

I didn't use the spectral radiance conversion (with M1 and M2) at all. I directly converted the luminance Y and the two chromaticities (x,y) to CIE XYZ tristimulus values, and then to RGB (standard rec 709 RGB primaries and white point) using a simple matrix multiply.

Using this XYZ->RGB conversion procedure (Rec. 709 RGB with D65 white point, see http://www.inforamp.net/%7Epoynton/notes/colour_and_gamma/ColorFAQ.html#RTFToC18 )

R 3.240479,-1.537150,-0.498535 X

G = -0.969256, 1.875992, 0.041556 Y

B 0.055648,-0.204043, 1.057311 Z

I obtained the following image (90 degree FOV, turbidity 3)

Here's my complete Renderman source code (only 100 lines, mostly initialization):

float PerezFunction(

float coeffs[5];

float theta, gamma;

float thetaS, lvz;

)

{

float A = coeffs[0], B = coeffs[1], C = coeffs[2], D = coeffs[3], E = coeffs[4];

float den = (1 + A * exp(B)) * (1 + C * exp(D * thetaS) + E * cos(thetaS)*cos(thetaS));

float num = (1 + A * exp(B / cos(theta))) * (1 + C * exp(D * gamma) + E * cos(gamma)*cos(gamma));

return (lvz * num / den);

}

float AngleBetween(float thetav, phiv, theta, phi) {

float cospsi = sin(thetav) * sin(theta) * cos(phi-phiv) + cos(thetav) * cos(theta);

if (cospsi > 1) return 0;

if (cospsi < -1) return PI;

return acos(cospsi);

}

surface

skycolor (

vector toSun = (1, -3, 1.5);

float turb = 3;

)

{

toSun = normalize(toSun);

uniform float thetaS = acos(zcomp(toSun));

uniform float phiS = atan(ycomp(toSun), xcomp(toSun));

uniform float theta2, theta3, T, T2, chi, zenith_Y, zenith_x, zenith_y;

uniform float perez_Y[5], perez_x[5], perez_y[5];

theta2 = thetaS * thetaS;

theta3 = theta2 * thetaS;

T = turb;

T2 = turb * turb;

chi = (4 / 9 - T / 120) * (PI - 2 * thetaS);

zenith_Y = (4.0453 * T - 4.9710) * tan(chi) - .2155 * T + 2.4192;

zenith_Y *= 1000; // conversion from kcd/m^2 to cd/m^2

zenith_x =

( 0.00165 * theta3 - 0.00375 * theta2 + 0.00209 * thetaS + 0.0) * T2 +

(-0.02903 * theta3 + 0.06377 * theta2 - 0.03202 * thetaS + 0.00394) * T +

( 0.11693 * theta3 - 0.21196 * theta2 + 0.06052 * thetaS + 0.25886);

zenith_y =

( 0.00275 * theta3 - 0.00610 * theta2 + 0.00317 * thetaS + 0.0) * T2 +

(-0.04214 * theta3 + 0.08970 * theta2 - 0.04153 * thetaS + 0.00516) * T +

( 0.15346 * theta3 - 0.26756 * theta2 + 0.06670 * thetaS + 0.26688);

perez_Y[0] = 0.17872 * T - 1.46303;

perez_Y[1] = -0.35540 * T + 0.42749;

perez_Y[2] = -0.02266 * T + 5.32505;

perez_Y[3] = 0.12064 * T - 2.57705;

perez_Y[4] = -0.06696 * T + 0.37027;

perez_x[0] = -0.01925 * T - 0.25922;

perez_x[1] = -0.06651 * T + 0.00081;

perez_x[2] = -0.00041 * T + 0.21247;

perez_x[3] = -0.06409 * T - 0.89887;

perez_x[4] = -0.00325 * T + 0.04517;

perez_y[0] = -0.01669 * T - 0.26078;

perez_y[1] = -0.09495 * T + 0.00921;

perez_y[2] = -0.00792 * T + 0.21023;

perez_y[3] = -0.04405 * T - 1.65369;

perez_y[4] = -0.01092 * T + 0.05291;

// shader code starts here

color skycolor = 0;

vector Iw = normalize(transform("world", I));

float theta = acos(zcomp(Iw));

float phi = atan(ycomp(Iw), xcomp(Iw));

float gamma = AngleBetween(theta, phi, thetaS, phiS);

// Compute xyY values

float x = PerezFunction(perez_x, theta, gamma, thetaS, zenith_x);

float y = PerezFunction(perez_y, theta, gamma, thetaS, zenith_y);

float Y = PerezFunction(perez_Y, theta, gamma, thetaS, zenith_Y);

float X = (x / y) * Y;

float Z = ((1 - x - y) / y) * Y;

skycolor = color (

1-exp(-(1/15000) * ( 3.240479 * X - 1.537150 * Y - 0.498535 * Z)),

1-exp(-(1/15000) * (-0.969256 * X + 1.875992 * Y + 0.041556 * Z)),

1-exp(-(1/15000) * ( 0.055648 * X - 0.204043 * Y + 1.057311 * Z))

);

if (zcomp(Iw) < 0) skycolor = 0;

Oi = 1;

Ci = skycolor;

}

As you can see, I'm using an exponential exposure function.

By the way, as you can see in the source code, X = x/y * Y, Z = (1-x-y)/y * Y

[Edit: Corrected my XYZ->RGB matrix, I had been using the 1958 NTSC matrix ]

[edited by - greeneggs on June 28, 2002 7:42:11 PM]

###
#16
Members - Reputation: **1802**

Posted 28 June 2002 - 12:33 PM

Now, IIRC, I had to alter the matrix aswell, in order to get similar shots, but it wasn't very much. I also added some kind of chromatic correction curve, in a (failed) attempt to make them more realistic. Since I don't have my source anymore, I can't give you the exact matrix values, but I'm almost sure I started by using the rec.709 HDTV matrix:

R 3.240479 -1.53715 -0.49853 X

G = -0.969256 1.875991 0.041556 Y

B 0.055648 -0.204043 1.057311 Z

Try those values, does it look better ? [edited answer: yes !]

As a reference, for Stimpson:

To convert from Yxy to CIE XYZ, you do:

Y = Y

X = x * (Y / y)

Z = (1 - x - y) * (Y / y)

Then, you go from XYZ to RGB by multiplying it with the above mentioned matrix.

/ Yann

[edited by - Yann L on June 28, 2002 7:51:05 PM]

###
#17
Members - Reputation: **130**

Posted 28 June 2002 - 12:51 PM

quote:Original post by Yann L

Hmm, your matrix seems way off, where did you get it from ?

Yes, sorry about that. I fixed my earlier post. Still, it doesn''t look quite like the Preetham et al paper. The sky has too much red and not enough orange. The images in the paper look better. Could the daylight curve stuff affect that?

###
#18
Members - Reputation: **1802**

Posted 28 June 2002 - 12:58 PM

quote:

Yes, sorry about that. I fixed my earlier post. Still, it doesn''t look quite like the Preetham et al paper. The sky has too much red and not enough orange. The images in the paper look better. Could the daylight curve stuff affect that?

Hmm, good question. As I mentioned, I did some additional tweaking to get some results, but it really wasn''t that much. My results were very near to the results in the paper (although not 100% the same).

The daylight curve does a direct modification on the x and y chroma components, so it could affect the colour. But what would be the point of doing that ? Why not directly doing a gamma correction on x,y and then go to XYZ/RGB ? The whole Yxy->spectral radiance, spectral radiance->Yxy->XYZ->RGB seems totally awkward to me, if it was just for a colour correction.

/ Yann

###
#19
Members - Reputation: **130**

Posted 29 June 2002 - 10:11 AM

First is a projection of the upper hemisphere of the view sphere down onto the equatorial disk. Times 1am, 2am, ... 12pm (the afternoon and evening are symmetrical). North is up, East is to the right (sort of upside down). The algorithm does not correctly compute the night sky, between about 9pm and 3am, as you can see.

Next sunrise, looking directly East. 15 minute intervals starting at 3am, ending at 7am. (90 degree FOV, images cropped to only show upper half (the sky))

Similarly sunset, looking East. 15 minute intervals starting at 5pm, ending at 9pm.

By symmetry you can think of the sunset pictures as sunrise pictures looking West, if you flip them left to right and look at them in reverse order.

Any comments? Effects which are missing? Suggestions for improving color matching? These images should be easily reproducible with the shader I posted earlier and BMRT.

###
#20
Members - Reputation: **1802**

Posted 29 June 2002 - 01:21 PM

quote:

Any comments? Effects which are missing? Suggestions for improving color matching?

Looks good to me. Did you modify the conversion matrix to get the shots ?

There are just two little things that crossed my mind: when looking directly into the sun (such as on your sunrise shots), the sun saturates too much. Although this could be easily solved using an exposure or similar high dynamic range function.

And the lighting at dusk/dawn, when the sun is below the horizon, isn''t right (as you already pointed out). It should be more red/orange, instead of the dominant yellow. But I think this is not so much a problem of your code, but more a limitation of this particular skylight model.

/ Yann