Given an original image and a convolved version of it, I need to calculate the convolution kernel. For example, given:
and
I would like to calculate the convolution kernel that generated the second image from the first (a gaussian in this case).
The images will be similar to the example, but I can change the configuration slightly to make the process easier (change the shape, have a number of shapes etc.)
While it can be done in Frequency Domain it requires delicate handling of the edges (Discrete Frequency Domain assumes Periodic Signals).
Hence I think the best approach is to build this problem as an optimization problem in spatial domain.
Given the objective function:
$$ \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} $$
Where $ h $ is the 2D convolution kernel and $ x $ is the 2D convolution image and $ y $ is a given 2D image.
What would be:
$$ \frac{\mathrm{d} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} }{\mathrm{d} h} $$
The easiest approach would be writing each case using Matrix Form of the convolution.
In this answer we assume the discrete convolution is applied only on valid support (Matching MATLAB's valid parameter for the convolution).
Namely, given $ x \in \mathbb{R}^{m \times n} $ and $ h \in \mathbb{R}^{k \times l} $ then $ h \ast x \in \mathbb{R}^{ \left( m - k + 1 \right) \times \left( n - l + 1 \right) } $. Needless to say $ m \geq k $ and $ n \geq l $ as otherwise the operation isn't well defined. Pay attention that this form of convolution isn't commutative.
The matrix form is given by:
$$ f \left( h \right) = \frac{1}{2} {\left\| X h - y \right\|}_{2}^{2} $$
Where $ X $ is the 2D Convolution Matrix Form of the image. Then:
$$ \frac{\mathrm{d} f \left( h \right) }{\mathrm{d} h} = {X}^{T} \left( X h - y \right) $$
The $ {X}^{T} $ forms a correlation (Versus Convolution) with full support of the operation (Equivalent of the full convolution shape in MATLAB syntax).
Hence we have:
$$ \frac{\mathrm{d} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} }{\mathrm{d} h} = x \star \left( h \ast x - y \right) $$
In MATLAB Code:
conv2(mX(end:-1:1, end:-1:1), (conv2(mX, mH, CONVOLUTION_MODE_VALID) - mY), CONVOLUTION_MODE_VALID);
I implemented vanilla Gradient Descent to solve the problem.
As can be seen the solver is converging to the correct solution.
I'd recommend you'd add Gradient Descent Acceleration step (Like FISTA) to make things much faster.
The full code is available on my StackExchange Signal Processing Q61043 GitHub Repository (Look at the SignalProcessing\Q61043 folder).
Take Fourier Transform of both original image and blurred image, then divide Fourier transform of blurred image by Fourier transform of original image. This will give you the Fourier transform of kernel (except at frequencies with zero value which lead to division by zero). Then take inverse Fourier transform and find the kernel in spatial domain.
