Unstable results using deconvlucy function

When I am using a deconvlucy function from Image Processing Toolbox recently to do deconvolution on a 3-dimensional cube, I observe instability within 30 iterations.

Here are the codes:

I = {delta}; %delta is a 3-D image (300*300*300) that is to be deconvolved, it is the convolution of a 3-D delta function and a 3-D PSF
track_rms = zeros(N_iter,1); %track the rms noise in the cube at every iteration
for i = 1:N_iter
    I = deconvlucy(I, psf, 1);
    %track rms
    power = I{2}(pos(1),pos(2),pos(3)); %power from the source
    track_rms(i) = sqrt((sum(I{2}(:).^2) - power^2)/noise_len); %remaining noise in the cube
end

Here is the plot for RMS scattered power of one particular location of the delta function:

Here is another plot of the power at one particular location from another test, the intensity at (121,121,121) is not supposed to be 0:

What I want to ask is about the instability of the above results. Could there be any bugs in the deconvlucy code, or is this normal for deconvlucy method? Am I missing some parameters here, maybe the DAMPAR parameter? I would like to see a converging, or at least stable results. How can I improve the results? Thank you.

Roy

Answer 1

[Bjorn Gustavsson]

Without knowing exactly what your PSF looks like or is there are added noise to your delta-3D image, it is a bit difficult to give explicit explanations...

...but before you settle on a single metric like your track_rms, I would look at how the deconvolution progresses, but putting in something like

subplot(1,3,1)
imagesc(squeeze(I(pos(1),:,:)))
subplot(1,3,2)
imagesc(squeeze(I(:,pos(2),:)))
title(i)
subplot(1,3,3)
imagesc(squeeze(I(:,:,pos(3))))
drawnow

(possibly with some hard-set caxis and colorbar decorations and whatnot) That way you'll see how the deconvolution progresses.

My hunch is that 40 iterations is a large number of iterations, and rather likely something goes numerically unstable - deconvolution is a high-frequency enhancement operation...

HTH

[roy image]

Thank you for reply. My PSF is in a strange shape shown below:

It is basically two cones with the top one upside down connected to each other. Everything is simulated in MATLAB and I didn't add any noise to the images or PSF.

I think even within the first 40 iterations, the results are not stable, shown in the two plots in my original post. Are you suggesting early stop, i.e. stop at 8? If that is case, how do I decide which iteration to stop? Thanks.

[Bjorn Gustavsson]

Is that your psf? Should I interpret that as a point in the centre of your volume gets spread out on voxels at the surface of this double-cone? Or are the point spread out into the volumes "inside" the cones?

If I had to deconvolve anything like this I'd do one iteration at a time and then looked at slices after each to see what's going on.

HTH

[roy image]

Yes, your understanding is correct. The PSF is the surface of this double cone.

1) One of the support staff from MATLAB suggests that doing iterations using a for loop is different from using NUMIT parameter of deconvlucy.

for i = 1:100
    I = deconvlucy(I, psf, 1);
end

it is different from

I = deconvlucy(I, psf, 100);

Do you have any idea about this?

2) Another thing that comes to me is the DAMPAR parameter of this function. It says "Iterations are suppressed for pixels that deviate beyond the DAMPAR value from their original value", but I don't quite understand what does "Iterations are suppressed" mean?

Roy

[Bjorn Gustavsson]

Then I think (speculation warning) that you'll run into aliasing-type problems? I don't know how (or why) 40 one-step iterations would be different than one 40-step iteration. But I suggest that you check what's happening between iterations by looking at suitable slices (or whatever plotting makes sense in your case). That would be the fastest way to see what is happening...

HTH