Change of initial estimate of lucy

[roflmaostc]

Hey,

with the current estimate, Lucy Richardson fails to remove the Poisson noise. Actually LR tries to generate details and therefore it won't touch the high frequencies (noise).
Therefore we change the initial estimate to something else namely the measured image convolved with the PSF.

A minimal example to convince yourself, is below.
You can also see the final image here. The left one is the current lucy, the right one is the fixed one.

using FFTW, Noise, Deconvolution, TestImages, Images

function lucy2(observed::AbstractArray, psf::AbstractArray; iterations::Integer = 1)
    @assert size(observed) == size(psf)
    @assert iterations >= 0

    psf_ft = fft(psf)
    psf_ft_conj = conj.(psf_ft)

    function lucystep(e)
        return e .* real(ifft(fft(observed ./ ifft(fft(e) .* psf_ft)) .* psf_ft_conj)) 
    end

    estimated = real(ifft(fft(observed) .* psf_ft))
    for i in 1:iterations
        estimated = lucystep(estimated)
    end

    return estimated
end

img = channelview(testimage("fabio_gray_256"))

x = LinRange(-40, 40, 256)
y = transpose(x)

# create PSF
psf = exp.(- (x.^2 .+ y.^2))
# normalize PSF to 1
psf ./= sum(psf)
psf = ifftshift(psf)

# blur
img_b = real(ifft(fft(img) .* fft(psf)))
# Poisson noise
img_noisy = poisson(img_b, 500)
# current from package
res_lr1 = lucy(img_noisy, psf, iterations=20);
# fixed one
res_lr2 = lucy2(img_noisy, psf, iterations=20);
colorview(Gray, [res_lr1[100:200, 100:200] res_lr2[100:200, 100:200]])

[jakubwro]

Awesome! Will review before the weekend. I have no access to my laptop now :(

[roflmaostc]

Thanks, but no hurry 😄 I'm already using the fixed local version.
The CI failed only a few times, maybe a rebuild fixes it.

[jakubwro]

The CI failed only a few times, maybe a rebuild fixes it.

Yes, rerun helped.

Sorry, I remember about this review, but this Poisson noise case is a bit tricky. I am not sure if initial blurring by PSF is the best way to cancel it. I am thinking about using different low pass filter related more to added Poisson noise, not PSF.
@roflmaostc, if this works for you, let's merge it and I will try to work on it later using my examples. If you have some some papers please share :)

[roflmaostc]

For the initial estimate you can even start with a constant image, RL will handle it.
I can't provide a proper source right now, but it is somehow "standard" to do so (according to my supervisor 😂).

Edit:
This paper mentions it on the second page: http://adsabs.harvard.edu/full/1994A%26AS..108..409B

[codecov-commenter]

Codecov Report

Merging #20 into master will not change coverage.
The diff coverage is 100.00%.

@@            Coverage Diff            @@
##            master       #20   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            3         3           
  Lines           35        35           
=========================================
  Hits            35        35           

Impacted Files	Coverage Δ
src/lucy.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e0c4a7c...e63127a. Read the comment docs.