spatially variant deconvolution (Richardson-Lucy)

I'm working on a deblurring problem, but need some help.

I have a spatially variant blurred image, thus multiple point spread functions are used to deblur the image. Now I have to recover the originial image by spatially variant deconvolution with the Richardson-Lucy deconvolution method. I would like to use the deconvlucy algorithm of MATLAB, however, this function is spatially invariant (only one psf can be used).

I read in several articles that it is possible to change the deconvlucy algorithm to a spatially variant one, but I haven't managed it myself.

Is someone familiar with this problem and able to help me?

My code so far:

if true
P = phantom('Modified Shepp-Logan',150);
subplot(1,3,1), imshow(P)
%%Generate space-variant Blurring
Sigma2 = 2;
Sigma3 = 4;
% Find pixels with  0 < value < 0,2 , give blurring 3
BW2 = (P <= 0.2);
BW2idx = find(BW2);
PBW2 = P.*BW2;
Gaus2 = imgaussfilt(PBW2, Sigma2);
% Find pixels with value = 0.3, give blurring 4
BW3 = (P > 0.2);
BW3idx = find(BW3);
PBW3 = P.*BW3;
Gaus3 = imgaussfilt(PBW3, Sigma3);
% Sum images
IM = imadd(Gaus2, Gaus3);
BlurredImage = IM;
subplot(1,3,2), imshow(BlurredImage)
%%Generate psf
% Place calculated sigmas at right place in the image, according to the segmentation of the tissues
hsize = 3*Sigma2+1;
Sigmas = 0.5*ones(150);
Sigmas(BW2idx)= Sigma2;
Sigmas(BW3idx)= Sigma3;
V = 150+hsize-1;
tsize = hsize-1;
rsize = tsize/2;
FinalKernel = zeros(V,1,'double');
for q = 1:150 % y-axis
  for w = 1:150 % x-axis
      if Sigmas(q,w) == 2
          psf = fspecial('gaussian',hsize,Sigma2);
      elseif Sigmas(q,w) == 4
          psf = fspecial('gaussian',hsize,Sigma3);
      end
      Matrix{q,w} = psf;
  end
end
BlurredImage = padarray(BlurredImage,[rsize rsize],0,'both')
NUMIT = 10
for q = 1:150 % y-axis
  for w = 1:150 % x-axis
      PSFMatrix = Matrix{q,w};
      % BlurredImage moet ook groter worden gemaakt. 
     Deblurred = deconvlucy2(BlurredImage(q+3,w+3), PSFMatrix,NUMIT);
     DeblurredImage{q,w} = Deblurred;
  end
end
subplot(1,3,3), imshow(Deblurred,[])
%colormap(gca, fire)
end

Answer 1

[Eric Brost]

Did you ever find a solution to this problem? I am starting to research the same problem, similar to the one solved in this paper: https://link.springer.com/chapter/10.1007/978-3-540-72823-8_46

I will let you know of what I find that works as I move forward.