function devolver=gndecon(noisy,gn);

%% Script for deconvolution of noisy images with BOUNDARY DATA
%% The images (NOISY) and 64x64 noisy convolver (GN=g+n) are in the file
%% noisy.mat
%% The main idea is that if one has the full convolved images including the
%% boundary data then one can ESTIMATE the deconvolver as a difference over
%% the averaged noisy images.
%% The principle is as follows. Matlab convolution conv(X,Y) returns a
%% vector U=conv(X,Y) whose k-th entry is the sum of {x_i y_j: i+j=k+1}
%% Then u_1=x_1y_1; u_2=x_1y_2+x_2y_1 ETC
%% suppose now that you want to estimate a PSF P from several noisy signals
%% NOISY =conv(P,ORIGINAL). Suppose that you have enough data so that the
%% hypothesis MEAN(ORIGINAL)=C holds where C is a constant. Then you have 
%% MEAN(NOISY)=(M_1,M_2,...) where M_1=C p_1, M_2= C(p_1+p_2), M_3=
%% C(p_1+p_2+p_3) ETC.
%% IN other words, 
%% DIFF(MEAN(NOISY))=C(p_2,p_3,...)
%% so if you adnoin C p_1 at the beginning then you get a constant multiple
%% of the convolver. Therefore, the deconvolved image will be a constant
%% multiple of the clean image.

%load noisy.mat       % load noisy images and convolver

%% Build deconvolver assuming that PSF is g'*g (tensor product)
%% if this is false the method will fail but there is a corresponding
%% method that works by averaging the UL,UR,LL,LR blocks of the noisy
%% images directly


s1=size(noisy,1);
s2=size(noisy,2);
s3=size(noisy,3);   % size of noisy images array

g1=size(gn,1);
g2=size(gn,2);      % size of true PSF matrix, both should be equal
g=floor(g2/2);
h=floor(g1/2);


    

edge=zeros(2*g+2,s3);       % extract border of the 'right' row
for i=1:s3
    el=noisy(1:g+1,h,i);
    er=noisy(s1-g:s1,h,i);
    edge(:,i)=cat(1,el,er);
end

edgeave=mean(edge,2);       % take average of border values
diff1=diff(edgeave(1:g));
diff2=flipud(diff(flipud(edgeave(g+3:2*g+2))));
volv1=cat(1,edgeave(1),diff1,diff2,edgeave(2*g+2));     % convolver from border average

%% in general you would take the corners of the images and define the
%% convolver matrix by taking suitable differences of the corner averages
%%


figure;
plot(gn(g,:)./max(gn(g,:)));          % normalized g-th row of convolver
hold on
plot(volv1./max(volv1),'r');        % normalized convolver

volv2=volv1*volv1';    % g1 x g1 convolution matrix

delta=zeros(s1-g1+1,s2-g2+1);
delta(1,1)=1;   
devolver=conv2(delta,volv2);