numarray Manual

21.3 Filter functions

The functions described in this section all perform some type of spatial filtering of the the input array: the elements in the output are some function of the values in the neighborhood of the corresponding input element. We refer to this neighborhood of elements as the filter kernel, which is often rectangular in shape but may also have an arbitrary footprint. Many of the functions described below allow you to define the footprint of the kernel, by passing a mask through the footprint parameter. For example a cross shaped kernel can be defined as follows:

>>> footprint = array([[0,1,0],[1,1,1],[0,1,0]])
>>> print footprint
[[0 1 0]
 [1 1 1]
 [0 1 0]]

Usually the origin of the kernel is at the center calculated by dividing the dimensions of the kernel shape by two. For instance, the origin of a one-dimensional kernel of length three is at the second element. Take for example the correlation of a one-dimensional array with a filter of length 3 consisting of ones:

>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> correlate1d(a, [1, 1, 1])
[0 0 1 1 1 0 0]

Sometimes it is convenient to choose a different origin for the kernel. For this reason most functions support the origin parameter which gives the origin of the filter relative to its center. For example:

>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> print correlate1d(a, [1, 1, 1], origin = -1)
[0 1 1 1 0 0 0]

The effect is a shift of the result towards the left. This feature will not be needed very often, but it may be useful especially for filters that have an even size. A good example is the calculation of backward and forward differences:

>>> a = [0, 0, 1, 1, 1, 0, 0]
>>> print correlate1d(a, [-1, 1])              ## backward difference
[ 0  0  1  0  0 -1  0]
>>> print correlate1d(a, [-1, 1], origin = -1) ## forward difference
[ 0  1  0  0 -1  0  0]

We could also have calculated the forward difference as follows:

>>> print correlate1d(a, [0, -1, 1])
[ 0  1  0  0 -1  0  0]

however, using the origin parameter instead of a larger kernel is more efficient. For multi-dimensional kernels origin can be a number, in which case the origin is assumed to be equal along all axes, or a sequence giving the origin along each axis.

Since the output elements are a function of elements in the neighborhood of the input elements, the borders of the array need to be dealt with appropriately by providing the values outside the borders. This is done by assuming that the arrays are extended beyond their boundaries according certain boundary conditions. In the functions described below, the boundary conditions can be selected using the mode parameter which must be a string with the name of the boundary condition. Following boundary conditions are currently supported:

Boundary condition	Description	Example
`"nearest"`	Use the value at the boundary	`[1 2 3]->[1 1 2 3 3]`
`"wrap"`	Periodically replicate the array	`[1 2 3]->[3 1 2 3 1]`
`"reflect"`	Reflect the array at the boundary	`[1 2 3]->[1 1 2 3 3]`
`"constant"`	Use a constant value, default value is 0.0	`[1 2 3]->[0 1 2 3 0]`

The "constant" mode is special since it needs an additional parameter to specify the constant value that should be used.

Note: The easiest way to implement such boundary conditions would be to copy the data to a larger array and extend the data at the borders according to the boundary conditions. For large arrays and large filter kernels, this would be very memory consuming, and the functions described below therefore use a different approach that does not require allocating large temporary buffers.

21.3.1 Correlation and convolution

correlate1d( input, weights, axis=-1, output=None, mode='reflect', cval=0.0, origin=0, output_type=None): The correlate1d function calculates a one-dimensional correlation along the given axis. The lines of the array along the given axis are correlated with the given weights. The weights parameter must be a one-dimensional sequences of numbers.

correlate( input, weights, output=None, mode='reflect', cval=0.0, origin=0, output_type=None): The function correlate implements multi-dimensional correlation of the input array with a given kernel.

convolve1d( input, weights, axis=-1, output=None, mode='reflect', cval=0.0, origin=0, output_type=None): The convolve1d function calculates a one-dimensional convolution along the given axis. The lines of the array along the given axis are convoluted with the given weights. The weights parameter must be a one-dimensional sequences of numbers.
Note: A convolution is essentially a correlation after mirroring the kernel. As a result, the origin parameter behaves differently than in the case of a correlation: the result is shifted in the opposite directions.

