Lucy-Richardson Deconvolution: So what is it?

[michael.h.f.wilkinson]

Lucy-Richardson deconvolution is a bit of a hot item at the moment, and I have found it very useful in sharpening solar images, producing clearly superior results to the wavelet sharpening I had been doing before.

After some playing with the algorithm, and implementing it in C, I thought it might be a good idea to make an attempt at a simple explanation of Lucy-Richardson deconvolution (also at the request of bigmakstutov). As all deconvolution methods, it tries to undo the blurring caused by optics and seeing. The effect of blurring is to replace every point source (star) in the image by a little blob called the point-spread function (PSF). While most obvious on stars, the effect is also present in extended sources. Here, the value of each pixel that it should have had without blur, is replaced by a weighted sum its original value and that of the neighbouring pixels. This process of replacing pixels by weighted sums of the neighbourhood is called convolution. The effect of convolution with a small blob is to lose fine detail and sharpness. In theory, we can perform a so-called inverse filter, which undoes the convolution, but this only works in the absence of noise. The key problem is that the noise in the image is not convolved with the PSF. This means that fine detail in the noise is not lost in the image. Inverse filtering boosts fine details indiscriminately, causing the fine detail in the noise to explode.

The LR deconvolution method is different in that it does not directly invert the process of convolution, but tries to estimate an image which when blurred by the PSF is as close as possible to the original input (or observed image). The idea behind this is simply that the observed image is a mix of the desired, unblurred image, convolved with the PSF, plus some Poisson noise (photon and dark current noise). Therefore, the aim is to find an image that, when blurred with the same PSF, gives a result consistent with the observed image. In the following we assume that some estimate of the PSF is known (a Gaussian function often works well).

The process starts with some first estimate (which can be the observed image, or even a flat, 50% grey image). At each iteration it computes a correction factor for each pixel, and multiplies each pixel of the current estimate with this correction factor. The correction factor approaches 1.00 as the process approaches the desired, deconvolved image, because multiplying by 1.00 does not change anything. The set of correction factors for each pixel forms an image in its own right, so computing the correction factors can be done by processing images. The first step is to convolve a copy of the current estimate with the PSF (we blur the current estimate). We then divide the observed image pixel by pixel by this blurred estimate. This yields a new image, which is then blurred by the PSF mirrored in the origin (i.e. you flip it horizontally and vertically). As most PSFs used in this context are symmetric this need often not be done. There are however cases in which you need to do it. After this blurring, you now have the desired correction factor. All that remains is to multiply the correction factor image with the current estimate, and go to the next iteration.

I did wonder about the behaviour of the correction factor in the case that the blurred current estimate contained zero values. This would lead to division-by-zero errors. A short analysis shows that this can be avoided by setting all correction-factor pixels to zero at those locations (or to any arbitrary factor). The reason for this is that these zero values in the blurred estimate occur only where the current estimate is zero itself. Whatever we use to multiply those pixels with, the result remains zero.

As you iterate the effect is to creep closer and closer to the final image. If the process converges (which is not guaranteed) adding more iterations does not sharpen the results further. What can happen is that more iterations become unstable due to noise, or a poor estimate of the PSF. In this case the pixel values can start to oscillate around the correct value. To prevent such oscillations some versions of the method introduce a strength parameter, which I think decreases the step size of the iterations. The idea is that more, smaller steps will not overshoot the desired solution. If the PSF estimate is wrong, this will not help much. Thus, when optimizing the parameters, focus on the PSF is more important than the number of iterations, unlike e.g. iterated unsharp masking. I do often enhance the image with a bit of unsharp masking after LR deconvolution

One key problem in deconvolution is knowledge of the PSF. In most cases it is supposed to be known, and often a Gaussian function with width tuned by the so-called standard deviation sigma can be used. There are methods called blind deconvolution which try to estimate the PSF as well. If you do have knowledge of the PSF, it is best not to use these.

Here are some results on prominences. The first image uses wavelets, the second LR deconvolution:

I will post a GUI-based program which allows batch processing later.

[Knight of Clear Skies]

Thanks for that, my first guess was 60s dance move.

[Mr Spock]

Excellent article.

You can see the difference too. The right image has more detail and looks more natural.

[Pete Presland]

the difference is very noticable, much finer detail visible.

[Davey-T]

Hi Michael,

