Deconvolution with Gauss Prior #5026
Conversation
|
Hello @alexjarosch! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found: There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻 Comment last updated at 2020-10-31 09:28:58 UTC |
|
I have created a test routine for my new deconvolution function in skimage/restoration/tests/test_restoration.py. It relies on a reference data file, similar to the test_richardson_lucy(): test. Even though I did register the data in skimage/data/_registry.py and the test working locally, I was not able to make it work on the continuous integration. Main problem is the test_download_all_with_pooch() test in skimage/data/tests/test_data.py, which breaks as my datafile is not yet on github. I am sure we can sort this out when my pull request is considered for merging. Until then I have commented out my respective test code in the respective files. |
There was a problem hiding this comment.
@alexjarosch thanks for this! The paper seems to be highly used so a priori I think this is a good addition to the library. The implementation is beautifully simple, too! I'd love to hear from other @scikit-image/core devs with more deconv experience (more than ~0 😂) what they think.
We're busy preparing our next release so I haven't done a complete review with reference to the papers, but I've done a first pass with some minor suggestions, as well as a link that I hope will get you unstuck with regards to the reference data.
Thanks again and I'm looking forward to the next pass!
skimage/data/tests/test_data.py
Outdated
| # https://github.com/scikit-image/scikit-image/pull/4666/files/26d5138b25b958da6e97ebf979e9bc36f32c3568#r422604863 | ||
| # https://github.com/scikit-image/scikit-image/pull/4666/files/ | ||
| # 26d5138b25b958da6e97ebf979e9bc36f32c3568#r422604863 |
There was a problem hiding this comment.
reverted in latest commit
skimage/data/_registry.py
Outdated
| # "restoration/tests/camera_dnp.npy": | ||
| # "05d6831f9c81f83fa820237052e894259f1a8f0ae2974449f60195439d369338", |
There was a problem hiding this comment.
To add sample data, please create a PR to gitlab.com/scikit-image/data (search for "pooch" in our contributing guide). You can see an example for doing this here:
https://github.com/scikit-image/scikit-image/pull/4922/files
There was a problem hiding this comment.
Just placed in a merge request on Gitlab. waiting for that one to merge and then I continue.
There was a problem hiding this comment.
Test still commented out, wait for merge on gitlab
skimage/restoration/deconvolution.py
Outdated
| Gx = np.fft.fft2(onesrow, s=(n, m)) | ||
| Gy = np.fft.fft2(onescol, s=(n, m)) | ||
| F = np.fft.fft2(psf, s=(n, m)) |
There was a problem hiding this comment.
Could we use fftn here in order to allow support for 3D+ images? Ditto for using np.flip(..., axis=i) instead of fliplr and flipud. I think with a small amount of effort this code can be generically n-dimensional.
There was a problem hiding this comment.
Please import the fft module from skimage._shared.fft as in the following file:
You can then call FFTs using that module:
fft.fftn(onescol, s=(n, m))This will result in using scipy.fft in place of numpy.fft whenver the user has SciPy >= 1.4. The SciPy implementation is preferable for three main reasons:
1.) it can perform transforms in single precision
2.) it can be configured to use multiple threads
3.) it can be overridden by alternative backends such as MKL or PyFFTW
There was a problem hiding this comment.
@grlee77 thanks, changed to skimage fft. done in latest commit.
@jni I still work on the conversion to N dimensions, which is not trivial. I will have to check if N dimensional makes sense for this type of de-convolution. As I saw the unsupervised wiener filter version is not N dimensional either, right?
Also would you then assume a N dimensional psf? Or a 2D psf applied to an N dimensional image?
There was a problem hiding this comment.
For the implementation proposed in alexjarosch#1, the psf has the same number of dimensions as the image. I can also add multichannel support if we want to allow e.g. 2d psf for a 2d + channels image, etc.
skimage/restoration/deconvolution.py
Outdated
| def DNP_Gauss_freq(image, psf, we=0.01, clip=True): | ||
| """ Deconvolution using natural image priors | ||
| with a Gaussian prior, solved in the frequency domain |
There was a problem hiding this comment.
| def DNP_Gauss_freq(image, psf, we=0.01, clip=True): | |
| """ Deconvolution using natural image priors | |
| with a Gaussian prior, solved in the frequency domain | |
| def gaussian_natural_prior(image, psf, we=0.01, clip=True): | |
| """Deconvolution using Gaussian natural image priors. |
(following PEP-8 and PEP-257)
I would also rename we to smoothness_weight.
Finally, I would like to see a bit of information in the argument description about what it does, rather than forcing users to go find a paper.
There was a problem hiding this comment.
done. thank you. done in latest commit
There was a problem hiding this comment.
The function name change to gaussian_natural_prior remains to be done (skimage functions should have snake case names)
|
Thank you for introducing me to this material, @alexjarosch. What are your thoughts on the sparse priors; shall we add those too? |
|
@jni thanks for the input and help, I will proceed and submit the test data formally and look into changing the code to 3D. So more commits to come. During that process I will answer your review suggestions. |
|
@stefanv, my pleasure, hope you enjoy this type of de-convolution as much as I do. Yes we should definitely add sparse priors as well, they are on my todo list. The solution routine works directly in the space domain and not in the frequency domain, so the solution need quite a bit more computational resources, but the results are stunning, even better control on "ringing". As my time allows I will look into coding the sparse priors as well. Help always appreciated ;) |
skimage/restoration/deconvolution.py
Outdated
| hs1 = np.int(np.floor((nf - 1) / 2)) | ||
| hs2 = np.int(np.floor((mf - 1) / 2)) |
There was a problem hiding this comment.
| hs1 = np.int(np.floor((nf - 1) / 2)) | |
| hs2 = np.int(np.floor((mf - 1) / 2)) | |
| hs1 = (nf - 1) // 2 | |
| hs2 = (mf - 1) // 2 |
There was a problem hiding this comment.
Thank you. Changed in latest commit
|
@stefanv , I have looked into the sparse prior codes and did some tests. For my applications (2D astronomical images and telescope psf's) the sparse prior methods do not create favourable results. So I placed an implementation in python lower on my priority list, however I still think for many applications it would be great to have. However computational cost will become rather unpleasantly large. |
alignment of the deconvolved result seems incorrect when this flip is present
Gradient terms are implemented via 1D FFTs and broadcasting for efficiency
|
Hi @alexjarosch. I was playing with this and went ahead and added n-dimensional support in PR alexjarosch#1 for your consideration. (If you merge that PR in your repository, the commits will show up here) I happened to be looking at deconvolution in relation to the following draft Dask blog post (dask/dask-blog#77) which is based on Richardson-Lucy, but clearly the algorithm proposed here would be faster than what is in the blog post because it does not require iteration. |
|
The data PR was recently merged, so things can proceed again here |
enhancements to Gauss prior deconvolution
| data += 0.1 * data.std() * np.random.standard_normal(data.shape) | ||
| deconvolved = restoration.gaussian_natural_prior(data, psf) | ||
|
|
||
| path = image_fetcher.fetch('restoration/tests/astronaut_dnp.npy') |
There was a problem hiding this comment.
How was this file generated? Was it saved using the implementation in this PR or was it created with a 3rd party implementation such as the Matlab code of the original authors?
If it is generated in the code in this PR, then this doesn't really test that the output here is valid, but will help ensure that any future changes keep output that is consistent with the current implementation.
There was a problem hiding this comment.
It is probably worth adding another test case similar to the above but where the original image is compared based on some metric to determine if it has improved. structural_similarity may be a good candidate for this. e.g. something like
from skimage.metrics import structural_similarity
ideal = rgb2gray(astronaut())
# CODE AS ABOVE TO DO THE DECONVOLUTION
ssim_blurred = structural_similarity(ideal, data)
ssim_deconv = structural_similarity(ideal, deconvoled)
assert ssim_deconv > ssim_blurredThere was a problem hiding this comment.
@grlee77 The file was created with a one to one translation from Matlab to Python based on the reference implementation from the original authors, referenced as [1] in the code description. I have checked the translated code on an individual line basis to ensure the results are the same in Matlab as well as in Python. So the image reflects the original code prior to our changes in the PR here.
There was a problem hiding this comment.
regarding the structural_similarity test, thats a great idea. However currently that test fails, probably due to either edge effects and/or the psf not being strong enough. Unfortunately I have no more time today to work on this, probably can take this on next week.
There was a problem hiding this comment.
@grlee77 The structural_similarity test works on data that has no noise added. I have implemented that test in the latest commit. Let me know if you are satisfied with this.
Description
I have implemented a new deconvolution method, based on natural image priors and solution done in the frequency domain. Only numpy is required. The main advantage to Richardson-Lucy is that artefacts on sharp edges, ie 'ringing' is strongly reduced.
Here the references
.. [1] http://groups.csail.mit.edu/graphics/CodedAperture/
.. [2] Levin, A., Fergus, R., Durand, F., & Freeman, W. T. (2007).
Deconvolution using natural image priors.
Massachusetts Institute of Technology,
Computer Science and Artificial Intelligence Laboratory, 3.
.. [3] Levin, A., Fergus, R., Durand, F., & Freeman, W. T. (2007).
Image and depth from a conventional camera with a coded aperture.
ACM transactions on graphics (TOG), 26(3), 70-es.
For reviewers
later.
__init__.py.doc/release/release_dev.rst.