%Chapter 9,  Fig 9.5, 9.6
%Eric Dubois, updated 2018-12-17
% Program to convert a color image to OZW, filter O2 and O3, and see the
% effect on the spectral density matrix.
%Required function:  WelchSDM_est
clear all, close all;
% load the input image, remove gamma correction and convert to OZW
IMG_RGB = im2double(imread('05.tif'));
% undo sRGB gamma correction on three components to get sRGB linear
% components
IMG_RGB_lin = sRGB_gamma(IMG_RGB);
% convert to OZW basis
%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];
%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_sRGBtoOZW = A_XYZtoOZW*inv(A_XYZtosRGB);
for C = 1:3
    IMG_OZW(:,:,C) = A_sRGBtoOZW(C,1)*IMG_RGB_lin(:,:,1) + ...
        A_sRGBtoOZW(C,2)*IMG_RGB_lin(:,:,2) + ...
        A_sRGBtoOZW(C,3)*IMG_RGB_lin(:,:,3);
end
%Apply a low pass filter to component O2 and O3
load('least_pth_filter_cof');
IMG_OZW_filt(:,:,1)=IMG_OZW(:,:,1);
IMG_OZW_filt(:,:,2) = imfilter(IMG_OZW(:,:,2),h,'symmetric','same');
IMG_OZW_filt(:,:,3) = imfilter(IMG_OZW(:,:,3),h,'symmetric','same');
%Convert back to sRGB
A_OZWtosRGB = inv(A_sRGBtoOZW);
for C=1:3
    IMG_RGB_filt(:,:,C)= A_OZWtosRGB(C,1)*IMG_OZW_filt(:,:,1) + ...
                         A_OZWtosRGB(C,2)*IMG_OZW_filt(:,:,2) + ... 
                         A_OZWtosRGB(C,3)*IMG_OZW_filt(:,:,3);
end
%apply gamma correction for display
IMG_RGB_filt_gc = sRGB_gamma_correction(IMG_RGB_filt);
%Display input and output images
%figure, imshow(IMG_RGB);
%figure, imshow(IMG_RGB_filt_gc);
%compute and plot the relevant spectral density matrix estimates
N=64;
WelchSDM_OZW = WelchSDM_est(IMG_OZW,N);
WelchSDM_OZW_filt = WelchSDM_est(IMG_OZW_filt,N);
WelchSDM_RGB_filt = WelchSDM_est(IMG_RGB_filt,N);
% Convert to a dB scale
WelchSDM_OZW_dB = 10*log10(abs(WelchSDM_OZW));
minplot = floor(min(min(min(min(WelchSDM_OZW_dB))))/10)*10; %Min over all the spectra
maxplot = ceil(max(max(max(max(WelchSDM_OZW_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_OZW_dB(1:2:N,1:2:N,C1,C2));
       axis([-.5 .5 -.5 .5 minplot maxplot]);
       view(-18,16);
       xlabel('u (c/px)'), ylabel('v (c/px)')
       zlabel('Spectral Density (dB)')
       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
close all
% For the filtered OZW components
WelchSDM_OZW_filt_dB = 10*log10(abs(WelchSDM_OZW_filt));
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_OZW_filt_dB(1:2:N,1:2:N,C1,C2));
       axis([-.5 .5 -.5 .5 minplot-10 maxplot]);
       view(-18,16);
       xlabel('/it u \rm (c/px)'), ylabel('\it v \rm (c/px)')
       zlabel('Spectral Density (dB)')
       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
% for the filtered RGB components
WelchSDM_RGB_filt_dB = 10*log10(abs(WelchSDM_RGB_filt));
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_RGB_filt_dB(1:2:N,1:2:N,C1,C2));
       axis([-.5 .5 -.5 .5 -50 20]);
       view(-18,16);
       xlabel('\it u \rm (c/px)'), ylabel('\it v \rm (c/px)')
       zlabel('Spectral Density (dB)')
       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