# Compute point-spread-function between original and blurred image

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=$$

• I don't that you can do it, in general. Let's consider a filter that averages all the values on the image, what would you do in that case? Sep 20 '19 at 19:01
• I tried some some blind deconvolution algorithms without success Sep 21 '19 at 16:25
• I agree with @nicoguaro that this seems unlikely to be solved without some assumptions or more data. If you had a lot of input-output image samples, you could perhaps find an operator to go from the output blurred image to the input one, but not sure if that's an option. Oct 4 '19 at 5:31
• I assume you have both input image f and output image g and want to compute (or estimate) the filter h. Blind deconvolution (en.wikipedia.org/wiki/Blind_deconvolution ) is for trying to reconstruct the input image f and is a much harder problem. Nov 11 '19 at 0:22

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.

• It works nice. Can we solve analytical the minimization problem : $\min_h \|Xh-b\|^2$ ? It seems it's a least square problem and we should have something like : $h=(X^TX)^{-1}X^Tb$ ? Same for the second optimization problem, so we could have a "analytical solution" (still it's an approximation of the original image) for $x_{est}$ Aug 27 '20 at 13:16
• Even though $X$ is sparse, I'd argue that noise can still be a problem if you go for the least square minimum. Rather iterate via the update $h_{k+1} = h_{k}-\alpha X^T(Xh_{k}-b)$ with $\alpha=\vert\vert X^TX \vert\vert^{-1}$ and stop at a given threshold. You can see if noise makes a problem when the norm of $h$ blows up. Either way, you can try it out.
– Ron
Aug 27 '20 at 13:22

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.