Pourquoi mon calcul de la couleur du ciel dans Mathematica est-il incorrect?

J'essaie d'implémenter un algorithme pour calculer la couleur du ciel sur la base de cet article (modèle de Perez). Avant de commencer à programmer le shader, je voulais tester le concept dans Mathematica. Il y a déjà des problèmes dont je ne peux pas me débarrasser. Peut-être que quelqu'un a déjà implémenté l'algorithme.

J'ai commencé avec des équations pour les luminances zénitales absolues Yz, xzet yzcomme proposé dans l'article (page 22). Les valeurs de Yzsemblent raisonnables. Le diagramme suivant montre Yzen fonction de la distance zénitale du soleil pour une turbidité Tde 5:

La fonction gamma (zénith, azimut, solarzenith, solarazimuth) calcule l'angle entre un point avec la distance zénitale donnée et l'azimut et le soleil à la position donnée. Cette fonction semble également fonctionner. Le diagramme suivant montre cet angle pour solarzenith=0.5et solarazimuth=0. zenithcroît de haut en bas (0 à Pi / 2), azimuthcroît de gauche à droite (-Pi à Pi). Vous pouvez voir clairement la position du soleil (la tache lumineuse, l'angle devient nul):

La fonction de Perez (F) et les coefficients ont été mis en œuvre comme indiqué dans l'article. Ensuite, les valeurs de couleur Yxy devraient être absolute value * F(z, gamma) / F(0, solarzenith). Je m'attends à ce que ces valeurs soient dans la plage [0,1]. Cependant, ce n'est pas le cas pour le composant Y (voir la mise à jour ci-dessous pour plus de détails). Voici quelques exemples de valeurs:

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

Voici le résultat actuel:

Le cahier Mathematica avec tous les calculs peut être trouvé ici et la version PDF ici .

Quelqu'un at-il une idée de ce que je dois changer pour obtenir les mêmes résultats que dans le document?

C comme le 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

Comme promis, j'ai écrit un article de blog sur le rendu du ciel. Vous pouvez le trouver ici .

Il y a deux erreurs dans la matrice utilisée pour xz: 1,00166 devrait être 0,00166 et 0,6052 devrait être 0,06052.

Il semble que ce soit un problème de mise à l'échelle des valeurs de couleur?