What Is the Difference between Difference of Gaussian, Laplace of Gaussian, and Mexican Hat Wavelet?

Question

There are three techniques used in CV that seem very similar to each other, but with subtle differences:

• Laplacian of Gaussian: $\nabla^2\left[g(x,y,t)\ast f(x,y)\right]$
• Difference of Gaussians: $ \left[g_1(x,y,t)\ast f(x,y)\right] - \left[g_2(x,y,t)\ast f(x,y)\right]$
• Convolution with Ricker wavelet: $\textrm{Ricker}(x,y,t)\ast f(x,y)$

As I understand it currently: DoG is an approximation of LoG. Both are used in blob detection, and both perform essentially as band-pass filters. Convolution with a Mexican Hat/Ricker wavelet seems to achieve very much the same effect.

I've applied all three techniques to a pulse signal (with requisite scaling to get the magnitudes similar) and the results are pretty darn close. In fact, LoG and Ricker look nearly identical. The only real difference I noticed is with DoG, I had 2 free parameters to tune ($\sigma_1$ and $\sigma_1$) vs 1 for LoG and Ricker. I also found the wavelet was the easiest/fastest, as it could be done with a single convolution (done via multiplication in Fourier space with FT of a kernel) vs 2 for DoG, and a convolution plus a Laplacian for LoG.

• What are the comparative advantages/disadvantages of each technique?
• Are there different use-cases where one outshines the other?

I also have the intuitive thought that on discrete samples, LoG and Ricker degenerate to the same operation, since $\nabla^2$ can be implemented as the kernel $$\begin{bmatrix}-1,& 2,& -1\end{bmatrix}\quad\text{or}\quad\begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \end{bmatrix}\quad\text{for 2D images}$$.

Applying that operation to a gaussian gives rise to the Ricker/Hat wavelet. Furthermore, since LoG and DoG are related to the heat diffusion equation, I reckon that I could get both to match with enough parameter fiddling.

