Here’s a post which combines my favourite bits of writing a blog – fairly mathematical, not too simple or difficult to implement, mostly based around pictures, not covered in my undergraduate education, and pretty damn useful in my job. Excited? You should be.
It’s important to be regular, in many aspects of life, but most importantly of all when dealing with ill posed problems. I’ll come onto a specific example later, but for now consider the general matrix problem
where the information contained in is known, and we would like to recover . The matrix is determined by the physical system under analysis, where the physical quantity is mapped under some linear transformation to a quantity which is measured.
Formally, this answer is correct and recovers exactly. However this is the real world and so we don’t know exactly, it is associated with some error . To see how the solution we recover varies as is altered, denote the error in as then
This is generically known as an ‘inverse problem’, where the ‘forward transform’ is known but the inverse transform problematic or impractical to calculate directly. Many problems can be translated into an inverse problem, a notable example being computed tomography which occupied a portion of my PhD.
Here let’s look at the simpler example of image deconvolution. When taking a picture, the ‘true image’ will be smeared across the image detector due to imperfections in the detector. This smearing can be modelled by a convolution
for some point-spread function (PSF) (also known as a kernel). The PSF acts to mix up the image such that a given measured pixel contains contributions from many ‘true’ pixels, and the resolution is worsened.
Lets do this. Below I construct a Gaussian PSF, generate its matrix, and apply it to an image to make a blurred image. I then invert the matrix, and apply that to the blurred image to recover the initial ‘true’ image.
Worked perfectly! What’s the problem then? Let’s add a tiny amount of noise to the blurred image, just 1 count on average (out of an 8-bit gray value range 0-255)
The blurred image looks, by eye, to be identical to the previous figure. The recovered image is utterly corrupted though, to the point where no semblance of image structure is maintained. What happened?
If I calculate the condition number of the matrix , I get . This is a very large condition number, and says that a 1% fractional change in the input image will generate a 100,000,000% change in the output image. This is, putting it mildly, unfortunate. To check I understand what’s going on, I computed recovered images for different levels of added noise and checked how far they were from the true image:
The blue line marks what should be expected given the condition number of , and clearly the recovered images follow this trend. To recover a recognisable image with error < 100, one requires an RMS measurement error of below 1 part in a million. This is very rarely achievable in practice, and demonstrates that this de-blurring technique will generally not work.
What is needed is a bit of regularisation.
The problem we are running into here is that, with some noise added, there can be many possible solutions for the recovered image. The large condition number of the matrix means that while , it is not necessarily true that .
i.e. is close to a solution in a least-squares sense, but which also isn’t too large. The parameter determines the relative weighting of each term, and is known as a regularising parameter. This regularisation scheme is known as Tikhonov regularisation.
Regularised image deconvolution
Let’s finally see how well this technique works to de-blur images. Here Gaussian noise is added of magnitude 1:
If the noise magnitude is lower, the deconvolution is more successful:
If the noise magnitude is larger, the deconvolution isn’t very successful in improving the resolution, but it does lower the noise level:
If, rather than the square of the image, the total variation is penalised, the regularisation process encourages image uniformity (or ‘gradient sparsity’). This a technique widely used in computerised tomography, where it increases the quality of the 3D reconstruction without increasing the radiation dose to the patient.
The code used in this post was written in Matlab, which deals well with large matrices. In particular the kernel matrix is large and sparse – for a 128×128 pixel image it is 16384×16384 and 99.5% sparse.
The scripts can be found here on Github. Keep the images small, even 128×128 takes a few tens of seconds to run.