Por que meu cálculo de cor do céu no Mathematica está incorreto?

17

Estou tentando implementar um algoritmo para calcular a cor do céu com base neste artigo (modelo de Perez). Antes de começar a programar o shader, eu queria testar o conceito no Mathematica. Já existem alguns problemas que não consigo me livrar. Talvez alguém já tenha implementado o algoritmo.

Comecei com equações para as luminâncias zenital absolutos Yz, xze yzcomo proposto no papel (página 22). Os valores para Yzparecem razoáveis. O diagrama a seguir mostra Yzcomo uma função da distância zenital do sol para uma turbidez Tde 5:

Yz (z)

A função gama (zenith, azimute, solarzenith, solarazimuth) calcula o ângulo entre um ponto com a distância zenital especificada e o azimute e o sol na posição especificada. Essa função parece funcionar também. O diagrama a seguir mostra esse ângulo para solarzenith=0.5e solarazimuth=0. zenithcresce de cima para baixo (0 a Pi / 2), azimuthcresce da esquerda para a direita (-Pi para Pi). Você pode ver claramente a posição do sol (o ponto brilhante, o ângulo se torna zero):

gama (zênite, azimute, 0,5,0)

A função de Perez (F) e os coeficientes foram implementados conforme apresentado no artigo. Os valores de cor Yxy devem ser absolute value * F(z, gamma) / F(0, solarzenith). Eu espero que esses valores estejam dentro do intervalo [0,1]. No entanto, esse não é o caso do componente Y (consulte a atualização abaixo para obter detalhes). Aqui estão alguns valores de amostra:

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

Aqui está o resultado atual:

Imagem RGB

O Notebook Mathematica com todos os cálculos pode ser encontrado aqui e a versão em PDF aqui .

Alguém tem uma idéia do que preciso alterar para obter os mesmos resultados que no artigo?

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

Solução

Como prometido, escrevi um artigo de blog sobre renderização do céu. Você pode encontrá-lo aqui .

Nico Schertler
fonte
Eu suspeito que mais pessoas aqui possam ajudá-lo se você tentar implementar o algoritmo no código real (shader ou não) em vez de no Mathematica.
Tétrada
2
Existe um Mathematica SE , embora você precise verificar as perguntas frequentes para ver se sua pergunta está sobre o assunto.
Michaelhouse
Bem, a questão não é realmente sobre o Mathematica, mas sobre o algoritmo. Eu adicionei a versão em PDF do notebook para que todos possam lê-lo. Tenho certeza de que a sintaxe é compreensível para um programador comum e provavelmente mais compreensível que o código shader.
Nico Schertler 27/03
@NicoSchertler: O problema é que acho que muitas pessoas aqui não entendem a sintaxe do Mathematica. Você provavelmente terá mais sorte se reescrever seu código em uma linguagem semelhante a C ou Python, pelo menos para os fins desta pergunta.
Panda Pyjama
2
A questão é realmente muito localizada e pode ser fechada, mas obrigado pelo link em papel, é interessante.
sam hocevar 29/03

Respostas:

4

Existem dois erros na matriz usados ​​para xz: 1,00166 deve ser 0,00166 e 0,6052 deve ser 0,06052.

sam hocevar
fonte
Obrigado pela correção. O resultado agora parece melhor, mas não pode estar correto. Eu apreciaria se você considerasse a pergunta atualizada.
Nico Schertler 29/03
-2

Parece que pode ser um problema de escala de valor de cor?

boobami
fonte
2
Embora sua suposição possa estar correta, seria mais útil fornecer explicações adicionais. Como você não pode responder a toda a pergunta, o que você escreveu deve ser um comentário.
31513 danijar
Isso não fornece uma resposta para a pergunta. Para criticar ou solicitar esclarecimentos a um autor, deixe um comentário abaixo da postagem - você sempre pode comentar em suas próprias postagens e, quando tiver reputação suficiente , poderá comentar em qualquer post .
Michaelhouse
1
Não entendo por que as sugestões não são toleradas aqui. Se você olhar para a solução acima, é uma questão de valor. Apontar as pessoas na direção certa é a melhor maneira de aprender do que fornecer soluções exatas o tempo todo, não é? E não, não posso comentar abaixo da pergunta dele porque não tenho permissão. Por isso, comentei aqui. Mas obrigado pelo rebaixamento. Muito gentil da sua parte e muito encorajador para novas pessoas como eu. Obrigado.
boobami