Warum ist meine Himmelsfarbenberechnung in Mathematica falsch?

17

Ich versuche einen Algorithmus zu implementieren, um die Himmelsfarbe basierend auf diesem Papier (Perez 'Modell) zu berechnen . Bevor ich den Shader programmiere, wollte ich das Konzept in Mathematica testen. Es gibt bereits einige Probleme, die ich nicht loswerden kann. Vielleicht hat jemand den Algorithmus schon implementiert.

Ich begann mit Gleichungen für den absoluten zenital Helligkeiten Yz, xzund yzwie in dem Papier (Seite 22) vorgeschlagen. Die Werte für Yzscheinen vernünftig zu sein. Das folgende Diagramm zeigt Yzin Abhängigkeit vom Zenitabstand der Sonne bei einer Trübung Tvon 5:

Yz (z)

Die Funktion Gamma (Zenit, Azimut, Solarzenit, Solarazimut) berechnet den Winkel zwischen einem Punkt mit dem angegebenen Zenitabstand und Azimut und der Sonne an der angegebenen Position. Diese Funktion scheint auch zu funktionieren. Das folgende Diagramm zeigt diesen Winkel für solarzenith=0.5und solarazimuth=0. zenithwächst von oben nach unten (0 bis Pi / 2), azimuthwächst von links nach rechts (-Pi bis Pi). Sie können den Sonnenstand deutlich sehen (der helle Punkt, der Winkel wird zu Null):

Gamma (Zenit, Azimut, 0,5,0)

Die Perez-Funktion (F) und die Koeffizienten wurden wie im Artikel angegeben implementiert. Dann sollten die Farbwerte Yxy sein absolute value * F(z, gamma) / F(0, solarzenith). Ich erwarte, dass diese Werte im Bereich [0,1] liegen. Dies ist jedoch bei der Y-Komponente nicht der Fall (Details siehe Update unten). Hier sind einige Beispielwerte:

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

Hier ist das aktuelle Ergebnis:

RGB-Bild

Das Mathematica-Notizbuch mit allen Berechnungen finden Sie hier und die PDF-Version hier .

Hat jemand eine Idee, was ich ändern muss, um die gleichen Ergebnisse wie in der Zeitung zu erzielen?

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

Lösung

Wie versprochen habe ich einen Blogartikel über das Rendern des Himmels geschrieben. Sie finden es hier .

Nico Schertler
quelle
Ich vermute, dass mehr Leute hier in der Lage wären, Ihnen zu helfen, wenn Sie versuchen würden, den Algorithmus in tatsächlichem Code (Shader oder auf andere Weise) anstatt in Mathematica zu implementieren.
Tetrad
2
Es gibt eine Mathematica SE , aber Sie müssten deren FAQ überprüfen, um festzustellen, ob Ihre Frage dort drüben zum Thema gehört.
MichaelHouse
Nun, die Frage bezieht sich nicht wirklich auf Mathematica, sondern auf den Algorithmus. Ich habe die PDF-Version des Notizbuchs hinzugefügt, damit jeder sie lesen kann. Ich bin sicher, dass die Syntax für einen gewöhnlichen Programmierer verständlich und wahrscheinlich verständlicher ist als Shader-Code.
Nico Schertler
@NicoSchertler: Das Problem ist, dass ich glaube, dass nicht viele Leute hier die Mathematica-Syntax verstehen. Sie werden wahrscheinlich mehr Glück haben, wenn Sie Ihren Code in einer C- oder Python-ähnlichen Sprache umschreiben, zumindest für die Zwecke dieser Frage.
Panda Pyjama
2
Die Frage ist wirklich zu lokalisiert und könnte geschlossen werden, aber danke für den Link auf dem Papier, es ist interessant.
Sam Hocevar

Antworten:

4

Die verwendete Matrix enthält zwei Fehler xz: 1.00166 sollte 0.00166 sein und 0.6052 sollte 0.06052 sein.

sam hocevar
quelle
Danke für die Korrektur. Das Ergebnis sieht jetzt besser aus, kann aber nicht richtig sein. Ich würde mich freuen, wenn Sie die aktualisierte Frage berücksichtigen.
Nico Schertler
-2

Es scheint, als könnte es sich um ein Problem mit der Farbwertskalierung handeln.

Boobami
quelle
2
Auch wenn Ihre Annahme richtig sein mag, wäre es hilfreicher, zusätzliche Erklärungen zu liefern. Da Sie nicht die ganze Frage beantworten können, sollte das, was Sie geschrieben haben, ein Kommentar unter der Frage sein.
Danijar
Dies gibt keine Antwort auf die Frage. Wenn Sie einen Autor kritisieren oder um Klarstellung bitten möchten, hinterlassen Sie einen Kommentar unter seinem Beitrag. Sie können jederzeit Ihre eigenen Beiträge kommentieren. Wenn Sie über eine ausreichende Reputation verfügen, können Sie jeden Beitrag kommentieren .
MichaelHouse
1
Ich verstehe nicht, warum hier Vorschläge überhaupt nicht toleriert werden. Wenn Sie sich die Lösung oben ansehen, handelt es sich um ein Wertproblem. Menschen in die richtige Richtung zu weisen, ist eine bessere Möglichkeit zu lernen, als ständig genaue Lösungen zu geben, nicht wahr? Und nein, ich kann seine Frage nicht kommentieren, weil ich das nicht darf. Deshalb habe ich hier kommentiert. Aber danke für die Herabstufung. Wirklich nett von dir und sehr ermutigend für neue Leute wie mich. Vielen Dank.
Boobami