%%%
%%% Show how an image is processed by different filters in a hypercolumn
%%%

%%% (c) Pietro Perona. California Institute of Technology. January 21, 2002

K = [1 2 1]/4; %% Elementary Gaussian
D1 = [-1 0 1]; %% 1st derivative Kernel
D2 = [1 -2 1]; %% 2nd derivative Kernel
NS = 6; %% Number of scales to be considered


%% Generate 1D Gaussians at decreasing scales of resolution
G1{1}=[0 1 0];
D1G1{1} = conv2(G1{1},D1);
D2G1{1} = conv2(G1{1},D2);

for s=2:NS+3,
    K = conv2(K,G1{s-1});
    G1{s} = K;
    D1G1{s} = conv2(G1{s},D1);
    D2G1{s} = conv2(G1{s},D2);
end;

%%% Display the kernels
for s=1:NS,
    figure(1); subplot(NS,1,s); plot(G1{s});
    figure(2); subplot(NS,1,s); plot(D1G1{s});
    figure(3); subplot(NS,1,s); plot(D2G1{s});
end;

%%%%
%%%% Impulse responses in 2D
%%%%

%% Create image containing impulse
L = 2^(NS-1);
impulse = zeros(L+1,L+1);
impulse(round(L/2),round(L/2))=1;

%% Calculate and display impulse responses
for s=1:NS,
    figure(4);
    subplot(4,NS,s); imagesc(conv2(conv2(impulse,G1{s}','same'),G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,1*NS+s); imagesc(conv2(conv2(impulse,G1{s}','same'),D1G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,2*NS+s); imagesc(conv2(conv2(impulse,G1{s}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,3*NS+s); imagesc(conv2(conv2(impulse,G1{s+3}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
end;

%% Calculate and display response to image with SOS figure
ima = double(imread('test_image_sos.jpg'));
for s=1:NS,
    figure(5);
    subplot(4,NS,s); imagesc(conv2(conv2(ima,G1{s}','same'),G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,1*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D1G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,2*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,3*NS+s); imagesc(conv2(conv2(ima,G1{s+3}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
end;

%% Calculate and display response to image with lines figure
ima = double(imread('test_image_lines.jpg'));
for s=1:NS,
    figure(6);
    subplot(4,NS,s); imagesc(conv2(conv2(ima,G1{s}','same'),G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,1*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D1G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,2*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,3*NS+s); imagesc(conv2(conv2(ima,G1{s+3}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
end;

%% Same with noise added
ima = ima + (max(ima(:))-min(ima(:)))*rand(size(ima));
for s=1:NS,
    figure(7);
    subplot(4,NS,s); imagesc(conv2(conv2(ima,G1{s}','same'),G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,1*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D1G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,2*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,3*NS+s); imagesc(conv2(conv2(ima,G1{s+3}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
end;

%% Calculate and display response to image with circles
ima = double(imread('test_image_circles.jpg'));
for s=1:NS,
    figure(8);
    subplot(4,NS,s); imagesc(conv2(conv2(ima,G1{s}','same'),G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,1*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D1G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,2*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,3*NS+s); imagesc(conv2(conv2(ima,G1{s+3}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
end;

%% Calculate and display response to image with circles
ima = double(imread('test_image_vernier.jpg'));
for s=1:NS,
    figure(8);
    subplot(4,NS,s); imagesc(conv2(conv2(ima,G1{s}','same'),G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,1*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D1G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,2*NS+s); imagesc(conv2(conv2(ima,G1{s}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
    subplot(4,NS,3*NS+s); imagesc(conv2(conv2(ima,G1{s+3}','same'),D2G1{s},'same')); colormap(gray); axis off; axis equal;
end;