Take the 2-minute tour ×
Game Development Stack Exchange is a question and answer site for professional and independent game developers. It's 100% free, no registration required.

I'm trying to implement an algorithm to calculate sky color based on this paper (Perez' model). Before I start programming the shader I wanted to test the concept in Mathematica. There are already some problems I cannot get rid off. Maybe someone already has implemented the algorithm.

I started with equations for the absolute zenital luminances Yz, xz and yz as proposed in the paper (page 22). The values for Yz seem to be reasonable. The following diagram shows Yz as a function of the zenital distance of the sun for a turbidity T of 5:

Yz(z)

The function gamma(zenith,azimuth,solarzenith,solarazimuth) calculates the angle between a point with the given zenital distance and azimuth and the sun at the given position. This function seems to work, too. The following diagram shows this angle for solarzenith=0.5 and solarazimuth=0. zenith grows from top down (0 to Pi/2), azimuth grows from left to right (-Pi to Pi). You can clearly see the position of the sun (the bright spot, angle becomes zero):

gamma(zenith,azimuth,0.5,0)

The Perez function (F) and coefficients have been implemented as given in the paper. Then the color values Yxy should be absolute value * F(z, gamma) / F(0, solarzenith). I expect those values to be within the range [0,1]. However, this is not the case for the Y component (see update below for details). Here are some sample values:

{Y, x, y}
{19.1548, 0.25984, 0.270379}
{10.1932, 0.248629, 0.267739]
{20.0393, 0.268119, 0.280024}

Here is the current result:

RGB image

The Mathematica Notebook with all calculations can be found here and the PDF version here.

Does anyone have an idea what I have to change to get the same results as in the paper?

C like code

// this function returns the zenital Y component for 
// a given solar zenital distance z and turbidity T
float Yz(float z, float T)
{
    return (4.0453 * T - 4.9710)*tan( (4.0f/9-T/120)*(Pi-2*z) ) - 0.2155 * T + 2.4192
}

// returns zenital x component
float xz(float z, float T)
{
    return //matrix calculation, see paper
}

// returns zenital y component
float yz(float z, float T)
{
    return //matrix calculation, see paper
}

// returns the rgb color of a Yxy color
Color RGB(float Y, float x, float y)
{
    Matrix m; //this is a CIE XYZ -> RGB conversion matrix
    Vector v;
    v.x = x/y*Y;
    v.y = Y;
    v.z = (1-x-y)/y*Y;
    v = M * v; //matrix-vector multiplication;
    return Color ( v.x, v.y, v.z );        
}

// returns the 5 coefficients (A-E) for the given turbidity T
float[5] CoeffY(float T)
{
    float[5] result;
    result[0] = 0.1787 * T - 1.4630;
    result[1] = -0.3554 * T + 0.4275;
    ...
    return result;
}

//same for Coeffx and Coeffy

// returns the angle between an observed point and the sun
float PerezGamma(float zenith, float azimuth, float solarzenith, float solarazimuth)
{
    return acos(sin(solarzenith)*sin(zenith)*cos(azimuth-solarazimuth)+cos(solarzenith)*cos(zenith));
}

// evalutes Perez' function F
// the last parameter is a function
float Perez(float zenith, float gamma, float T, t->float[5] coeffs)
{
    return (1+coeffs(T)[0] * exp(coeffs(T)[1]/cos(zenith)) *
           (1+coeffs(T)[2] * exp(coeffs(T)[3]*gamma) + 
            coeffs(T)[4]*pow(cos(gamma),2))
}

// calculates the color for a given point
YxyColor calculateColor(float zenith, float azimuth, float solarzenith, float solarazimuth, float T)
{
    YxyColor c;
    float gamma = PerezGamma(zenith, azimuth, solarzenith, solarazimuth);
    c.Y = Yz(solarzenith, T) * Perez(zenith, gamma, T, CoeffY) / Perez(0, solarzenith, T, CoeffY);
    c.x = xz(solarzenith, T) * Perez(zenith, gamma, T, Coeffx) / Perez(0, solarzenith, T, Coeffx);
    c.y = yz(solarzenith, T) * Perez(zenith, gamma, T, Coeffy) / Perez(0, solarzenith, T, Coeffy); 
    return c;
}

// draws an image of the sky
void DrawImage()
{
    for(float z from 0 to Pi/2) //zenithal distance
    {
        for(float a from -Pi to Pi) //azimuth
        {
            YxyColor c = calculateColor(zenith, azimuth, 1, 0, 5);
            Color rgb = RGB(c.Y, c.x, c.y);
            setNextColor(rgb);
        }
        newline();
    }
}

Solution

As promised I wrote a blog article about rendering the sky. You can find it here.

share|improve this question
    
I suspect that more people here would be able to help you if you were to try to implement the algorithm in actual code (shader or otherwise) instead of in Mathematica. –  Tetrad Mar 27 '13 at 22:16
2  
There is a Mathematica SE, though you would have to check their FAQ to see if your question is on topic over there. –  Byte56 Mar 27 '13 at 22:20
    
Well, the question is not really about Mathematica, but about the algorithm. I added the PDF version of the notebook, so everyone can read it. I'm sure that the syntax is comprehensible for a common programmer and probably more comprehensible than shader code. –  Nico Schertler Mar 27 '13 at 22:58
    
@NicoSchertler: The problem is that I don't think many people in here understand Mathematica syntax. You will probably have more luck if you rewrite your code in a C-like or Python-like language, at least for the purposes of this question. –  Panda Pajama Mar 29 '13 at 9:44
2  
The question is really too localised and might get closed, but thanks for the paper link, it's interesting. –  Sam Hocevar Mar 29 '13 at 10:01
show 3 more comments

2 Answers 2

up vote 4 down vote accepted

There are two errors in the matrix used for xz: 1.00166 should be 0.00166, and 0.6052 should be 0.06052.

share|improve this answer
    
Thanks for the correction. The result now looks better, but can't be correct. I would appreciate if you consider the updated question. –  Nico Schertler Mar 29 '13 at 13:43
add comment

It seems like it might be color value scaling issue?

share|improve this answer
2  
Though your assumption may be correct it would be more helpful to provide additional explanations. Since you can't answer the whole question what you wrote should be a comment under the question. –  danijar Mar 29 '13 at 9:56
    
This does not provide an answer to the question. To critique or request clarification from an author, leave a comment below their post - you can always comment on your own posts, and once you have sufficient reputation you will be able to comment on any post. –  Byte56 Apr 4 '13 at 22:27
    
I don't understand why suggestions are not tolerated here at all. If you look at the solution above, it is a value issue. Pointing folks in the right direction is better way to learn than giving exact solutions all the time isn't it? And no, I can't comment below his question because I am not allowed to. Which is why I commented here. But thanks for the demote. Really nice of you and very encouraging to new folks like me. Thank you. –  boobami Apr 6 '13 at 3:56
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.