function im_deblur = LucyMartin(im, len, angle, iterations)

% Restores Image Using Lucy - Richardson Method
%
% im_deblur = LucyMartin(im, len, angle, iterations)
%
% Inputs: im, len, angle, iterations
% Returns: im_deblur
%
% im:    Input Image
% angle: Blur Angle (Degrees) 
% len:   Blur length (Number of Pixels)
% iterations: Number of Iterations Through L-R Algorithm
%
% im_deblur: Restored Image

% Create psf of Blur
psf = fspecial('motion',len,angle);

% Performing Median Filtering 
im = medfilt2(abs(im));

% Convert psf to OTF of Desired Size
% OTF is Optical Transfer Function
OTF = psf2otf(psf,size(im));

% Initialize Estimate to Recorded Image
est = im;

i = 1;
while i<=iterations
    
    % Convert Estimate to Frequency Domain
    fest = fft2(est);
    
    % Multiply OTF with Estimate in Frequency Domain
    fblest = OTF.*fest;
    
    % Convert Blurred Image Estimate to Spatial Domain
    imest = ifft2(fblest);
    
    % Ratio of Blurred Image and Estimate of Deblurred Image
    iratio = im./imest;
    
    % Convert Ratio to Frequency Domain
    firatio = fft2(iratio);
    
    % Create Correction Vector
    corrvec = OTF .* firatio;
    
    % Convert Correction Vector to Spatial Domain
    icorrvec = ifft2(corrvec);
    
    % Multiply Correction Vector & Estimate of Deblurred Image, Create Next Estimate
    aftercorr = icorrvec.*est;
    est = aftercorr;

    i = i+1;
end

% Restored Image
im_deblur = abs(est);