2
$\begingroup$

So, I acquired some images by scanning a radiation source with a square detector like in the following gif.

enter image description here

Where the dashed grid represents reality, the 3x3 square my detector, and the 4x4 my acquired data. Obviously, because I'm moving the detector in steps smaller than its size, my image is convoluted. Assuming my detector has uniform sensitivity, how do I go about deconvolving it?

The point spread function should be simple, but I can't seem to find examples anywhere.

(I can work with Python or Matlab, so feel free to give any examples you want)

The data I'm working with consists of gaussian-ish peaks surrounded by noise like so: enter image description here

$\endgroup$
  • $\begingroup$ Do we know anything about the area in dashed area which isn't marked as blue? $\endgroup$ – Royi Jan 26 at 17:21
  • $\begingroup$ Ok I didn't explain myself correctly. At the center of the image there should be a gaussian-ish bump (from the source). Around it I'm acquiring mostly noise, so the edges don't really matter . $\endgroup$ – Frederico Costa Jan 26 at 17:24
  • $\begingroup$ How big is the grid in the real world problem? $\endgroup$ – Royi Jan 26 at 17:34
  • $\begingroup$ Right now, most of mine are 200x200 $\endgroup$ – Frederico Costa Jan 26 at 17:36
  • 1
    $\begingroup$ I've been using it and it works like a charm. I also noticed that for higher radii, lowering the lambda parameter bears quite impressive results. $\endgroup$ – Frederico Costa Jan 27 at 15:46
2
$\begingroup$

Your model is exactly a Convolution with Uniform Kernel where the output is what is called the Valid Part of the Convolution.

In MATLAB lingo it will be using conv2(mA, mK, 'valid').

So the way to solve it will be using a matrix form of the convolution and solving the linear system of equations.

Let's use the Lenna Image as input (Size was reduced for faster calculations):

enter image description here

We have a uniform kernel for the sensor model.
The output of the convolution with uniform kernel is given by:

enter image description here

The output from the sensor is both blurred and smaller (Less 2 rows and 2 columns as it is 3x3 kernel) just as in your model. This is the model of Valid Convolution.

In Matrix form what we have is:

$$ \boldsymbol{b} = K \boldsymbol{a} $$

Where $ \boldsymbol{b} $ is the column stack vector of the output image, $ \boldsymbol{a} $ is the column stack vector of the input image and $ K $ is the convolution operator (Valid Convolution) in matrix form. In the code it is done in the function CreateConvMtx2D().

So, now all we need is to restore the image by solving the Matrix Equation.
Yet the issue is the equation is Underdetermined System and the matrix has high condition number which suggest not to solve this equation directly.

The solution is to use some kind of regularization of the least squares form of the problem:

$$ \arg \min_{\boldsymbol{a}} \frac{1}{2} {\left\| K \boldsymbol{a} - \boldsymbol{b} \right\|}_{2}^{2} + \lambda r \left( \boldsymbol{a} \right) $$

Where $ r \left( \boldsymbol{a} \right) $ is the regularization term. In the optimal case the regularization should match the prior knowledge on the problem. For instance, in Image Processing we can assume a Piece Wise Smooth / Constant Model which matches the Total Variation regularization.

Since we have no knowledge here, we will use the classic regularization to handle the Condition Number - Tikhonov Regularization:

$$ \arg \min_{\boldsymbol{a}} \frac{1}{2} {\left\| K \boldsymbol{a} - \boldsymbol{b} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| \boldsymbol{a} \right\|}_{2}^{2} = {\left( {K}^{T} K + \lambda I \right)}^{-1} {K}^{T} \boldsymbol{b} $$

The output is given by (For $ \lambda = 0.005 $):

enter image description here

We can see that near the edge we have some artifacts which are due to the fact the system is Underdetermined and we have less equations to describe those pixels.
One can use the $ \lambda $ parameter to balance between how sharp the output is (Yet with artifacts) to how smooth it is, basically governing the level inversion of the system.

I advise playing with the parameter to find the best balance for your case but more than that, find a better regularization. Since the information you're after looks smooth you can use something in that direction.

The full MATLAB code is available on my StackExchange Signal Processing Q63449 GitHub Repository.

Enjoy...

$\endgroup$
  • $\begingroup$ I was looking closer at your code and I noticed that in the final script, in "Solution By linear algebra" you are using mA in the reshape instead of vA. Replacing it, ruins the restored image $\endgroup$ – Frederico Costa Jan 26 at 15:49
  • $\begingroup$ Please notify me when you do (btw thank you very much for your help) $\endgroup$ – Frederico Costa Jan 26 at 17:21
  • 1
    $\begingroup$ I update my answer with some explanation of what I did. I have many answer to similar problems, you should look at them as well. $\endgroup$ – Royi Jan 28 at 7:59
1
$\begingroup$

Below is an attempt to do what you're asking in Python.

First, the dashed item:

Thing to be measured.

Then the sensor. It's uniform,so just comes out as black.

Sensor

Then the output of the sensor (convolve the thing to be measured with the sensor).

Measurements.

Finally, the output of the deconvolution.

Output

Note that the output is not precisely the same as the input, but it's pretty close.


Code Only Below

#Import all libraries we will use
from matplotlib import pyplot as plt, rcParams, rc
from scipy import ndimage
import random
import numpy as np
import cv2
from skimage import color, data, restoration

N = 6

img = np.zeros((N,N),np.uint8)

for x in range(N):
    for y in range(N):
        #We use "0" for black color (do nothing) and "1" for white color (change pixel value to [255,255,255])
        if (x == 2 or x == 3) and (y == 2 or y == 3):
            img[x,y] = 1

cv2.imwrite("img.png",img)
plt.figure()
plt.imshow(img)

M = 4
sensor = np.zeros((M,M), np.uint8)

for x in range(M):
    for y in range(M):
        sensor[x,y] = 1


plt.figure()
plt.imshow(sensor)        
cv2.imwrite("sensor.png",sensor)

measurement = ndimage.convolve(img, sensor, mode='reflect', cval=0.0)
plt.figure();
plt.imshow(measurement)
cv2.imwrite("measurement.png",measurement)


deconvolved_measurement = restoration.richardson_lucy(measurement, sensor, iterations=30)
plt.figure();
plt.imshow(deconvolved_measurement)
cv2.imwrite("deconvolved_measurement.png",deconvolved_measurement)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.