Convert light frequency to RGB?

Question

Does anyone know of any formula for converting a light frequency to an RGB value?

Solution 1

Here's a detailed explanation of the entire conversion process: http://www.fourmilab.ch/documents/specrend/. Source code included!

Solution 2

For lazy guys (like me), here is an implementation in java of the code found in @user151323 's answer (that is, just a simple translation from pascal code found in Spectra Lab Report):

static private final double Gamma = 0.80;
static private final double IntensityMax = 255;

/**
 * Taken from Earl F. Glynn's web page:
 * <a href="http://www.efg2.com/Lab/ScienceAndEngineering/Spectra.htm">Spectra Lab Report</a>
 */
public static int[] waveLengthToRGB(double Wavelength) {
    double factor;
    double Red, Green, Blue;

    if((Wavelength >= 380) && (Wavelength < 440)) {
        Red = -(Wavelength - 440) / (440 - 380);
        Green = 0.0;
        Blue = 1.0;
    } else if((Wavelength >= 440) && (Wavelength < 490)) {
        Red = 0.0;
        Green = (Wavelength - 440) / (490 - 440);
        Blue = 1.0;
    } else if((Wavelength >= 490) && (Wavelength < 510)) {
        Red = 0.0;
        Green = 1.0;
        Blue = -(Wavelength - 510) / (510 - 490);
    } else if((Wavelength >= 510) && (Wavelength < 580)) {
        Red = (Wavelength - 510) / (580 - 510);
        Green = 1.0;
        Blue = 0.0;
    } else if((Wavelength >= 580) && (Wavelength < 645)) {
        Red = 1.0;
        Green = -(Wavelength - 645) / (645 - 580);
        Blue = 0.0;
    } else if((Wavelength >= 645) && (Wavelength < 781)) {
        Red = 1.0;
        Green = 0.0;
        Blue = 0.0;
    } else {
        Red = 0.0;
        Green = 0.0;
        Blue = 0.0;
    }

    // Let the intensity fall off near the vision limits

    if((Wavelength >= 380) && (Wavelength < 420)) {
        factor = 0.3 + 0.7 * (Wavelength - 380) / (420 - 380);
    } else if((Wavelength >= 420) && (Wavelength < 701)) {
        factor = 1.0;
    } else if((Wavelength >= 701) && (Wavelength < 781)) {
        factor = 0.3 + 0.7 * (780 - Wavelength) / (780 - 700);
    } else {
        factor = 0.0;
    }


    int[] rgb = new int[3];

    // Don't want 0^x = 1 for x <> 0
    rgb[0] = Red == 0.0 ? 0 : (int)Math.round(IntensityMax * Math.pow(Red * factor, Gamma));
    rgb[1] = Green == 0.0 ? 0 : (int)Math.round(IntensityMax * Math.pow(Green * factor, Gamma));
    rgb[2] = Blue == 0.0 ? 0 : (int)Math.round(IntensityMax * Math.pow(Blue * factor, Gamma));

    return rgb;
}

Solution 3

General idea:

1. Use CEI color matching functions to convert wavelength to XYZ color.
2. Convert XYZ to RGB
3. Clip components to [0..1] and multiply by 255 to fit in the unsigned byte range.

Steps 1 and 2 may vary.

There are several color matching functions, available as tables or as analytic approximations (suggested by @Tarc and @Haochen Xie). Tables are best if you need a smooth preсise result.

There is no single RGB color space. Multiple transformation matrices and different kinds of gamma correction may be used.

Below is the C# code I came up with recently. It uses linear interpolation over the "CIE 1964 standard observer" table and sRGB matrix + gamma correction.

static class RgbCalculator {

    const int
         LEN_MIN = 380,
         LEN_MAX = 780,
         LEN_STEP = 5;

