Page MenuHome
Paste P222

BlackbodyApproxTest.cpp
ActivePublic

Authored by Sv. Lockal (lockal) on May 4 2015, 12:32 PM.
#include <iostream>
#include <iomanip>
#include <cmath>
#include <cstdio>
using namespace std;
struct float3 { float x, y, z; };
float3 make_float3(float x, float y, float z)
{
float3 a = {x, y, z};
return a;
}
float3 max(float3 a, float3 b)
{
return make_float3(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z));
}
float3 operator*(const float3 a, float f)
{
return make_float3(a.x*f, a.y*f, a.z*f);
}
float3 operator/=(float3& a, float f)
{
float invf = 1.0f/f;
return a = a * invf;
}
float3 xyz_to_rgb(float x, float y, float z)
{
return make_float3(3.240479f * x + -1.537150f * y + -0.498535f * z,
-0.969256f * x + 1.875991f * y + 0.041556f * z,
0.055648f * x + -0.204043f * y + 1.057311f * z);
}
float linear_rgb_to_gray(float3 c)
{
return c.x*0.2126f + c.y*0.7152f + c.z*0.0722f;
}
#define BB_DRAPER 800.0f
#define BB_MAX_TABLE_RANGE 12000.0f
#define BB_TABLE_XPOWER 1.5f
#define BB_TABLE_YPOWER 5.0f
#define BB_TABLE_SPACING 2.0f
float3 bb(float temperature) {
const float cie_colour_match[81][3] = {
{0.0014f,0.0000f,0.0065f}, {0.0022f,0.0001f,0.0105f}, {0.0042f,0.0001f,0.0201f},
{0.0076f,0.0002f,0.0362f}, {0.0143f,0.0004f,0.0679f}, {0.0232f,0.0006f,0.1102f},
{0.0435f,0.0012f,0.2074f}, {0.0776f,0.0022f,0.3713f}, {0.1344f,0.0040f,0.6456f},
{0.2148f,0.0073f,1.0391f}, {0.2839f,0.0116f,1.3856f}, {0.3285f,0.0168f,1.6230f},
{0.3483f,0.0230f,1.7471f}, {0.3481f,0.0298f,1.7826f}, {0.3362f,0.0380f,1.7721f},
{0.3187f,0.0480f,1.7441f}, {0.2908f,0.0600f,1.6692f}, {0.2511f,0.0739f,1.5281f},
{0.1954f,0.0910f,1.2876f}, {0.1421f,0.1126f,1.0419f}, {0.0956f,0.1390f,0.8130f},
{0.0580f,0.1693f,0.6162f}, {0.0320f,0.2080f,0.4652f}, {0.0147f,0.2586f,0.3533f},
{0.0049f,0.3230f,0.2720f}, {0.0024f,0.4073f,0.2123f}, {0.0093f,0.5030f,0.1582f},
{0.0291f,0.6082f,0.1117f}, {0.0633f,0.7100f,0.0782f}, {0.1096f,0.7932f,0.0573f},
{0.1655f,0.8620f,0.0422f}, {0.2257f,0.9149f,0.0298f}, {0.2904f,0.9540f,0.0203f},
{0.3597f,0.9803f,0.0134f}, {0.4334f,0.9950f,0.0087f}, {0.5121f,1.0000f,0.0057f},
{0.5945f,0.9950f,0.0039f}, {0.6784f,0.9786f,0.0027f}, {0.7621f,0.9520f,0.0021f},
{0.8425f,0.9154f,0.0018f}, {0.9163f,0.8700f,0.0017f}, {0.9786f,0.8163f,0.0014f},
{1.0263f,0.7570f,0.0011f}, {1.0567f,0.6949f,0.0010f}, {1.0622f,0.6310f,0.0008f},
{1.0456f,0.5668f,0.0006f}, {1.0026f,0.5030f,0.0003f}, {0.9384f,0.4412f,0.0002f},
{0.8544f,0.3810f,0.0002f}, {0.7514f,0.3210f,0.0001f}, {0.6424f,0.2650f,0.0000f},
{0.5419f,0.2170f,0.0000f}, {0.4479f,0.1750f,0.0000f}, {0.3608f,0.1382f,0.0000f},
{0.2835f,0.1070f,0.0000f}, {0.2187f,0.0816f,0.0000f}, {0.1649f,0.0610f,0.0000f},
{0.1212f,0.0446f,0.0000f}, {0.0874f,0.0320f,0.0000f}, {0.0636f,0.0232f,0.0000f},
{0.0468f,0.0170f,0.0000f}, {0.0329f,0.0119f,0.0000f}, {0.0227f,0.0082f,0.0000f},
{0.0158f,0.0057f,0.0000f}, {0.0114f,0.0041f,0.0000f}, {0.0081f,0.0029f,0.0000f},
{0.0058f,0.0021f,0.0000f}, {0.0041f,0.0015f,0.0000f}, {0.0029f,0.0010f,0.0000f},
{0.0020f,0.0007f,0.0000f}, {0.0014f,0.0005f,0.0000f}, {0.0010f,0.0004f,0.0000f},
{0.0007f,0.0002f,0.0000f}, {0.0005f,0.0002f,0.0000f}, {0.0003f,0.0001f,0.0000f},
{0.0002f,0.0001f,0.0000f}, {0.0002f,0.0001f,0.0000f}, {0.0001f,0.0000f,0.0000f},
{0.0001f,0.0000f,0.0000f}, {0.0001f,0.0000f,0.0000f}, {0.0000f,0.0000f,0.0000f}
};
float X = 0;
float Y = 0;
float Z = 0;
const double c1 = 3.74183e-16; // 2*pi*h*c^2, W*m^2
const double c2 = 1.4388e-2; // h*c/k, m*K
// h is Planck's const, k is Boltzmann's
const float dlambda = 5.0f * 1e-9f; // in meters
/* from OSL "spectrum_to_XYZ" */
for(int n = 0; n < 81; ++n) {
float lambda = 380.0f + 5.0f * n;
double wlm = lambda * 1e-9f; // Wavelength in meters
// N.B. spec_intens returns result in W/m^2 but it's a differential,
// needs to be scaled by dlambda!
float spec_intens = float((c1 * pow(wlm, -5.0)) / (exp(c2 / (wlm * temperature)) -1.0));
float Me = spec_intens * dlambda;
X += Me * cie_colour_match[n][0];
Y += Me * cie_colour_match[n][1];
Z += Me * cie_colour_match[n][2];
}
/* Convert from xyz color space */
float3 col = xyz_to_rgb(X, Y, Z);
/* Clamp to zero if values are smaller */
col = max(col, make_float3(0.0f, 0.0f, 0.0f));
float l = linear_rgb_to_gray(col);
if(l != 0.0f)
col /= l;
return col;
}
float3 bbx(float t) {
/* Calculate color in range 800..12000 using an approximation
* a/x+bx+c for R and G and ((at + b)t + c)t + d) for B
* Max absolute error for RGB is (0.00095, 0.00077, 0.00057),
* which is enough to get the same 8 bit/channel color.
*/
const float rc[6][3] = {
{ 2.52432244e+03f, -1.06185848e-03f, 3.11067539e+00f },
{ 3.37763626e+03f, -4.34581697e-04f, 1.64843306e+00f },
{ 4.10671449e+03f, -8.61949938e-05f, 6.41423749e-01f },
{ 4.66849800e+03f, 2.85655028e-05f, 1.29075375e-01f },
{ 4.60124770e+03f, 2.89727618e-05f, 1.48001316e-01f },
{ 3.78765709e+03f, 9.36026367e-06f, 3.98995841e-01f },
};
const float gc[6][3] = {
{ -7.50343014e+02f, 3.15679613e-04f, 4.73464526e-01f },
{ -1.00402363e+03f, 1.29189794e-04f, 9.08181524e-01f },
{ -1.22075471e+03f, 2.56245413e-05f, 1.20753416e+00f },
{ -1.42546105e+03f, -4.01730887e-05f, 1.44002695e+00f },
{ -1.18134453e+03f, -2.18913373e-05f, 1.30656109e+00f },
{ -5.00279505e+02f, -4.59745390e-06f, 1.09090465e+00f },
};
const float bc[6][4] = {
{ 0.0f, 0.0f, 0.0f, 0.0f }, /* zeros should be optimized by compiler */
{ 0.0f, 0.0f, 0.0f, 0.0f },
{ 0.0f, 0.0f, 0.0f, 0.0f },
{ -2.02524603e-11f, 1.79435860e-07f, -2.60561875e-04f, -1.41761141e-02f },
{ -2.22463426e-13f, -1.55078698e-08f, 3.81675160e-04f, -7.30646033e-01f },
{ 6.72595954e-13f, -2.73059993e-08f, 4.24068546e-04f, -7.52204323e-01f },
};
if (t >= 12000.0f) return make_float3(0.826270103f, 0.994478524f, 1.56626022);
/* Define a macro to reduce stack usage for nvcc */
#define MAKE_BB_RGB(i) make_float3(\
rc[i][0] / t + rc[i][1] * t + rc[i][2],\
gc[i][0] / t + gc[i][1] * t + gc[i][2],\
((bc[i][0] * t + bc[i][1]) * t + bc[i][2]) * t + bc[i][3])
if (t >= 6365.0f) return MAKE_BB_RGB(5);
if (t >= 3315.0f) return MAKE_BB_RGB(4);
if (t >= 1902.0f) return MAKE_BB_RGB(3);
if (t >= 1449.0f) return MAKE_BB_RGB(2);
if (t >= 1167.0f) return MAKE_BB_RGB(1);
if (t >= 965.0f) return MAKE_BB_RGB(0);
/* For 800 <= t < 965 color does not change in OSL, so keep color the same */
return make_float3(4.70366907f, 0.0f, 0.0f);
}
int main()
{
for (int i = 800; i <= 12000; i++) {
float3 c1 = bb(i);
float3 c2 = bbx(i);
cout << std::setprecision(9);
cout << i << "\t" << c1.x << "\t" << c1.y << "\t" << c1.z;
cout << "\t" << c2.x << "\t" << c2.y << "\t" << c2.z;
cout << "\t" << c2.x - c1.x << "\t" << c2.y - c1.y << "\t" << c2.z - c1.z;
cout << endl;
}
return 0;
}

Event Timeline

Sv. Lockal (lockal) changed the title of this paste from untitled to BlackbodyApproxTest.cpp.
Sv. Lockal (lockal) updated the paste's language from autodetect to cpp.
Sv. Lockal (lockal) added a project: Cycles.