Converting a spectrum to a colour

(32 comments)

This article presents a Python script to map a spectrum of wavelengths to a representation of a colour. There is no unique way to do this, but the formulation used here is based on the CIE colour matching functions, $\bar{x}(\lambda)$, $\bar{y}(\lambda)$ and $\bar{z}(\lambda)$. These model the chromatic response of a "standard observer" by mapping a power spectrum of wavelengths, $P(\lambda)$, to a set of tristimulus values, $X$, $Y$ and $Z$, analogous to the actual response of the three types of cone cell in the human eye.

\begin{align*} X &= \int P(\lambda)\bar{x}(\lambda)\mathrm{d}\lambda,\\ Y &= \int P(\lambda)\bar{y}(\lambda)\mathrm{d}\lambda,\\ Z &= \int P(\lambda)\bar{z}(\lambda)\mathrm{d}\lambda, \end{align*}

CIE colour-matching functions

$X$, $Y$ and $Z$ can be normalized by dividing by their sum (at the expense of losing information about the brightness of the light):

$$ x = \frac{X}{X+Y+Z}, \quad y = \frac{Y}{X+Y+Z}, \quad z = \frac{Z}{X+Y+Z} = 1 - x - y $$

In this way, only two parameters, $x$ and $y$ are needed to describe the colour (more accurately, the chromaticity) of the light. The CIE standard chromaticity diagram is shown below.

CIE standard chromaticity diagram

Original image by user Spigget, licensed under Creative Commons Attribution-Share Alike 3.0 Unported.

Further conversion of $(x, y)$ to RGB values for output by a display device requires transformation by the appropriate chromaticity matrix. Geometrically, this maps points in the above colour "tongue" onto the subset of points within the RGB "gamut", the indicated triangular region. A colour system may be defined by a matrix of three primary colour chromaticities (the vertices of the triangle) and a white point: a set of chromaticiy coordinates defining the "colour" white for some purpose.

$$ \left( \begin{array}{lll} x_r & x_g & x_b \\ y_r & y_g & y_b \\ z_r & z_g & z_b \end{array} \right) \left( \begin{array}{l} r \\ g \\ b \\ \end{array} \right) = \left( \begin{array}{l} x \\ y \\ z \\ \end{array} \right). $$

Multiplication of the vector of $(x,y,z)$ values by the inverse of this matrix therefore gives the RGB values describing the corresponding colour within the system being used.

Not all $(x, y)$ pairs map to points within the RGB gamut (they would give negative values for one or more component): one way to deal with this is to "desaturate" by raising the values of all components equally until they are all non-negative.

The code below defines a class, ColourSystem, for representing and using colour systems, and instantiates a few particular examples. The CIE matching function is read in from the file cie-cmf.txt.

# colour_system.py
import numpy as np

def xyz_from_xy(x, y):
    """Return the vector (x, y, 1-x-y)."""
    return np.array((x, y, 1-x-y))

class ColourSystem:
    """A class representing a colour system.

    A colour system defined by the CIE x, y and z=1-x-y coordinates of
    its three primary illuminants and its "white point".

    TODO: Implement gamma correction

    """

    # The CIE colour matching function for 380 - 780 nm in 5 nm intervals
    cmf = np.loadtxt('cie-cmf.txt', usecols=(1,2,3))

    def __init__(self, red, green, blue, white):
        """Initialise the ColourSystem object.

        Pass vectors (ie NumPy arrays of shape (3,)) for each of the
        red, green, blue  chromaticities and the white illuminant
        defining the colour system.

        """

        # Chromaticities
        self.red, self.green, self.blue = red, green, blue
        self.white = white
        # The chromaticity matrix (rgb -> xyz) and its inverse
        self.M = np.vstack((self.red, self.green, self.blue)).T 
        self.MI = np.linalg.inv(self.M)
        # White scaling array
        self.wscale = self.MI.dot(self.white)
        # xyz -> rgb transformation matrix
        self.T = self.MI / self.wscale[:, np.newaxis]

    def xyz_to_rgb(self, xyz, out_fmt=None):
        """Transform from xyz to rgb representation of colour.

        The output rgb components are normalized on their maximum
        value. If xyz is out the rgb gamut, it is desaturated until it
        comes into gamut.

        By default, fractional rgb components are returned; if
        out_fmt='html', the HTML hex string '#rrggbb' is returned.

        """

        rgb = self.T.dot(xyz)
        if np.any(rgb < 0):
            # We're not in the RGB gamut: approximate by desaturating
            w = - np.min(rgb)
            rgb += w
        if not np.all(rgb==0):
            # Normalize the rgb vector
            rgb /= np.max(rgb)

        if out_fmt == 'html':
            return self.rgb_to_hex(rgb)
        return rgb

    def rgb_to_hex(self, rgb):
        """Convert from fractional rgb values to HTML-style hex string."""

        hex_rgb = (255 * rgb).astype(int)
        return '#{:02x}{:02x}{:02x}'.format(*hex_rgb)

    def spec_to_xyz(self, spec):
        """Convert a spectrum to an xyz point.

        The spectrum must be on the same grid of points as the colour-matching
        function, self.cmf: 380-780 nm in 5 nm steps.

        """

        XYZ = np.sum(spec[:, np.newaxis] * self.cmf, axis=0)
        den = np.sum(XYZ)
        if den == 0.:
            return XYZ
        return XYZ / den

    def spec_to_rgb(self, spec, out_fmt=None):
        """Convert a spectrum to an rgb value."""

        xyz = self.spec_to_xyz(spec)
        return self.xyz_to_rgb(xyz, out_fmt)