I have a copy of AIP4WIN that came with the Berry / Burnell book, which has Lucy / Richadson deconvolution as an option, it has several variables that can be set but seems to involve some trial and error. Have you tried it on Solar ?

Dave

[Zakalwe]

Excellent post and description!

[JB80]

Hi Michael,

I have a copy of AIP4WIN that came with the Berry / Burnell book, which has Lucy / Richadson deconvolution as an option, it has several variables that can be set but seems to involve some trial and error. Have you tried it on Solar ?

Dave

Hi Dave,

I have done a brief test with AIP4WIN but nothing thorough, there is a gif in this thread that does some comparing....

http://stargazerslounge.com/topic/227612-decon-madness/?hl=%2Baip4win#entry2453631

Like you say there are a number of variables and I was/am still a bit confused by it so maybe there is a good workflow in it to be discovered.

[callump]

You can also to L-R deconvolution using IRIS.

I am not sure if STSCI originated it, but it came to the fore to correct the early Hubble Space Telescope images before COSTAR was installed.

Callum

[NickK]

Here's my LR GPU implementation with 20 iterations - single sub on the left, two frame stacked and 20 iteration LR on the stacked image on the right:

[riklaunim]

I use LR deconvolution on planetary stacks and it works very well (in older Astra Image I have).

[michael.h.f.wilkinson]

Here's my LR GPU implementation with 20 iterations - single sub on the left, two frame stacked and 20 iteration LR on the stacked image on the right:

Interesting stuff. Is your code CUDA based, or do you use OpenCL? Is it a plugin for Pixinsight?

My main reason for writing my simple piece of code was the need for batch processing. I might take a look at Pixinsight to see what it can do.

[philj]

Excellent write up Michael, Ive been recommended this for lunar in the pastbut never got round to it yet, I shall give it a try now I reckon

[NickK]

Interesting stuff. Is your code CUDA based, or do you use OpenCL? Is it a plugin for Pixinsight?

My main reason for writing my simple piece of code was the need for batch processing. I might take a look at Pixinsight to see what it can do.

The code is OpenCL/Grand Central Dispatch/Objective-C++, the PixInisight was just to put the two side by side. This is the screenshot in total:

If you're thinking of PixInsight - it's probably worth checking the threading model otherwise the entire UI will grind to a halt whilst processing..

[michael.h.f.wilkinson]

Thanks, I assume that you mean Objective C (it being Mac OS-X), or is it a mix of Objective C with C++. I have just added OpenMP pragmas in my code and get a useful speed-up of 4.95 (hyperthreading, I assume) on my quad core Core i5: 100 iterations in 5.51s on a 3.44 Mpixel image (all double precision math). Not in GPU range, so I might add CUDA or OpenCL support later. A higher priority is a simple GUI wrapper

[NickK]

Thanks, I assume that you mean Objective C (it being Mac OS-X), or is it a mix of Objective C with C++. I have just added OpenMP pragmas in my code and get a useful speed-up of 4.95 (hyperthreading, I assume) on my quad core Core i5: 100 iterations in 5.51s on a 3.44 Mpixel image (all double precision math). Not in GPU range, so I might add CUDA or OpenCL support later. A higher priority is a simple GUI wrapper

Impressive is that FFT based?

The next step for me is to move the code to C++/OpenCL because interworking Apple's GCD and OpenCL "sucks" (the polite version).

[michael.h.f.wilkinson]

Impressive is that FFT based?

The next step for me is to move the code to C++/OpenCL because interworking Apple's GCD and OpenCL "sucks" (the polite version).

Mine is based on a recursive IIR filter algorithm for Gaussian convolution by Ted Young and Lucas van Vliet. FFT methods are slower. I might add more kernels later.

[NickK]

Mine is based on a recursive IIR filter algorithm for Gaussian convolution by Ted Young and Lucas van Vliet. FFT methods are slower. I might add more kernels later.

That's what I thought (as FFT is slower)

[michael.h.f.wilkinson]

That's what I thought (as FFT is slower)

FFT is never a good idea for kernels with small support. The nice thing about the recursive filter is that it just needs six multiply and add operations per dimension used, regardless of the value of sigma. When I add other kernels, I will first go for finite-support functions, using a direct implementation of convolution. That is usually faster.

[NickK]

So per pixel in each dimension that sounds about right if the gaussian is pre-calculated (as it would be). My initial figures for the GPU (0.4 sec for 100 iteration IIRC - is more maths that would be required) but I'm more familiar with the FFT so do that first With FFT based alignment anyway it makes it simple to extend without time to research.

