Compute point-spread-function between original and blurred image

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} vertvert Xh-bvertvert^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} vertvert Hx_{est}-bvertvert^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.