illuminant_D65 = xyz_from_xy(0.3127, 0.3291)
cs_hdtv = ColourSystem(red=xyz_from_xy(0.67, 0.33),
                       green=xyz_from_xy(0.21, 0.71),
                       blue=xyz_from_xy(0.15, 0.06),
                       white=illuminant_D65)

cs_smpte = ColourSystem(red=xyz_from_xy(0.63, 0.34),
                        green=xyz_from_xy(0.31, 0.595),
                        blue=xyz_from_xy(0.155, 0.070),
                        white=illuminant_D65)

cs_srgb = ColourSystem(red=xyz_from_xy(0.64, 0.33),
                       green=xyz_from_xy(0.30, 0.60),
                       blue=xyz_from_xy(0.15, 0.06),
                       white=illuminant_D65)

Here's one application of the ColourSystem class: to visualize the colour of a black body at a given temperature, and example given on this excellent page on the colour rendering of spectra. The spectral radiance of a black body is given by the Planck function:

$$ B(\lambda; T) = \frac{2hc^2}{\lambda^5}\frac{1}{\exp\left(\frac{hc}{\lambda k_\mathrm{B}T}\right) - 1} $$

Feeding $B(\lambda; T)$ as spec to the function ColourSystem.spec_to_rgb returns the RGB components corresponding to the colour of a black body at temperature $T$. This is visualized for different temperatures below.

import numpy as np
from scipy.constants import h, c, k
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

from colour_system import cs_hdtv
cs = cs_hdtv

def planck(lam, T):
    """ Returns the spectral radiance of a black body at temperature T.

    Returns the spectral radiance, B(lam, T), in W.sr-1.m-2 of a black body
    at temperature T (in K) at a wavelength lam (in nm), using Planck's law.

    """

    lam_m = lam / 1.e9
    fac = h*c/lam_m/k/T
    B = 2*h*c**2/lam_m**5 / (np.exp(fac) - 1)
    return B

fig, ax = plt.subplots()

# The grid of visible wavelengths corresponding to the grid of colour-matching
# functions used by the ColourSystem instance.
lam = np.arange(380., 781., 5)

