%% Code for experimenting with model of hyperacuity Spillman-Werner '90
%% chapter 10, fig 18

% P. Perona, Caltech, Jan21, 2008

clear

PAUSE_LEN = 0.1;
%% Parameters of the stimuli experiment

% Vernier stimuli
LINE_LEN = 32;
LINE_WID = 2; 
LINE_GAP = 2; %% Half vertical distance between line ends
LINE_OFFSET = 6;%% Half horizontal distance between lines
ROWS = 2*LINE_LEN+1; % Rows and columns of the image
COLS = 2*LINE_LEN+1;

ima0 = zeros(ROWS,COLS);
ima0(1:2*LINE_LEN+1,LINE_LEN+1+[-LINE_WID:LINE_WID])=1;

ima1 = zeros(ROWS,COLS);
ima1(1:LINE_LEN-LINE_GAP,LINE_LEN + 1+ [-LINE_WID:LINE_WID]-LINE_OFFSET) = 1.0;
ima1(LINE_LEN+1+[1+LINE_GAP:LINE_LEN],LINE_LEN +1+[-LINE_WID:LINE_WID] + LINE_OFFSET)=1.0;

ima2 = ima1(:,end:-1:1);

figure(1);
subplot(2,3,1); imagesc(ima0); axis off; title('Vertical line');
subplot(2,3,2); imagesc(ima1); axis off; title('Vernier Left');
subplot(2,3,3); imagesc(ima2); axis off; title('Vernier Right');


%% Parameters of the simple cell

GAUS_X = 12;
GAUS_Y =  4;

k = [1 2 1]/4; %% Basic kernel

kx = k;
ky = k;
dyy = 10*[-1 2 -1];

% Blur kernels to obtain desired sigma
for s=1:round(2*GAUS_X^2),
    kx = conv2(kx,k,'full'); 
end;

for s=1:round(2*GAUS_Y^2),
    ky = conv2(ky,k,'full'); 
end;
ky = conv2(ky,dyy,'valid');

kx_halflen = floor(length(kx)/2);
ky_halflen = floor(length(ky)/2);


figure(2);
plot(-kx_halflen:kx_halflen,kx,'g'); hold on; plot(-ky_halflen:ky_halflen,ky,'r'); hold off
title('x and y kernels');

k2D = ky' * kx;
figure(3);
imagesc(k2D)

THETA = -0:5:180;
COLS_PLOT = ceil(sqrt(length(THETA)));
ROWS_PLOT = ceil(length(THETA/COLS_PLOT));

for theta=THETA,
    i_theta = find(theta == THETA);
    k_rot{i_theta} = imrotate(k2D,theta,'bilinear');
    k_rot{i_theta} = k_rot{i_theta};
    [rr,cc] = size(k_rot{i_theta});
    ci = ceil(rr/2);
    cj = ceil(cc/2);

    k_rot{i_theta} = k_rot{i_theta}(ci+[-LINE_LEN:LINE_LEN],cj+[-LINE_LEN:LINE_LEN]);
    figure(1); subplot(2,3,4); axis off
    imagesc(k_rot{i_theta});
    axis equal
    title(['\theta = ' num2str(theta)]);
    pause(PAUSE_LEN);
end;



%% Compute the response curve of the simple cell to the two stimuli

for i_theta = 1:length(THETA),
    pixelwise_r0 = k_rot{i_theta}.*ima0;
    pixelwise_r1 = k_rot{i_theta}.*ima1;
    pixelwise_r2 = k_rot{i_theta}.*ima2;
    
    r0(i_theta) = sum(pixelwise_r0(:));
    r1(i_theta) = sum(pixelwise_r1(:));
    r2(i_theta) = sum(pixelwise_r2(:));
    
    figure(4); 
    subplot(1,3,1); imagesc(pixelwise_r0);
    subplot(1,3,2); imagesc(pixelwise_r1);
    subplot(1,3,3); imagesc(pixelwise_r2);
    pause(PAUSE_LEN);
end;

figure(1); 
subplot(2,3,5); plot(THETA,r1,'r',THETA,r2,'g');
xlabel('orientation (degrees)');

title('Response to Vernier');
legend('Left','Right','Location','Best');
subplot(2,3,6); 
plot(THETA,r2-r1,'b','LineWidth',3);

hold on;
plot(THETA,r0,'m','LineWidth',3);
hold off;
xlabel('orientation (degrees)');
title('Response to Vernier and to Line');
legend('Right-Left','Line','Location','Best');

print -depsc vernier_figure