I'm trying to implement an algorithm of de-convolution in Haskell and couldn't find a simpler one than Richardson Lucy. I looked up at the existing matlab/python implementation but am unable to understand from where to start or how exactly to implement. The library I want to use is https://github.com/lehins/hip. If someone can provide an outline of some implementation or some general idea about the functions with some code snippets, that would be very helpful to me.
Thanks in advance!
The algorithm is actually pretty straightforward. Using the notation on the Wikipedia page for Richardson-Lucy deconvolution, if an underlying image u0 was convolved by a kernel p to produce an observed image d, then you can iterate the function:
deconvolve p d u = u * conv (transpose p) (d / conv p u)
over u with an initial starting estimate (of d, for example) to get a progressively better estimate of u0.
In HIP, the actual one-step deconvolve function might look like:
deconvolve :: Image VS X Double
-> Image VS RGB Double
-> Image VS RGB Double
-> Image VS RGB Double
deconvolve p d u
= u * conv (transpose p) (d / conv p u)
where conv = convolve Edge
and you could use something like this:
let us = iterate (deconvolve p d) d
u10 = us !! 10 -- ten iterations
An example of a full program is:
import Graphics.Image as I
import Graphics.Image.Interface as I
import Prelude as P
blur :: Image VS X Double
blur = blur' / scalar (I.sum blur')
where blur' = fromLists [[0,0,4,3,2]
,[0,1,3,4,3]
,[1,2,3,3,4]
,[0,1,2,1,0]
,[0,0,1,0,0]]
deconvolve :: Image VS X Double
-> Image VS RGB Double
-> Image VS RGB Double
-> Image VS RGB Double
deconvolve p d u
= u * conv (transpose p) (d / conv p u)
where conv = convolve Edge
main :: IO ()
main = do
-- original underlying image
u0 <- readImage' "images/frog.jpg" :: IO (Image VS RGB Double)
-- the kernel
let p = blur
-- blurred imaged
let d = convolve Edge p u0
-- iterative deconvolution
let us = iterate (deconvolve p d) d
u1 = us !! 1 -- one iteration
u2 = us !! 20 -- twenty iterations
let output = makeImage (rows u0, cols u0 * 4)
(\(r,c) ->
let (i, c') = c `quotRem` cols u0
in index ([u0,d,u1,u2] !! i) (r,c'))
:: Image VS RGB Double
writeImage "output.jpg" output
which generates the following image of (left-to-right) the original frog, the blurred frog, a one-fold deconvolution, and a twenty-fold deconvolution.