for i in range(24):
    # T = 500 to 12000 K
    T = 500*i + 500

    # Calculate the black body spectrum and the HTML hex RGB colour string
    # it looks like
    spec = planck(lam, T)
    html_rgb = cs.spec_to_rgb(spec, out_fmt='html')

    # Place and label a circle with the colour of a black body at temperature T
    x, y = i % 6, -(i // 6)
    circle = Circle(xy=(x, y*1.2), radius=0.4, fc=html_rgb)
    ax.add_patch(circle)
    ax.annotate('{:4d} K'.format(T), xy=(x, y*1.2-0.5), va='center',
                ha='center', color=html_rgb)

# Set the limits and background colour; remove the ticks
ax.set_xlim(-0.5,5.5)
ax.set_ylim(-4.35, 0.5)
ax.set_xticks([])
ax.set_yticks([])
ax.set_facecolor('k')
# Make sure our circles are circular!
ax.set_aspect("equal")
plt.show()

Black body colours

Current rating: 4.3

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

Axel Schindler 7 years, 8 months ago

Thanks for sharing. The link for the cie-cmf.txt file is not working, however

Link | Reply
Current rating: 5

christian 7 years, 8 months ago

Thanks, Axel -- I've fixed this broken link now.

Link | Reply
Current rating: 4

Mukyu 5 years, 5 months ago

Are there any reference about the xyz2rgb part ... ? The matlab seems to be different from it. Matlab only implement matrix transform and gamma correction, but the code you proposed included "desaturating" and "normalize a rgb vector" ...

Link | Reply
Current rating: 5

asad 5 years, 2 months ago

i have problem with ''from colour_system import cs_hdtv''
is colour_system a module?
i can't find something like it in internet?and python gives error

Link | Reply
Currently unrated

christian 5 years, 2 months ago

colour_system.py is the filename you should give the code above the program you're running (the second block of code). I've added a comment line to make it a bit clearer.

Link | Reply
Currently unrated

Claude Falbriard 5 years, 1 month ago

Under Anaconda Juypyter Python 3.6 I was able to instantiate the class with the following statement:
cs = ColourSystem(red=xyz_from_xy(0.67, 0.33),
green=xyz_from_xy(0.21, 0.71),
blue=xyz_from_xy(0.15, 0.06),
white=illuminant_D65)
No need for class import as the python code is present in the previous notebook cell.
There is also an outdated matplotlib call, please use:
#ax.set_axis_bgcolor('k') # outdated
ax.set_facecolor('k')
Hope this also works for you. Regards, Claude

Link | Reply
Currently unrated

christian 5 years, 1 month ago

Yup, if you run both blocks in a Jupyter Notebook you don't need to import anything.

Thanks for pointing out the updated matplotlib call: I wrote this before set_axis_bgcolor was deprecated. I will update.

Link | Reply
Currently unrated

blackle 4 years, 10 months ago

for gamma correction I made this function:

def linear_srgb_to_rgb(rgb):
nonlinearity = np.vectorize(lambda x: 12.92*x if x < 0.0031308 else 1.055*(x**(1.0/2.4))-0.055)
return nonlinearity(rgb)

see https://en.wikipedia.org/wiki/SRGB#The_forward_transformation_(CIE_XYZ_to_sRGB) for the formula

Link | Reply
Currently unrated

MD 4 years, 1 month ago

I find that if i insert a spectrum containing every wavelength at equal intensity, it returns light orange, instead of white. How does that happen? Equal intensities for every wavelength should appear as white light (255, 255, 255).

Link | Reply
Currently unrated

fook 3 years, 11 months ago

Because the white point of sRGB displays is D65, not illuminant E.

Link | Reply
Current rating: 5

David 3 years, 10 months ago

Hi, I want to convert a reflectance spectrum into a colour. I guess that this can be done with that script.
Perhaps, a silly question, but I am inexperienced with coding.
How can I include my spectrum into the script, for this to be converted in a colour and RGC colour values?
Many thanks! Cheer, David

Link | Reply
Currently unrated

christian 3 years, 10 months ago

You should be OK to provide your spectrum as spec in a call to spec_to_rgb(), but it should be on the same grid of points as the colour-matching function, self.cmf: 380-780 nm in 5 nm steps. I hope it goes OK!

Link | Reply
Currently unrated

David 3 years, 10 months ago

Hi Christian, many thanks for replying!
Ok, my spectrum is from 380 to 780 nm with 5 nm in step
Shall I call it in this line?

cmf = np.loadtxt('cie-cmf.txt', usecols=(1,2,3))

Obviously, changing the name to that of my spectrum?

Thanks and cheers,

David

Link | Reply
Currently unrated

christian 3 years, 10 months ago

That's the colour matching function. You'll want to replace spec in this line: html_rgb = cs.spec_to_rgb(spec, out_fmt='html') with your spectrum, spec, as a suitable array.

Link | Reply
Current rating: 5

David 3 years, 9 months ago

Many thanks, Christian!

Link | Reply
Currently unrated

razan 8 months, 2 weeks ago

Hello David,
I am also trying to calculate the colour of a spectrum and have no knowledge with coding. Did you get the colour calculation correctly ? if so, do you mind sharing it please

Link | Reply
Currently unrated

Yannis 3 years, 7 months ago

Hi, thanks for providing the code for spectrum to RGB conversion. What I'd like to ask is how can I convert a spectrum in a different spectral range or with a different interval. Thank you in advance!

Link | Reply
Currently unrated

christian 3 years, 7 months ago

I think in this case it might be best to interpolate / bin your spectrum to the same interval and wavelength spacing as the colour-matching function.

Link | Reply
Current rating: 5

Hakanai 2 years ago

The integral at the top of the article has a dλ in it, but the code doesn't have an equivalent to that, as currently written.

If you throw in some code to multiply by Δλ as you're summing the values, it should make it work the same irrespective of how many samples you have.

Link | Reply
Currently unrated

Clay 3 years, 3 months ago

How can I avoid the intensity normalization? I made normalization optional in spec_to_xyz() and rgb_to_xyz(), but it's not working, and I still can't get brown, gray, or any dark colors. I'm trying to approximate mixing three monochromatic light sources of adjustable frequency and amplitude. Thanks

Link | Reply
Currently unrated

xiao7 2 years, 8 months ago

whether a limited 2D plane can Reproduce all visible colors in theory?

Link | Reply
Currently unrated

christian 2 years, 8 months ago

Well, you have three different types of cone receptor each giving a different response to a given light stimulus, and if you take out the total brightness you might approximate that you are left with two independent variables...

Link | Reply
Currently unrated

Hakanai 2 years ago

The one flaw I see in the logic proposed here is:

return XYZ / den

By dividing the by the denominator you will always get a point on the X+Y+Z=1 plane, which is possibly the intent here, but in general it will give the wrong result, because a spectrum with double all the values should give an XYZ result with double all the values.

Link | Reply
Currently unrated

christian 2 years ago

You do indeed lose the brightness information by normalizing in this way: the X + Y + Z = 1 plane contains the chromaticity information only.

Link | Reply
Currently unrated

Ben 1 year, 8 months ago

This is great! I want to include it in some notebook based learning material I am writing to support ocean colour applications. How should I credit you?

Link | Reply
Currently unrated

BF 1 year, 4 months ago

My understanding is that the x,y (chromaticity) returned by spec_to_xyz can be used to look-up a color on the CIE color-space "tongue", but that the same x,y cannot be used to look-up a color on the RGB "triangle". Is that correct? Thank you for the great tool, explanations, and comments Q&A

Link | Reply
Currently unrated

christian 1 year, 4 months ago

I'm glad you find it useful! The (x, y) returned by spec_to_xyz correspond to the CIE colour-space tongue; the RGB representation does not cover the whole of this colour-space: only colours within the indicated triangle can be represented. To obtain the actual values of (R, G, B), you need the transformation matrix (T) given.

Link | Reply
Currently unrated

Jaeyun Ko 10 months ago

Hello. Thanks for the spectrum-to-color code. How can I get some input of spectrum, to run this code? I know the valu of wavelength and intensity.
Can you help me writing spectrum input code?

Link | Reply
Currently unrated

Antoine 8 months, 3 weeks ago

Hello, I believe that cie-cmf.txt file you posted is valid for D65 lighting at 10 degrees, is that correct? How do I adjust it for different lightings from different angles?

Link | Reply
Currently unrated

Bruno S 5 months, 3 weeks ago

Thanks for the great insights. Building on the code shown, I have built a small online demo regarding this topic (runs python in the browser via pyodide). If anyone wants to try it out, it can be found at https://bluemi.github.io/color-spectrum-showcase.html.

Link | Reply
Currently unrated

christian 5 months, 3 weeks ago

Nice – and in the browser, too!

Link | Reply
Currently unrated

Дарья 1 week, 1 day ago

Hi! To study, you really need to figure out this code. I copied it, downloaded this text file and ran the code, but ended up with a runfile('C:/Users/дарья/Desktop/untitled2.py ', wdir='C:/Users/дарья/Desktop '). And how to get exactly the coordinates of the color? I really need help. thank you in advance

Link | Reply
Currently unrated

New Comment

required

required (not published)

optional

required