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

Question

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.

Answer 1

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));
mB = ApplyGaussianBlur(mA, gaussianKernelStd, STD_TO_RADIUS_FACTOR);

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

Answer 2

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.