(I'm still getting my feet wet with this stuff to feel free to correct/clarify any of this!)

Answer 1

Laplace of Gaussian

The Laplace of Gaussian (LoG) of image $f$ can be written as

$$ \nabla^2 (f * g) = f * \nabla^2 g $$

with $g$ the Gaussian kernel and $*$ the convolution. That is, the Laplace of the image smoothed by a Gaussian kernel is identical to the image convolved with the Laplace of the Gaussian kernel. This convolution can be further expanded, in the 2D case, as

$$ f * \nabla^2 g = f * \left(\frac{\partial^2}{\partial x^2}g+\frac{\partial^2}{\partial y^2}g\right) = f * \frac{\partial^2}{\partial x^2}g + f * \frac{\partial^2}{\partial y^2}g $$

Thus, it is possible to compute it as the addition of two convolutions of the input image with second derivatives of the Gaussian kernel (in 3D this is 3 convolutions, etc.). This is interesting because the Gaussian kernel is separable, as are its derivatives. That is,

$$ f(x,y) * g(x,y) = f(x,y) * \left( g(x) * g(y) \right) = \left( f(x,y) * g(x) \right) * g(y) $$

meaning that instead of a 2D convolution, we can compute the same thing using two 1D convolutions. This saves a lot of computations. For the smallest thinkable Gaussian kernel you'd have 5 samples along each dimension. A 2D convolution requires 25 multiplications and additions, two 1D convolutions require 10. The larger the kernel, or the more dimensions in the image, the more significant these computational savings are.

Thus, the LoG can be computed using four 1D convolutions. The LoG kernel itself, though, is not separable.

There is an approximation where the image is first convolved with a Gaussian kernel and then $\nabla^2$ is implemented using finite differences, leading to the 3x3 kernel with -4 in the middle and 1 in its four edge neighbors.

The Ricker wavelet or Mexican hat operator are identical to the LoG, up to scaling and normalization.

Difference of Gaussians

The difference of Gaussians (DoG) of image $f$ can be written as

$$ f * g_{(1)} - f * g_{(2)} = f * (g_{(1)} - g_{(2)}) $$

So, just as with the LoG, the DoG can be seen as a single non-separable 2D convolution or the sum (difference in this case) of two separable convolutions. Seeing it this way, it looks like there is no computational advantage to using the DoG over the LoG. However, the DoG is a tunable band-pass filter, the LoG is not tunable in that same way, and should be seen as the derivative operator it is. The DoG also appears naturally in the scale-space setting, where the image is filtered at many scales (Gaussians with different sigmas), the difference between subsequent scales is a DoG.

There is an approximation to the DoG kernel that is separable, reducing computational cost by half, though that approximation is not isotropic, leading to rotational dependence of the filter.

I once showed (for myself) the equivalence of the LoG and DoG, for a DoG where the difference in sigma between the two Gaussian kernels is infinitesimally small (up to scaling). I don't have records of this, but it was not difficult to show.

Other forms of computing these filters

Laurent's answer mentions recursive filtering, and the OP mentions computation in the Fourier domain. These concepts apply to both the LoG and the DoG.

The Gaussian and its derivatives can be computed using a causal and anti-causal IIR filter. So all 1D convolutions mentioned above can be applied in constant time w.r.t. the sigma. Note that this is only efficient for larger sigmas.

Likewise, any convolution can be computed in the Fourier domain, so both the DoG and LoG 2D kernels can be transformed to the Fourier domain (or rather computed there) and applied by multiplication.

In conclusion

There are no significant differences in the computational complexity of these two approaches. I have yet to find a good reason to approximate the LoG using the DoG.

Answer 2

The Ricker wavelet, the (isotropic) Marr wavelet, the Mexican hat or the Laplacian of Gaussians belong to be the same concept: continuous admissible wavelets (satisfying certain conditions). Traditionally, the Ricker wavelet is the 1D version. The Marr wavelet or the Mexican hat are names given in the context of 2D image decompositions, you can consider for instance Section 2.2 of A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity, Signal Processing, 2011, L. Jacques et al. The Laplacian of Gaussian is the multidimensional generalization.

However, in practice, people accepts different types of discretizations, at different levels.

I tend to believe (unless given more details) that a $3\times 3$ discrete gradient kernel applied to a Gaussian is not the original Ricker, but a simplification, that explains subtle differences in the graph. I am interested in references. Indeed, you can have at least two natural discretizations of the $3\times 3$ Laplacian operator (4- and 8-neighbors):

$$\begin{pmatrix} 0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \\ \end{pmatrix}$$

or $$\begin{pmatrix} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \\ \end{pmatrix}$$ There are also other approximations, for instance with a $5\times 5$ kernel, or other avatars of the Laplacian/Laplacian of Gaussian.

With a proper choice in their variance ratios $\sigma_1$ and $\sigma_2$ (usually around 1.6), a difference of Gaussians provides a nice separable approximation to the LoG (see for instance Fast Almost-Gaussian Filtering, P. Kovesi). Those Gaussians can in turn be approximated by recursive approximate Gaussians.

But other ratios have been used, in some Laplacian pyramids for instance, that turn DoG more into more generic bandpass filters or edge detectors.

Last reference: Image Matching Using Generalized Scale-Space Interest Points, T. Lindeberg, 2015.

Answer 3

Let's see how DoG approximates LoG for the 2D case (for an image, e.g.). By derivative theorem of convolution (by associativity and commutativity),

$$\nabla^2[f(x, y) \ast G(x, y)] = \nabla^2 G(x, y) \ast I (x, y)$$

where the 2D Gaussian (separable) pdf $G(.)$ and the corresponding first and second order derivatives are defined as follows:

$$G(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{y^2}{2\sigma^2}}$$

$$G_x(x,y)=\frac{1}{2\pi\sigma^4}(-x)e^{-\frac{x^2+y^2}{2\sigma^2}}$$

$$G_y(x,y)=\frac{1}{2\pi\sigma^4}(-y)e^{-\frac{x^2+y^2}{2\sigma^2}}$$

$$G_{xx}(x,y)=\frac{1}{2\pi\sigma^4}e^{-\frac{x^2+y^2}{2\sigma^2}}\left(\frac{x^2}{\sigma^2}-1\right)$$ $$G_{xy}(x,y)=\frac{1}{2\pi\sigma^4}e^{-\frac{x^2+y^2}{2\sigma^2}}\left(\frac{xy}{\sigma^2}\right)$$ $$G_{yy}(x,y)=\frac{1}{2\pi\sigma^4}e^{-\frac{x^2+y^2}{2\sigma^2}}\left(\frac{y^2}{\sigma^2}-1\right)$$

Hence, $$\nabla^2 G(x,y)=G_{xx}+G_{yy}=\frac{1}{2\pi\sigma^4}e^{-\frac{x^2+y^2}{2\sigma^2}}\left(\frac{x^2+y^2}{\sigma^2}-2\right)$$.

Now, let's prove that the Laplacian of a Gaussian (LoG) is (a scaled version of) the derivative with respect to $2\sigma^2$ of a Gaussian. Observe that

$$\ln G(x, y, \sigma^2) = -\frac{x^2+y^2}{2\sigma^2} - \ln{\pi} - \ln{2\sigma^2}$$

$$\implies \frac{\partial }{\partial (2\sigma^2)}\ln G = \frac{1}{G} \frac{\partial G}{\partial (2\sigma^2)} = \frac{x^2+y^2}{\left(2\sigma^2\right)^2} - \frac{1}{2\sigma^2}$$

$$\begin{align}\implies \frac{\partial }{\partial \left(2\sigma^2\right)} G(x,y,\sigma^2) = & \frac{G(x,y,\sigma^2)}{4\sigma^2}\left(\frac{x^2+y^2}{\sigma^2} - 2\right) \\ = & \frac{1}{4}.\frac{1}{2\pi\sigma^4}e^{-\frac{x^2+y^2} {2\sigma^2}}\left(\frac{x^2+y^2}{\sigma^2} - 2\right) \\ = & \frac{1}{4}LoG(x, y, \sigma^2) \end{align}$$

$$LoG(x, y, \sigma^2) \propto \frac{\partial }{\partial \left(\sigma^2\right)}G(x, y, \sigma^2) = \lim\limits_{h \to 0}\frac{G(x, y, \sigma^2+h)-G(x, y, \sigma^2)}{h}$$

That is, it is (upto scaling) the limit of one Gaussian minus a just smaller Gaussian. For this reason, we can approximate it as the difference of two Gaussians at two different scales (s.d.s):

$$Lo𝐺_𝜎=\nabla^2 G_𝜎≈𝐺_{𝜎_1}−𝐺_{𝜎_2}=DoG_{𝜎_1, 𝜎_2} \text{ , with } 𝜎_1>𝜎_2$$

The following figure shows how DoG closely approximates LoG: