function z = Blur (mode,x,psf)
% The blur operator.
% Applies a blurring operator to the image x,  implemented by convolution
% with kernel passed in the parameter struct 'par'.
% input:	'mode' admits values 'fwd', for the forward transform, and
%                     'trn' for the adjoint transform.
%                     'x' - the input column vector, which usually
%                     represents a 2D image reordered as a column.
%                     'par' - the struct of parameters generated by
%                     setup_Blur.m
% output:   'z' is the result z=Hx (in case of forward transform) or z=H'x
%                       (in case of adjoint transform).
% blurring kernel:      a kernel of size 15x15 is used, which (i,j)
%                     element, -7<i,j<7 , is computed by kernel(i,j)=(1+ii^2+jj^2)^(-1), 
%                     with following normalization of the sum of the elements in
%                     kernel to 1.
% edges handling: we add margins of width 7 to the input 2D umage,
%                     generated as mirror reflection of the image edges.
%                     After the convolution with the blurring kernel, the
%                     14-pixel wide boundaries of resulting image are
%                     removed.

ws   =  size(psf);
t    =  (ws-1)/2;

if mode == 'fwd',
    s  = x;
    se = [s(:,t:-1:1,:), s, s(:,end:-1:end-t+1,:)];
    se = [se(t:-1:1,:,:); se; se(end:-1:end-t+1,:,:)];
    
    if size(x,3)==3
        z(:,:,1) = conv2(se(:,:,1),psf,'valid');
        z(:,:,2) = conv2(se(:,:,2),psf,'valid');
        z(:,:,3) = conv2(se(:,:,3),psf,'valid');
    else
        z = conv2(se,psf,'valid');
    end

elseif mode == 'trn', 
    y  = x;
    ye = [y(:,t:-1:1,:), y, y(:,end:-1:end-t+1,:)];
    ye = [ye(t:-1:1,:,:); ye; ye(end:-1:end-t+1,:,:)];
    
    if size(ye,3)==3
        z(:,:,1) = conv2(ye(:,:,1),psf,'valid');
        z(:,:,2) = conv2(ye(:,:,2),psf,'valid');
        z(:,:,3) = conv2(ye(:,:,3),psf,'valid');
    else
        z = conv2(ye,psf,'valid');
    end
end
