%Chapter 8,  Fig 8.3
%Eric Dubois, updated 2018-12-14
%Only the two-dimensional frequency responses are plotted in the figure
%Program to implement aspects of the Zhang and Wandell S-CIELAB model
%Parameters from Wandell

clear all; close all;
w1 = [0.921 0.105 -0.108];
w2 = [0.531 0.330];
w3 = [0.488 0.371];
s1 = [ 0.0283 0.133 4.336];
s2 = [0.0392 0.494];
s3 = [0.0536 0.386];
%Normalize weights to sum to 1 as per Johnson and Fairchild
%These values appear in Table 8.1
w1 = w1/sum(w1);
w2 = w2/sum(w2);
w3 = w3/sum(w3);
%Calculate k values so that each Gaussian filter has a DC gain of 1.
%These values appear in Table 8.1
k1 = ones(1,3)./(s1.^2*pi);
k2 = ones(1,2)./(s2.^2*pi);
k3 = ones(1,2)./(s3.^2*pi);
%Plot all the frequency responses
%frequencies are in cycles per degree
%Plot the 1D profiles (v=0) and the 2D perspective plots
umax = 30;   %maximum frequency to plot in cycles per degree
usp1 = 0.05; %spacing for 1D plot
usp2 = 2.;  %spacing for 2D plot
u1 = -umax:usp1:umax;
u = -umax:usp2:umax; v = u;
[U,V]=meshgrid(u,v);
%channel 1
for m = 1:3
    Z1_1D(:,m)=exp(-(pi^2)*((u1).^2)*(s1(m)^2));
    Z1(:,:,m)=exp(-(pi^2)*((U).^2+(V).^2)*(s1(m)^2));
end
Z1t_1D = w1(1) * Z1_1D(:,1) + w1(2) * Z1_1D(:,2) + w1(3) * Z1_1D(:,3);
Z1t = w1(1) * Z1(:,:,1) + w1(2) * Z1(:,:,2) + w1(3) * Z1(:,:,3);
%plot the 1D response
figure;
plot(u1,Z1t_1D,'k');
xlabel('spatial frequency (c/deg)');
ylabel('response');
set(gcf,'Color',[1 1 1]);
%plot the 2D response
figure;
mesh(U,V,Z1t)	% generate the perspective plot
colormap([0 0 0]);	% use black only
xlabel('u (c/deg)'), ylabel('v (c/deg)');
set(gca,'ydir','reverse');
axis([-30 30 -30 30 0 1.2]);
set(gcf,'Color',[1 1 1]);
%channel 2
for m = 1:2
    Z2_1D(:,m)=exp(-(pi^2)*((u1).^2)*(s2(m)^2));
    Z2(:,:,m)=exp(-(pi^2)*((U).^2+(V).^2)*(s2(m)^2));
end
Z2t_1D = w2(1) * Z2_1D(:,1) + w2(2) * Z2_1D(:,2) ;
Z2t = w2(1) * Z2(:,:,1) + w2(2) * Z2(:,:,2) ;
%plot the 1D response
figure;
plot(u1,Z2t_1D,'k');
xlabel('spatial frequency (c/deg)');
ylabel('response');
set(gcf,'Color',[1 1 1]);
set(gcf, 'PaperSize', [4.8 3.7]);
figure;
mesh(U,V,Z2t)	% generate the perspective plot
colormap([0 0 0]);	% use black only
xlabel('u (c/deg)'), ylabel('v (c/deg)');
set(gca,'ydir','reverse');
axis([-30 30 -30 30 0 1.2]);
set(gcf,'Color',[1 1 1]);
%channel 3
for m = 1:2
    Z3_1D(:,m)=exp(-(pi^2)*((u1).^2)*(s3(m)^2));
    Z3(:,:,m)=exp(-(pi^2)*((U).^2+(V).^2)*(s3(m)^2));
end
Z3t_1D = w3(1) * Z3_1D(:,1) + w3(2) * Z3_1D(:,2) ;
Z3t = w3(1) * Z3(:,:,1) + w3(2) * Z3(:,:,2) ;
%plot the 1D response
figure;
plot(u1,Z3t_1D,'k');
xlabel('spatial frequency (c/deg)');
ylabel('response');
set(gcf,'Color',[1 1 1]);

figure;
mesh(U,V,Z3t)	% generate the perspective plot
colormap([0 0 0]);	% use black only
xlabel('u (c/deg)'), ylabel('v (c/deg)');
set(gca,'ydir','reverse');
axis([-30 30 -30 30 0 1.2]);
set(gcf,'Color',[1 1 1]);

%Plot all the impulse responses
%Spatial units are degrees 
xmax = 10; %maximum spatial extent in degrees for channel 1
xsp1 = .04; % spacing between samples in degrees
x1 = -xmax:xsp1:xmax;
%Channel 1
for m=1:3
    H1_1D(:,m) = w1(m)*k1(m)*exp(-x1.^2/(s1(m)^2));
end
H1t_1D = H1_1D(:,1) + H1_1D(:,2) + H1_1D(:,3);
%plot central portion from -1 to 1 deg at full amplitude
figure;
plot(x1(226:276),H1t_1D(226:276),'k')
xlabel('spatial extent (deg)');
ylabel('impulse response');
axis([-1, 1, -50, 450]);
set(gcf,'Color',[1 1 1]);
set(gcf, 'PaperSize', [4.8 3.7]);
%plot full extent but vertically magnified to show the small negative lobe
figure;
plot(x1,H1t_1D,'k')
xlabel('spatial extent (deg)');
ylabel('impulse response');
set(gcf,'Color',[1 1 1]);
axis([-xmax,xmax,-.005,.005])
%Channel 2
xmax2 = 1;
xsp2 = .01; % spacing between samples in degrees for channel 2
x2 = -xmax2:xsp2:xmax2;
for m=1:2
     H2_1D(:,m) = w2(m)*k2(m)*exp(-x2.^2/(s2(m)^2));
end
H2t_1D = H2_1D(:,1) + H2_1D(:,2) ;
figure,
plot(x2,H2t_1D,'k')
xlabel('spatial extent (deg)');
ylabel('impulse response');
set(gcf,'Color',[1 1 1]);
%Channel 3
for m=1:2
     H3_1D(:,m) = w3(m)*k3(m)*exp(-x2.^2/(s3(m)^2));
end
H3t_1D = H3_1D(:,1) + H3_1D(:,2) ;
figure,
plot(x2,H3t_1D,'k')
xlabel('spatial extent (deg)');
ylabel('impulse response');
set(gcf,'Color',[1 1 1]);