¿Por qué mi cálculo del color del cielo en Mathematica es incorrecto?

17

Estoy tratando de implementar un algoritmo para calcular el color del cielo basado en este documento (modelo de Pérez). Antes de comenzar a programar el sombreador, quería probar el concepto en Mathematica. Ya hay algunos problemas que no puedo solucionar. Quizás alguien ya haya implementado el algoritmo.

Comencé con ecuaciones para las luminancias zenitales absolutas Yz, xzy yzcomo se propone en el documento (página 22). Los valores para Yzparecen ser razonables. El siguiente diagrama muestra Yzen función de la distancia cenital del sol para una turbidez Tde 5:

Yz (z)

La función gamma (cenit, azimut, solarzenith, solarazimuth) calcula el ángulo entre un punto con la distancia cenital dada y el azimut y el sol en la posición dada. Esta función también parece funcionar. El siguiente diagrama muestra este ángulo para solarzenith=0.5y solarazimuth=0. zenithcrece de arriba hacia abajo (0 a Pi / 2), azimuthcrece de izquierda a derecha (-Pi a Pi). Puede ver claramente la posición del sol (el punto brillante, el ángulo se convierte en cero):

gamma (cenit, azimut, 0.5,0)

La función de Pérez (F) y los coeficientes se han implementado como se indica en el documento. Entonces los valores de color Yxy deberían ser absolute value * F(z, gamma) / F(0, solarzenith). Espero que esos valores estén dentro del rango [0,1]. Sin embargo, este no es el caso para el componente Y (consulte la actualización a continuación para obtener más detalles). Aquí hay algunos valores de muestra:

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

Aquí está el resultado actual:

Imagen RGB

El Cuaderno de Mathematica con todos los cálculos se puede encontrar aquí y la versión en PDF aquí .

¿Alguien tiene una idea de lo que tengo que cambiar para obtener los mismos resultados que en el documento?

C como código

// 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();
    }
}

Solución

Como prometí, escribí un artículo de blog sobre cómo representar el cielo. Lo puedes encontrar aquí .

Nico Schertler
fuente
Sospecho que más personas aquí podrían ayudarlo si tratara de implementar el algoritmo en código real (sombreador u otro) en lugar de en Mathematica.
Tetrad
2
Hay un SE de Mathematica , aunque tendrías que consultar sus preguntas frecuentes para ver si tu pregunta es sobre el tema allí.
MichaelHouse
Bueno, la pregunta no es realmente sobre Mathematica, sino sobre el algoritmo. Agregué la versión PDF del cuaderno, para que todos puedan leerlo. Estoy seguro de que la sintaxis es comprensible para un programador común y probablemente más comprensible que el código de sombreador.
Nico Schertler
@NicoSchertler: El problema es que no creo que mucha gente aquí entienda la sintaxis de Mathematica. Probablemente tendrá más suerte si reescribe su código en un lenguaje similar a C o Python, al menos a los fines de esta pregunta.
Panda Pyjama
2
La pregunta es realmente demasiado localizada y podría cerrarse, pero gracias por el enlace en papel, es interesante.
sam hocevar

Respuestas:

4

Hay dos errores en la matriz utilizada para xz: 1.00166 debería ser 0.00166 y 0.6052 debería ser 0.06052.

sam hocevar
fuente
Gracias por la corrección. El resultado ahora se ve mejor, pero no puede ser correcto. Le agradecería si considera la pregunta actualizada.
Nico Schertler
-2

Parece que podría ser un problema de escala de valor de color?

boobami
fuente
2
Aunque su suposición puede ser correcta, sería más útil proporcionar explicaciones adicionales. Como no puede responder a toda la pregunta, lo que escribió debe ser un comentario debajo de la pregunta.
danijar
Esto no proporciona una respuesta a la pregunta. Para criticar o solicitar una aclaración de un autor, deje un comentario debajo de su publicación; siempre puede comentar sus propias publicaciones, y una vez que tenga suficiente reputación podrá comentar cualquier publicación .
MichaelHouse
1
No entiendo por qué las sugerencias no son toleradas aquí en absoluto. Si nos fijamos en la solución anterior, es un problema de valor. Señalar a la gente en la dirección correcta es una mejor manera de aprender que dar soluciones exactas todo el tiempo, ¿no? Y no, no puedo comentar debajo de su pregunta porque no tengo permitido hacerlo. Es por eso que comenté aquí. Pero gracias por el descenso de categoría. Muy amable de su parte y muy alentador para gente nueva como yo. Gracias.
boobami