Add blind Richardson-Lucy deconvolution (supersedes #3524)

[bioimage-analysis]

Description

Follow up on @Anti-Xyz pull request #3524, added most of the modifications.

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[pep8speaks]

Hello @bioimage-analysis! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file skimage/data/_registry.py:

Line 121:80: E501 line too long (110 > 79 characters)
Line 122:80: E501 line too long (120 > 79 characters)

• In the file skimage/restoration/deconvolution.py:

Line 349:38: E226 missing whitespace around arithmetic operator
Line 349:41: E226 missing whitespace around arithmetic operator
Line 349:44: E226 missing whitespace around arithmetic operator

• In the file skimage/restoration/tests/test_restoration.py:

Line 112:5: E265 block comment should start with '# '

Comment last updated at 2020-05-26 16:29:03 UTC

[sciunto]

I made a quick reading to solve the CI issue.

[sciunto]

There is still a failure in our CI. See https://travis-ci.org/github/scikit-image/scikit-image/jobs/687078292

I encourage you to build and run the test on your machine. Please, refer to https://scikit-image.org/docs/stable/contribute.html#build-environment-setup

[bioimage-analysis]

@sciunto sorry for the delay, things are reopening and I am running a little out of time. To test blind richardson, we need to fetch "reconstruction_blind_RL.npy", but the file is not on the scikit-image server so I am not sure how to deal with that with the function fetch?

[sciunto]

No problem! You need to install pooch

[bioimage-analysis]

@sciunto sorry, I was not clear, everything work well on my computer when I use a local path, the probably is if I try to fetch the file on scikit-image, it doesn't work. Here is the error I get:

KeyError                                  Traceback (most recent call last)
<ipython-input-20-0de14c355e12> in <module>
----> 1 test_blind_richardson_lucy()

<ipython-input-6-cc4e3b6eb27d> in test_blind_richardson_lucy()
     30                                                  iterations=iterations)
     31 
---> 32     path = fetch('restoration/tests/reconstruction_blind_RL.npy')
     33     #path = 'skimage/restoration/tests/reconstruction_blind_RL.npy'
     34 

~/Desktop/scikit-image/skimage/_shared/testing.py in fetch(data_filename)
    214     """Attempt to fetch data, but if unavailable, skip the tests."""
    215     try:
--> 216         return data._fetch(data_filename)
    217     except ConnectionError:
    218         pytest.skip(f'Unable to download {data_filename}')

~/Desktop/scikit-image/skimage/data/__init__.py in _fetch(data_filename)
    125     """
    126     resolved_path = osp.join(data_dir, '..', data_filename)
--> 127     expected_hash = registry[data_filename]
    128 
    129     # Case 1:

KeyError: 'restoration/tests/reconstruction_blind_RL.npy'

[jni]

@bioimage-analysis you need to add a registry entry in skimage/data/registry/_registry.py, that maps the path of the file to the SHA256 (I think) hash of the file.

[jni]

@bioimage-analysis it's looking good! =D

[sciunto]

@jni should we host the file in https://gitlab.com/scikit-image/data ?

[jni]

I would add blank lines between paragraphs here.

[jni]

@sciunto

Yes I think so!

@bioimage-analysis

It looks like the documentation build is failing:

https://dev.azure.com/scikit-image/scikit-image/_build/results?buildId=3464&view=logs&j=31eb6c88-23c5-50ae-eda4-a7f72508e1eb&t=f7291549-81f8-5748-7064-2be81963872c&l=570

d:\a\1\s\doc\source\auto_examples\filters\plot_blind_deconvolution.rst:24:Unexpected indentation.

Have you tried making the docs locally? I can't say that I understand the sphinx error because the lines don't correspond to the lines in the source .py file. I guess the best way is to attempt the build and then inspect the generated .rst file.

[jni]

@sciunto

Yes I think so!

Having said this, 156KB is not a huge deal for this repo, and this PR has taken a long time already, so I wouldn't hold it up for this.

[bioimage-analysis]

Getting closer, thank you for the advice @jni

[anki-xyz]

@bioimage-analysis @jni @sciunto thanks for taking care of this PR!

[emmanuelle]

Suggested change

	We present here a usage example of a blind deconvolution algorithms
	We present here a usage example of a blind deconvolution algorithm

[emmanuelle]

the names are misleading here, if they mean width and height, since it's the other way around.

[emmanuelle]

I think you can just invert them :-)

[emmanuelle]

Suggested change

[emmanuelle]

Thanks a lot @bioimage-analysis ! I left a couple of small comments about the tutorial. About the algorithm itself, I had a quick look at the paper by Fish et al. A comment in the code says that the implementation is not exactly what's in the paper because it does not work, could you please elaborate on this? Also, what is the message of the tutorial about the ideal number of iterations? I tried reducing the intensity of the noise and even with a very small noise the loss function was never monotonous (whereas it is in the paper for a low noise intensity), what do you think about this?

[anki-xyz]

I revisited the code and honestly, I think the code is not fully correct. And this is why "the hack" is needed. Everything actually refers to this part of the paper:

Where g(x) is the PSF estimation, and f(x) the reconstruction, c(x) is the corrupted image.

However, there should be k iterations, i.e. cycles, of m RL iterations (see also Fig. 1).

Here I programed it from scratch, reconstructing a cross similar to the paper:

# Create a cross
cross = np.zeros((128, 128))
cross[40:80, 55:65] = 1
cross[55:65, 40:80] = 1

# Add Gaussian and Poisson noise
# c is corrupted image
c = gaussian(cross, sigma=5) * np.random.poisson(lam=20., size=cross.shape)

# Initial guess for object distribution
f = np.full(c.shape, 0.5)
# f = c.copy() # not sure if this makes sense at all, but I keep it here for reference
# Initial guess for PSF
g = np.full(c.shape, 0.5)

# Hyperparameter
cycles = 50
m = 5 # RL iterations
epsilon = 1e-9 # No division by 0

# main cycles
for k in range(cycles):
    for i in range(m):  # m RL iterations, refining PSF
        g = convolve((c / (convolve(g, f, mode='same')+epsilon)), f[::-1,::-1], mode='same') *  g

    for i in range(m): # m RL iterations, refining reconstruction
        f = convolve((c / (convolve(f, g, mode='same')+epsilon)), g[::-1,::-1], mode='same') *  f

And this is what is now (!) happening:

Thus, I think the code is

• not correct
• misses the split between RL iterations of refining PSF and reconstruction, respectively
• does not use main cycles
• and thus, not producing what @emmanuelle is expecting.