    static readonly double[]
        X = {
                0.000160, 0.000662, 0.002362, 0.007242, 0.019110, 0.043400, 0.084736, 0.140638, 0.204492, 0.264737,
                0.314679, 0.357719, 0.383734, 0.386726, 0.370702, 0.342957, 0.302273, 0.254085, 0.195618, 0.132349,
                0.080507, 0.041072, 0.016172, 0.005132, 0.003816, 0.015444, 0.037465, 0.071358, 0.117749, 0.172953,
                0.236491, 0.304213, 0.376772, 0.451584, 0.529826, 0.616053, 0.705224, 0.793832, 0.878655, 0.951162,
                1.014160, 1.074300, 1.118520, 1.134300, 1.123990, 1.089100, 1.030480, 0.950740, 0.856297, 0.754930,
                0.647467, 0.535110, 0.431567, 0.343690, 0.268329, 0.204300, 0.152568, 0.112210, 0.081261, 0.057930,
                0.040851, 0.028623, 0.019941, 0.013842, 0.009577, 0.006605, 0.004553, 0.003145, 0.002175, 0.001506,
                0.001045, 0.000727, 0.000508, 0.000356, 0.000251, 0.000178, 0.000126, 0.000090, 0.000065, 0.000046,
                0.000033
            },

        Y = {
                0.000017, 0.000072, 0.000253, 0.000769, 0.002004, 0.004509, 0.008756, 0.014456, 0.021391, 0.029497,
                0.038676, 0.049602, 0.062077, 0.074704, 0.089456, 0.106256, 0.128201, 0.152761, 0.185190, 0.219940,
                0.253589, 0.297665, 0.339133, 0.395379, 0.460777, 0.531360, 0.606741, 0.685660, 0.761757, 0.823330,
                0.875211, 0.923810, 0.961988, 0.982200, 0.991761, 0.999110, 0.997340, 0.982380, 0.955552, 0.915175,
                0.868934, 0.825623, 0.777405, 0.720353, 0.658341, 0.593878, 0.527963, 0.461834, 0.398057, 0.339554,
                0.283493, 0.228254, 0.179828, 0.140211, 0.107633, 0.081187, 0.060281, 0.044096, 0.031800, 0.022602,
                0.015905, 0.011130, 0.007749, 0.005375, 0.003718, 0.002565, 0.001768, 0.001222, 0.000846, 0.000586,
                0.000407, 0.000284, 0.000199, 0.000140, 0.000098, 0.000070, 0.000050, 0.000036, 0.000025, 0.000018,
                0.000013
            },

        Z = {
                0.000705, 0.002928, 0.010482, 0.032344, 0.086011, 0.197120, 0.389366, 0.656760, 0.972542, 1.282500,
                1.553480, 1.798500, 1.967280, 2.027300, 1.994800, 1.900700, 1.745370, 1.554900, 1.317560, 1.030200,
                0.772125, 0.570060, 0.415254, 0.302356, 0.218502, 0.159249, 0.112044, 0.082248, 0.060709, 0.043050,
                0.030451, 0.020584, 0.013676, 0.007918, 0.003988, 0.001091, 0.000000, 0.000000, 0.000000, 0.000000,
                0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
                0.000000
            };

    static readonly double[]
        MATRIX_SRGB_D65 = {
             3.2404542, -1.5371385, -0.4985314,
            -0.9692660,  1.8760108,  0.0415560,
             0.0556434, -0.2040259,  1.0572252
        };

    public static byte[] Calc(double len) {
        if(len < LEN_MIN || len > LEN_MAX)
            return new byte[3];

        len -= LEN_MIN;
        var index = (int)Math.Floor(len / LEN_STEP);
        var offset = len - LEN_STEP * index;

        var x = Interpolate(X, index, offset);
        var y = Interpolate(Y, index, offset);
        var z = Interpolate(Z, index, offset);

        var m = MATRIX_SRGB_D65;

        var r = m[0] * x + m[1] * y + m[2] * z;
        var g = m[3] * x + m[4] * y + m[5] * z;
        var b = m[6] * x + m[7] * y + m[8] * z;

        r = Clip(GammaCorrect_sRGB(r));
        g = Clip(GammaCorrect_sRGB(g));
        b = Clip(GammaCorrect_sRGB(b));

        return new[] { 
            (byte)(255 * r),
            (byte)(255 * g),
            (byte)(255 * b)
        };
    }

    static double Interpolate(double[] values, int index, double offset) {
        if(offset == 0)
            return values[index];

        var x0 = index * LEN_STEP;
        var x1 = x0 + LEN_STEP;
        var y0 = values[index];
        var y1 = values[1 + index];

        return y0 + offset * (y1 - y0) / (x1 - x0);
    }

