Laplacian Variance Filtering Techniques

• Laplacian Variance Filtering is a set of techniques that compute the variance of a Laplacian response to quantify high-frequency activity and image sharpness.
• It leverages both global and local window statistics, using isotropic, anisotropic, and quarter Laplacian operators to adapt to various processing needs.
• Applications range from blur detection and noise estimation to edge-preserving smoothing and optimization in machine learning.

Laplacian variance filtering denotes a family of signal and image processing procedures in which a Laplacian or Laplacian-like response is first computed and then aggregated, weighted, or regularized through a variance-related quantity. In its canonical form, the method computes the variance of a discrete Laplacian response and uses that scalar as a blur or sharpness indicator. In broader usage, the same principle extends to local-window variance maps, multiscale Laplacian-band energies, Laplacian-based noise variance estimators, and Laplacian-derived variance reduction operators in optimization. The common motif is that second-order structure, or a Laplacian surrogate, is treated as a proxy for high-frequency activity, with variance quantifying its spatial, spectral, or stochastic dispersion (So et al., 2024, Sumiya et al., 2022, Khmou et al., 2013, Curtin et al., 2018).

1. Canonical definition and statistical interpretation

The standard formulation begins with a grayscale image $I(x,y)$, a discrete Laplacian operator $\mathcal{L}$, and the Laplacian response

$L(x,y) = (\mathcal{L} * I)(x,y).$

If the image contains $N$ pixels, the mean Laplacian value is

$\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$

and the Variance of Laplacian is

$$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$

The solar image quality study describes this procedure as “convolving the image with a 3×3 Laplacian operator and calculating the variance, which is the squared standard deviation,” and uses the resulting scalar as a per-image sharpness score. In that setting, higher variance corresponds to sharper imagery because narrow, high-contrast edges produce a wider distribution of Laplacian values, whereas blur attenuates those responses and narrows the distribution (So et al., 2024).

A local variant replaces the global average over all pixels by a windowed computation. The quarter Laplacian study states the standard local formulation as

$$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$

and interprets it as a local measure for focus, texture, or edge activity. In this form, Laplacian variance filtering is not restricted to global blur detection; it becomes a spatially varying activity measure that can guide smoothing, enhancement, or segmentation decisions (Gong et al., 2021).

The statistical meaning of the construction is tied to the Laplacian as a second-derivative operator. Smooth regions tend to yield $L(x,y)\approx 0$, while edges, fine textures, and other high-spatial-frequency structures produce large positive or negative responses. The variance therefore measures the spread of these responses rather than their sign. A common misconception is that Laplacian variance is intrinsically a per-pixel feature map; in the solar study it is used purely as a scalar per image, whereas in edge-aware and multiscale settings it is naturally extended to local maps or bandwise energies (So et al., 2024).

2. Operator design: isotropic, anisotropic, and quarter Laplacians

The behavior of Laplacian variance filtering depends strongly on the underlying Laplacian operator. The quarter Laplacian filter was introduced as a corner- and edge-preserving alternative to the standard isotropic discrete Laplacian. The reference isotropic kernel in that work is

$$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$

which is selected because it has the most isotropic Fourier spectrum among the listed standard kernels. The same study explicitly notes that such an isotropic Laplacian “does not preserve edges”: in diffusion, it blurs across them (Gong et al., 2021).

The quarter Laplacian replaces the full $3\times 3$ support by four oriented quarter windows. Each oriented kernel uses only a $\mathcal{L}$0 effective support: three neighbors with equal weight $\mathcal{L}$1 and a central weight $\mathcal{L}$2. For an image $\mathcal{L}$3, the per-orientation responses are

$\mathcal{L}$4

and the filter selects the orientation with the smallest-magnitude response,

$\mathcal{L}$5

This selection makes the overall operator nonlinear and anisotropic. The stated intuition is that the chosen quarter window is the one that least crosses a strong edge, so the filter diffuses along local structure rather than across it. The paper attributes the edge- and corner-preserving behavior specifically to this argmin selection and to the smaller support region (Gong et al., 2021).