convolve( input, weights, output=None, mode='reflect', cval=0.0, origin=0, output_type=None): The function convolve implements multi-dimensional convolution of the input array with a given kernel.
Note: A convolution is essentially a correlation after mirroring the kernel. As a result, the origin parameter behaves differently than in the case of a correlation: the results is shifted in the opposite direction.

21.3.2 Smoothing filters

gaussian_filter1d( input, sigma, axis=-1, order=0, output=None, mode='reflect', cval=0.0, output_type=None): The gaussian_filter1d function implements a one-dimensional Gaussian filter. The standard-deviation of the Gaussian filter is passed through the parameter sigma. Setting order=0 corresponds to convolution with a Gaussian kernel. An order of 1, 2, or 3 corresponds to convolution with the first, second or third derivatives of a Gaussian. Higher order derivatives are not implemented.

gaussian_filter( input, sigma, order=0, output=None, mode='reflect', cval=0.0, output_type=None): The gaussian_filter function implements a multi-dimensional Gaussian filter. The standard-deviations of the Gaussian filter along each axis are passed through the parameter sigma as a sequence or numbers. If sigma is not a sequence but a single number, the standard deviation of the filter is equal along all directions. The order of the filter can be specified separately for each axis. An order of 0 corresponds to convolution with a Gaussian kernel. An order of 1, 2, or 3 corresponds to convolution with the first, second or third derivatives of a Gaussian. Higher order derivatives are not implemented. The order parameter must be a number, to specify the same order for all axes, or a sequence of numbers to specify a different order for each axis.
Note: The multi-dimensional filter is implemented as a sequence of one-dimensional Gaussian filters. The intermediate arrays are stored in the same data type as the output. Therefore, for output types with a lower precision, the results may be imprecise because intermediate results may be stored with insufficient precision. This can be prevented by specifying a more precise output type.

uniform_filter1d( input, size, axis=-1, output=None, mode='reflect', cval=0.0, origin=0, output_type=None): The uniform_filter1d function calculates a one-dimensional uniform filter of the given size along the given axis.

uniform_filter( input, size, output=None, mode='reflect', cval=0.0, origin=0, output_type=None): The uniform_filter implements a multi-dimensional uniform filter. The sizes of the uniform filter are given for each axis as a sequence of integers by the size parameter. If size is not a sequence, but a single number, the sizes along all axis are assumed to be equal.
Note: The multi-dimensional filter is implemented as a sequence of one-dimensional uniform filters. The intermediate arrays are stored in the same data type as the output. Therefore, for output types with a lower precision, the results may be imprecise because intermediate results may be stored with insufficient precision. This can be prevented by specifying a more precise output type.

21.3.3 Filters based on order statistics

minimum_filter1d( input, size, axis=-1, output=None, mode='reflect', cval=0.0, origin=0): The minimum_filter1d function calculates a one-dimensional minimum filter of given size along the given axis.

maximum_filter1d( input, size, axis=-1, output=None, mode='reflect', cval=0.0, origin=0): The maximum_filter1d function calculates a one-dimensional maximum filter of given size along the given axis.

minimum_filter( input, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0): The minimum_filter function calculates a multi-dimensional minimum filter. Either the sizes of a rectangular kernel or the footprint of the kernel must be provided. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint, if provided, must be an array that defines the shape of the kernel by its non-zero elements.

maximum_filter( input, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0): The maximum_filter function calculates a multi-dimensional maximum filter. Either the sizes of a rectangular kernel or the footprint of the kernel must be provided. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint, if provided, must be an array that defines the shape of the kernel by its non-zero elements.

rank_filter( input, rank, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0): The rank_filter function calculates a multi-dimensional rank filter. The rank may be less then zero, i.e., rank=-1 indicates the largest element. Either the sizes of a rectangular kernel or the footprint of the kernel must be provided. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint, if provided, must be an array that defines the shape of the kernel by its non-zero elements.

