Gaussian filter

From SubSurfWiki
Jump to navigation Jump to search
Kernels for various linear filters. (a) Mean with uniform weights and square support. (b) Mean with uniform weights and circular support. (c) Mean with triangular weights and circular support; the mean is computed with n = 25, the total of the weight factors. (d) A discrete approximation to a Gaussian with σ = 1, for which the mean is computed with n = 273.

There are many other linear smoothing filters, but the most important one is the Gaussian filter, which applies weights according to the Gaussian distribution (d in the figure)[1].

The key parameter is σ, which controls the extent of the kernel and consequently the degree of smoothing (and how long the algorithm takes to execute). More aggressive than the mean filter, the Gaussian filter deals with random noise more effectively (Figures 1d and 2d). But it still simply mixes the noise into the result and smooths indiscriminately across edges. Again, it is imperative to remove spikes before applying this filter.

Optimizing

It's important that the Gaussian kernel is not too small — it must be large enough to 'fade out' without sharp edges. Edges have large bandwidth, so they introduce aliasing artifacts into the result after convolution.[2]

The optimal value for σ is between about 0.8 and 1.0.[2]

Implementation in Python

We can make a Gaussian kernel in Python:[3]

def gaussian_kernel(size, size_y=None):
    size = int(size)
    if not size_y:
        size_y = size
    else:
        size_y = int(size_y)
    x, y = numpy.mgrid[-size:size+1, -size_y:size_y+1]
    g = numpy.exp(-(x**2/float(size)+y**2/float(size_y)))
    return g / g.sum()

# Make the Gaussian by calling the function
gaussian_kernel_array = gaussian_kernel(5)

plt.imshow(gaussian_kernel_array, cmap=plt.get_cmap('jet'), interpolation='nearest')
plt.colorbar()
plt.show()

Gaussian kernel.png

But we don't really need to — we can just use the SciPy signal processing library:

import scipy.ndimage
gaussian = scipy.ndimage.filters.gaussian_filter(noisy_horizon, 3.0)

plt.imshow(gaussian, aspect=0.5, vmin=vmin, vmax=vmax)
plt.show()

See also

External links

References

  1. Hall, M (2007). Smooth operator: smoothing seismic horizons and attributes. The Leading Edge 26 (1), January 2007, p16-20. doi:10.1190/1.2431821
  2. 2.0 2.1 Blog post by Cris Luengo.
  3. Filtering horizons — IPython Notebook on smoothing horizons.