For Laplacian variance filtering, the quarter Laplacian provides an operator replacement rather than a different statistic. The same variance formula can be applied to $\mathcal{L}$6, yielding a quarter-Laplacian variance

$\mathcal{L}$7

with

$\mathcal{L}$8

The paper characterizes this substitution as producing a response that is more local, orientation-adaptive, and edge-aware than the standard Laplacian. It also emphasizes that the quarter Laplacian can be implemented efficiently via a sliding $\mathcal{L}$9 box sum, so the runtime is almost the same as that of the original Laplacian filter despite the nonlinear selection step (Gong et al., 2021).

3. Multiscale and nonlinear extensions

Multiscale Laplacian variance filtering arises naturally when Laplacian responses are organized into pyramids or scale spaces. In fast multi-layer Laplacian enhancement, an edge-aware smoother $L(x,y) = (\mathcal{L} * I)(x,y).$0 is used to build multiple smoothed versions $L(x,y) = (\mathcal{L} * I)(x,y).$1, and the corresponding detail bands are formed as differences such as $L(x,y) = (\mathcal{L} * I)(x,y).$2 and $L(x,y) = (\mathcal{L} * I)(x,y).$3. The paper writes the filtered reconstruction as a sum of a base layer and detail layers, with optional nonlinear mappings $L(x,y) = (\mathcal{L} * I)(x,y).$4, and shows that the result can be rewritten as a linear combination of Laplacians derived from affinity matrices. It further introduces a normalization-free approximation

$L(x,y) = (\mathcal{L} * I)(x,y).$5

so that the filter itself is interpreted as identity plus a graph-Laplacian term. A structure mask derived from the sum of affinities suppresses detail amplification in flat or noisy regions and preserves it in structured regions, which makes the method a practical multiscale Laplacian-band gating scheme (Talebi et al., 2016).

Local Laplacian filtering extends this multiscale viewpoint by making the remapping intensity-aware and edge-aware at each scale. In the Gaussian Fourier Pyramid method for Local Laplacian Filter, the remap has the form

$L(x,y) = (\mathcal{L} * I)(x,y).$6

and the expensive remap evaluation is approximated by Fourier series expansion. The resulting Fourier LLF retains the local Laplacian design while requiring only precomputed Gaussian Fourier pyramids. The paper reports that Fourier LLF has higher accuracy than fast LLF for the same number of pyramids and that it supports parameter-adaptive filtering because the pyramids depend only on the input image, whereas remap dependence remains in scalar coefficients. The same study also presents a pixel-by-pixel enhancement scheme in which “the degree of enhancement is determined according to the variance of the pixels in a local window; the low variance is a small enhancement,” directly coupling local variance to Laplacian-band manipulation (Sumiya et al., 2022).

A more radical generalization replaces the linear Laplacian by the nonlinear $L(x,y) = (\mathcal{L} * I)(x,y).$7-Laplacian

$L(x,y) = (\mathcal{L} * I)(x,y).$8

and defines filtering through the associated flow $L(x,y) = (\mathcal{L} * I)(x,y).$9. The $N$0-Laplacian spectra framework rigorously defines decomposition, reconstruction, filtering, and spectrum for such nonlinear flows, including finite-time extinction for $N$1 and a spectral function $N$2 constructed by fractional derivatives in time. Filtering is then expressed as

$N$3

This suggests a nonlinear analogue of multiscale Laplacian variance filtering in which the relevant quantity is spectral energy over extinction times rather than variance of a linear second-derivative response (Cohen et al., 2019).

Another nonlinear extension is the fourth-order Laplacian-based denoising model

$N$4

with evolution equation

$N$5

