Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Smoothing a 2D array along only one axis

I wonder if anyone could help me extend the smoothing example in the SciPy cookbook to a 2D problem.

This script works great for smoothing a 1D function, and they also give code for a 2D smoothing in both axis (i.e. blurring an image).

However, I'd like to apply this function to a 2D dataset, but only along one axis (x direction). I could do this in a loop, by inspecting each slice in y, applying the 1D convolution, then rebuilding the array. But this seems poor coding technique.

Therefore, I wonder how to do it in 2D? I imagine I need to make a 2D kernel with weights changing along one direction only, but I'm not sure how to do this, or which convolve function to use (numpy.convolve, scipy.signal.convolve, scipy.ndimage.filters.convolve1d etc.)

like image 401
IanRoberts Avatar asked Jan 07 '23 17:01

IanRoberts


2 Answers

Perhaps the simplest option is to use one of the 1D filters in scipy.ndimage.filters:

from scipy import ndimage
from scipy.misc import lena

img = lena()

# a uniform (boxcar) filter with a width of 50
boxcar = ndimage.uniform_filter1d(img, 50, 1)

# a Gaussian filter with a standard deviation of 10
gauss = ndimage.gaussian_filter1d(img, 10, 1)

You could also use the non-1D versions of the filters like this: ndimage.gaussian_filter(img, (0, 10)) (i.e. set the filter width to 0 for the axes that you don't want to smooth along).

To smooth with an arbitrary kernel, you could use scipy.ndimage.convolve1d:

import numpy as np

kern = np.hanning(50)   # a Hanning window with width 50
kern /= kern.sum()      # normalize the kernel weights to sum to 1

hanning = ndimage.convolve1d(img, kern, 1)

Here's what the various outputs look like:

from matplotlib import pyplot as plt

fig, ax = plt.subplots(2, 2, figsize=(8, 8))
ax[0, 0].imshow(img)
ax[0, 0].set_title('Original')
ax[0, 1].imshow(boxcar)
ax[0, 1].set_title('Boxcar filter (width = 50)')
ax[1, 0].imshow(gauss)
ax[1, 0].set_title(r'Gaussian filter ($\sigma$ = 10)')
ax[1, 1].imshow(hanning)
ax[1, 1].set_title(r'Hanning window (width = 50)')

for aa in ax.flat:
    aa.set_axis_off()

fig.tight_layout()
plt.show()

enter image description here

like image 151
ali_m Avatar answered Jan 10 '23 23:01

ali_m


Introduction

It seems like you ought to be able to do ny=1 to perform 1D convolution of a 2D image, but that reveals that the cookbooks functions are actually using length 2 * n + 1 kernels. This made me think that you could use ny=0, but that creates a 0/0 in the definition of the kernel. So, no luck there either. :( Based on this, I'm of the opinion that cookbook is not great for your purposes, so I've provided an alternate method of doing what you're asking.

To perform smoothing of a 2D array by convolution along 1 dimension only, all you need to do is make a 2D array (kernel) that has a shape of 1 along one of the dimensions,

import numpy as np

kern = np.ones((11, 1))  # This will smooth along columns

And normalize it so that it sums to one,

kern /= kern.sum()

Then convolve it with your signal,

import scipy.signal as signal

X, Y = np.mgrid[-70:70, -70:70]
Z = np.cos((X**2+Y**2)/200.) + np.random.normal(size=X.shape)

Z_smooth = signal.convolve(Z, kern)

This ought to give you something like this, Boxcar window.

A Nicer Kernel

Above I've used the 'boxcar' kernel (constant values), which many people consider to be somewhat crude. People often prefer to use sharper or smoother filters (e.g. 'hanning', or 'gaussian' as in the cookbook).

kern_hanning = signal.hanning(11)[:, None]
kern_hanning /= kern_hanning.sum()
kern_gauss7 = signal.gaussian(11, 7)[:, None]
kern_gauss7 /= kern_gauss7.sum()
kern_gauss3 = signal.gaussian(11, 3)[:, None]
kern_gauss3 /= kern_gauss3.sum()

These different windows look like this,

Kernels

After applying these filters you'll get something like,

Hanning Kernel. Gauss3 Kernel Gauss7 Kernel

Note that the 'Gauss7' kernel is pretty much the same as the boxcar, and so it produces very similar results in the output. The hanning window, on the other hand, is much finer and so it produces a sharper filter of the data (much less smearing across rings).

like image 38
farenorth Avatar answered Jan 10 '23 22:01

farenorth