Fix wiener filter for real value for the noise

[roflmaostc]

Hey,

I think I found a bug in the Wiener implementation. There is the possibility to use a single value for the noise instead of a noise array. This option is actually quite useful.

However, the code transforms this noise parameter into a array where the value of noise is at each entry. Then, in the wiener implementation we are performing a Fourier transform of this constant array. This leads to the problem, that the noise power spectrum is zero everywhere except at (1, 1) where it has the value of noise (Fourier transform of a constant gives a delta peak).

In my solution, instead of passing a constant array to the wiener function we just call the respective _wiener* function ourself but directly with a constant array representing the noise power spectrum.

I compared my code and the existing one. In the current implementation it can occur that the final image contains NaNs because sometimes it happens that we divide by (a complex) zero. This won't happen with a constant noise power spectrum.
Also I discovered that for the case wiener(noisy_img, noisy_img, 10) the current code only changes the value of fft(noisy_img)[1,1] therefore reducing only the global brightness of the noisy_img. This is reasonable if we look at the Wiener filter definition.

I didn't update the test code because I don't know how it was generated.

Thanks,

Felix

[giordano]

To be honest, I'm not an expert in digital signal processing, this was mostly a fun project 🙂

What you say makes sense to me, but I'd like also a review by @jakubwro, if he feels like to.

I compared my code and the existing one. In the current implementation it can occur that the final image contains NaNs because sometimes it happens that we divide by (a complex) zero. This won't happen with a constant noise power spectrum.

Could you please add a regression test for this?

I didn't update the test code because I don't know how it was generated.

In the worst way possible: run the function, get the result, use it as reference value for the test. If you know some reference to use for testing, please do use it by all means!

[roflmaostc]

Ah ok, I see 😄

Yeah, I'm also working on another issue right now.
In the Lucy Richardson, by derivation the psf is flipped and convolved. You are using the complex conjugate in Fourier Space which is analytically the same but causes sometimes a few problems.

I'll try to get into the test cases in the next few days and I'll keep you updated if I find a proper solution for it 👍

[jakubwro]

Hi!
I didn't run that but after reading the code and this PR discussion my understanding is that previously we treated noise real parameter as brightness level, now we are interpreting it as white noise level. The new approach makes more sense for me and is matching what the doc is saying: https://juliadsp.github.io/Deconvolution.jl/dev/#wiener-function-1

In the Lucy Richardson, by derivation the psf is flipped and convolved. You are using the complex conjugate in Fourier Space which is analytically the same but causes sometimes a few problems.

@roflmaostc, the fact is that I was not really sure if should I flip PSF or use FT conjugate. I think code was shorter and ringing was less visible with approach I used. Would love to see what are your thoughts. Please be careful about cases when PSF is not symmetrical 😄

[jakubwro]

Code looks good, but tests need to be green. We need also to make sure that scripts from doc/ are still working.

[giordano]

Thanks @jakubwro for the review! Once we have the new tests and documentation is up-to-date I'll be happy to merge!

[jakubwro]

I have some suggestion and I am not insisting to include it. There is some code repetition. I would trade some performance for readability.

function whitenoise(noiselevel::Real, size)
    result = zeros(size)
    result[1] = noiselevel
    return result
end

function wiener(input::AbstractArray, signal::AbstractArray, noiselevel::Real)
    return wiener(input, signal, whitenoise(noiselevel, size(input)))
end

function wiener(input::AbstractArray, signal::AbstractArray,
                noiselevel::AbstractArray, blurring::AbstractArray)
    return wiener(input, signal, whitenoise(noiselevel, size(input)), blurring)
end

[roflmaostc]

Yeah, I was also not sure how to handle this repetition.

You mean as performance drop this extra fft call on the whitenoise array?
Furthermore I was wondering if we should rewrite the wiener function itself because we can also save the whitenoise array completely and instead just add this noise value point wise in the denominator.

return real(ifft(Y .* conj(H) ./ (abs2.(H) .+ noise ./ S)))

Btw, as you might have seen I'm also updating the documentation. I recognized that the example in the documentation is outdated but working code is in time-series.jl. So I transferred the code from this file to the documentation.
Also I would like to add an example of a deconvolution of a time series with the Wiener filter.

[jakubwro]

You mean as performance drop this extra fft call on the whitenoise array?

Yes, that was my concern.

Furthermore I was wondering if we should rewrite the wiener function itself because we can also save the whitenoise array completely and instead just add this noise value point wise in the denominator.

I am not sure if white noise is so common to optimise for this case. If you and @giordano agree, let's do that.

Btw, as you might have seen I'm also updating the documentation

Thanks!

Also I would like to add an example of a deconvolution of a time series with the Wiener filter.

This is awesome, but I think some kind of STFT should be used for that. I would like to see such example :)

[roflmaostc]

I pushed the noisy time series and also added a example of using the constant noise parameter.

In this case it doesn't work so well because the image is affected heavily by noise. With less noise and more blurring, the constant noise parameter is actually more powerful.

[jakubwro]

Docs look good now, let's ignore my comment and after making tests pass I will merge

[roflmaostc]

I didn't have a proper idea for fixing the testcases, so I just copied the print output of the function.

[codecov-io]

Codecov Report

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

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

Impacted Files	Coverage Δ
src/wiener.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 6acd8ed...b357e36. Read the comment docs.

[jakubwro]

I didn't have a proper idea for fixing the testcases, so I just copied the print output of the function.

That's fine. For lucy I am using different approach, not relying on fixed data, but expecting that difference of restored to original is smaller than input to original. I'm not sure if this is way to go. Due to code duplication code coverage is smaller now, let's make it 100% again :)

[roflmaostc]

Should I fix something or what do you suggest?

[jakubwro]

Should I fix something or what do you suggest?

The perfect solution will be to avoid duplication at all 😄 I think we can just add missing test cases. If you can please do, otherwise I will do it tomorrow, I have no access to my computer with Julia environment. Thanks for fixing the bug!

[roflmaostc]

You mean extending the previous test cases for wiener or just fixing the failing ones?
Right now no test cases fail anymore

[jakubwro]

You mean extending the previous test cases for wiener or just fixing the failing ones?
Right now no test cases fail anymore

I meant to add missing test cases. Previously there were code reuse between scalar and array noise. After your changes some code was left not covered by tests.

[roflmaostc]

Ah year, I see. I will do them now

[giordano]

Thank you both for the work on this pull requests!!

[roflmaostc]

I think the readthedocs is still deprecated, could you maybe update it using the current master?

[jakubwro]

I think the readthedocs is still deprecated, could you maybe update it using the current master?

Right, now we should use only Documenter version, will change URL