percentile_filter( input, percentile, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0): The percentile_filter function calculates a multi-dimensional percentile filter. The percentile may be less then zero, i.e., percentile=-20 equals percentile=80. Either the sizes of a rectangular kernel or the footprint of the kernel must be provided. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint, if provided, must be an array that defines the shape of the kernel by its non-zero elements.

median_filter( input, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0): The median_filter function calculates a multi-dimensional median filter. Either the sizes of a rectangular kernel or the footprint of the kernel must be provided. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint if provided, must be an array that defines the shape of the kernel by its non-zero elements.

21.3.4 Derivatives

Derivative filters can be constructed in several ways. The function gaussian_filter1d described in section 21.3.2 can be used to calculate derivatives along a given axis using the order parameter. Other derivative filters are the Prewitt and Sobel filters:

prewitt( input, axis=-1, output=None, mode='reflect', cval=0.0): The prewitt function calculates a derivative along the given axis.

sobel( input, axis=-1, output=None, mode='reflect', cval=0.0): The sobel function calculates a derivative along the given axis.

The Laplace filter is calculated by the sum of the second derivatives along all axes. Thus, different Laplace filters can be constructed using different second derivative functions. Therefore we provide a general function that takes a function argument to calculate the second derivative along a given direction and to construct the Laplace filter:

generic_laplace( input, derivative2, output=None, mode='reflect', cval=0.0, output_type=None, extra_arguments = (), extra_keywords =)

The function generic_laplace calculates a laplace filter using the function passed through derivative2 to calculate second derivatives. The function derivative2 should have the following signature:

derivative2(input, axis, output, mode, cval, *extra_arguments, **extra_keywords)

It should calculate the second derivative along the dimension axis. If output is not None it should use that for the output and return None, otherwise it should return the result. mode, cval have the usual meaning.

The extra_arguments and extra_keywords arguments can be used to pass a tuple of extra arguments and a dictionary of named arguments that are passed to derivative2 at each call.

For example:

>>> def d2(input, axis, output, mode, cval):
...     return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)
... 
>>> a = zeros((5, 5))
>>> a[2, 2] = 1
>>> print generic_laplace(a, d2)
[[ 0  0  0  0  0]
 [ 0  0  1  0  0]
 [ 0  1 -4  1  0]
 [ 0  0  1  0  0]
 [ 0  0  0  0  0]]

To demonstrate the use of the extra_arguments argument we could do:

>>> def d2(input, axis, output, mode, cval, weights):
...     return correlate1d(input, weights, axis, output, mode, cval, 0,)
... 
>>> a = zeros((5, 5))
>>> a[2, 2] = 1
>>> print generic_laplace(a, d2, extra_arguments = ([1, -2, 1],))
[[ 0  0  0  0  0]
 [ 0  0  1  0  0]
 [ 0  1 -4  1  0]
 [ 0  0  1  0  0]
 [ 0  0  0  0  0]]

or:

>>> print generic_laplace(a, d2, extra_keywords = {'weights': [1, -2, 1]})
[[ 0  0  0  0  0]
 [ 0  0  1  0  0]
 [ 0  1 -4  1  0]
 [ 0  0  1  0  0]
 [ 0  0  0  0  0]]

The following two functions are implemented using generic_laplace by providing appropriate functions for the second derivative function:

laplace( input, output=None, mode='reflect', cval=0.0, output_type=None): The function laplace calculates the Laplace using discrete differentiation for the second derivative (i.e. convolution with [1, -2, 1]).

gaussian_laplace( input, sigma, output=None, mode='reflect', cval=0.0, output_type=None): The function gaussian_laplace calculates the Laplace using gaussian_filter to calculate the second derivatives. The standard-deviations of the Gaussian filter along each axis are passed through the parameter sigma as a sequence or numbers. If sigma is not a sequence but a single number, the standard deviation of the filter is equal along all directions.

The gradient magnitude is defined as the square root of the sum of the squares of the gradients in all directions. Similar to the generic Laplace function there is a generic_gradient_magnitude function that calculated the gradient magnitude of an array:

generic_gradient_magnitude( input, derivative, output=None, mode='reflect', cval=0.0, output_type=None, extra_arguments = (), extra_keywords =)

