Untitled

%%
% It was a while ago now, but on
% <https://blogs.mathworks.com/steve/files/illuminant-d65.png April 27>
% I started explaining how I made this plot, which is from
% <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *DIPUM3E*> (_Digital Image Processing Using MATLAB_, 3rd ed.):
%
% <<https://blogs.mathworks.com/steve/files/illuminant-d65.png>>
%
% 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
% <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *DIPUM3E*> functions. These are available to you in _MATLAB Color
% Tools_ on the
% <https://www.mathworks.com/matlabcentral/fileexchange/64161-matlab-color-tools
% File Exchange> and on <https://github.com/mathworks/matlab-color-tools
% GitHub>. The entire set of
% <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *DIPUM3E*> functions is also on
% <https://github.com/dipum/dipum-toolbox GitHub>.
%
%% 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
% <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *DIPUM3E*> function |lambda2xyz|. But it is really simple: just
% interpolate into the CIE XYZ color matching functions, which are
% returned by the <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *DIPUM3E*> function |colorMatchingFunctions|.

T = colorMatchingFunctions;
head(T)

%%
% 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)

%%
% 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
% <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *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
% <http://www.imageprocessingplace.com/DIPUM-3E/dipum3e_main_page.htm
% *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.
a=5;

%% 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


%%
% _Copyright 2020 The MathWorks, Inc._