%% MATLAB code for implementing Lucy-Richardson-Rosen algorithm
A=double(imread('File location')); %Read the image file of PSF and convert it into double precision arrays
A=A(:,:,1); %Extract the channel of interest - RGB, A(:,:,1) extracts Red A(:,:,2) extracts Green and A(:,:,3) extracts Blue
A=A/max(max(A));% Normalise the matrix between 0 and 1
imagesc(A) % Display image
B=double(imread('File location'));%Read the image file of object intensity and convert it into double precision arrays
B=B(:,:,1);
B=B/max(max(B));
figure;imagesc(B)


%%
PSF=A; % Assign A to PSF
OI=B; % Assign B to Object intensity
S1 = OI; % Initial guess solution is set as the object intensity pattern
OTF = psf2otf(PSF,size(OI)); % Convert PSF to OTF
iterations = 10;% Enter the iteration number - start with say 5
figure;colormap turbo % Open new figure and set colormap turbo
for i=1:iterations %start for loop
    i % show iteration number 
    FC = (ifft2(fft2(S1).*OTF)); % Forward Convolution
    ratio = OI./FC; % Calculate ratio
    ratio_f = fft2(ratio); % Fourier transform of ratio 
    alpha = 0.4; % Enter alpha value between 0 and 1; when alpha is 1 and beta is 1, it is Lucy-Richardson Algorithm
    beta = 1; % Enter beta value between 0 and 1
    residue = ifft2(conj((abs(OTF).^alpha).*exp(1i*angle(OTF))).*((abs(ratio_f).^beta).*exp(1i*angle(ratio_f)))); % Apply Non-linear reconstruction 
    S1 = residue.*S1; % Calculate the next solution
    imagesc(abs(S1).^1); % Display the solution
    pause(0.1) % pause to see the evolution of the solution
end % end for loop
result = abs(S1); % Final solution
imagesc(result) % Display solution