# Estimating Convolution Kernel from Input and Output Images

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.)

• Could you share the code to generate those images? It will make things easier to share working code with you. – Royi Oct 3 '19 at 22:53
• if you could give any input, you simply could use a delta function as input and the output will be your kernel. – Mohammad M Oct 4 '19 at 7:36
• I think you better share the images in a non loss format (PNG). As it changes teh output image and might hurt the process. – David Oct 5 '19 at 13:11

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.

### The Graident with Respect to the Convolution Kernel

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

### Gradient with Respect to Convolution Kernel $$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);


## Numeric Solution

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.

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.

• This will partially work. Pay attention that this solution assumes the output image is a result of Periodic Convolution. In this specific image it won't be a problem as it is all zeros near the edges. Yet for real world images this, if not done correctly with padding, will result in artifacts in the edge. – Royi Oct 4 '19 at 0:53
• @Royi Thanks, I'm aware of these details, but i was trying to keep my answer as short as possible. – Mohammad M Oct 4 '19 at 7:43
• But for many users it will just won't work. As most images has non zero pixels near the edges. – Royi Oct 4 '19 at 7:44