    static double GammaCorrect_sRGB(double c) {
        if(c <= 0.0031308)
            return 12.92 * c;

        var a = 0.055;
        return (1 + a) * Math.Pow(c, 1 / 2.4) - a;
    }

    static double Clip(double c) {
        if(c < 0)
            return 0;
        if(c > 1)
            return 1;
        return c;
    }
}

Result for the 400-700 nm range:

Solution 4

Although this is an old question and already gets a handful good answers, when I tried to implement such conversion functionality in my application I was not satisfied with the algorithms already listed here and did my own research, which gave me some good result. So I'm going to post a new answer.

After some researchs I came across this paper, Simple Analytic Approximations to the CIE XYZ Color Matching Functions, and tried to adopt the introduced multi-lobe piecewise Gaussian fit algorithm in my application. The paper only described the functions to convert a wavelength to the corresponding XYZ values, so I implemented a function to convert XYZ to RGB in the sRGB color space and combined them. The result is fantastic and worth sharing:

/**
 * Convert a wavelength in the visible light spectrum to a RGB color value that is suitable to be displayed on a
 * monitor
 *
 * @param wavelength wavelength in nm
 * @return RGB color encoded in int. each color is represented with 8 bits and has a layout of
 * 00000000RRRRRRRRGGGGGGGGBBBBBBBB where MSB is at the leftmost
 */
public static int wavelengthToRGB(double wavelength){
    double[] xyz = cie1931WavelengthToXYZFit(wavelength);
    double[] rgb = srgbXYZ2RGB(xyz);

    int c = 0;
    c |= (((int) (rgb[0] * 0xFF)) & 0xFF) << 16;
    c |= (((int) (rgb[1] * 0xFF)) & 0xFF) << 8;
    c |= (((int) (rgb[2] * 0xFF)) & 0xFF) << 0;

    return c;
}

/**
 * Convert XYZ to RGB in the sRGB color space
 * <p>
 * The conversion matrix and color component transfer function is taken from http://www.color.org/srgb.pdf, which
 * follows the International Electrotechnical Commission standard IEC 61966-2-1 "Multimedia systems and equipment -
 * Colour measurement and management - Part 2-1: Colour management - Default RGB colour space - sRGB"
 *
 * @param xyz XYZ values in a double array in the order of X, Y, Z. each value in the range of [0.0, 1.0]
 * @return RGB values in a double array, in the order of R, G, B. each value in the range of [0.0, 1.0]
 */
public static double[] srgbXYZ2RGB(double[] xyz) {
    double x = xyz[0];
    double y = xyz[1];
    double z = xyz[2];

    double rl =  3.2406255 * x + -1.537208  * y + -0.4986286 * z;
    double gl = -0.9689307 * x +  1.8757561 * y +  0.0415175 * z;
    double bl =  0.0557101 * x + -0.2040211 * y +  1.0569959 * z;

    return new double[] {
            srgbXYZ2RGBPostprocess(rl),
            srgbXYZ2RGBPostprocess(gl),
            srgbXYZ2RGBPostprocess(bl)
    };
}

/**
 * helper function for {@link #srgbXYZ2RGB(double[])}
 */
private static double srgbXYZ2RGBPostprocess(double c) {
    // clip if c is out of range
    c = c > 1 ? 1 : (c < 0 ? 0 : c);

    // apply the color component transfer function
    c = c <= 0.0031308 ? c * 12.92 : 1.055 * Math.pow(c, 1. / 2.4) - 0.055;

    return c;
}

/**
 * A multi-lobe, piecewise Gaussian fit of CIE 1931 XYZ Color Matching Functions by Wyman el al. from Nvidia. The
 * code here is adopted from the Listing 1 of the paper authored by Wyman et al.
 * <p>
 * Reference: Chris Wyman, Peter-Pike Sloan, and Peter Shirley, Simple Analytic Approximations to the CIE XYZ Color
 * Matching Functions, Journal of Computer Graphics Techniques (JCGT), vol. 2, no. 2, 1-11, 2013.
 *
 * @param wavelength wavelength in nm
 * @return XYZ in a double array in the order of X, Y, Z. each value in the range of [0.0, 1.0]
 */
