From f4a19753515517e39cb61bd2447bae288d6e6a40 Mon Sep 17 00:00:00 2001 From: Mikko Rasa Date: Sun, 13 Jul 2014 16:48:14 +0300 Subject: [PATCH] Improve lighting color computations This produces a bit too bluish skylight and a bit too reddish sunlight, but it's a step in the right direction. --- source/3d/layout.cpp | 182 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 178 insertions(+), 4 deletions(-) diff --git a/source/3d/layout.cpp b/source/3d/layout.cpp index 29e6bb4..f30e427 100644 --- a/source/3d/layout.cpp +++ b/source/3d/layout.cpp @@ -1,3 +1,4 @@ +#include #include "beamgate.h" #include "layout.h" #include "signal.h" @@ -9,6 +10,128 @@ using namespace std; using namespace Msp; +namespace { + +double optical_thickness(double zenith_cosine) +{ + double r = 635; // Radius of Earth divided by atmospheric thickness at zenith + double rc = r*zenith_cosine; + return sqrt(2*r+1+rc*rc)-rc; +} + +GL::Color spectrum_to_rgb(const double *spectrum) +{ + // Data from http://cvrl.ioo.ucl.ac.uk/cmfs.htm + static double color_matching[] = + { + 3.769647E-03, 4.146161E-04, 1.847260E-02, + 9.382967E-03, 1.059646E-03, 4.609784E-02, + 2.214302E-02, 2.452194E-03, 1.096090E-01, + 4.742986E-02, 4.971717E-03, 2.369246E-01, + 8.953803E-02, 9.079860E-03, 4.508369E-01, + 1.446214E-01, 1.429377E-02, 7.378822E-01, + 2.035729E-01, 2.027369E-02, 1.051821E+00, + 2.488523E-01, 2.612106E-02, 1.305008E+00, + 2.918246E-01, 3.319038E-02, 1.552826E+00, + 3.227087E-01, 4.157940E-02, 1.748280E+00, + 3.482554E-01, 5.033657E-02, 1.917479E+00, + 3.418483E-01, 5.743393E-02, 1.918437E+00, + 3.224637E-01, 6.472352E-02, 1.848545E+00, + 2.826646E-01, 7.238339E-02, 1.664439E+00, + 2.485254E-01, 8.514816E-02, 1.522157E+00, + 2.219781E-01, 1.060145E-01, 1.428440E+00, + 1.806905E-01, 1.298957E-01, 1.250610E+00, + 1.291920E-01, 1.535066E-01, 9.991789E-01, + 8.182895E-02, 1.788048E-01, 7.552379E-01, + 4.600865E-02, 2.064828E-01, 5.617313E-01, + 2.083981E-02, 2.379160E-01, 4.099313E-01, + 7.097731E-03, 2.850680E-01, 3.105939E-01, + 2.461588E-03, 3.483536E-01, 2.376753E-01, + 3.649178E-03, 4.277595E-01, 1.720018E-01, + 1.556989E-02, 5.204972E-01, 1.176796E-01, + 4.315171E-02, 6.206256E-01, 8.283548E-02, + 7.962917E-02, 7.180890E-01, 5.650407E-02, + 1.268468E-01, 7.946448E-01, 3.751912E-02, + 1.818026E-01, 8.575799E-01, 2.438164E-02, + 2.405015E-01, 9.071347E-01, 1.566174E-02, + 3.098117E-01, 9.544675E-01, 9.846470E-03, + 3.804244E-01, 9.814106E-01, 6.131421E-03, + 4.494206E-01, 9.890228E-01, 3.790291E-03, + 5.280233E-01, 9.994608E-01, 2.327186E-03, + 6.133784E-01, 9.967737E-01, 1.432128E-03, + 7.016774E-01, 9.902549E-01, 8.822531E-04, + 7.967750E-01, 9.732611E-01, 5.452416E-04, + 8.853376E-01, 9.424569E-01, 3.386739E-04, + 9.638388E-01, 8.963613E-01, 2.117772E-04, + 1.051011E+00, 8.587203E-01, 1.335031E-04, + 1.109767E+00, 8.115868E-01, 8.494468E-05, + 1.143620E+00, 7.544785E-01, 5.460706E-05, + 1.151033E+00, 6.918553E-01, 3.549661E-05, + 1.134757E+00, 6.270066E-01, 2.334738E-05, + 1.083928E+00, 5.583746E-01, 1.554631E-05, + 1.007344E+00, 4.895950E-01, 1.048387E-05, + 9.142877E-01, 4.229897E-01, 0.000000E+00, + 8.135565E-01, 3.609245E-01, 0.000000E+00, + 6.924717E-01, 2.980865E-01, 0.000000E+00, + 5.755410E-01, 2.416902E-01, 0.000000E+00, + 4.731224E-01, 1.943124E-01, 0.000000E+00, + 3.844986E-01, 1.547397E-01, 0.000000E+00, + 2.997374E-01, 1.193120E-01, 0.000000E+00, + 2.277792E-01, 8.979594E-02, 0.000000E+00, + 1.707914E-01, 6.671045E-02, 0.000000E+00, + 1.263808E-01, 4.899699E-02, 0.000000E+00, + 9.224597E-02, 3.559982E-02, 0.000000E+00, + 6.639960E-02, 2.554223E-02, 0.000000E+00, + 4.710606E-02, 1.807939E-02, 0.000000E+00, + 3.292138E-02, 1.261573E-02, 0.000000E+00, + 2.262306E-02, 8.661284E-03, 0.000000E+00, + 1.575417E-02, 6.027677E-03, 0.000000E+00, + 1.096778E-02, 4.195941E-03, 0.000000E+00, + 7.608750E-03, 2.910864E-03, 0.000000E+00, + 5.214608E-03, 1.995557E-03, 0.000000E+00, + 3.569452E-03, 1.367022E-03, 0.000000E+00, + 2.464821E-03, 9.447269E-04, 0.000000E+00, + 1.703876E-03, 6.537050E-04, 0.000000E+00, + 1.186238E-03, 4.555970E-04, 0.000000E+00, + 8.269535E-04, 3.179738E-04, 0.000000E+00, + 5.758303E-04, 2.217445E-04, 0.000000E+00, + 4.058303E-04, 1.565566E-04, 0.000000E+00, + 2.856577E-04, 1.103928E-04, 0.000000E+00, + 2.021853E-04, 7.827442E-05, 0.000000E+00, + 1.438270E-04, 5.578862E-05, 0.000000E+00, + 1.024685E-04, 3.981884E-05, 0.000000E+00, + 7.347551E-05, 2.860175E-05, 0.000000E+00, + 5.259870E-05, 2.051259E-05, 0.000000E+00, + 3.806114E-05, 1.487243E-05, 0.000000E+00, + 2.758222E-05, 1.080001E-05, 0.000000E+00, + 2.004122E-05, 7.863920E-06, 0.000000E+00, + 1.458792E-05, 5.736935E-06, 0.000000E+00, + 1.068141E-05, 4.211597E-06, 0.000000E+00, + 7.857521E-06, 3.106561E-06, 0.000000E+00, + 5.768284E-06, 2.286786E-06, 0.000000E+00, + 4.259166E-06, 1.693147E-06, 0.000000E+00, + 3.167765E-06, 1.262556E-06, 0.000000E+00, + 2.358723E-06, 9.422514E-07, 0.000000E+00, + 1.762465E-06, 7.053860E-07, 0.000000E+00, + }; + + // Compute CIE XYZ tristimulus values + double x = 0; + double y = 0; + double z = 0; + for(unsigned i=0; i<89; ++i) + { + x += spectrum[i]*color_matching[i*3]; + y += spectrum[i]*color_matching[i*3+1]; + z += spectrum[i]*color_matching[i*3+2]; + } + + // Compute linear RGB color + return GL::Color(3.2406*x-1.5372*y-0.4986*z, -0.9689*x+1.8758*y+0.0415*z, 0.0557*x-0.2040*y+1.0570*z); +} + +} + namespace R2C2 { Layout3D::Layout3D(Layout &l): @@ -17,7 +140,8 @@ Layout3D::Layout3D(Layout &l): { // South, 15° from zenith sun.set_position(GL::Vector4(0, -0.259, 0.966, 0)); - sun.set_diffuse(GL::Color(0.7)); + sun.set_diffuse(GL::Color(0.0)); + sun.set_specular(GL::Color(0.0)); lighting.set_ambient(GL::Color(0.2)); lighting.attach(0, sun); @@ -75,9 +199,59 @@ void Layout3D::tick() Vector diff = Vector(sun.get_position())-sun_dir; if(diff.norm()>0.0025f) sun.set_position(GL::Vector4(sun_dir, 0.0f)); - // TODO replace these very unscientific guesses with correct formulas - lighting.set_ambient(GL::Color(pow(0.2f*(sun_dir.z+1.0f), 2.0f))); - sun.set_diffuse(GL::Color(0.7f*max(sun_dir.z, 0.0f))); + + double T = 5777; + double h = 6.62606957e-34; // Planck constant + double k = 1.3806488e-23; // Boltzmann constant + double c = 299792458; // Speed of light in vacuum + double wavelengths[89]; + for(unsigned i=0; i<89; ++i) + wavelengths[i] = (390+i*5)/1e9; + + // Compute the spectrum of sunlight in space + double incoming_spectrum[89]; + for(unsigned i=0; i<89; ++i) + incoming_spectrum[i] = (2*h*pow(c, 2))/pow(wavelengths[i], 5)/(exp(h*c/(wavelengths[i]*k*T))-1); + + // Apply rayleigh scattering to get the spectrum at the surface + double sun_thickness = optical_thickness(sun_dir.z); + double transmitted_spectrum[89]; + for(unsigned i=0; i<89; ++i) + { + double f = 7.9e-27/pow(wavelengths[i], 4); + transmitted_spectrum[i] = incoming_spectrum[i]*exp(-f*sun_thickness); + } + + GL::Color sunlight_color = spectrum_to_rgb(transmitted_spectrum)*2e-15; + /* When the sun is setting, the amount of direct illlumination is + proportional to the fraction of sun still over the horizon. This formula + assumes a square sun, but the difference is hardly noticeable. */ + double cutoff = sun_dir.z>0.005 ? 1 : sun_dir.z>-0.005 ? (sun_dir.z+0.005)*100 : 0; + sun.set_diffuse(sunlight_color*cutoff); + + double skylight_spectrum[89] = {}; + double total_weight = 0; + for(int x=-2; x<=2; ++x) + for(int y=-2; y<=2; ++y) + for(int z=0; z<=2; ++z) + { + if(abs(x)!=2 && abs(y)!=2 && abs(z)!=2) + continue; + + double l_sq = x*x+y*y+z*z; + double l = sqrt(l_sq); + double view_thickness = optical_thickness(z/l); + double w = z/l/l_sq; + for(unsigned i=0; i<89; ++i) + { + double f = 7.9e-27/pow(wavelengths[i], 4); + skylight_spectrum[i] += w*incoming_spectrum[i]*(exp(-f*sun_thickness)-exp(-f*view_thickness))/(1-sun_thickness/view_thickness); + } + total_weight += w; + } + + GL::Color skylight_color = spectrum_to_rgb(skylight_spectrum)*(2e-15/total_weight); + lighting.set_ambient(skylight_color*0.5); } void Layout3D::object_added(Object &o) -- 2.43.0