• Create Account

# Having trouble implementing a sky colouring algorithm...

Old topic!

Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

38 replies to this topic

### #1Stimpson  Members   -  Reputation: 122

Like
Likes
Like

Posted 19 June 2002 - 03:08 AM

I've been trying to implement Nishita's paper "A Fast Display Method of Sky Color Using Basis Functions". If anyone has implemented this paper let me know! But anyway, the problem I have is computing the weights: To calculate a weight (eq. 10) you have to sum the intensities around the circle. But to calculate a given intensity, (Eq. 5) you need to sum the weights!! You shouldn't need weight information to calculate an intensity value! Argh. It's late, and i'm frustrated. Dave. [edited by - Stimpson on June 19, 2002 10:17:40 AM]

### #2greeneggs  Members   -  Reputation: 130

Like
Likes
Like

Posted 19 June 2002 - 04:32 AM

I think you don't quite understand the DNKY paper. DNKY's paper only deals with storing the skylight distribution. The methodology assumes you have a method for calculating the skylight everywhere (e.g., some complicated model, or photographs) and suggests using basis functions to store the data efficiently. So,

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]

### #3Yann L  Members   -  Reputation: 1802

Like
Likes
Like

Posted 19 June 2002 - 05:26 AM

Greeneggs is right. This paper presents a kind of compression system for sky gradients. The method should work rather well, but be aware of the fact, that you''ll need a *very* high amount of samples to get a good result.

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

### #4Stimpson  Members   -  Reputation: 122

Like
Likes
Like

Posted 19 June 2002 - 05:40 AM

Thanks guys I knew i was missing something major! I was originally going to implement "A practical analytic model for daylight" but i didn''t actually like the look of the examples given.

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

### #5bishop_pass  Members   -  Reputation: 109

Like
Likes
Like

Posted 19 June 2002 - 06:08 AM

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.

### #6Yann L  Members   -  Reputation: 1802

Like
Likes
Like

Posted 19 June 2002 - 07:02 AM

OK, I'll try to dig it out. As I said, it's been a while I tried that. It was only a quick'n'dirty hack, just to see if it suited my needs, more or less wrapped around Brian Smits sources. All in all, the results were similar to the images in the paper, but the colours were quite un-natural, IMO.

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

### #7okonomiyaki  Members   -  Reputation: 548

Like
Likes
Like

Posted 19 June 2002 - 07:08 AM

If you guys don''t know about it-
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:

Likes

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.

### #9Stimpson  Members   -  Reputation: 122

Like
Likes
Like

Posted 26 June 2002 - 04:31 AM

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.

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

### #10Yann L  Members   -  Reputation: 1802

Like
Likes
Like

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:

Likes

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 normalinline 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?

### #12greeneggs  Members   -  Reputation: 130

Like
Likes
Like

Posted 26 June 2002 - 08:48 AM

The paper "Trichromatic model of daylight variation" goes into a little more detail on the correct way of scaling the spectral data, and converting to RGB. It looks like they convert the spectral distribution to eight temperatures 4000, 4500, 5000, 5500, 8000, 8500, 9000, 9500 degrees K (this conversion depends on the daylight curve used, so you should use the "trichromatic.." expressions for M1 and M2 rather than the "practical.." expressions -- they are very different!). Then table 5 gives the scaling conversions from temperatures to XYZ, which you convert to RGB by a matrix multiplication.

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]

### #13Yann L  Members   -  Reputation: 1802

Like
Likes
Like

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]

### #14Stimpson  Members   -  Reputation: 122

Like
Likes
Like

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.

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

### #15greeneggs  Members   -  Reputation: 130

Like
Likes
Like

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   XG = -0.969256, 1.875992, 0.041556   YB    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);}surfaceskycolor (  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]

### #16Yann L  Members   -  Reputation: 1802

Like
Likes
Like

Posted 28 June 2002 - 12:33 PM

Hmm, your matrix seems way off, where did you get it from ? [edited answer: from 1958, that explains it *grin*]

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  XG = -0.969256  1.875991  0.041556 YB    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]

### #17greeneggs  Members   -  Reputation: 130

Like
Likes
Like

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?

### #18Yann L  Members   -  Reputation: 1802

Like
Likes
Like

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

### #19greeneggs  Members   -  Reputation: 130

Like
Likes
Like

Posted 29 June 2002 - 10:11 AM

Here are some images from the skylight algorithm (turbidity 2, latitude 20; longitude 0, early April; linear scaling the RGB components by 1/15000 but no exponential exposure function).

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.

### #20Yann L  Members   -  Reputation: 1802

Like
Likes
Like

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

Old topic!

Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

PARTNERS