How to use a total variation regularisation with richardson lucy deconvolution ?

Algorithms
1

Hello,
i would like to use a Richardson Lucy deconvolution algorithm with TV regularisation with ITK.
There is already a RL algorithm implemented: https://itk.org/Doxygen44/html/classitk_1_1RichardsonLucyDeconvolutionImageFilter.html
This is ok it work only with a kernel in uint16 but it work.
How to implement a regularisation like in Dey et al 2006 than use for iteration n+1
divergence [ (gradient Image(n) ) / gradient norm Image(n) ) ]
Which filters can i use to calulate the divergence of the gradient divide by gradient norm of the image?
Thank you,
Cyril Jaudet, Phd
CLCC Balcesse, Caen

1 Like
2

I think i find a solution which i am not 100% sure.
divergence [(gradient image (n)/gradient norm (Image(n))]=laplacian(image(n))/normgradient(image(n)
I use itk LaplacianImageFilter() and GradientMagnitudeImageFilter(). It seems ok.
So example in python:
#####initialize filter ###############
laplacian= sitk.LaplacianImageFilter()
normgradient=sitk.GradientMagnitudeImageFilter()
divide=sitk.DivideImageFilter()
multiply=sitk.MultiplyImageFilter()
Substract=sitk.SubtractImageFilter()
Cast=sitk.CastImageFilter()
############TVregulation terme calculation##########
##########terme regularisation##############
image_cast=sitk.Cast(image,sitk.sitkFloat64)
L=laplacian.Execute(image_cast)
NG=normgradient.Execute(image_cast)
NG=sitk.Cast(NG,sitk.sitkFloat64)
i1=divide.Execute(L, NG )
i2=multiply.Execute( i1, alpha) #alphaTV=0.02 for this example
i3=Substract.Execute(1,i2)
i4=divide.Execute(1,i3)
Have a nice day,
Cyril

1 Like
3

That you for sharing your solution.

A tip with using SimpleITK, is that we have overloaded image operators. So there is no need to go through the trouble of using the Filters directly (unless you really want to). You could more simply do something like this:

1.0/(((L/NG)*alpha)-1.0)

For many of the deconvolution filters, the primary input and the kernel input are expected to be the same image type and pixel type. So if you want to use float, for the kernel cast you image to float too.

HTH,
Brad

4

Thank you, it will be much more simple like this.
Best regards,
Cyril