public static double[] cie1931WavelengthToXYZFit(double wavelength) {
    double wave = wavelength;

    double x;
    {
        double t1 = (wave - 442.0) * ((wave < 442.0) ? 0.0624 : 0.0374);
        double t2 = (wave - 599.8) * ((wave < 599.8) ? 0.0264 : 0.0323);
        double t3 = (wave - 501.1) * ((wave < 501.1) ? 0.0490 : 0.0382);

        x =   0.362 * Math.exp(-0.5 * t1 * t1)
            + 1.056 * Math.exp(-0.5 * t2 * t2)
            - 0.065 * Math.exp(-0.5 * t3 * t3);
    }

    double y;
    {
        double t1 = (wave - 568.8) * ((wave < 568.8) ? 0.0213 : 0.0247);
        double t2 = (wave - 530.9) * ((wave < 530.9) ? 0.0613 : 0.0322);

        y =   0.821 * Math.exp(-0.5 * t1 * t1)
            + 0.286 * Math.exp(-0.5 * t2 * t2);
    }

    double z;
    {
        double t1 = (wave - 437.0) * ((wave < 437.0) ? 0.0845 : 0.0278);
        double t2 = (wave - 459.0) * ((wave < 459.0) ? 0.0385 : 0.0725);

        z =   1.217 * Math.exp(-0.5 * t1 * t1)
            + 0.681 * Math.exp(-0.5 * t2 * t2);
    }

    return new double[] { x, y, z };
}

my code is written in Java 8, but it shouldn't be hard to port it to lower versions of Java and other languages.

Solution 5

I guess I might as well follow up my comment with a formal answer. The best option is to use the HSV colour space - though the hue represents the wavelength it is not a one-to-one comparison.

View more solutions

Comments

[Francis Davey]

And the Fourmilab article makes the important point that some colours are not representable in RGB (bright oranges being a good example) because you cannot "make" arbitrary colours of light by adding three primary colours together, whatever our physics teachers may have told us (well mine did). Too bad, but in practice not usually fatal.

[Hassedev]