The paper states that this nonlinear filter preserves discontinuities while filtering out noise, produces piecewise linear restored data, and avoids the staircase effect commonly observed with total variation denoising methods. Although the paper does not use the phrase “variance of Laplacian,” it explicitly treats $N$6 as the operative second-order quantity within a robust curvature penalty, which places it within the broader Laplacian variance filtering family (Hoang, 2024).

4. Noise estimation, manifold filtering, and stochastic variance reduction

Laplacian variance filtering also appears as an estimator rather than as a direct image transform. In 3D video denoising with a Kalman filter, a Laplacian-based fast AWGN estimator is used to infer the noise standard deviation: $N$7 with

$N$8

The estimated variance then sets the Kalman covariances, with the study choosing $N$9 after PSNR-based tuning. The paper reports a root mean square error of $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$0 for the estimator on a $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$1 test image across 20 noise levels, and gives video denoising examples in which PSNR increases from $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$2 dB to $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$3 dB for one frame of “gmissa.avi” and from $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$4 dB to $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$5 dB for one frame of “gflower” (Khmou et al., 2013).

On manifolds with known symmetry, the graph Laplacian itself has an explicit variance term in its convergence rate. The $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$6-invariant graph Laplacian constructs affinities not only between sampled points $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$7 and $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$8, but between all pairs generated by the action of a compact unitary group $\mu_L = \frac{1}{N}\sum_{x,y} L(x,y),$9. For uniformly sampled data, the paper shows

$$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$0

where the standard graph Laplacian has a variance term depending on $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$1 rather than $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$2. The same study demonstrates denoising on noisy samples from an $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$3-invariant manifold: for $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$4, the noisy MSE is $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$5, standard GL denoised MSE is $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$6, and $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$7-GL denoised MSE is $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$8; for $$\sigma^2_{\mathrm{Lap}} = \frac{1}{N}\sum_{x,y}\big(L(x,y)-\mu_L\big)^2.$$9, the corresponding values are $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$0, $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$1, and $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$2 (Rosen et al., 2023).

In optimization, the object being filtered is the stochastic gradient rather than an image. Laplacian Smoothing Gradient Descent replaces the raw gradient by

$$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$3

where $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$4 is the one-dimensional discrete Laplacian with periodic boundary conditions. Because $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$5 is circulant, its inverse is computed efficiently by FFT. The paper states that the method can dramatically reduce the variance, allow a larger step size, and improve generalization accuracy. It also proves that the smoothing preserves the mean, increases the smallest component, and decreases the largest component of the gradient vector. Empirically, on logistic regression for MNIST, the maximum variance for batch size $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$6 drops from $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$7 with $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$8 to $$\sigma^2_L(x,y) = \frac{1}{|W|}\sum_{(u,v)\in W}\left(L(u,v)-\mu_L(x,y)\right)^2, \qquad \mu_L(x,y)=\frac{1}{|W|}\sum_{(u,v)\in W}L(u,v),$$9 with $L(x,y)\approx 0$0, and to $L(x,y)\approx 0$1 with $L(x,y)\approx 0$2 (Curtin et al., 2018).

5. Established applications

The most direct application is blur and sharpness assessment. In solar image quality assessment, the Variance of Laplacian method was applied to over 65,000 individual solar images acquired on more than 160 days. After image selection and filtering, 65,172 out of 288,756 images remained, covering 168 days. The study defined a daily-normalized sharpness measure

$L(x,y)\approx 0$3

and found a clear daily trend: DVL is relatively low in the morning, peaks around local noon, and decreases again toward sunset. When plotted against solar altitude, DVL increases with altitude and plateaus above roughly $L(x,y)\approx 0$4. The paper then models the dependence on zenith angle $L(x,y)\approx 0$5 through Image Sharpness Functions analogous to airmass formulae, including

$L(x,y)\approx 0$6

with $L(x,y)\approx 0$7, $L(x,y)\approx 0$8, $L(x,y)\approx 0$9, $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$0, and

$$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$1

