# Finding the Best Gaussian Smoothing Kernel to Minimize the Discrepancy Between Two Images

Suppose we have two grayscale images, $A$ and $B$. $A$ and $B$ very strongly resemble each other, such that the mean of the absolute difference $\lvert A - B\rvert$ is fairly small. Suppose further that $B$ actually appears to be a blurred version of $A$, although it was not necessarily obtained through a blurring process. I'd like to find the ideal Gaussian kernel $g$ to smooth $A$ with to maximize its correspondence to $B$, in other words, select $g$ s.t.

$$\underset{g}{\arg\min}\sum{\lvert(g*A) - B\rvert}$$

where $g$, $A$, and $B$ are 2D arrays and $*$ is the convolution operator.

• I don't understand the setup, is $|A-B|$ "fairly small" before or after the blurring of B? – A_A May 11 '18 at 10:02
• Working on this... – Royi May 11 '18 at 14:28
• @A_A the key is that $B$ wasn't constructed by blurring A. It was actually generated by a neural network; but it very strongly resembles a slightly blurry version of $A$. I'd like to estimate a gaussian kernel that most closely replicates the effect. – eriophora May 11 '18 at 16:55
• can you tell us why you want to estimate the kernel (by the way, assuming it's zero-mean and circular, that kernel only has one relevant parameter: its variance)? Are you maybe just after a "degree of blurring", which might be easy to capture as a frequency domain property? – Marcus Müller May 12 '18 at 0:58

This is closely related to Blind Deconvolution.
The only difference is we limit our self to a very specific type of blur kernels.

The nice thing about the Gaussian Kernel is being defined by single parameter - The Standard Deviation of the kernel.
The less nice thing is the connection isn't linear.

## Optimization Problem

Let's define a classic non linear model for this problem:

$$\arg \min_{\sigma} \frac{1}{2} {\left\| A \left( \sigma \right) x - b \right\|}_{2}^{2}$$

Where $A \left( \sigma \right)$ is the convolution matrix generated by a Gaussian Kernel parameterized by $\sigma$, $x$ is the original image ($A$ in your question) in a vector shape, and $b$ is the blurred image ($B$ in your question).
This is a classical Non Linear Least Squares problem which can be solved by MATLAB using lsqnonlin().

## Code Sample

This is the main part of the code:

mA = rand([numRows, numCols]);

gaussianKernelStd = kernelStdLowerBound + ((kernelStdUpperBound - kernelStdLowerBound) * rand(1));

% Objective Functions
hObjFun = @(kernelStd) reshape(ApplyGaussianBlur(mA, kernelStd, STD_TO_RADIUS_FACTOR) - mB, [numPx, 1]);
estKernelEst    = lsqnonlin(hObjFun, initKernelStd, kernelStdLowerBound, kernelStdUpperBound, sSolverOptions);


## Results As can be seen, the estimation is almost perfect.
On large images it might take time (Using some tricks of the Gaussian Filter in the Fourier Domain one could do that there with major speed up), but still, it is not free.

The full code is available on my StackExchange Signal Processing Q49121 GitHub Repository.

• This could be precisely what I'm looking for. The method should also be robust to the fact that $B$ wasn't generated by gaussian smoothing in the first place but merely looks like it was. I'll write out some code right now and see if it works, thanks so much! – eriophora May 11 '18 at 16:57
• One question, though, when you say the original image "in vector shape" is this equivalent to a reshaping (in the matlab [or my precisely, python--which is my language of choice]) of the image into a column vector? – eriophora May 11 '18 at 16:58
• I am accepting this as the answer, although I didn't end up using it exactly. I implemented it in Tensorflow and took an iterative approach to optimize the Gaussian's sigma over many many images. However, framing the problem as what is effectively a traditional optimization task was the key I needed. Thanks, Royi! – eriophora May 11 '18 at 22:14
• Don't use the pre existing functions. You will basically not learn anything if you do. – mathreadler May 13 '18 at 9:14

Considering that the convolution is a multiplication in the Fourier domain, this problem can be converted to a very simple fitting problem.

$B$ is a blurred version of $A$. Thus we have $\hat{B} = \hat{A} G$, with $\hat B$ the DFT of $B$. This can be rewritten into an expression $G = \hat{B} / \hat{A}$, with a point-wise division. Next, simply fit a Gaussian to $G$. You can do this while ignoring the pixels where $\hat{A}$ is close to zero.

• This is general classic Deconvolution. Yet it doesn't take into account The Gaussian Blur model @eriophora asked for. – Royi May 12 '18 at 6:39
• @Royi: Why do you su that? After fitting a Gaussian you have the parameters to a Gaussian kernel that best approximates the linear filter that was applied. – Cris Luengo May 12 '18 at 13:21
• The division could create peaks in places different then the origin. Then fitting Gaussian spatially to the filter might put its origin on a different spot. Yet, since blurring doesn't change the DC and Gaussian filter is LPF it means the fit won't make sense. You can enforce origin to be in a specific place yet I think in practice it will yield poor results. – Royi May 12 '18 at 13:25
• By the way, you could do that in Spatial domain. Solving linear system of blurring (Forcing Circulant Matrix) and then do the same Gaussian Kernel fit. I still think it won't work very well. – Royi May 12 '18 at 13:27
• @Royi: in the spatial domain it’s a more complex problem. The circulant matrix is a lot larger than the image itself. I don’t know if the DFT method would work better than yours or not, but it’s easy to try out, and well worth the effort, because if it works well, it’s very simple and efficient. Of course, you fit a Gaussian centered at 0 and with unit weight at zero, so you only need to fit two or three parameters, depending on how much you want to constrain it. – Cris Luengo May 12 '18 at 13:40