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

4 files changed, 148 insertions, 20 deletions

diff --git a/source/blender/blenkernel/intern/smoke.c b/source/blender/blenkernel/intern/smoke.c
index c3e9d7d793c..62874939a6d 100644
--- a/source/blender/blenkernel/intern/smoke.c
+++ b/source/blender/blenkernel/intern/smoke.c
@@ -128,7 +128,6 @@ void smoke_initBlenderRNA(struct FLUID_3D *UNUSED(fluid), float *UNUSED(alpha),
                           float *UNUSED(flame_vorticity), float *UNUSED(flame_ignition_temp), float *UNUSED(flame_max_temp)) {}
 struct DerivedMesh *smokeModifier_do(SmokeModifierData *UNUSED(smd), Scene *UNUSED(scene), Object *UNUSED(ob), DerivedMesh *UNUSED(dm), bool UNUSED(for_render)) { return NULL; }
 float smoke_get_velocity_at(struct Object *UNUSED(ob), float UNUSED(position[3]), float UNUSED(velocity[3])) { return 0.0f; }
-void flame_get_spectrum(unsigned char *UNUSED(spec), int UNUSED(width), float UNUSED(t1), float UNUSED(t2)) {}
 
 #endif /* WITH_SMOKE */
 
diff --git a/source/blender/blenlib/BLI_math_color.h b/source/blender/blenlib/BLI_math_color.h
index d70dfcd9e58..5e6b1256d30 100644
--- a/source/blender/blenlib/BLI_math_color.h
+++ b/source/blender/blenlib/BLI_math_color.h
@@ -145,6 +145,8 @@ MINLINE void rgba_char_args_set(char col[4], const char r, const char g, const c
 MINLINE void rgba_char_args_test_set(char col[4], const char r, const char g, const char b, const char a);
 MINLINE void cpack_cpy_3ub(unsigned char r_col[3], const unsigned int pack);
 
+void blackbody_temperature_to_rgb_table(float *r_table, int width, float min, float max);
+
 /********* lift/gamma/gain / ASC-CDL conversion ***********/
 
 void lift_gamma_gain_to_asc_cdl(float *lift, float *gamma, float *gain, float *offset, float *slope, float *power);
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;
+	}
+}
diff --git a/source/blender/editors/space_view3d/drawvolume.c b/source/blender/editors/space_view3d/drawvolume.c
index 102bfc55a2c..72d84b4be04 100644
--- a/source/blender/editors/space_view3d/drawvolume.c
+++ b/source/blender/editors/space_view3d/drawvolume.c
@@ -98,19 +98,19 @@ static GPUTexture *create_flame_spectrum_texture(void)
 
 	GPUTexture *tex;
 	int i, j, k;
-	unsigned char *spec_data = malloc(SPEC_WIDTH * 4 * sizeof(unsigned char));
-	float *spec_pixels = malloc(SPEC_WIDTH * 4 * 16 * 16 * sizeof(float));
+	float *spec_data = MEM_mallocN(SPEC_WIDTH * 4 * sizeof(float), "spec_data");
+	float *spec_pixels = MEM_mallocN(SPEC_WIDTH * 4 * 16 * 16 * sizeof(float), "spec_pixels");
 
-	flame_get_spectrum(spec_data, SPEC_WIDTH, 1500, 3000);
+	blackbody_temperature_to_rgb_table(spec_data, SPEC_WIDTH, 1500, 3000);
 
 	for (i = 0; i < 16; i++) {
 		for (j = 0; j < 16; j++) {
 			for (k = 0; k < SPEC_WIDTH; k++) {
 				int index = (j * SPEC_WIDTH * 16 + i * SPEC_WIDTH + k) * 4;
 				if (k >= FIRE_THRESH) {
-					spec_pixels[index] = ((float)spec_data[k * 4]) / 255.0f;
-					spec_pixels[index + 1] = ((float)spec_data[k * 4 + 1]) / 255.0f;
-					spec_pixels[index + 2] = ((float)spec_data[k * 4 + 2]) / 255.0f;
+					spec_pixels[index] = (spec_data[k * 4]);
+					spec_pixels[index + 1] = (spec_data[k * 4 + 1]);
+					spec_pixels[index + 2] = (spec_data[k * 4 + 2]);
 					spec_pixels[index + 3] = MAX_FIRE_ALPHA * (
 					        (k > FULL_ON_FIRE) ? 1.0f : (k - FIRE_THRESH) / ((float)FULL_ON_FIRE - FIRE_THRESH));
 				}
@@ -123,8 +123,8 @@ static GPUTexture *create_flame_spectrum_texture(void)
 
 	tex = GPU_texture_create_1D(SPEC_WIDTH, spec_pixels, NULL);
 
-	free(spec_data);
-	free(spec_pixels);
+	MEM_freeN(spec_data);
+	MEM_freeN(spec_pixels);
 
 #undef SPEC_WIDTH
 #undef FIRE_THRESH