with $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$2, $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$3, $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$4, $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$5, $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$6, $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$7. The study further proposes the same methodology for the Moon, bright planets, and defocused stars (So et al., 2024).

Edge-aware image processing is another established domain. The quarter Laplacian filter is used inside a diffusion-type smoothing pipeline

$$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$8

with the paper observing that after about 10 iterations the result stabilizes and fixing $$K_\mathrm{iso}= \begin{pmatrix} \frac{1}{12} & \frac{1}{6} & \frac{1}{12}\ \frac{1}{6} & -1 & \frac{1}{6}\ \frac{1}{12} & \frac{1}{6} & \frac{1}{12} \end{pmatrix},$$9 in its applications. The reported tasks are edge-preserving image smoothing, texture enhancement through an unsharp-mask-style decomposition into smoothed base and detail layer, and low-light image enhancement through multiscale illumination decomposition. The paper compares the smoothing application to Domain Transform filter, Guided Filter, Relative Total Variation, and Static/Dynamic filter, and attributes fewer halo or gradient reversal artifacts to the quarter Laplacian formulation (Gong et al., 2021).

Multiscale detail editing and tone manipulation constitute a further application area. Fast multi-layer Laplacian enhancement is described as supporting image smoothing, sharpening, and tone manipulation with low computational and memory requirements, while the structure mask is used specifically to avoid noise or artifact magnification and detail loss. The local Laplacian literature addresses the same application space from a different direction: LLF and Fourier LLF aim to preserve edge-aware decompositions while reducing halos and supporting parameter-adaptive filtering (Talebi et al., 2016, Sumiya et al., 2022).

6. Assumptions, sensitivities, and recurrent limitations

The canonical VL statistic is simple, but its interpretation is conditional. The solar study emphasizes that the Laplacian is a high-pass operator, so image noise and compression artifacts can increase high-frequency content and thus increase VL. It also notes that absolute VL values depend on aperture and focal length, pixel size relative to the seeing disk, wavelength band, camera response, and JPG quality. For this reason the paper introduces DVL as a relative, within-day normalization and states that any Image Sharpness Function is unique to each location and sensitive to sky conditions (So et al., 2024).

Edge-aware operator design introduces a different set of trade-offs. The quarter Laplacian achieves locality and edge preservation through a $3\times 3$0 effective support and an $3\times 3$1 selection rule, but the same paper makes clear that the operator is still local and linear per orientation, that strong pixel-scale noise may require additional robustification, and that the orientation selection is based solely on minimizing $3\times 3$2. It also notes that no explicit multiscale or anisotropic diffusion PDE is developed; the method simply iterates the discrete update (Gong et al., 2021).

Multiscale methods incur approximation and representation constraints. Fast LLF depends on sampled remap intensities, so sparse sampling reduces fidelity; Fourier LLF addresses this but is formulated for grayscale, with color handled by applying LLF to $3\times 3$3 in YUV and leaving $3\times 3$4 unchanged. The paper states that fully multichannel distance metrics would require high-dimensional Gaussian kernel approximations. This suggests that the variance term in Laplacian variance filtering is often not the limiting difficulty; the more delicate issue is how faithfully the underlying Laplacian-like transform represents edge structure, scale interactions, and color dependencies (Sumiya et al., 2022).

A broader misconception is that Laplacian variance filtering is a single algorithm. The literature represented here shows instead a spectrum of formulations: scalar global sharpness measures, local window variances, multiscale Laplacian-band decompositions, Laplacian-based noise variance estimators, nonlinear $3\times 3$5-Laplacian spectra, fourth-order curvature flows, group-invariant graph Laplacians, and Laplacian smoothing operators for stochastic gradients. What unifies them is not a fixed implementation, but the use of Laplacian responses as second-order structure surrogates and the use of variance, energy, or variance-like dispersion to control inference, filtering, or regularization (Cohen et al., 2019, Rosen et al., 2023, Curtin et al., 2018).