Smoke (fire): Move spectrum code from C++ (intern/) to C code (BLI)

This change is for a few reasons:
- it works with color, and (therefore) will need to be color managed, at
some point. This will be much easier to do if the code is closer to the
actual color management code (in Blender's core, so to speak).
- it has nothing to do with the actual fire simulation, as it is just
used to create a lookup table
- it can be reused for other purposes (i.e. in Blender internal
renderer, if people are interrested in a blackbody node à la Cycles)
- cleanup: some functions (`contrain_rgb`, `xyz_to_rgb`) already exist
in BLI

Reviewers: brecht

Reviewed By: brecht

Subscribers: brecht

Differential Revision: https://developer.blender.org/D1719

1 files changed, 138 insertions, 11 deletions

diff --git a/source/blender/blenlib/intern/math_color.c b/source/blender/blenlib/intern/math_color.c
index e149fdb2d26..6d7d24c79c3 100644
--- a/source/blender/blenlib/intern/math_color.c
+++ b/source/blender/blenlib/intern/math_color.c
@@ -455,7 +455,7 @@ void minmax_rgb(short c[3])
 	else if (c[2] < 0) c[2] = 0;
 }
 
-/*If the requested RGB shade contains a negative weight for
+/* If the requested RGB shade contains a negative weight for
  * one of the primaries, it lies outside the color gamut
  * accessible from the given triple of primaries.  Desaturate
  * it by adding white, equal quantities of R, G, and B, enough
@@ -463,21 +463,15 @@ void minmax_rgb(short c[3])
  * components were modified, zero otherwise.*/
 int constrain_rgb(float *r, float *g, float *b)
 {
-	float w;
-
-	/* Amount of white needed is w = - min(0, *r, *g, *b) */
-
-	w = (0 < *r) ? 0 : *r;
-	w = (w < *g) ? w : *g;
-	w = (w < *b) ? w : *b;
-	w = -w;
+	/* Amount of white needed */
+	const float w = -min_ffff(0.0f, *r, *g, *b);
 
 	/* Add just enough white to make r, g, b all positive. */
-
-	if (w > 0) {
+	if (w > 0.0f) {
 		*r += w;
 		*g += w;
 		*b += w;
+
 		return 1; /* Color modified to fit RGB gamut */
 	}
 
@@ -659,3 +653,136 @@ void rgb_to_lab(float r, float g, float b, float *ll, float *la, float *lb)
 	rgb_to_xyz(r, g, b, &x, &y, &z);
 	xyz_to_lab(x, y, z, ll, la, lb);
 }
+
+static void xyz_to_lms(float x, float y, float z, float *l, float *m, float *s)
+{
+	*l =  0.3897f * x + 0.6890f * y - 0.0787f * z;
+	*m = -0.2298f * x + 1.1834f * y + 0.0464f * z;
+	*s = z;
+}
+
+static void lms_to_xyz(float l, float m, float s, float *x, float *y, float *z)
+{
+	*x =  1.9102f * l - 1.1121f * m + 0.2019f * s;
+	*y =  0.3709f * l + 0.6290f * m + 0.0000f * s;
+	*z = s;
+}
+
+static void normalize_rgb(float rgb[3])
+{
+	const float max = max_fff(rgb[0], rgb[1], rgb[2]);
+
+	if (max > 0.0f) {
+		mul_v3_fl(rgb, 1.0f / max);
+	}
+}
+
+/* Color rendering of spectra, adapted from public domain code by John Walker,
+ * http://www.fourmilab.ch/
+ */
+
+static void spectrum_to_xyz(float temperature, float xyz[3])
+{
+	int i;
+	float lambda, x = 0.0f, y = 0.0f, z = 0.0f, xyz_sum;
+
+	/* CIE colour matching functions xBar, yBar, and zBar for wavelengths from
+	 * 380 through 780 nanometers, every 5 nanometers.
+	 * For a wavelength lambda in this range:
+	 *
+	 * cie_colour_match[(lambda - 380) / 5][0] = xBar
+	 * cie_colour_match[(lambda - 380) / 5][1] = yBar
+	 * cie_colour_match[(lambda - 380) / 5][2] = zBar
+	 */
+
+	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}
+	};
+
+	for (i = 0, lambda = 380.0f; lambda < 780.1f; i++, lambda += 5.0f) {
+		/* wavelength in meter */
+		const float wlm = lambda * 1e-9f;
+		const float Me = (3.74183e-16f * powf(wlm, -5.0f)) / (expf(1.4388e-2f / (wlm * temperature)) - 1.0f);
+
+		x += Me * cie_colour_match[i][0];
+		y += Me * cie_colour_match[i][1];
+		z += Me * cie_colour_match[i][2];
+	}
+
+	xyz_sum = (x + y + z);
+
+	xyz[0] = x / xyz_sum;
+	xyz[1] = y / xyz_sum;
+	xyz[2] = z / xyz_sum;
+}
+
+void blackbody_temperature_to_rgb_table(float *r_table, int width, float min, float max)
+{
+	int i, j = 0, dj = 1;
+	float rgb[3], xyz[3], lms[3], lms_w[3];
+	float bb_temp;
+
+	if (min < max) {
+		SWAP(float, min, max);
+		j = width - 1;
+		dj = -1;
+	}
+
+	for (i = 0; i < width; i++, j += dj) {
+		bb_temp = min + (max - min) / (float)width * (float)i;
+
+		/* integrate blackbody radiation spectrum to XYZ */
+		spectrum_to_xyz(bb_temp, xyz);
+
+		/* normalize highest temperature to white (in LMS system) */
+		xyz_to_lms(xyz[0], xyz[1], xyz[2], &lms[0], &lms[1], &lms[2]);
+
+		if (i == 0) {
+			lms_w[0] = 1.0f / lms[0];
+			lms_w[1] = 1.0f / lms[1];
+			lms_w[2] = 1.0f / lms[2];
+		}
+
+		mul_v3_v3(lms, lms_w);
+
+		lms_to_xyz(lms[0], lms[1], lms[2], &xyz[0], &xyz[1], &xyz[2]);
+
+		/* convert to RGB */
+		xyz_to_rgb(xyz[0], xyz[1], xyz[2], &rgb[0], &rgb[1], &rgb[2], BLI_XYZ_CIE);
+		constrain_rgb(&rgb[0], &rgb[1], &rgb[2]);
+		normalize_rgb(rgb);
+
+		copy_v3_v3(&r_table[(j << 2)], rgb);
+
+		if (rgb[2] > 0.1f)
+			r_table[(j << 2) + 3] = rgb[2];
+		else
+			r_table[(j << 2) + 3] = 0.0f;
+	}
+}