Compute point-spread-function between original and blurred image

Question

Take an image $f$ with some characters on it (below, hjFu3). Let's apply a filter $h$ on it to obtain a second image $g$ where the text is not visible.

Is there a way to compute what kind of filter $h$ we have applied ? The aim being that this filter should be applied to decode a new blurred image that was blurred with the same filter $h$.

---

I know that I should use a deconvolution algorithm. The point is I don't know the PSF function, so a blind algorithm could the solution. However I don't have any success using typical blind deconvolution algorithm of matlab.

I would like to know if we can compute the psf function given that we have the orignal and final image.

---

Data for the problem :

$f=$

$\qquad \quad \large \downarrow PSF ?$

$g=$

Answer 1

We construct an operator based on the assumption that the system is a linear space invariant system. The blurred image is denoted $b$ and the input is denoted $x$.

Since the convolution is commutative, we can write

$ \begin{align} b &= h*x\\ &=x*h \end{align} $

So we can have two equal representations using the matrices $H$ and $X$ corresponding to the integral operations involving $x$ and $h$, i.e.

$ \begin{align} b &= Hx\\ &=Xh \end{align} $

Since we know about $b$ and $x$ exactly, we would like to know about $h$. Our task is thus to construct the operator $X$ based on our knowledge on $x$ and the properties of the imaging system.

Here, I take old code I have lying around to construct the matrix. Forgive me if I forgot some details, but the details for the handling of convolution matrices that have certain structures is nicely explained in the book:

Deblurring Images: Matrices, Spectra, and Filtering

Authors: Per Christian Hansen, James G. Nagy, and Dianne P. O’Leary

Either way, we construct the matrix $X$ based on the following code. I think this corresponds to periodic boundary conditions, but take it with a grain of salt.

X = toeplitz(x,circshift(flipud(x),1));

We then solve for the optimal $h$ by e.g. gradient descent for the objective function:

$ \begin{align} \vert\vert Xh-b\vert\vert^2_2 \end{align} $

We thus know about the filter coefficients in this case. Conversely, based on the computed coefficients $h$ we can construct the matrix $H$ via the same formula as in the above code. Then, given $b$ and the operator $H$, we would like to compute the optimal $x_{est}$ via solving the objective function:

$ \begin{align} \vert\vert Hx_{est}-b\vert\vert^2_2 \end{align} $

Since the original $x$ has sharp edges, using normal gradient descent to solve for the optimal $x_{est}$ would be rather smooth. Since we know that most parts of the image is zero, we choose an additional projection step onto the scaled standard simplex for each iteration (i.e. nonnegative iterative soft-thresholding). The results are shown below.

Note: To plot the filter coefficients $h$ as a PSF, we need to shift them a little bit via (s being the size of the images):

hs = reshape(circshift(h(:),(numel(h)+s(1))/2),s)

Furthermore, I only took a certain color type of your images and subtracted the background such that the background was set to zero.

Answer 2

Assuming that the filter is a linear FIR filter and you know (or choose) the size of the filter, this can be solved the following way.

We have a known image f and an output image g, and we want to determine what filter h takes f->g.

For concreteness, let's assume that the filter is 5x5 filter, although the method applies to any other size. The problem is then to determine the 25 filter coefficients. Each point in the output image is thus a linear combination of 25 input points (known). Thus we have 25 unknowns and a number of equations equal to the number of pixels in g (which should be same as f). There will be uncertainty in pixels near the edge, where you will have to make decisions about what algorithm was used (assume zeros beyond the edge, copy pixels, etc.).

If the images are sufficiently large, this is in general solvable (degenerate cases obviously exist where you do not have enough independent data, such as an image with constant values).

More specifically, you can set up a linear equation to solve

$Ax=b$

where the unknowns $x \in \mathbb{R}^{25}$ are the 25 values in the filter h, $b \in \mathbb{R}^n$ are the values of pixels in the output image arranged in a one-dimensional vector, and $A \in \mathbb{R}^{n \times 25}$ has the values of the input input image arranged in rows. Each row has the 25 pixels that correspond to the input to the filter for outputting that row. That is, the ith row has the 25 pixels that used in computing the ith output pixel (in vector b).

Then you solve that with $x=A\setminus b$ in MATLAB or octave, or numpy.linalg.lstsq in numpy.

You can see how accurately your filter works by comparing $Ax$ with $b$ and perhaps displaying in an image format.