There seems to be a bug in your code. If the wavelength is for example 439.5, your function returns black. The original code on the site was working with integers, I believe (I don't know pascal at all). I suggest to change Wavelength<=439 to Wavelength<440.

[Tarc]

You're right! Thank you for pointing this out to me :) Already corrected.

[GrayFace]

In addition to it: en.wikipedia.org/wiki/Srgb The article was written before sRGB standard was widely adopted. Also note the "Calculations assume the 2° standard colorimetric observer" phrase, which means CIE 1931 table found in accompanying source to the paper should be used and not CIE 1964.

[Tomáš Zato]

It would be nice to provide some example how to use the code. It requires function as an argument, uses temperature to calculate colors and such things. One would be happy to know what to delete and change to get it to work.

[Andrew]

This looks fancy, but why did you encode it into a single int? Wouldn't a client just have to pull the values back out again?

[Haochen Xie]

Encoding an (A)RGB color value in a single int is a common convention, and the UI toolkit I was working on uses that convention. You're free to change the interface to suit your application for sure.

[Baddack]

DoubleUnaryOperator transfer = c -> c <= 0.0031308 ? c * 12.92 : 1.055 * Math.pow(c, 1. / 2.4) - 0.055; transfer = transfer.compose(c -> c > 1 ? 1 : (c < 0 ? 0 : c)); This line is throwing me errors. C# doesn't know what a Double Unary Operator is. I don't know either, Ive googled it and seems to be a Java thing. Is there a way to do this line in C#?

[Baddack]

I think I've figured it out after thinking about it more. I think all this function is doing is deciding if the RGB is within range, if it is it returns the color otherwise it returns 0 for that color. That's a fancy way to do it.

[Haochen Xie]

@Baddack, you're right: it's just a fancy way to make some further transformation on the calculated values. I cannot remember exactly, but I think it first applies a gamma correction, then cuts off out of range values. Maybe I should have it done in a separate method, but I wasn't actually thinking about sharing the code while writing it, and it was a toy project in which I needed this conversion.

[Baddack]

@HaochenXie, thanks for sharing, nice code. I was able to get something working. Cheers

[Haochen Xie]

@Baddack I've dug out the project that I needed this conversion, and rewrote this part without using java 8 lambda so the code is more clear. I actually remembered incorrectly about what the transfer DoubleUnaryOperator was doing (thus the explanation in my previous comment are not correct), so please check the new code.

[Baddack]

Thanks for the update, this code is very useful! I appreciate your help.

[Haochen Xie]

@Baddack i'm glad that the code helps you. and if you don't mind, could you please upvote it so it may potentially help more people?

[Baddack]

What does Math.pow(c, 1. / 2.4) do? I can't do this in C#. It's doing c to the power of ? What does the notation 1. / x mean in Math.pow?

[Haochen Xie]

@Baddack Math.pow(c, 1. / 2.4) = c^(1/2.4), i.e. raise c to the power of 1/2.4; 1. is just 1 but the type will be double instead of int

[Haochen Xie]

@Baddack it seems that in C#, Math.pow is just Math.Pow. See msdn.microsoft.com/en-us/library/system.math.pow.aspx

[Baddack]

Thanks, I just didn't know what the 1. was, but I'm on the same page now. Cheers

[Ruslan]

While answer by Tarc looks more or less convincing, yours has some anomaly in the violet zone. Notice the strange bluish line there.

[Ruslan]

Your link is dead.

[Haochen Xie]

@Ruslan since this algorithm is a analytical fit of the CIE standard observer (which could be considered to be the "precise" model), there are errors. But from the paper, if you look at the Figure 1 on page 7 (compare (d) with (f)), this method provides quite a close approximation. Especially if you look at (f), you could see that there is also a bluish line even in the standard model. Also, color perception of pure light source varies personally, so this level of error is probably negligible.

[Ruslan]

Thanks, I see now. One more thing: I do understand why you clamp the resulting values from negative to zero, but why do you clamp from large (i.e. >1) to 1 instead of rescaling the whole visible range by a constant of about 0.4 to avoid saturation (which does change hue, e.g. in a range between green and red)?

[Haochen Xie]

@Ruslan I think you're talking about the function #srgbXYZ2RGBPostprocess(..). Please take a look at the reference I've given in the comments of function #srgbXYZ2RGB(..), I was simply following the algorithm the reference gives.

[Haochen Xie]

@Ruslan you could for sure develop your own XYZ to RGB conversion function as the best conversion matrix (as well as clipping method) varies from application to application

[Ruslan]

Ah, maybe clipping to one is really needed when one clips to zero. I guess it's some sort of compensation. Thanks for your answers.

[phorgan1]

This is really interesting to me. I have an idea to use something like this to give a normal response, but use a WXYZ response to mimic the response of tetrachromats who have a fourth cone which responds to a frequency far enough from any of the other normal three types of cones. That might let me take source images and infer the differences they see. N.B. they don't see new colors, it's that lights that blend, (sum), for example, to a particular yellow seems identical to a yellow of a particular frequency for most of us, but for them, the light wouldn't blend to that yellow at all.

[phorgan1]

Of course, for a particular RGB color, it could have been arrived at in a lot of ways. The green of a leaf could come from the filtering out of everything but green, or the green could have been filtered out, but nano characteristics could cause blue and yellow to reflect and look identical to the green. Given an image rather than the light, is there any way that I can differentiate?

[Violet Giraffe]

It is worth noting that only a small subset of all possible visible wavelengths can be exactly represented in the RGB color space. The conversion process is quite intricate and ambiguous. See physics.stackexchange.com/a/94446/5089 and physics.stackexchange.com/a/419628/5089

[Rodrigo Borba]

Is it expected to have repeated RFB to some frequencies? (RED): 652 - rgb(255, 0, 0) | 660 - rgb(255, 0, 0) | 692 - rgb(255, 0, 0) | 700 - rgb(255, 0, 0) | ...

[fuweichin]

This algorithm has very similar results with @shital-shah 's "Method 1", I trust this one more, but maybe Clip(GammaCorrect_sRGB(r)) should be replaced with GammaCorrect_sRGB(Clip(r)) although it doesn't seem to affect final result.