Added Richardson-Lucy deconvolution

.gitignore

@@ -1,4 +1,4 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
-/doc/_build/
+/doc/build/

README.md

@@ -72,6 +72,29 @@ case a constant noise with that value will be assumed (this is a good
 approximation in the case of
 [white noise](https://en.wikipedia.org/wiki/White_noise)).
 
+### `lucy`
+
+```julia
+lucy(observed, psf[, iterations])
+```
+
+The [Richardson-Lucy deconvolution](https://en.wikipedia.org/wiki/Richardson-Lucy_deconvolution)
+is an iterative method based on Bayesian inference for restoration of signal
+that is convolved with a point spread function.
+
+The `lucy` function can be used to apply the Richardson-Lucy deconvolution
+method to a digital signal. The arguments are:
+
+* `observed`: the observed blurred signal
+* `psf`: the point spread function (the blurring kernel)
+* `iterations` (optional argument): the number of iterations
+
+First two arguments must be arrays, all with the same size, and all of them
+in the time/space domain (they will be converted to the frequency domain
+internally using `fft` function).  Argument `iterations` is an integer number.
+The more iterations is specified the better result should be if the solution
+converges and it is going to converge if PSF is estimated well.
+
 Examples
 --------
 
@@ -117,6 +140,30 @@ polished2 = wiener(blurred_img, img2, noise2, blurring)
 # imshow(polished) # ...the polished image
 # imshow(polished2) # ...the second polished image
 ```
+### Richardson-Lucy deconvolution
+
+Here is an example of use of `lucy` function to perform the Richardson-Lucy
+deconvolution of an image blurred by kernel that models spherical lens aberration.
+
+``` julia
+using Images, TestImages, Deconvolution, FFTW, ZernikePolynomials, ImageView
+
+img = channelview(testimage("cameraman"))
+
+# model of lens aberration
+blurring = evaluateZernike(LinRange(-16,16,512), [12, 4, 0], [1.0, -1.0, 2.0], index=:OSA)
+blurring = fftshift(blurring)
+blurring = blurring ./ sum(blurring)
+
+blurred_img = fft(img) .* fft(blurring) |> ifft |> real
+
+@time restored_img = lucy(blurred_img, blurring, iterations=1000)
+
+imshow(img)
+imshow(blurring)
+imshow(blurred_img)
+imshow(restored_img)
+```
 
 Development
 -----------

doc/src/index.md

@@ -27,7 +27,7 @@ Older versions are also available for Julia 0.4-0.7.
 
 ## Usage
 
-Currently `Deconvolution.jl` provides only one methd, but others will
+Currently `Deconvolution.jl` provides only two methods, but others will
 hopefully come in the future.
 
 ### `wiener` function
@@ -44,7 +44,7 @@ Theoretically, the Wiener deconvolution method requires the knowledge of
 the original signal, the blurring function, and the noise. However,
 these conditions are difficult to met (and, of course, if you know the
 original signal you do not need to perform a deconvolution in order to
-recover the signal itself), but a strenght of the Wiener deconvolution
+recover the signal itself), but a strength of the Wiener deconvolution
 is that it works in the frequency domain, so you only need to know with
 good precision the power spectra of the signal and the noise. In
 addition, most signals of the same class have fairly similar power
@@ -74,6 +74,25 @@ number, in which case a constant noise with that value will be assumed
 (this is a good approximation in the case of [white
 noise](https://en.wikipedia.org/wiki/White_noise)).
 
+### `lucy` function
+
+The [Richardson-Lucy deconvolution](https://en.wikipedia.org/wiki/Richardson-Lucy_deconvolution)
+is an iterative method based on Bayesian inference for restoration of signal
+that is convolved with a point spread function.
+
+The `lucy` function can be used to apply the Richardson-Lucy deconvolution
+method to a digital signal. The arguments are:
+
+-   `observed`: the digital signal
+-   `psf`: the point spread function
+-   `iterations` (optional argument): the number of iterations
+
+First two arguments must be arrays, all with the same size, and all of them
+in the time/space domain (they will be converted to the frequency domain
+internally using `fft` function). Argument `iterations` is an integer number.
+The more iterations is specified the better result should be if the solution
+converges and it is going to converge if PSF is estimated well.
+
 ## Examples
 
 ### Wiener deconvolution
@@ -192,6 +211,37 @@ view(polished2) # ...the second polished image
 
 ![image](wiener-cameraman.jpg)
 
+### Richardson-Lucy deconvolution
+
+#### Blurred image
+
+Here is an example of use of `lucy` function to perform the Richardson-Lucy
+deconvolution of an image convolved with point spread function that models lens
+aberration.
+
+``` {.sourceCode .julia}
+    using Images, TestImages, Deconvolution, FFTW, ZernikePolynomials, ImageView
+
+    img = channelview(testimage("cameraman"))
+
+    # model of lens aberration
+    blurring = evaluateZernike(LinRange(-16,16,512), [12, 4, 0], [1.0, -1.0, 2.0], index=:OSA)
+    blurring = fftshift(blurring)
+    blurring = blurring ./ sum(blurring)
+
+    blurred_img = fft(img) .* fft(blurring) |> ifft |> real
+
+    @time restored_img_200 = lucy(blurred_img, blurring, iterations=200)
+    @time restored_img_2000 = lucy(blurred_img, blurring, iterations=2000)
+
+    imshow(img)
+    imshow(blurred_img)
+    imshow(restored_img_200)
+    imshow(restored_img_2000)
+```
+
+![image](lucy-cameraman.jpg)
+
 ## Development
 
 The package is developed at

doc/src/lucy-cameraman.jl

@@ -0,0 +1,23 @@
+using Images, TestImages
+using Deconvolution
+using FFTW
+using ZernikePolynomials #https://github.com/rdoelman/ZernikePolynomials.jl
+
+img = channelview(testimage("cameraman"))
+
+blurring = evaluateZernike(LinRange(-16,16,512), [12, 4, 0], [1.0, -1.0, 2.0], index=:OSA)
+blurring = fftshift(blurring)
+blurring = blurring ./ sum(blurring)
+
+blurred_img = fft(img) .* fft(blurring) |> ifft |> real
+
+@time result_200_iterations = lucy(blurred_img, blurring, iterations=200)
+@time result_2000_iterations = lucy(blurred_img, blurring, iterations=2000)
+
+(N, M) = size(blurred_img)
+
+out = ((vcat(hcat(img, ones(M, 20), blurred_img),
+             ones(20, 2M+20),
+             hcat(result_200_iterations, ones(M, 20), result_2000_iterations))))
+
+save("lucy-cameraman.jpg", Images.clamp01nan.(out))

src/Deconvolution.jl

@@ -16,5 +16,6 @@ module Deconvolution
 using FFTW
 
 include("wiener.jl")
+include("lucy.jl")
 
 end # module

src/lucy.jl

@@ -0,0 +1,49 @@
+### lucy.jl ---  Lucy deconvolution
+#
+# Copyright (C) 2019 Jakub Wronowski.
+#
+# Maintainer: Mosè Giordano <mose AT gnu DOT org>
+# Keywords: richadson, lucy, deconvolution, signal processing
+#
+# This file is a part of Deconvolution.jl.
+#
+# License is MIT "Expat".
+#
+### Commentary:
+#
+# This file provides the `lucy' function to perform Richardson-Lucy deconvolution.
+#
+### Code:
+
+export lucy
+
+function lucy(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 = observed
+    for i in 1:iterations
+        estimated = lucystep(estimated)
+    end
+
+    return estimated
+end
+
+"""
+    lucy(observed, psf[, iterations])
+
+Returns the Richardson-Lucy deconvolution of `observed`, using the point spread
+function `psf`.
+
+All arguments must be arrays in the time/space domain and all of the same size,
+they will be converted into the frequency domain internally using the `fft`
+function.
+"""
+lucy

test/runtests.jl

@@ -28,3 +28,25 @@ end
          -1.6455701546601862,-0.4926325667795202,-1.0323127293402468,
          -0.8380228416032021,1.8214493186523952,0.7447075584686396,1.0566048110361044]
 end
+
+##### Richardson-Lucy deconvolution
+
+@testset "Richardson-Lucy deconvolution tests" begin
+    k=-0.001
+    blurring_ft  = exp.(-k * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
+    blurred_s_ft = fftshift(blurring_ft) .* fft(s)
+    blurred_s    = real(ifft(blurred_s_ft))
+    blurring     = ifft(fftshift(blurring_ft))
+
+    estimated = lucy(blurred_s, blurring, iterations=1)
+
+    @test sum(abs.(s .- estimated)) < sum(abs.(s .- blurred_s))
+
+    k =-0.0015 # should improve also when blur params are not perfectly estimated 
+    blurring_ft  = exp.(-k * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
+    blurring     = ifft(fftshift(blurring_ft))
+
+    estimated = lucy(blurred_s, blurring, iterations=1)
+
+    @test sum(abs.(s .- estimated)) < sum(abs.(s .- blurred_s))
+end