It was a while ago now, but on April 27 I started explaining how I made this plot, which is from DIPUM3E (Digital Image Processing Using MATLAB, 3rd ed.):

In today's follow-up, I'll discuss how I computed the colors of the spectrum to display below the x-axis. I will use and refer to several DIPUM3E functions. These are available to you in MATLAB Color Tools on the File Exchange and on GitHub. The entire set of DIPUM3E functions is also on GitHub.

Contents

Visible wavelength

What should the x-axis limits be on this plot? In other words, what is the visible wavelength that we are interested in? You will see different limits being used in different places. The limits used here, 380 nm to 780 nm, are those given in Berns, R. S. (2000). Billmeyer and Saltzman's Principals of Color Technology, 3rd ed., John Wiley & Sons, NJ.

lambda = 380:780;

Find XYZ values for the spectral colors

The first computational step is to find the XYZ values for each value of lambda. This computation can be found in the DIPUM3E function lambda2xyz. But it is really simple: just interpolate into the CIE XYZ color matching functions, which are returned by the DIPUM3E function colorMatchingFunctions.

T = colorMatchingFunctions;
head(T)

ans =

  8×4 table

    lambda        x             y             z     
    ______    __________    __________    __________

     360       0.0001299     3.917e-06     0.0006061
     361      0.00014585    4.3936e-06    0.00068088
     362       0.0001638    4.9296e-06    0.00076515
     363        0.000184    5.5321e-06    0.00086001
     364      0.00020669    6.2082e-06    0.00096659
     365       0.0002321     6.965e-06      0.001086
     366      0.00026073    7.8132e-06     0.0012206
     367      0.00029308    8.7673e-06     0.0013727

To find the XYZ values for a specific lambda, such as 500, we can use interp1:

interp1(T.lambda,[T.x,T.y,T.z],500)

ans =

    0.0049    0.3230    0.2720

Or we can find the XYZ values for all of the wavelengths we are interested in.

XYZ = interp1(T.lambda,[T.x,T.y,T.z],lambda(:));

plot(lambda(:),XYZ)
title("XYZ values for spectral wavelengths")
legend("X","Y","Z")

Try a simple XYZ to RGB conversion

Let's try simply converting these XYZ values directly to RGB using the Image Processing Toolbox function xyz2rgb.

RGB = xyz2rgb(XYZ);
drawColorbar(RGB)

(The function drawColorbar is below. It uses the DIPUM3E function colorSwatches.)

That doesn't look like a very good spectral color plot to me. It seems uneven, with several patches that seem to be mostly one color. What's going on with this? We can see the problem if we draw a line plot of the RGB values. (The plotrgb and shadeGamutRegion functions are down below.)

close
plotrgb(lambda(:),RGB)
shadeGamutRegion

The gray shaded region shows the range between 0 and 1; this is the displayable range of colors. Everything outside that range (negative values, or values greater than 1), can't be displayed exactly. These out-of-range values are being clipped to the displayable range, and that leads to a bad results.

I'm going to show you the method used by the DIPUM3E function spectrumColors. The method is a variation of the one described in: Andrew Young (2012). Rendering Spectra (https://aty.sdsu.edu/explain/optics/rendering.html). Retrieved July 16, 2020.

Work with linear RGB values

First, let's convert the XYZ values to "linear" RGB values. The typical RGB values we see for image pixels are related nonlinearly to light intensity, and linear RGB values are more appropriate for the following averaging and scaling steps. The Image Processing Toolbox function xyz2rgb can optionally convert to linear values.

RGB_lin = xyz2rgb(XYZ,'ColorSpace','linear-rgb');
plotrgb(lambda(:),RGB_lin)
title("Linear RGB values for spectral wavelengths")

Heuristic scaling of linear RGB values

We want to modify those curves so that they fall within the range [0,1] and produce a reasonably accurate, smoothly varying, and attractive representation of the spectral colors.

The next thing we'll do is scale so that the maximum linear RGB value is 1.0. (Note: Young (2012) divides by a fixed value of 2.34.)

RGB_lin = RGB_lin / max(RGB_lin(:));

plotrgb(lambda(:),RGB_lin)

Now, one component at a time, and for each spectral color, mix in a sufficient amount neutral gray with the same Y to bring negative component values up to 0.

Y = XYZ(:,2);
for k = 1:3
   C = RGB_lin(:,k);
   F = Y ./ (Y - C);

   % No scaling is needed for component values that are already
   % nonnegative.
   F(C >= 0) = 1;

   RGB_lin = Y + F.*(RGB_lin - Y);
end

plotrgb(lambda(:),RGB_lin)

Next, to get brighter spectral colors, including a good yellow, scale up the linear RGB values, allowing them to get higher than 1.0. Then, for each color, scale all components back down, if necessary, so that the maximum component value is 1.0. Note: [Young 2012] uses a scale factor of 1.85.

RGB_lin = RGB_lin * 2.5;
S = max(RGB_lin,[],2);
S = max(S,1);
RGB_lin = RGB_lin ./ S;

plotrgb(lambda(:),RGB_lin)

Smooth the curves

Smooth out the linear RGB curves to eliminate discontinuous first derivatives. This helps the spectrum to appear smoother, reducing sharp transition points. Note: This step is not in [Young 2012].

RGB_lin = conv2(RGB_lin,ones(21,1)/21,'same');
plotrgb(lambda(:),RGB_lin)

Eliminate small negative numbers and numbers slightly greater than 1 that have been introduced through floating-point round-off.

RGB_lin = min(max(RGB_lin,0),1);

Convert to nonlinear RGB values for the final result

Convert to nonlinear sRGB values suitable for display on a computer monitor.

RGB = lin2rgb(RGB_lin);

plotrgb(lambda(:),RGB)

drawColorbar(RGB)

Next time, I'll talk about how to draw the spectral color scale underneath the plot.

Utility functions

function drawColorbar(rgb_colors)
f = figure;
f.Position(4) = f.Position(4) / 5;
colorSwatches(rgb_colors,0)
daspect([40 1 1])
axis off
end

function plotrgb(x,rgb_colors)
% Pick the colors we want to use from the normal line color order.
c = lines(7);
blue = c(1,:);
red = c(2,:);
green = c(5,:);
clf
plot(x,rgb_colors(:,1),'Color',red);
hold on
plot(x,rgb_colors(:,2),'Color',green);
plot(x,rgb_colors(:,3),'Color',blue)
hold off
grid on
axis tight
yl = ylim;
ylim([min(yl(1)-0.05,-0.05) max(yl(2)+0.05,1.05)])

legend("R","G","B")
xlabel("wavelength (nm)")

end

function shadeGamutRegion
xl = xlim;
xx = [xl(1) xl(2) xl(2) xl(1) xl(1)];
yy = [0 0 1 1 0];
p = patch(xx,yy,[0.5 0.5 0.5],"FaceAlpha",0.1,...
    "EdgeAlpha",0.1,"HandleVisibility","off");
end

Get the MATLAB code

Published with MATLAB® R2020a