Smoke (fire): Move spectrum code from C++ (intern/) to C code (BLI)
authorKévin Dietrich <kevin.dietrich@mailoo.org>
Sat, 9 Jan 2016 03:48:31 +0000 (04:48 +0100)
committerKévin Dietrich <kevin.dietrich@mailoo.org>
Sat, 9 Jan 2016 03:48:31 +0000 (04:48 +0100)
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

intern/smoke/CMakeLists.txt
intern/smoke/extern/smoke_API.h
intern/smoke/intern/smoke_API.cpp
intern/smoke/intern/spectrum.cpp [deleted file]
intern/smoke/intern/spectrum.h [deleted file]
source/blender/blenkernel/intern/smoke.c
source/blender/blenlib/BLI_math_color.h
source/blender/blenlib/intern/math_color.c
source/blender/editors/space_view3d/drawvolume.c

index 8cda0fd..57678ec 100644 (file)
@@ -40,7 +40,6 @@ set(SRC
        intern/FLUID_3D_SOLVERS.cpp
        intern/FLUID_3D_STATIC.cpp
        intern/LU_HELPER.cpp
-       intern/spectrum.cpp
        intern/SPHERE.cpp
        intern/WTURBULENCE.cpp
        intern/smoke_API.cpp
@@ -54,7 +53,6 @@ set(SRC
        intern/LU_HELPER.h
        intern/MERSENNETWISTER.h
        intern/OBSTACLE.h
-       intern/spectrum.h
        intern/SPHERE.h
        intern/VEC3.h
        intern/WAVELET_NOISE.h
index 08dbded..a68a587 100644 (file)
@@ -99,9 +99,6 @@ void smoke_export(struct FLUID_3D *fluid, float *dt, float *dx, float **dens, fl
 void smoke_turbulence_export(struct WTURBULENCE *wt, float **dens, float **react, float **flame, float **fuel,
                                                         float **r, float **g, float **b, float **tcu, float **tcv, float **tcw);
 
-/* flame spectrum */
-void flame_get_spectrum(unsigned char *spec, int width, float t1, float t2);
-
 /* data fields */
 int smoke_has_heat(struct FLUID_3D *fluid);
 int smoke_has_fuel(struct FLUID_3D *fluid);
index d79aaf7..5cab0b3 100644 (file)
@@ -30,7 +30,6 @@
 
 #include "FLUID_3D.h"
 #include "WTURBULENCE.h"
-#include "spectrum.h"
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -449,11 +448,6 @@ extern "C" void smoke_turbulence_set_noise(WTURBULENCE *wt, int type, const char
        wt->setNoise(type, noisefile_path);
 }
 
-extern "C" void flame_get_spectrum(unsigned char *spec, int width, float t1, float t2)
-{
-       spectrum(t1, t2, width, spec);
-}
-
 extern "C" int smoke_has_heat(FLUID_3D *fluid)
 {
        return (fluid->_heat) ? 1 : 0;
diff --git a/intern/smoke/intern/spectrum.cpp b/intern/smoke/intern/spectrum.cpp
deleted file mode 100644 (file)
index 15e4d0c..0000000
+++ /dev/null
@@ -1,428 +0,0 @@
-/*
-                Colour Rendering of Spectra
-
-                       by John Walker
-                  http://www.fourmilab.ch/
-                 
-                Last updated: March 9, 2003
-
-           This program is in the public domain.
-
-    For complete information about the techniques employed in
-    this program, see the World-Wide Web document:
-
-             http://www.fourmilab.ch/documents/specrend/
-            
-    The xyz_to_rgb() function, which was wrong in the original
-    version of this program, was corrected by:
-    
-           Andrew J. S. Hamilton 21 May 1999
-           Andrew.Hamilton@Colorado.EDU
-           http://casa.colorado.edu/~ajsh/
-
-    who also added the gamma correction facilities and
-    modified constrain_rgb() to work by desaturating the
-    colour by adding white.
-    
-    A program which uses these functions to plot CIE
-    "tongue" diagrams called "ppmcie" is included in
-    the Netpbm graphics toolkit:
-       http://netpbm.sourceforge.net/
-    (The program was called cietoppm in earlier
-    versions of Netpbm.)
-
-*/
-
-#include <stdio.h>
-#include <math.h>
-#include "spectrum.h"
-
-/* A colour system is defined by the CIE x and y coordinates of
-   its three primary illuminants and the x and y coordinates of
-   the white point. */
-
-struct colourSystem {
-    const char *name;              /* Colour system name */
-    double xRed, yRed,             /* Red x, y */
-           xGreen, yGreen,         /* Green x, y */
-           xBlue, yBlue,           /* Blue x, y */
-           xWhite, yWhite,         /* White point x, y */
-          gamma;                   /* Gamma correction for system */
-};
-
-/* White point chromaticities. */
-
-#if 0
-#define IlluminantC     0.3101, 0.3162         /* For NTSC television */
-#define IlluminantD65   0.3127, 0.3291         /* For EBU and SMPTE */
-#endif
-#define IlluminantE    0.33333333, 0.33333333  /* CIE equal-energy illuminant */
-
-/*  Gamma of nonlinear correction.
-
-    See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at:
-    
-       http://www.poynton.com/ColorFAQ.html
-       http://www.poynton.com/GammaFAQ.html
-*/
-
-#define GAMMA_REC709   0               /* Rec. 709 */
-
-static struct colourSystem
-                  /* Name                  xRed    yRed    xGreen  yGreen  xBlue  yBlue    White point        Gamma   */
-#if 0 /* UNUSED */
-    NTSCsystem  =  { "NTSC",               0.67,   0.33,   0.21,   0.71,   0.14,   0.08,   IlluminantC,    GAMMA_REC709 },
-    EBUsystem   =  { "EBU (PAL/SECAM)",    0.64,   0.33,   0.29,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 },
-    SMPTEsystem =  { "SMPTE",              0.630,  0.340,  0.310,  0.595,  0.155,  0.070,  IlluminantD65,  GAMMA_REC709 },
-    HDTVsystem  =  { "HDTV",               0.670,  0.330,  0.210,  0.710,  0.150,  0.060,  IlluminantD65,  GAMMA_REC709 },
-#endif
-
-    CIEsystem   =  { "CIE",                0.7355, 0.2645, 0.2658, 0.7243, 0.1669, 0.0085, IlluminantE,    GAMMA_REC709 };
-
-#if 0 /* UNUSED */
-    Rec709system = { "CIE REC 709",        0.64,   0.33,   0.30,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 };
-#endif
-
-/*                         UPVP_TO_XY
-
-    Given 1976 coordinates u', v', determine 1931 chromaticities x, y
-    
-*/
-
-#if 0 /* UNUSED */
-static void upvp_to_xy(double up, double vp, double *xc, double *yc)
-{
-    *xc = (9 * up) / ((6 * up) - (16 * vp) + 12);
-    *yc = (4 * vp) / ((6 * up) - (16 * vp) + 12);
-}
-#endif
-
-/*                         XY_TO_UPVP
-
-    Given 1931 chromaticities x, y, determine 1976 coordinates u', v'
-    
-*/
-
-#if 0 /* UNUSED */
-static void xy_to_upvp(double xc, double yc, double *up, double *vp)
-{
-    *up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3);
-    *vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3);
-}
-#endif
-
-/*                             XYZ_TO_RGB
-
-    Given an additive tricolour system CS, defined by the CIE x
-    and y chromaticities of its three primaries (z is derived
-    trivially as 1-(x+y)), and a desired chromaticity (XC, YC,
-    ZC) in CIE space, determine the contribution of each
-    primary in a linear combination which sums to the desired
-    chromaticity.  If the  requested chromaticity falls outside
-    the Maxwell  triangle (colour gamut) formed by the three
-    primaries, one of the r, g, or b weights will be negative. 
-
-    Caller can use constrain_rgb() to desaturate an
-    outside-gamut colour to the closest representation within
-    the available gamut and/or norm_rgb to normalise the RGB
-    components so the largest nonzero component has value 1.
-    
-*/
-
-static void xyz_to_rgb(struct colourSystem *cs,
-                       double xc, double yc, double zc,
-                       double *r, double *g, double *b)
-{
-    double xr, yr, zr, xg, yg, zg, xb, yb, zb;
-    double xw, yw, zw;
-    double rx, ry, rz, gx, gy, gz, bx, by, bz;
-    double rw, gw, bw;
-
-    xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
-    xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
-    xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);
-
-    xw = cs->xWhite;  yw = cs->yWhite;  zw = 1 - (xw + yw);
-
-    /* xyz -> rgb matrix, before scaling to white. */
-    
-    rx = (yg * zb) - (yb * zg);  ry = (xb * zg) - (xg * zb);  rz = (xg * yb) - (xb * yg);
-    gx = (yb * zr) - (yr * zb);  gy = (xr * zb) - (xb * zr);  gz = (xb * yr) - (xr * yb);
-    bx = (yr * zg) - (yg * zr);  by = (xg * zr) - (xr * zg);  bz = (xr * yg) - (xg * yr);
-
-    /* White scaling factors.
-       Dividing by yw scales the white luminance to unity, as conventional. */
-       
-    rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw;
-    gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw;
-    bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw;
-
-    /* xyz -> rgb matrix, correctly scaled to white. */
-    
-    rx = rx / rw;  ry = ry / rw;  rz = rz / rw;
-    gx = gx / gw;  gy = gy / gw;  gz = gz / gw;
-    bx = bx / bw;  by = by / bw;  bz = bz / bw;
-
-    /* rgb of the desired point */
-
-    *r = (rx * xc) + (ry * yc) + (rz * zc);
-    *g = (gx * xc) + (gy * yc) + (gz * zc);
-    *b = (bx * xc) + (by * yc) + (bz * zc);
-}
-
-/*                            INSIDE_GAMUT
-
-     Test whether a requested colour is within the gamut
-     achievable with the primaries of the current colour
-     system.  This amounts simply to testing whether all the
-     primary weights are non-negative. */
-
-#if 0 /* UNUSED */
-static int inside_gamut(double r, double g, double b)
-{
-    return (r >= 0) && (g >= 0) && (b >= 0);
-}
-#endif
-
-/*                          CONSTRAIN_RGB
-
-    If the requested RGB shade contains a negative weight for
-    one of the primaries, it lies outside the colour gamut 
-    accessible from the given triple of primaries.  Desaturate
-    it by adding white, equal quantities of R, G, and B, enough
-    to make RGB all positive.  The function returns 1 if the
-    components were modified, zero otherwise.
-    
-*/
-
-static int constrain_rgb(double *r, double *g, double *b)
-{
-    double 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;
-
-    /* Add just enough white to make r, g, b all positive. */
-    
-    if (w > 0) {
-        *r += w;  *g += w; *b += w;
-        return 1;                     /* Colour modified to fit RGB gamut */
-    }
-
-    return 0;                         /* Colour within RGB gamut */
-}
-
-/*                          GAMMA_CORRECT_RGB
-
-    Transform linear RGB values to nonlinear RGB values. Rec.
-    709 is ITU-R Recommendation BT. 709 (1990) ``Basic
-    Parameter Values for the HDTV Standard for the Studio and
-    for International Programme Exchange'', formerly CCIR Rec.
-    709. For details see
-    
-       http://www.poynton.com/ColorFAQ.html
-       http://www.poynton.com/GammaFAQ.html
-*/
-
-#if 0 /* UNUSED */
-static void gamma_correct(const struct colourSystem *cs, double *c)
-{
-    double gamma;
-
-    gamma = cs->gamma;
-
-    if (gamma == GAMMA_REC709) {
-       /* Rec. 709 gamma correction. */
-       double cc = 0.018;
-       
-       if (*c < cc) {
-           *c *= ((1.099 * pow(cc, 0.45)) - 0.099) / cc;
-       } else {
-           *c = (1.099 * pow(*c, 0.45)) - 0.099;
-       }
-    } else {
-       /* Nonlinear colour = (Linear colour)^(1/gamma) */
-       *c = pow(*c, 1.0 / gamma);
-    }
-}
-
-static void gamma_correct_rgb(const struct colourSystem *cs, double *r, double *g, double *b)
-{
-    gamma_correct(cs, r);
-    gamma_correct(cs, g);
-    gamma_correct(cs, b);
-}
-#endif
-
-/*                         NORM_RGB
-
-    Normalise RGB components so the most intense (unless all
-    are zero) has a value of 1.
-    
-*/
-
-static void norm_rgb(double *r, double *g, double *b)
-{
-#define Max(a, b)   (((a) > (b)) ? (a) : (b))
-    double greatest = Max(*r, Max(*g, *b));
-
-    if (greatest > 0) {
-       *r /= greatest;
-       *g /= greatest;
-       *b /= greatest;
-    }
-#undef Max
-}
-
-/*                          SPECTRUM_TO_XYZ
-
-    Calculate the CIE X, Y, and Z coordinates corresponding to
-    a light source with spectral distribution given by  the
-    function SPEC_INTENS, which is called with a series of
-    wavelengths between 380 and 780 nm (the argument is 
-    expressed in meters), which returns emittance at  that
-    wavelength in arbitrary units.  The chromaticity
-    coordinates of the spectrum are returned in the x, y, and z
-    arguments which respect the identity:
-
-            x + y + z = 1.
-*/
-
-static void spectrum_to_xyz(double (*spec_intens)(double wavelength),
-                            double *x, double *y, double *z)
-{
-    int i;
-    double lambda, X = 0, Y = 0, Z = 0, XYZ;
-
-    /* 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
-
-       To save memory, this table can be declared as floats
-       rather than doubles; (IEEE) float has enough 
-       significant bits to represent the values. It's declared
-       as a double here to avoid warnings about "conversion
-       between floating-point types" from certain persnickety
-       compilers. */
-
-    static double cie_colour_match[81][3] = {
-        {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
-        {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
-        {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
-        {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
-        {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
-        {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
-        {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
-        {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
-        {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
-        {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
-        {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
-        {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
-        {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
-        {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
-        {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
-        {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
-        {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
-        {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
-        {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
-        {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
-        {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
-        {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
-        {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
-        {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
-        {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
-        {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
-        {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
-    };
-
-    for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
-        double Me;
-
-        Me = (*spec_intens)(lambda);
-        X += Me * cie_colour_match[i][0];
-        Y += Me * cie_colour_match[i][1];
-        Z += Me * cie_colour_match[i][2];
-    }
-    XYZ = (X + Y + Z);
-    *x = X / XYZ;
-    *y = Y / XYZ;
-    *z = Z / XYZ;
-}
-
-/*                            BB_SPECTRUM
-
-    Calculate, by Planck's radiation law, the emittance of a black body
-    of temperature bbTemp at the given wavelength (in metres).  */
-
-static double bbTemp = 5000;          /* Hidden temperature argument
-                                         to BB_SPECTRUM. */
-static double bb_spectrum(double wavelength)
-{
-    double wlm = wavelength * 1e-9;   /* Wavelength in meters */
-
-    return (3.74183e-16 * pow(wlm, -5.0)) /
-           (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0);
-}
-
-static void xyz_to_lms(double x, double y, double z, double* l, double* m, double* s)
-{
-       *l =  0.3897*x + 0.6890*y - 0.0787*z;
-       *m = -0.2298*x + 1.1834*y + 0.0464*z;
-       *s = z;
-}
-
-static void lms_to_xyz(double l, double m, double s, double* x, double *y, double* z)
-{
-       *x =  1.9102*l - 1.1121*m + 0.2019*s;
-       *y =  0.3709*l + 0.6290*m + 0.0000*s;
-       *z = s;
-}
-
-void spectrum(double t1, double t2, int N, unsigned char *d)
-{
-       int i,j,dj;
-       double X,Y,Z,R,G,B,L,M,S, Lw, Mw, Sw;
-       struct colourSystem *cs = &CIEsystem;
-
-       j = 0; dj = 1;
-       if (t1<t2) {
-               double t = t1;
-               t1 = t2;
-               t2 = t;
-               j = N-1; dj=-1;
-       }
-
-       for (i=0; i<N; i++) {
-               bbTemp = t1 + (t2-t1)/N*i;
-
-                       // integrate blackbody radiation spectrum to XYZ
-               spectrum_to_xyz(bb_spectrum, &X, &Y, &Z);
-
-                       // normalize highest temperature to white (in LMS system)
-               xyz_to_lms(X,Y,Z,&L,&M,&S);
-               if (i==0) {
-                       Lw=1/L; Mw=1/M; Sw=1/S;
-               }
-               L *= Lw; M *= Mw; S *= Sw;
-               lms_to_xyz(L,M,S,&X,&Y,&Z);
-
-                       // convert to RGB
-               xyz_to_rgb(cs, X, Y, Z, &R, &G, &B);
-               constrain_rgb(&R, &G, &B);
-               norm_rgb(&R, &G, &B);
-               d[(j<<2)] = (unsigned char) ((double)R*255);
-               d[(j<<2)+1] = (unsigned char) ((double)G*255);
-               d[(j<<2)+2] = (unsigned char) ((double)B*255);
-               d[(j<<2)+3] = (B>0.1)? B*255 : 0;
-               j += dj;
-       }
-}
diff --git a/intern/smoke/intern/spectrum.h b/intern/smoke/intern/spectrum.h
deleted file mode 100644 (file)
index 3ffd41f..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef __SPECTRUM_H
-#define __SPECTRUM_H
-
-void spectrum(double t1, double t2, int n, unsigned char *d);
-
-#endif
index c3e9d7d..6287493 100644 (file)
@@ -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 */
 
index d70dfcd..5e6b125 100644 (file)
@@ -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);
index e149fdb..6d7d24c 100644 (file)
@@ -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;
+       }
+}
index 102bfc5..72d84b4 100644 (file)
@@ -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