There's a couple of options you have for GPU performance enhancement - use of an interpolated local 1D gaussian table means removing one of the global data reads, but that does remove the capability for non-symetric unless you start doing distortions depending on the coordinates and rotation of elongated motion blur, localised varying convolution over the image etc. Also for realtime it's also possible to use less stringent optimisations - possibly trading quality (i.e. some ringing) for the speed optimisation - using boxes etc to simulate rather than have a large extern gaussian data table means it can be coded rather than data driven.

Personally I will keep it data driven for flexibility for the input PSF.. but I would be really tempted for 100fps video for solar to make a fast-and-furious mode using the above optimisations. The FFT based LR that I'm doing - does full image sized PSF at the moment - so there's a further set of optimisations I could do that will make it a lot faster.. but all takes time.. and that's what I've been lacking over the last few weeks.

Double precision is only really just starting to be supported by GPUs.. users will not understand this limitation.. especially apple users. 

I'm not an expert as you do this for as part of your university position - half the fun for me is implementing it from the maths

[NickK]

From the PixInsight website

GPU AccelerationIn preparation: Acceleration with Graphics Processing Units.

We plan on implementing GPU acceleration using the CUDA development platform. GPU acceleration will be integrated at a low level with our PixInsight Class Library C++ development framework, so it will be available to the PixInsight core application—including native implementations of JavaScript runtime objects—and most existing tools. Hopefully the first GPU-accelerated versions of PixInsight will become available during 2014. Stay tuned!

A bit dumb going down the CUDA route.. as it's not going to be portable to Intel and AMD graphics units..

However I suspect they'll have an update at some point with the GPU and fingers crossed they may even do some realtime processing

It may be worth you getting in contact with them and discussing integration...

[michael.h.f.wilkinson]

So per pixel in each dimension that sounds about right if the gaussian is pre-calculated (as it would be). My initial figures for the GPU (0.4 sec for 100 iteration IIRC - is more maths that would be required) but I'm more familiar with the FFT so do that first With FFT based alignment anyway it makes it simple to extend without time to research.

There's a couple of options you have for GPU performance enhancement - use of an interpolated local 1D gaussian table means removing one of the global data reads, but that does remove the capability for non-symetric unless you start doing distortions depending on the coordinates and rotation of elongated motion blur, localised varying convolution over the image etc. Also for realtime it's also possible to use less stringent optimisations - possibly trading quality (i.e. some ringing) for the speed optimisation - using boxes etc to simulate rather than have a large extern gaussian data table means it can be coded rather than data driven.

Personally I will keep it data driven for flexibility for the input PSF.. but I would be really tempted for 100fps video for solar to make a fast-and-furious mode using the above optimisations. The FFT based LR that I'm doing - does full image sized PSF at the moment - so there's a further set of optimisations I could do that will make it a lot faster.. but all takes time.. and that's what I've been lacking over the last few weeks.

Double precision is only really just starting to be supported by GPUs.. users will not understand this limitation.. especially apple users. 

I'm not an expert as you do this for as part of your university position - half the fun for me is implementing it from the maths

From the PixInsight website

A bit dumb going down the CUDA route.. as it's not going to be portable to Intel and AMD graphics units..

However I suspect they'll have an update at some point with the GPU and fingers crossed they may even do some realtime processing

It may be worth you getting in contact with them and discussing integration...

I do not use Gaussian tables or precomputed values. A recursive filter works on a totally different basis than a direct implementation of convolution. Six adds and muls per dimension is impossible for sigma=50 in the latter case. In the recursive car it is possible.

[NickK]

I'd initially assumed you'd have used the gaussian signal as a point of ideal reference for constructing the IIR filter - I've started researching it.. so I may have some more dumb comments/questions

[NickK]

Actually have you thought about using IIR in stacking? Utilising the time portion?

[sharkmelley]

Hi Michael,

I've been reading a bit about Moffat functions recently.  Apparently they are a better approximation to the point spread function of stars than Gaussians.  In other words, the shape of a star caused by seeing related atmospheric turbulence tends to be more like a Moffat rather than a Gaussian.  If this is true, it might be worth experimenting with a Moffat function kernel for your LR deconvolution.  If it's true for stars then it I guess it would also be true for the PSF of solar images.

Mark

[NickK]

PixInight has guassian and moffat. Another option is to perform motion analysis then use that as the PSF.