%Chapter 9,  Fig 9.3, 9.4
%Eric Dubois, updated 2018-12-17
% Generate six plots of components of the spectral density matrix for a
% sample color image
% Compute the power spectra of a color image. Consider removing gamma
% correction.
%Required function:  WelchSDM_est
clear all; close all;
N = 64;
% Specify the basis B to use by providing the matrix A_sRGB2B. Uncomment
% the one to use, including the text label for use in output file name.
% A_sRGB2B = eye(3);  % the sRGB basis
% bas = 'sRGB';
%1931 CIE XYZ to sRGB -- matrix from Poynton
A_XYZtosRGB = [3.240479 -1.537150 -0.498535; ...
              -0.969256  1.875992  0.041556; ...
               0.055648 -0.204043  1.057311];
%W0 to XYZ where W0 is the basis defined on pages 49-50 of my color book.
A_W0toXYZ = [ 0.9505  -0.1661 0.1397; ...
              1.0  0.0  0.0; ...
              1.0891 0.4013 0.4809];
%A_sRGB2B = inv(A_W0toXYZ)*inv(A_XYZtosRGB);  %WO from the color book.
%bas = 'WO';
%OZW of Zhang and Wandell
A_XYZtoOZW = [  0.279  0.720 -0.107;  ...
              -0.449 0.290 -0.077;  ...
              0.086 -0.590 0.501];
A_sRGB2B = A_XYZtoOZW*inv(A_XYZtosRGB); 
bas = 'OZW';
%
IMG = im2double(imread('05.tif'));
% undo sRGB gamma correction on three components to get sRGB linear
% components
IMG_lin = sRGB_gamma(IMG);
WelchSDM = WelchSDM_est(IMG_lin,N);
%Convert to the basis B
I = size(WelchSDM,1);
J = size(WelchSDM,2);
WelchSDM_B=zeros(size(WelchSDM));
for i=1:I
    for j=1:J
       WelchSDM_B(i,j,:,:)= A_sRGB2B * squeeze(WelchSDM(i,j,:,:)) *transpose(A_sRGB2B);
    end
end
% Convert to a dB scale
WelchSDM_dB = 10*log10(abs(WelchSDM_B));
minplot = floor(min(min(min(min(WelchSDM_dB))))/10)*10; %Min over all the spectra
maxplot = ceil(max(max(max(max(WelchSDM_dB))))/10)*10; %Max over all the spectra
% frequency grid
ux=-1/2:1/N:1/2-(1/N);
uy=ux;
[u, v]=meshgrid(ux,uy);
% Plot the magnitude of the six unique spectra
for C1 = 1:3
    for C2 = C1:3
       figure;
       mesh(u(1:2:N,1:2:N),v(1:2:N,1:2:N),WelchSDM_dB(1:2:N,1:2:N,C1,C2));
       axis([-.5 .5 -.5 .5 minplot maxplot]);
       view(-18,16);
       xlabel('\it u \rm (c/px)'), ylabel('\it v \rm (c/px)')
       zlabel('Spectral Density (dB)')
       title(strcat('Basis: ',bas, ' Component:',int2str(C1),int2str(C2)));
       colormap([0 0 0]);     % use black only
       set(gca,'ydir','reverse');  %make v axis point downward
       set(gcf,'Color',[1 1 1]);
       set(gca,'fontname','times')          %set the plot font to Times
    end
end