The function generic_gradient_magnitude calculates a gradient magnitude using the function passed through derivative to calculate first derivatives. The function derivative should have the following signature:

derivative(input, axis, output, mode, cval, *extra_arguments, **extra_keywords)

It should calculate the derivative along the dimension axis. If output is not None it should use that for the output and return None, otherwise it should return the result. mode, cval have the usual meaning.

The extra_arguments and extra_keywords arguments can be used to pass a tuple of extra arguments and a dictionary of named arguments that are passed to derivative at each call.

For example, the sobel function fits the required signature:

>>> a = zeros((5, 5))
>>> a[2, 2] = 1
>>> print generic_gradient_magnitude(a, sobel)
[[0 0 0 0 0]
 [0 1 2 1 0]
 [0 2 0 2 0]
 [0 1 2 1 0]
 [0 0 0 0 0]]

See the documentation of generic_laplace for examples of using the extra_arguments and extra_keywords arguments.

The sobel and prewitt functions fit the required signature and can therefore directly be used with generic_gradient_magnitude. The following function implements the gradient magnitude using Gaussian derivatives:

gaussian_gradient_magnitude( input, sigma, output=None, mode='reflect', cval=0.0, output_type=None): The function gaussian_gradient_magnitude calculates the gradient magnitude using gaussian_filter to calculate the first derivatives. The standard-deviations of the Gaussian filter along each axis are passed through the parameter sigma as a sequence or numbers. If sigma is not a sequence but a single number, the standard deviation of the filter is equal along all directions.

21.3.5 Generic filter functions

To implement filter functions, generic functions can be used that accept a callable object that implements the filtering operation. The iteration over the input and output arrays is handled by these generic functions, along with such details as the implementation of the boundary conditions. Only a callable object implementing a callback function that does the actual filtering work must be provided. The callback function can also be written in C and passed using a CObject (see 21.11 for more information).

generic_filter1d( input, function, filter_size, axis=-1, output=None, mode="reflect", cval=0.0, origin=0, output_type=None, extra_arguments = (), extra_keywords =)

The generic_filter1d function implements a generic one-dimensional filter function, where the actual filtering operation must be supplied as a python function (or other callable object). The generic_filter1d function iterates over the lines of an array and calls function at each line. The arguments that are passed to function are one-dimensional arrays of the tFloat64 type. The first contains the values of the current line. It is extended at the beginning end the end, according to the filter_size and origin arguments. The second array should be modified in-place to provide the output values of the line. For example consider a correlation along one dimension:

>>> a = arange(12, shape = (3,4))
>>> print correlate1d(a, [1, 2, 3])
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]

The same operation can be implemented using generic_filter1d as follows:

 
>>> def fnc(iline, oline):
...     oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:]
... 
>>> print generic_filter1d(a, fnc, 3)
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]

Here the origin of the kernel was (by default) assumed to be in the middle of the filter of length 3. Therefore, each input line was extended by one value at the beginning and at the end, before the function was called.

Optionally extra arguments can be defined and passed to the filter function. The extra_arguments and extra_keywords arguments can be used to pass a tuple of extra arguments and/or a dictionary of named arguments that are passed to derivative at each call. For example, we can pass the parameters of our filter as an argument:

 
>>> def fnc(iline, oline, a, b):
...     oline[...] = iline[:-2] + a * iline[1:-1] + b * iline[2:]
... 
>>> print generic_filter1d(a, fnc, 3, extra_arguments = (2, 3))
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]

 
>>> print generic_filter1d(a, fnc, 3, extra_keywords = {'a':2, 'b':3})
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]

generic_filter( input, function, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0, output_type=None, extra_arguments = (), extra_keywords =)

The generic_filter function implements a generic filter function, where the actual filtering operation must be supplied as a python function (or other callable object). The generic_filter function iterates over the array and calls function at each element. The argument of function is a one-dimensional array of the tFloat64 type, that contains the values around the current element that are within the footprint of the filter. The function should return a single value that can be converted to a double precision number. For example consider a correlation:

>>> a = arange(12, shape = (3,4))
>>> print correlate(a, [[1, 0], [0, 3]])
[[ 0  3  7 11]
 [12 15 19 23]
 [28 31 35 39]]

The same operation can be implemented using generic_filter as follows:

 
>>> def fnc(buffer): 
...     return (buffer * array([1, 3])).sum()
... 
>>> print generic_filter(a, fnc, footprint = [[1, 0], [0, 1]])
[[ 0  3  7 11]
 [12 15 19 23]
 [28 31 35 39]]

Here a kernel footprint was specified that contains only two elements. Therefore the filter function receives a buffer of length equal to two, which was multiplied with the proper weights and the result summed.

When calling generic_filter, either the sizes of a rectangular kernel or the footprint of the kernel must be provided. The size parameter, if provided, must be a sequence of sizes or a single number in which case the size of the filter is assumed to be equal along each axis. The footprint, if provided, must be an array that defines the shape of the kernel by its non-zero elements.

 
>>> def fnc(buffer, weights): 
...     weights = asarray(weights)
...     return (buffer * weights).sum()
... 
>>> print generic_filter(a, fnc, footprint = [[1, 0], [0, 1]], extra_arguments = ([1, 3],))
[[ 0  3  7 11]
 [12 15 19 23]
 [28 31 35 39]]

 
>>> print generic_filter(a, fnc, footprint = [[1, 0], [0, 1]], extra_keywords= {'weights': [1, 3]})
[[ 0  3  7 11]
 [12 15 19 23]
 [28 31 35 39]]

These functions iterate over the lines or elements starting at the last axis, i.e. the last index changest the fastest. This order of iteration is garantueed for the case that it is important to adapt the filter dependening on spatial location. Here is an example of using a class that implements the filter and keeps track of the current coordinates while iterating. It performs the same filter operation as described above for generic_filter, but additionally prints the current coordinates:

>>> a = arange(12, shape = (3,4))
>>> 
>>> class fnc_class:
...     def __init__(self, shape):
...         # store the shape:
...         self.shape = shape
...         # initialize the coordinates:
...         self.coordinates = [0] * len(shape)
...         
...     def filter(self, buffer):
...         result = (buffer * array([1, 3])).sum()
...         print self.coordinates
...         # calculate the next coordinates:
...         axes = range(len(self.shape))
...         axes.reverse()
...         for jj in axes:
...             if self.coordinates[jj] < self.shape[jj] - 1:
...                 self.coordinates[jj] += 1
...                 break
...             else:
...                 self.coordinates[jj] = 0
...         return result
... 
>>> fnc = fnc_class(shape = (3,4))
>>> print generic_filter(a, fnc.filter, footprint = [[1, 0], [0, 1]]) 
[0, 0]
[0, 1]
[0, 2]
[0, 3]
[1, 0]
[1, 1]
[1, 2]
[1, 3]
[2, 0]
[2, 1]
[2, 2]
[2, 3]
[[ 0  3  7 11]
 [12 15 19 23]
 [28 31 35 39]]

For the generic_filter1d function the same approach works, except that this function does not iterate over the axis that is being filtered. The example for generic_filte1d then becomes this:

>>> a = arange(12, shape = (3,4))
>>> 
>>> class fnc1d_class:
...     def __init__(self, shape, axis = -1):
...         # store the filter axis:
...         self.axis = axis
...         # store the shape:
...         self.shape = shape
...         # initialize the coordinates:
...         self.coordinates = [0] * len(shape)
...         
...     def filter(self, iline, oline):
...         oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:]
...         print self.coordinates
...         # calculate the next coordinates:
...         axes = range(len(self.shape))
...         # skip the filter axis:
...         del axes[self.axis]
...         axes.reverse()
...         for jj in axes:
...             if self.coordinates[jj] < self.shape[jj] - 1:
...                 self.coordinates[jj] += 1
...                 break
...             else:
...                 self.coordinates[jj] = 0
... 
>>> fnc = fnc1d_class(shape = (3,4))
>>> print generic_filter1d(a, fnc.filter, 3)
[0, 0]
[1, 0]
[2, 0]
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]

Send comments to the NumArray community.