\chapter{Multi-dimensional image processing}
\label{cha:ndimage}

%begin{latexonly}
\makeatletter
\py@reset
\makeatother
%end{latexonly}
\declaremodule[numarray.ndimage]{extension}{numarray.nd_image}
\moduleauthor{Peter Verveer}{verveer@users.sourceforge.net}
\modulesynopsis{Multidimensional image analysis functions}

\begin{quote}
  The \module{numarray.nd\_image} module provides functions for
  multidimensional image analysis.
\end{quote}

\section{Introduction}
Image processing and analysis are generally seen as operations on
two-dimensional arrays of values. There are however a number of fields
where images of higher dimensionality must be analyzed. Good examples
of these are medical imaging and biological imaging. \module{numarray}
is suited very well for this type of applications due its inherent
multi-dimensional nature. The \module{numarray.nd\_image} packages
provides a number of general image processing and analysis functions
that are designed to operate with arrays of arbitrary dimensionality.
The packages currently includes functions for linear and non-linear
filtering, binary morphology, B-spline interpolation, and object
measurements.

\section{Properties shared by all functions}
All functions share some common properties. Notably, all functions allow the 
specification of an output array with the \var{output} argument. With this 
argument you can specify an array that will be changed in-place with the 
result with the operation. In this case the result is not returned. Usually, 
using the \var{output} argument is more efficient, since an existing array 
is used to store the result.

The type of arrays returned is dependent on the type of operation, but it is  in most cases equal to the type of the input. If, however, the \var{output} argument is used, the type of the result is equal to the type of the specified output argument. If no output argument is given, it is still possible to specify what the result of the output should be. This is done by simply assigning the desired numarray type object to the output argument. For example:
\begin{verbatim}
>>> print correlate(arange(10), [1, 2.5])
[ 0  2  6  9 13 16 20 23 27 30]
>>> print correlate(arange(10), [1, 2.5], output = Float64)
[  0.    2.5   6.    9.5  13.   16.5  20.   23.5  27.   30.5]                   
\end{verbatim}
\note{In previous versions of \module{numarray.nd\_image}, some functions accepted the \var{output_type} argument to achieve the same effect. This argument is still supported, but its use will generate an deprecation warning. In a future version all instances of this argument will be removed. The preferred way to specify an output type, is by using the \var{output} argument, either by specifying an output array of the desired type, or by specifying the type of the output that is to be returned.}
%
\section{Filter functions}
\label{sec:ndimage: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 \var{footprint} 
parameter. For example a cross shaped kernel can be defined as follows:
\begin{verbatim}
>>> footprint = array([[0,1,0],[1,1,1],[0,1,0]])
>>> print footprint
[[0 1 0]
 [1 1 1]
 [0 1 0]]
\end{verbatim}
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:
\begin{verbatim}
>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> correlate1d(a, [1, 1, 1])
[0 0 1 1 1 0 0]
\end{verbatim}
Sometimes it is convenient to choose a different origin for the kernel. For
this reason most functions support the \var{origin} parameter which gives 
the origin of the filter relative to its center. For example:
\begin{verbatim}
>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> print correlate1d(a, [1, 1, 1], origin = -1)
[0 1 1 1 0 0 0]
\end{verbatim}
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:
\begin{verbatim}
>>> 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]
\end{verbatim}
We could also have calculated the  forward difference as follows:
\begin{verbatim}
>>> print correlate1d(a, [0, -1, 1])
[ 0  1  0  0 -1  0  0]
\end{verbatim}
however, using the origin parameter instead of a larger kernel is more
efficient. For multi-dimensional kernels \var{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 \var{mode} parameter which must be a 
string with the name of the boundary condition.  Following boundary 
conditions are currently supported:
\begin{tableiii}{l|l|l}{constant}{Boundary condition}{Description}{Example}
  \lineiii{"nearest"}{Use the value at the boundary}
  {\constant{[1 2 3]->[1 1 2 3 3]}}
  \lineiii{"wrap"}{Periodically replicate the array}
  {\constant{[1 2 3]->[3 1 2 3 1]}}
  \lineiii{"reflect"}{Reflect the array at the boundary}
  {\constant{[1 2 3]->[1 1 2 3 3]}}
  \lineiii{"constant"}{Use a constant value, default value is 0.0}
  {\constant{[1 2 3]->[0 1 2 3 0]}}
\end{tableiii}
The \constant{"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.}

\subsection{Correlation and convolution}

\begin{funcdesc}{correlate1d}{input, weights, axis=-1, output=None, 
    mode='reflect', cval=0.0, origin=0, output_type=None} The
  \function{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 \var{weights}. The \var{weights} parameter must 
  be a one-dimensional sequences of numbers.
\end{funcdesc}

\begin{funcdesc}{correlate}{input, weights, output=None, mode='reflect', 
    cval=0.0, origin=0, output_type=None} The function \function{correlate}
  implements multi-dimensional correlation of the input array with a given
  kernel.
\end{funcdesc}

\begin{funcdesc}{convolve1d}{input, weights, axis=-1, output=None, 
    mode='reflect', cval=0.0, origin=0, output_type=None} The
  \function{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 \var{weights}. The \var{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 \var{origin} parameter behaves differently than 
  in the case of a correlation: the result is shifted in the opposite 
  directions.}
\end{funcdesc}

\begin{funcdesc}{convolve}{input, weights, output=None, mode='reflect', 
    cval=0.0, origin=0, output_type=None} The function \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 \var{origin} parameter behaves differently than 
  in the case of a correlation: the results is shifted in the opposite 
  direction.}
\end{funcdesc}

\subsection{Smoothing filters}
\label{sec:ndimage:filter-functions:smoothing}

\begin{funcdesc}{gaussian_filter1d}{input, sigma, axis=-1, order=0, 
    output=None, mode='reflect', cval=0.0, output_type=None} The
  \function{gaussian_filter1d} function implements a one-dimensional 
  Gaussian
  filter. The standard-deviation of the Gaussian filter is passed through 
  the parameter \var{sigma}. Setting \var{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.
\end{funcdesc}

\begin{funcdesc}{gaussian_filter}{input, sigma, order=0, output=None, 
  mode='reflect', cval=0.0, output_type=None} The 
  \function{gaussian_filter} function implements a multi-dimensional 
  Gaussian filter. The standard-deviations of the Gaussian filter along 
  each axis are passed through the parameter \var{sigma} as a sequence or 
  numbers.  If \var{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 \var{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.}
\end{funcdesc}

\begin{funcdesc}{uniform_filter1d}{input, size, axis=-1, output=None,
    mode='reflect', cval=0.0, origin=0, output_type=None} The
  \function{uniform_filter1d} function calculates a one-dimensional uniform
  filter of the given \var{size} along the given axis.
\end{funcdesc}

\begin{funcdesc}{uniform_filter}{input, size, output=None, mode='reflect', 
    cval=0.0, origin=0, output_type=None} The \function{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 
  \var{size} parameter. If \var{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.}
  \end{funcdesc}

\subsection{Filters based on order statistics}

\begin{funcdesc}{minimum_filter1d}{input, size, axis=-1, output=None, 
    mode='reflect', cval=0.0, origin=0} The \function{minimum_filter1d}
  function calculates a one-dimensional minimum filter of given \var{size}
  along the given axis.
\end{funcdesc}

\begin{funcdesc}{maximum_filter1d}{input, size, axis=-1, output=None, 
    mode='reflect', cval=0.0, origin=0} The \function{maximum_filter1d}
  function calculates a one-dimensional maximum filter of given \var{size}
  along the given axis.
\end{funcdesc}

\begin{funcdesc}{minimum_filter}{input,  size=None, footprint=None, 
    output=None, mode='reflect', cval=0.0, origin=0} The
  \function{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 \var{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 \var{footprint}, if provided, 
  must be an array that defines the shape of the kernel by its non-zero 
  elements.
\end{funcdesc}

\begin{funcdesc}{maximum_filter}{input,  size=None, footprint=None, 
    output=None, mode='reflect', cval=0.0, origin=0} The
  \function{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 \var{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 \var{footprint}, if provided, 
  must be an array that defines the shape of the kernel by its non-zero 
  elements.
\end{funcdesc}

\begin{funcdesc}{rank_filter}{input, rank, size=None, footprint=None,
  output=None, mode='reflect', cval=0.0, origin=0} The 
  \function{rank_filter}
  function calculates a multi-dimensional rank filter.  The \var{rank} may 
  be less then zero, i.e., \var{rank}=-1 indicates the largest element. 
  Either the sizes of a rectangular kernel or the footprint of the kernel 
  must be provided. The \var{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 \var{footprint}, if provided, 
  must be an array that defines the shape of the kernel by its non-zero 
  elements.
\end{funcdesc}

\begin{funcdesc}{percentile_filter}{input, percentile, size=None, 
  footprint=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{percentile_filter} function calculates a multi-dimensional
  percentile filter.  The \var{percentile} may be less then zero, i.e.,
  \var{percentile}=-20 equals \var{percentile}=80. Either the sizes of a 
  rectangular kernel or the footprint of the kernel must be provided. The 
  \var{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 \var{footprint}, if provided, must be an array that 
  defines the shape of the kernel by its non-zero elements.
\end{funcdesc}

\begin{funcdesc}{median_filter}{input, size=None, footprint=None, 
  output=None, mode='reflect', cval=0.0, origin=0} The 
  \function{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 \var{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 \var{footprint} if provided, 
  must be an array that defines the shape of the kernel by its non-zero 
  elements.
\end{funcdesc}

\subsection{Derivatives}

Derivative filters can be constructed in several ways. The function
\function{gaussian_filter1d} described in section
\ref{sec:ndimage:filter-functions:smoothing} can be used to calculate
derivatives along a given axis using the \var{order} parameter. Other
derivative filters are the Prewitt and Sobel filters:

\begin{funcdesc}{prewitt}{input, axis=-1, output=None, mode='reflect', 
  cval=0.0} The \function{prewitt} function calculates a derivative along 
  the given axis.
\end{funcdesc}

\begin{funcdesc}{sobel}{input, axis=-1, output=None, mode='reflect', 
  cval=0.0} The \function{sobel} function calculates a derivative along 
  the given axis.
\end{funcdesc}

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:

\begin{funcdesc}{generic_laplace}{input, derivative2, output=None,
  mode='reflect', cval=0.0, output_type=None, extra_arguments = (), 
  extra_keywords = {}} The function 
  \function{generic_laplace} calculates a laplace filter using the
  function passed through \var{derivative2} to calculate second 
  derivatives. The function \function{derivative2} should have the 
  following signature:

  \function{derivative2(input, axis, output, mode, cval, *extra_arguments, **extra_keywords)}
  
  It should calculate the second derivative along the dimension \var{axis}. 
  If \var{output} is not \constant{None} it should use that for the output 
  and return \constant{None}, otherwise it should return the result. 
  \var{mode}, \var{cval} have the usual meaning.
  
  The \var{extra_arguments} and \var{extra_keywords} arguments can be used 
  to pass a tuple of extra arguments and a dictionary of named 
  arguments that are passed to \function{derivative2} at each call.

  For example:
\begin{verbatim}
>>> 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]]
\end{verbatim}
To demonstrate the use of the \var{extra_arguments} argument we could do:
\begin{verbatim}
>>> 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]]
\end{verbatim}
or:
\begin{verbatim}
>>> 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]]
\end{verbatim}
\end{funcdesc}

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

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

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

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 \function{generic_gradient_magnitude} function that 
calculated the gradient magnitude of an array:

\begin{funcdesc}{generic_gradient_magnitude}{input, derivative,
  output=None, mode='reflect', cval=0.0, output_type=None, 
  extra_arguments = (), extra_keywords = {}} The 
  function \function{generic_gradient_magnitude} calculates a gradient 
  magnitude using the function passed through \var{derivative} to calculate 
  first derivatives. The function \function{derivative} should have the 
  following signature:

  \function{derivative(input, axis, output, mode, cval, *extra_arguments, **extra_keywords)}
  
  It should calculate the derivative along the dimension \var{axis}. If
  \var{output} is not \constant{None} it should use that for the output and
  return \constant{None}, otherwise it should return the result. 
  \var{mode}, \var{cval} have the usual meaning.
  
  The \var{extra_arguments} and \var{extra_keywords} arguments can be used 
  to pass a tuple of extra arguments and a dictionary of named 
  arguments that are passed to \function{derivative} at each call.

  For example, the \function{sobel} function fits the required signature:
\begin{verbatim}
>>> 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]]
\end{verbatim}
See the documentation of \function{generic_laplace} for examples of using the \var{extra_arguments} and \var{extra_keywords} arguments.
\end{funcdesc}

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

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

\subsection{Generic filter functions}
\label{sec:ndimage:genericfilters}
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 \ref{sec:ndimage:ccallbacks} for more 
information).

\begin{funcdesc}{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 \function{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 
  \function{generic_filter1d} function iterates over the lines of an array 
  and calls \var{function} at each line. The arguments that are passed to 
  \var{function} are one-dimensional arrays of the \constant{tFloat64} 
  type. The first contains the values of the current line. It is extended 
  at the beginning end the end, according to the \var{filter_size} and 
  \var{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:

\begin{verbatim}
>>> a = arange(12, shape = (3,4))
>>> print correlate1d(a, [1, 2, 3])
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]
\end{verbatim}
The same operation can be implemented using \function{generic_filter1d} as 
follows:
\begin{verbatim} 
>>> 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]]
\end{verbatim}
  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 \var{extra_arguments} and \var{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:
\begin{verbatim} 
>>> 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]]
\end{verbatim}
or
\begin{verbatim} 
>>> print generic_filter1d(a, fnc, 3, extra_keywords = {'a':2, 'b':3})
[[ 3  8 14 17]
 [27 32 38 41]
 [51 56 62 65]]
\end{verbatim}
\end{funcdesc}

\begin{funcdesc}{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 \function{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 \function{generic_filter} 
  function iterates over the array and calls \var{function} at each 
  element. The argument of \var{function} is a one-dimensional array of the 
  \constant{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:

\begin{verbatim}
>>> 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]]
\end{verbatim}
The same operation can be implemented using \function{generic_filter} as 
follows:
\begin{verbatim} 
>>> 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]]
\end{verbatim}
  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 \function{generic_filter}, either the sizes of a rectangular 
  kernel or the footprint of the kernel must be provided. The \var{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 \var{footprint}, if provided, must be an array that defines the shape 
  of the kernel by its non-zero elements.

  Optionally extra arguments can be defined and passed to the filter 
  function. The \var{extra_arguments} and \var{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:
\begin{verbatim} 
>>> 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]]
\end{verbatim}
or
\begin{verbatim} 
>>> 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]]
\end{verbatim}
\end{funcdesc}

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 
\function{generic_filter}, but additionally prints the current coordinates:
\begin{verbatim}
>>> 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]]
\end{verbatim}

For the \function{generic_filter1d} function the same approach works, except that this function does not iterate over the axis that is being filtered. The example for \function{generic_filte1d} then becomes this:
\begin{verbatim}
>>> 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]]
\end{verbatim}

\section{Fourier domain filters}
The functions described in this section perform filtering operations in the
Fourier domain. Thus, the input array of such a function should be 
compatible with an inverse Fourier transform function, such as the 
functions from the \module{numarray.fft} module. We therefore have to deal 
with arrays that may be the result of a real or a complex Fourier 
transform. In the case of a real Fourier transform only half of the of the 
symmetric complex transform is stored. Additionally, it needs to be known 
what the length of the axis was that was transformed by the real fft.  The 
functions described here provide a parameter \var{n} that in the case of a 
real transform must be equal to the length of the real transform axis 
before transformation. If this parameter is less than zero, it is assumed 
that the input array was the result of a complex Fourier transform. The 
parameter \var{axis} can be used to indicate along which axis the real 
transform was executed.

\begin{funcdesc}{fourier_shift}{input, shift, n=-1, axis=-1, output=None}
  The \function{fourier_shift} function multiplies the input array with the
  multi-dimensional Fourier transform of a shift operation for the given 
  shift. The \var{shift} parameter is a sequences of shifts for each 
  dimension, or a single value for all dimensions.
\end{funcdesc}

\begin{funcdesc}{fourier_gaussian}{input, sigma, n=-1, axis=-1, output=None}
  The \function{fourier_gaussian} function multiplies the input array with 
  the multi-dimensional Fourier transform of a Gaussian filter with given
  standard-deviations \var{sigma}. The \var{sigma} parameter is a sequences 
  of values for each dimension, or a single value for all dimensions.
\end{funcdesc}

\begin{funcdesc}{fourier_uniform}{input, size, n=-1, axis=-1, output=None}
  The \function{fourier_uniform} function multiplies the input array with 
  the multi-dimensional Fourier transform of a uniform filter with given
  sizes \var{size}. The \var{size} parameter is a sequences of
  values for each dimension, or a single value for all dimensions.
\end{funcdesc}

\begin{funcdesc}{fourier_ellipsoid}{input, size, n=-1, axis=-1, 
  output=None}
  The \function{fourier_ellipsoid} function multiplies the input array with 
  the multi-dimensional Fourier transform of a elliptically shaped filter 
  with given sizes \var{size}. The \var{size} parameter is a sequences of 
  values for each dimension, or a single value for all dimensions.  
  \note{This function is
    only implemented for dimensions 1, 2, and 3.}
\end{funcdesc}

\section{Interpolation functions}
This section describes various interpolation functions that are based on
B-spline theory. A good introduction to B-splines can be found in: M. 
Unser, "Splines: A Perfect Fit for Signal and Image Processing," IEEE 
Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, November 1999.
\subsection{Spline pre-filters}
Interpolation using splines of an order larger than 1 requires a pre-
filtering step. The interpolation functions described in section
\ref{sec:ndimage:interpolation} apply pre-filtering by calling
\function{spline_filter}, but they can be instructed not to do this by 
setting the \var{prefilter} keyword equal to \constant{False}.  This is 
useful if more than one interpolation operation is done on the same array. 
In this case it is more efficient to do the pre-filtering only once and use 
a prefiltered array as the input of the interpolation functions. The 
following two functions implement the pre-filtering:

\begin{funcdesc}{spline_filter1d}{input, order=3, axis=-1, output=None,
    output_type=numarray.Float64} The \function{spline_filter1d} function
  calculates a one-dimensional spline filter along the given axis. An 
  output array can optionally be provided. The order of the spline must be 
  larger then 1 and less than 6.
\end{funcdesc}

\begin{funcdesc}{spline_filter}{input, order=3, output=None, 
    output_type=numarray.Float64} The \function{spline_filter} function
  calculates a multi-dimensional spline filter.
  
  \note{The multi-dimensional filter is implemented as a sequence of
    one-dimensional spline filters. The intermediate arrays are stored in 
    the same data type as the output. Therefore, if an output 
    with a limited precision is requested, the results may be imprecise 
    because intermediate results may be stored with insufficient precision. 
    This can be prevented by specifying a output type of high precision.}
\end{funcdesc}

\subsection{Interpolation functions}
\label{sec:ndimage:interpolation}
Following functions all employ spline interpolation to effect some type of
geometric transformation of the input array. This requires a mapping of the
output coordinates to the input coordinates, and therefore the possibility
arises that input values outside the boundaries are needed. This problem 
is solved in the same way as described in section
\ref{sec:ndimage:filter-functions} for the multi-dimensional filter 
functions. Therefore these functions all support a \var{mode} parameter 
that determines how the boundaries are handled, and a \var{cval} parameter 
that gives a constant value in case that the \constant{'constant'} mode is 
used.

\begin{funcdesc}{geometric_transform}{input, mapping, output_shape=None,
    output_type=None, output=None, order=3, mode='constant', cval=0.0,
    prefilter=True, extra_arguments = (), extra_keywords = {}} The \function{geometric_transform} function applies an
  arbitrary geometric transform to the input. The given \var{mapping} 
  function is called at each point in the output to find the corresponding 
  coordinates in the input.  \var{mapping} must be a callable object that 
  accepts a tuple of length equal to the output array rank and returns the 
  corresponding input coordinates as a tuple of length equal to the input 
  array rank. The output shape and output type can optionally be provided. 
  If not given they are equal to the input shape and type.
  
  For example:
\begin{verbatim}
>>> a = arange(12, shape=(4,3), type = Float64)
>>> def shift_func(output_coordinates):
...     return (output_coordinates[0] - 0.5, output_coordinates[1] - 0.5)
... 
>>> print geometric_transform(a, shift_func)
[[ 0.      0.      0.    ]
 [ 0.      1.3625  2.7375]
 [ 0.      4.8125  6.1875]
 [ 0.      8.2625  9.6375]]  
\end{verbatim}  

  Optionally extra arguments can be defined and passed to the filter 
  function. The \var{extra_arguments} and \var{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 shifts in our example as arguments:

\begin{verbatim}
>>> def shift_func(output_coordinates, s0, s1):
...     return (output_coordinates[0] - s0, output_coordinates[1] - s1)
... 
>>> print geometric_transform(a, shift_func, extra_arguments = (0.5, 0.5))
[[ 0.      0.      0.    ]
 [ 0.      1.3625  2.7375]
 [ 0.      4.8125  6.1875]
 [ 0.      8.2625  9.6375]]  
\end{verbatim}  
or
\begin{verbatim}
>>> print geometric_transform(a, shift_func, extra_keywords = {'s0': 0.5, 's1': 0.5})
[[ 0.      0.      0.    ]
 [ 0.      1.3625  2.7375]
 [ 0.      4.8125  6.1875]
 [ 0.      8.2625  9.6375]]  
\end{verbatim}  

\note{The mapping function can also be written in C and passed using a CObject. See \ref{sec:ndimage:ccallbacks} for more information.}
\end{funcdesc}

\begin{funcdesc}{map_coordinates}{input, coordinates, output_type=None, 
    output=None, order=3, mode='constant', cval=0.0, prefilter=True} 
  The function \function{map_coordinates} applies an arbitrary coordinate
  transformation using the given array of coordinates. The shape of the 
  output is derived from that of the coordinate array by dropping the first 
  axis. The parameter \var{coordinates} is used to find for each point in 
  the output the corresponding coordinates in the input. The values of 
  \var{coordinates} along the first axis are the coordinates in the input 
  array at which the output value is found. (See also the numarray 
  \function{coordinates} function.) Since the coordinates may be non-
  integer coordinates, the value of the input at these coordinates is 
  determined by spline interpolation of the requested order. Here is an 
  example that interpolates a 2D array at (0.5, 0.5) and (1, 2):
\begin{verbatim}
>>> a = arange(12, shape=(4,3), type = numarray.Float64)
>>> print a
[[  0.   1.   2.]
 [  3.   4.   5.]
 [  6.   7.   8.]
 [  9.  10.  11.]]
>>> print map_coordinates(a, [[0.5, 2], [0.5, 1]])
[ 1.3625  7.    ]
\end{verbatim}
\end{funcdesc}

\begin{funcdesc}{affine_transform}{input, matrix, offset=0.0, 
  output_shape=None, output_type=None, output=None, order=3, 
  mode='constant', cval=0.0, prefilter=True} The 
  \function{affine_transform} function applies an affine transformation to 
  the input array. The given transformation \var{matrix} and \var{offset} 
  are used to find for each point in the output the corresponding 
  coordinates in the input.  The value of the input at the
  calculated coordinates is determined by spline interpolation of the 
  requested order. The transformation \var{matrix} must be two-dimensional 
  or can also be given as a one-dimensional sequence or array.  In the 
  latter case, it is assumed that the matrix is diagonal. A more efficient 
  interpolation algorithm is then applied that exploits the separability of 
  the problem.  The output shape and output type can optionally be 
  provided. If not given they are equal to the input shape and type.
\end{funcdesc}

\begin{funcdesc}{shift}{input, shift, output_type=None, output=None, 
  order=3, mode='constant', cval=0.0, prefilter=True} The \function{shift} 
  function returns a shifted version of the input, using spline 
  interpolation of the requested \var{order}.
\end{funcdesc}

\begin{funcdesc}{zoom}{input, zoom, output_type=None, output=None, order=3, 
    mode='constant', cval=0.0, prefilter=True} The \function{zoom} function
  returns a rescaled version of the input, using spline interpolation of 
  the requested \var{order}.
\end{funcdesc}

\begin{funcdesc}{rotate}{input, angle, axes=(-1, -2), reshape=1,
    output_type=None, output=None, order=3, mode='constant', cval=0.0,
    prefilter=True} The \function{rotate} function returns the input array
  rotated in the plane defined by the two axes given by the parameter
  \var{axes}, using spline interpolation of the requested \var{order}. The
  angle must be given in degrees. If \var{reshape} is true, then the size 
  of the output array is adapted to contain the rotated input.
\end{funcdesc}

\section{Binary morphology}
\label{sec:ndimage:binary-morphology}

\begin{funcdesc}{generate_binary_structure}{rank, connectivity}
  The \function{generate_binary_structure} functions generates a binary
  structuring element for use in binary morphology operations. The 
  \var{rank} of the structure must be provided. The size of the structure 
  that is returned is equal to three in each direction. The value of each 
  element is equal to one if the square of the Euclidean distance from the 
  element to the center is less or equal to \var{connectivity}. For 
  instance, two dimensional 4-connected and 8-connected structures are 
  generated as follows:
\begin{verbatim}
>>> print generate_binary_structure(2, 1)
[[0 1 0]
 [1 1 1]
 [0 1 0]]
>>> print generate_binary_structure(2, 2)
[[1 1 1]
 [1 1 1]
 [1 1 1]]
\end{verbatim}
\end{funcdesc}

Most binary morphology functions can be expressed in terms of the basic
operations erosion and dilation:

\begin{funcdesc}{binary_erosion}{input, structure=None, iterations=1,
    mask=None, output=None, border_value=0, origin=0} The
  \function{binary_erosion} function implements binary erosion of arrays of
  arbitrary rank with the given structuring element. The origin parameter
  controls the placement of the structuring element as described in section
  \ref{sec:ndimage:filter-functions}. If no structuring element is 
  provided, an element with connectivity equal to one is generated using
  \function{generate_binary_structure}. The \var{border_value} parameter 
  gives the value of the array outside boundaries. The erosion is repeated
  \var{iterations} times. If \var{iterations} is less than one, the erosion 
  is repeated until the result does not change anymore. If a \var{mask} 
  array is given, only those elements with a true value at the 
  corresponding mask element are modified at each iteration.
\end{funcdesc}

\begin{funcdesc}{binary_dilation}{input, structure=None, iterations=1,
    mask=None, output=None, border_value=0, origin=0} The
  \function{binary_dilation} function implements binary dilation of arrays 
  of arbitrary rank with the given structuring element. The origin 
  parameter controls the placement of the structuring element as described 
  in section \ref{sec:ndimage:filter-functions}. If no structuring element 
  is provided, an element with connectivity equal to one is generated using
  \function{generate_binary_structure}. The \var{border_value} parameter 
  gives the value of the array outside boundaries. The dilation is repeated
  \var{iterations} times.  If \var{iterations} is less than one, the 
  dilation is repeated until the result does not change anymore. If a 
  \var{mask} array is given, only those elements with a true value at the 
  corresponding mask element are modified at each iteration.

  Here is an example of using \function{binary_dilation} to find all 
  elements that touch the border, by repeatedly dilating an empty array 
  from the border using the data array as the mask:
\begin{verbatim}
>>> struct = array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
>>> a = array([[1,0,0,0,0], [1,1,0,1,0], [0,0,1,1,0], [0,0,0,0,0]])
>>> print a
[[1 0 0 0 0]
 [1 1 0 1 0]
 [0 0 1 1 0]
 [0 0 0 0 0]]
>>> print binary_dilation(zeros(a.shape), struct, -1, a, border_value=1)
[[1 0 0 0 0]
 [1 1 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]]
\end{verbatim}
\end{funcdesc}

The \function{binary_erosion} and \function{binary_dilation} functions both
have an \var{iterations} parameter which allows the erosion or dilation to 
be repeated a number of times. Repeating an erosion or a dilation with a 
given structure \constant{n} times is equivalent to an erosion or a 
dilation with a structure that is \constant{n-1} times dilated with itself. 
A function is provided that allows the calculation of a structure that is 
dilated a number of times with itself:

\begin{funcdesc}{iterate_structure}{structure, iterations, origin=None} 
  The \function{iterate_structure} function returns a structure by dilation 
  of the input structure \var{iteration} - 1 times with itself. For 
  instance:
  \begin{verbatim}
>>> struct = generate_binary_structure(2, 1)
>>> print struct
[[0 1 0]
 [1 1 1]
 [0 1 0]]
>>> print iterate_structure(struct, 2)
[[0 0 1 0 0]
 [0 1 1 1 0]
 [1 1 1 1 1]
 [0 1 1 1 0]
 [0 0 1 0 0]]
\end{verbatim}
  If the origin of the original structure is equal to 0, then it is also 
  equal to 0 for the iterated structure. If not, the origin must also be 
  adapted if the equivalent of the \var{iterations} erosions or dilations 
  must be achieved with the iterated structure. The adapted origin is 
  simply obtained by multiplying with the number of iterations. For 
  convenience the
  \function{iterate_structure} also returns the adapted origin if the
  \var{origin} parameter is not \constant{None}:
\begin{verbatim}
>>> print iterate_structure(struct, 2, -1)
(array([[0, 0, 1, 0, 0],
       [0, 1, 1, 1, 0],
       [1, 1, 1, 1, 1],
       [0, 1, 1, 1, 0],
       [0, 0, 1, 0, 0]], type=Bool), [-2, -2])
\end{verbatim}
\end{funcdesc}

Other morphology operations can be defined in terms of erosion and d
dilation. Following functions provide a few of these operations for 
convenience:

\begin{funcdesc}{binary_opening}{input, structure=None, iterations=1,
  output=None, origin=0} The \function{binary_opening} function implements
  binary opening of arrays of arbitrary rank with the given structuring
  element. Binary opening is equivalent to a binary erosion followed by a
  binary dilation with the same structuring element. The origin parameter
  controls the placement of the structuring element as described in section
  \ref{sec:ndimage:filter-functions}. If no structuring element is 
  provided, an element with connectivity equal to one is generated using
  \function{generate_binary_structure}. The \var{iterations} parameter 
  gives the number of erosions that is performed followed by the same 
  number of dilations.
\end{funcdesc}

\begin{funcdesc}{binary_closing}{input, structure=None, iterations=1,
  output=None, origin=0} The \function{binary_closing} function implements
  binary closing of arrays of arbitrary rank with the given structuring
  element. Binary closing is equivalent to a binary dilation followed by a
  binary erosion with the same structuring element. The origin parameter
  controls the placement of the structuring element as described in section
  \ref{sec:ndimage:filter-functions}. If no structuring element is 
  provided, an element with connectivity equal to one is generated using
  \function{generate_binary_structure}. The \var{iterations} parameter   
  gives the number of dilations that is performed followed by the same 
  number of erosions.
\end{funcdesc}

\begin{funcdesc}{binary_fill_holes}{input, structure = None, output = None, 
origin = 0} The \function{binary_fill_holes} function is used to close 
holes in objects in a binary image, where the structure defines the 
connectivity of the holes. The origin parameter controls the placement of 
the structuring element as described in section \ref{sec:ndimage:filter-
functions}. If no structuring element is provided, an element with 
connectivity equal to one is generated using 
\function{generate_binary_structure}. 
\end{funcdesc}

\begin{funcdesc}{binary_hit_or_miss}{input, structure1=None, 
  structure2=None, output=None, origin1=0, origin2=None} The 
  \function{binary_hit_or_miss}
  function implements a binary hit-or-miss transform of arrays of arbitrary
  rank with the given structuring elements.  The hit-or-miss transform is
  calculated by erosion of the input with the first structure, erosion of   
  the logical \emph{not} of the input with the second structure, followed 
  by the logical \emph{and} of these two erosions.  The origin parameters 
  control the placement of the structuring elements as described in section
  \ref{sec:ndimage:filter-functions}. If \var{origin2} equals 
  \constant{None} it is set equal to the \var{origin1} parameter. If the 
  first structuring element is not provided, a structuring element with 
  connectivity equal to one is generated using 
  \function{generate_binary_structure}, if \var{structure2} is not 
  provided, it is set equal to the logical \emph{not} of \var{structure1}.
\end{funcdesc}

\section{Grey-scale morphology}
\label{sec:ndimage:grey-morphology}

Grey-scale morphology operations are the equivalents of binary morphology
operations that operate on arrays with arbitrary values. Below we describe 
the grey-scale equivalents of erosion, dilation, opening and closing. These
operations are implemented in a similar fashion as the filters described in
section \ref{sec:ndimage:filter-functions}, and we refer to this section 
for the description of filter kernels and footprints, and the handling of 
array borders. The grey-scale morphology operations optionally take a 
\var{structure} parameter that gives the values of the structuring element. 
If this parameter is not given the structuring element is assumed to be 
flat with a value equal to zero. The shape of the structure can optionally 
be defined by the \var{footprint} parameter. If this parameter is not 
given, the structure is assumed to be rectangular, with sizes equal to the 
dimensions of the \var{structure} array, or by the \var{size} parameter if 
\var{structure} is not given. The \var{size} parameter is only used if both 
\var{structure} and \var{footprint} are not given, in which case the 
structuring element is assumed to be rectangular and flat with the 
dimensions given by \var{size}. The \var{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 \var{footprint} 
parameter, if provided, must be an array that defines the shape of the 
kernel by its non-zero elements.

Similar to binary erosion and dilation there are operations for grey-scale
erosion and dilation:

\begin{funcdesc}{grey_erosion}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{grey_erosion} function calculates a multi-dimensional grey-
  scale erosion.
\end{funcdesc}

\begin{funcdesc}{grey_dilation}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{grey_dilation} function calculates a multi-dimensional grey-
  scale dilation.
\end{funcdesc}

Grey-scale opening and closing operations can be defined similar to their
binary counterparts:

\begin{funcdesc}{grey_opening}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{grey_opening} function implements grey-scale opening of arrays 
  of arbitrary rank. Grey-scale opening is equivalent to a grey-scale 
  erosion followed by a grey-scale dilation.
\end{funcdesc}

\begin{funcdesc}{grey_closing}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{grey_closing} function implements grey-scale closing of arrays 
  of arbitrary rank. Grey-scale opening is equivalent to a grey-scale 
  dilation followed by a grey-scale erosion.
\end{funcdesc}

\begin{funcdesc}{morphological_gradient}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{morphological_gradient} function implements a grey-scale
  morphological gradient of arrays of arbitrary rank. The grey-scale
  morphological gradient is equal to the difference of a grey-scale 
  dilation and a grey-scale erosion.
\end{funcdesc}

\begin{funcdesc}{morphological_laplace}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{morphological_laplace} function implements a grey-scale
  morphological laplace of arrays of arbitrary rank. The grey-scale
  morphological laplace is equal to the sum of a grey-scale dilation and a
  grey-scale erosion minus twice the input.
\end{funcdesc}

\begin{funcdesc}{white_tophat}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{white_tophat} function implements a white top-hat filter of 
  arrays of arbitrary rank. The white top-hat is equal to the difference of 
  the input and a grey-scale opening.
\end{funcdesc}

\begin{funcdesc}{black_tophat}{input, size=None, footprint=None, 
    structure=None, output=None, mode='reflect', cval=0.0, origin=0} The
  \function{black_tophat} function implements a black top-hat filter of 
  arrays of arbitrary rank. The black top-hat is equal to the difference of 
  the a grey-scale closing and the input.
\end{funcdesc}

\section{Distance transforms}
\label{sec:ndimage:grey-morphology}
Distance transforms are used to calculate the minimum distance from each
element of an object to the background. The following functions implement
distance transforms for three different distance metrics: Euclidean, City
Block, and Chessboard distances.

\begin{funcdesc}{distance_transform_cdt}{input, structure="chessboard",
  return_distances=True, return_indices=False, distances=None, 
  indices=None} The function \function{distance_transform_cdt} uses a 
  chamfer type algorithm to calculate the distance transform of the input, 
  by replacing each object element (defined by values larger than zero) 
  with the shortest distance to the background (all non-object elements). 
  The structure determines the type of chamfering that is done. If the 
  structure is equal to 'cityblock' a structure is generated using 
  \function{generate_binary_structure} with a squared distance equal to 1. 
  If the structure is equal to 'chessboard', a structure is generated using 
  \function{generate_binary_structure} with a squared distance equal to the 
  rank of the array. These choices correspond to the common interpretations 
  of the cityblock and the chessboard distancemetrics in two dimensions.
  
  In addition to the distance transform, the feature transform can be
  calculated. In this case the index of the closest background element is
  returned along the first axis of the result.  The \var{return_distances}, 
  and \var{return_indices} flags can be used to indicate if the distance 
  transform, the feature transform, or both must be returned.
  
  The \var{distances} and \var{indices} arguments can be used to give 
  optional output arrays that must be of the correct size and type (both
  \constant{Int32}).

  The basics of the algorithm used to implement this function is described
  in: G. Borgefors, "Distance transformations in arbitrary dimensions.",
  Computer Vision, Graphics, and Image Processing, 27:321--345, 1984.
\end{funcdesc}

\begin{funcdesc}{distance_transform_edt}{input, sampling=None,
  return_distances=True, return_indices=False, distances=None, 
  indices=None} The function \function{distance_transform_edt} calculates 
  the exact euclidean distance transform of the input, by replacing each 
  object element (defined by values larger than zero) with the shortest 
  euclidean distance to the background (all non-object elements).
  
  In addition to the distance transform, the feature transform can be
  calculated. In this case the index of the closest background element is
  returned along the first axis of the result.  The \var{return_distances}, 
  and \var{return_indices} flags can be used to indicate if the distance 
  transform, the feature transform, or both must be returned.
  
  Optionally the sampling along each axis can be given by the 
  \var{sampling} parameter which should be a sequence of length equal to 
  the input rank, or a single number in which the sampling is assumed to be 
  equal along all axes.

  The \var{distances} and \var{indices} arguments can be used to give 
  optional output arrays that must be of the correct size and type 
  (\constant{Float64} and \constant{Int32}).
  
  The algorithm used to implement this function is described in: C. R. 
  Maurer, Jr., R. Qi, and V. Raghavan, "A linear time algorithm for 
  computing exact euclidean distance transforms of binary images in 
  arbitrary dimensions. IEEE Trans. PAMI 25, 265-270, 2003.
\end{funcdesc}

\begin{funcdesc}{distance_transform_bf}{input, metric="euclidean",
  sampling=None, return_distances=True, return_indices=False, 
  distances=None, indices=None} The function  
  \function{distance_transform_bf} uses a brute-force algorithm to 
  calculate the distance transform of the input, by replacing each object 
  element (defined by values larger than zero) with the shortest distance 
  to the background (all non-object elements).  The metric must be one of 
  \constant{"euclidean"}, \constant{"cityblock"}, or 
  \constant{"chessboard"}.
  
  In addition to the distance transform, the feature transform can be
  calculated. In this case the index of the closest background element is
  returned along the first axis of the result.  The \var{return_distances}, 
  and \var{return_indices} flags can be used to indicate if the distance 
  transform, the feature transform, or both must be returned.
  
  Optionally the sampling along each axis can be given by the 
  \var{sampling} parameter which should be a sequence of length equal to 
  the input rank, or a single number in which the sampling is assumed to be 
  equal along all axes. This parameter is only used in the case of the 
  euclidean distance transform.

  The \var{distances} and \var{indices} arguments can be used to give 
  optional output arrays that must be of the correct size and type 
  (\constant{Float64} and \constant{Int32}).

  \note{This function uses a slow brute-force algorithm, the function
    \function{distance_transform_cdt} can be used to more efficiently 
    calculate cityblock and chessboard distance transforms. The function
    \function{distance_transform_edt} can be used to more efficiently 
    calculate the exact euclidean distance transform.}
\end{funcdesc}

\section{Segmentation and labeling}
Segmentation is the process of separating objects of interest from the
background. The most simple approach is probably intensity thresholding, 
which is easily done with \module{numarray} functions:
\begin{verbatim}
>>> a = array([[1,2,2,1,1,0],
...            [0,2,3,1,2,0],
...            [1,1,1,3,3,2],
...            [1,1,1,1,2,1]])
>>> print where(a > 1, 1, 0)
[[0 1 1 0 0 0]
 [0 1 1 0 1 0]
 [0 0 0 1 1 1]
 [0 0 0 0 1 0]]
\end{verbatim}

The result is a binary image, in which the individual objects still need to 
be identified and labeled.  The function \function{label} generates an 
array where each object is assigned a unique number:

\begin{funcdesc}{label}{input, structure=None, output=None}
  The \function{label} function generates an array where the objects in the
  input are labeled with an integer index. It returns a tuple consisting of 
  the array of object labels and the number of objects found, unless the
  \var{output} parameter is given, in which case only the number of objects 
  is returned. The connectivity of the objects is defined by a structuring
  element. For instance, in two dimensions using a four-connected 
  structuring element gives:
\begin{verbatim}
>>> a = array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
>>> s = [[0, 1, 0], [1,1,1], [0,1,0]]
>>> print label(a, s)
(array([[0, 1, 1, 0, 0, 0],
       [0, 1, 1, 0, 2, 0],
       [0, 0, 0, 2, 2, 2],
       [0, 0, 0, 0, 2, 0]]), 2)
\end{verbatim}
These two objects are not connected because there is no way in which we can
place the structuring element such that it overlaps with both objects. 
However, an 8-connected structuring element results in only a single 
object:
\begin{verbatim}
>>> a = array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
>>> s = [[1,1,1], [1,1,1], [1,1,1]]
>>> print label(a, s)[0]
[[0 1 1 0 0 0]
 [0 1 1 0 1 0]
 [0 0 0 1 1 1]
 [0 0 0 0 1 0]]
\end{verbatim}
If no structuring element is provided, one is generated by calling
\function{generate_binary_structure} (see section \ref{sec:ndimage:
morphology}) using a connectivity of one (which in 2D is the 4-connected 
structure of the first example).  The input can be of any type, any value 
not equal to zero is taken to be part of an object. This is useful if you 
need to 're-label' an array of object indices, for instance after removing 
unwanted objects. Just apply the label function again to the index array. 
For instance:
\begin{verbatim}
>>> l, n = label([1, 0, 1, 0, 1])
>>> print l
[1 0 2 0 3]
>>> l = where(l != 2, l, 0)
>>> print l
[1 0 0 0 3]
>>> print label(l)[0]
[1 0 0 0 2]
\end{verbatim}

\note{The structuring element used by \function{label} is assumed to be
  symmetric.}
\end{funcdesc}

There is a large number of other approaches for segmentation, for instance 
from an estimation of the borders of the objects that can be obtained for 
instance by derivative filters. One such an approach is watershed 
segmentation.  The function \function{watershed_ift} generates an array 
where each object is assigned a unique label, from an array that localizes 
the object borders, generated for instance by a gradient magnitude filter. 
It uses an array containing initial markers for the objects:
\begin{funcdesc}{watershed_ift}{input, markers, structure=None, 
  output=None} The \function{watershed_ift} function applies a watershed 
  from markers algorithm, using an Iterative Forest Transform, as described 
  in: P. Felkel, R.  Wegenkittl, and M. Bruckschwaiger, "Implementation and 
  Complexity of the Watershed-from-Markers Algorithm Computed as a Minimal 
  Cost Forest.", Eurographics 2001, pp. C:26-35.
  
  The inputs of this function are the array to which the transform is 
  applied, and an array of markers that designate the objects by a unique 
  label, where any non-zero value is a marker. For instance:
\begin{verbatim}
>>> input = array([[0, 0, 0, 0, 0, 0, 0],
...                [0, 1, 1, 1, 1, 1, 0],
...                [0, 1, 0, 0, 0, 1, 0],
...                [0, 1, 0, 0, 0, 1, 0],
...                [0, 1, 0, 0, 0, 1, 0],
...                [0, 1, 1, 1, 1, 1, 0],
...                [0, 0, 0, 0, 0, 0, 0]], numarray.UInt8)
>>> markers = array([[1, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 2, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0]], numarray.Int8)
>>> print watershed_ift(input, markers)
[[1 1 1 1 1 1 1]
 [1 1 2 2 2 1 1]
 [1 2 2 2 2 2 1]
 [1 2 2 2 2 2 1]
 [1 2 2 2 2 2 1]
 [1 1 2 2 2 1 1]
 [1 1 1 1 1 1 1]]
\end{verbatim}
  
  Here two markers were used to designate an object (marker=2) and the
  background (marker=1).  The order in which these are processed is 
  arbitrary: moving the marker for the background to the lower right corner 
  of the array yields a different result:
\begin{verbatim}
>>> markers = array([[0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 2, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 1]], numarray.Int8)
>>> print watershed_ift(input, markers)
[[1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1]
 [1 1 2 2 2 1 1]
 [1 1 2 2 2 1 1]
 [1 1 2 2 2 1 1]
 [1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1]]
\end{verbatim}
  The result is that the object (marker=2) is smaller because the second 
  marker was processed earlier. This may not be the desired effect if the 
  first marker was supposed to designate a background object. Therefore
  \function{watershed_ift} treats markers with a negative value explicitly 
  as background markers and processes them after the normal markers. For 
  instance, replacing the first marker by a negative marker gives a result 
  similar to the first example:
\begin{verbatim}
>>> markers = array([[0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 2, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, 0],
...                  [0, 0, 0, 0, 0, 0, -1]], numarray.Int8)
>>> print watershed_ift(input, markers)
[[-1 -1 -1 -1 -1 -1 -1]
 [-1 -1  2  2  2 -1 -1]
 [-1  2  2  2  2  2 -1]
 [-1  2  2  2  2  2 -1]
 [-1  2  2  2  2  2 -1]
 [-1 -1  2  2  2 -1 -1]
 [-1 -1 -1 -1 -1 -1 -1]]
\end{verbatim}
  
  The connectivity of the objects is defined by a structuring element. If 
  no structuring element is provided, one is generated by calling
  \function{generate_binary_structure} (see section
  \ref{sec:ndimage:morphology}) using a connectivity of one (which in 2D is 
  a 4-connected structure.) For example, using an 8-connected structure 
  with the last example yields a different object:
\begin{verbatim}
>>> print watershed_ift(input, markers,
...                     structure = [[1,1,1], [1,1,1], [1,1,1]])
[[-1 -1 -1 -1 -1 -1 -1]
 [-1  2  2  2  2  2 -1]
 [-1  2  2  2  2  2 -1]
 [-1  2  2  2  2  2 -1]
 [-1  2  2  2  2  2 -1]
 [-1  2  2  2  2  2 -1]
 [-1 -1 -1 -1 -1 -1 -1]]
\end{verbatim}

\note{The implementation of \function{watershed_ift} limits the data types 
of the input to \constant{UInt8} and \constant{UInt16}.}
\end{funcdesc}

\section{Object measurements}
Given an array of labeled objects, the properties of the individual objects 
can be measured. The \function{find_objects} function can be used to 
generate a list of slices that for each object, give the smallest sub-array 
that fully contains the object:

\begin{funcdesc}{find_objects}{input, max_label=0}
  The \function{find_objects} finds all objects in a labeled array and 
  returns a list of slices that correspond to the smallest regions in the 
  array that contains the object. For instance:
\begin{verbatim}
>>> a = array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
>>> l, n = label(a)
>>> f = find_objects(l)
>>> print a[f[0]]
[[1 1]
 [1 1]]
>>> print a[f[1]]
[[0 1 0]
 [1 1 1]
 [0 1 0]]
\end{verbatim}
\function{find_objects} returns slices for all objects, unless the
\var{max_label} parameter is larger then zero, in which case only the first
\var{max_label} objects are returned. If an index is missing in the 
\var{label} array, \constant{None} is return instead of a slice. For 
example:
\begin{verbatim}
>>> print find_objects([1, 0, 3, 4], max_label = 3)
[(slice(0, 1, None),), None, (slice(2, 3, None),)]
\end{verbatim}
\end{funcdesc}

The list of slices generated by \function{find_objects} is useful to find 
the position and dimensions of the objects in the array, but can also be 
used to perform measurements on the individual objects. Say we want to find 
the sum of the intensities of an object in image:
\begin{verbatim}
>>> image = arange(4*6,shape=(4,6))
>>> mask = array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
>>> labels = label(mask)[0]
>>> slices = find_objects(labels)
\end{verbatim}
Then we can calculate the sum of the elements in the second object:
\begin{verbatim}
>>> print where(labels[slices[1]] == 2, image[slices[1]], 0).sum()
80
\end{verbatim}
That is however not particularly efficient, and may also be more 
complicated for other types of measurements. Therefore a few measurements 
functions are defined that accept the array of object labels and the index 
of the object to be measured. For instance calculating the sum of the 
intensities can be done by:
\begin{verbatim}
>>> print sum(image, labels, 2)
80.0
\end{verbatim}
For large arrays and small objects it is more efficient to call the 
measurement functions after slicing the array:
\begin{verbatim}
>>> print sum(image[slices[1]], labels[slices[1]], 2)
80.0
\end{verbatim}
Alternatively, we can do the measurements for a number of labels with a 
single function call, returning a list of results. For instance, to measure 
the sum of the values of the background and the second object in our 
example we give a list of labels:
\begin{verbatim}
>>> print sum(image, labels, [0, 2])
[178.0, 80.0]
\end{verbatim}

The measurement functions described below all support the \var{index} 
parameter to indicate which object(s) should be measured. The default value 
of \var{index} is \constant{None}. This indicates that all elements where 
the label is larger than zero should be treated as a single object and 
measured. Thus, in this case the \var{labels} array is treated as a mask 
defined by the elements that are larger than zero. If \var{index} is a 
number or a sequence of numbers it gives the labels of the objects that are 
measured. If \var{index} is a sequence, a list of the results is returned. 
Functions that return more than one result, return their result as a tuple 
if \var{index} is a single number, or as a tuple of lists, if \var{index} 
is a sequence.

\begin{funcdesc}{sum}{input, labels=None, index=None}
  The \function{sum} function calculates the sum of the elements of the 
  object with label(s) given by \var{index}, using the \var{labels} array 
  for the object labels. If \var{index} is \constant{None}, all elements 
  with a non-zero label value are treated as a single object. If 
  \var{label} is \constant{None}, all elements of \var{input} are used in 
  the calculation.
\end{funcdesc}

\begin{funcdesc}{mean}{input, labels=None, index=None}
  The \function{mean} function calculates the mean of the elements of the
  object with label(s) given by \var{index}, using the \var{labels} array 
  for the object labels. If \var{index} is \constant{None}, all elements 
  with a non-zero label value are treated as a single object. If 
  \var{label} is \constant{None}, all elements of \var{input} are used in 
  the calculation.
\end{funcdesc}

\begin{funcdesc}{variance}{input, labels=None, index=None}
  The \function{variance} function calculates the variance of the elements 
  of the object with label(s) given by \var{index}, using the \var{labels} 
  array for the object labels. If \var{index} is \constant{None}, all 
  elements with a non-zero label value are treated as a single object. If 
  \var{label} is \constant{None}, all elements of \var{input} are used in 
  the calculation.
\end{funcdesc}

\begin{funcdesc}{standard_deviation}{input, labels=None, index=None}
  The \function{standard_deviation} function calculates the standard 
  deviation of the elements of the object with label(s) given by 
  \var{index}, using the \var{labels} array for the object labels. If 
  \var{index} is \constant{None}, all elements with a non-zero label value 
  are treated as a single object. If \var{label} is \constant{None}, all 
  elements of \var{input} are used in the calculation.
\end{funcdesc}

\begin{funcdesc}{minimum}{input, labels=None, index=None}
  The \function{minimum} function calculates the minimum of the elements of 
  the object with label(s) given by \var{index}, using the \var{labels} 
  array for the object labels. If \var{index} is \constant{None}, all 
  elements with a non-zero label value are treated as a single object. If 
  \var{label} is \constant{None}, all elements of \var{input} are used in 
  the calculation.
\end{funcdesc}

\begin{funcdesc}{maximum}{input, labels=None, index=None}
  The \function{maximum} function calculates the maximum of the elements of 
  the object with label(s) given by \var{index}, using the \var{labels} 
  array for the object labels. If \var{index} is \constant{None}, all 
  elements with a non-zero label value are treated as a single object. If 
  \var{label} is \constant{None}, all elements of \var{input} are used in 
  the calculation.
\end{funcdesc}

\begin{funcdesc}{minimum_position}{input, labels=None, index=None}
  The \function{minimum_position} function calculates the position of the
  minimum of the elements of the object with label(s) given by \var{index},
  using the \var{labels} array for the object labels. If \var{index} is
  \constant{None}, all elements with a non-zero label value are treated as 
  a single object. If \var{label} is \constant{None}, all elements of 
  \var{input} are used in the calculation.
\end{funcdesc}

\begin{funcdesc}{maximum_position}{input, labels=None, index=None}
  The \function{maximum_position} function calculates the position of the
  maximum of the elements of the object with label(s) given by \var{index},
  using the \var{labels} array for the object labels. If \var{index} is
  \constant{None}, all elements with a non-zero label value are treated as 
  a single object. If \var{label} is \constant{None}, all elements of 
  \var{input} are used in the calculation.
\end{funcdesc}

\begin{funcdesc}{extrema}{input, labels=None, index=None}
  The \function{extrema} function calculates the minimum, the maximum, and 
  their positions, of the elements of the object with label(s) given by 
  \var{index}, using the \var{labels} array for the object labels. If 
  \var{index} is \constant{None}, all elements with a non-zero label value 
  are treated as a single object. If \var{label} is \constant{None}, all 
  elements of \var{input} are used in the calculation. The result is a 
  tuple giving the minimum, the maximum, the position of the mininum and 
  the postition of the maximum. The result is the same as a tuple formed by 
  the results of the functions \function{minimum}, \function{maximum}, 
  \function{minimum_position}, and \function{maximum_position} that are 
  described above.
\end{funcdesc}

\begin{funcdesc}{center_of_mass}{input, labels=None, index=None}
  The \function{center_of_mass} function calculates the center of mass of 
  the of the object with label(s) given by \var{index}, using the 
  \var{labels} array for the object labels. If \var{index} is 
  \constant{None}, all elements with a non-zero label value are treated as 
  a single object. If \var{label} is \constant{None}, all elements of 
  \var{input} are used in the calculation.
\end{funcdesc}

\begin{funcdesc}{histogram}{input, min, max, bins, labels=None, index=None}
  The \function{histogram} function calculates a histogram of 
  the of the object with label(s) given by \var{index}, using the 
  \var{labels} array for the object labels. If \var{index} is 
  \constant{None}, all elements with a non-zero label value are treated as 
  a single object. If \var{label} is \constant{None}, all elements of 
  \var{input} are used in the calculation. Histograms are defined by their 
  minimum (\var{min}), maximum (\var{max}) and the number of bins 
  (\var{bins}). They are returned as one-dimensional arrays of type Int32. 
\end{funcdesc}

\section{Extending \module{nd\_image} in C}
\label{sec:ndimage:ccallbacks}
\subsection{C callback functions}
A few functions in the \module{numarray.nd\_image} take a call-back 
argument. This can be a python function, but also a CObject containing a 
pointer to a C function. To use this feature, you must write your own C 
extension that defines the function, and define a python function that 
returns a CObject containing a pointer to this function.

An example of a function that supports this is 
\function{geometric_transform} (see section \ref{sec:ndimage:
interpolation}). You can pass it a python callable object that defines a 
mapping from all output coordinates to corresponding coordinates in the 
input array. This mapping function can also be a C function, which 
generally will be much more efficient, since the overhead of calling a
python function at each element is avoided.

For example to implement a simple shift function we define the following 
function:
\begin{verbatim}
static int 
_shift_function(int *output_coordinates, double* input_coordinates,
                int output_rank, int input_rank, void *callback_data)
{
  int ii;
  /* get the shift from the callback data pointer: */
  double shift = *(double*)callback_data;
  /* calculate the coordinates: */
  for(ii = 0; ii < irank; ii++)
    icoor[ii] = ocoor[ii] - shift;
  /* return OK status: */
  return 1;
}
\end{verbatim}
This function is called at every element of the output array, passing the 
current coordinates in the \var{output_coordinates} array. On return, the 
\var{input_coordinates} array must contain the coordinates at which the 
input is interpolated. The ranks of the input and output array are passed 
through \var{output_rank} and \var{input_rank}. The value of the shift is 
passed through the \var{callback_data} argument, which is a pointer to 
void. The function returns an error status, in this case always 1, since no 
error can occur.

A pointer to this function and a pointer to the shift value must be passed 
to \function{geometric_transform}. Both are passed by a single CObject 
which is created by the following python extension function:
\begin{verbatim}
static PyObject *
py_shift_function(PyObject *obj, PyObject *args)
{
  double shift = 0.0;
  if (!PyArg_ParseTuple(args, "d", &shift)) {
    PyErr_SetString(PyExc_RuntimeError, "invalid parameters");
    return NULL;
  } else {
    /* assign the shift to a dynamically allocated location: */
    double *cdata = (double*)malloc(sizeof(double));
    *cdata = shift;
    /* wrap function and callback_data in a CObject: */
    return PyCObject_FromVoidPtrAndDesc(_shift_function, cdata,
                                        _destructor);
  }
}
\end{verbatim}
The value of the shift is obtained and then assigned to a dynamically 
allocated memory location. Both this data pointer and the function pointer 
are then wrapped in a CObject, which is returned. Additionally, a pointer 
to a destructor function is given, that will free the memory we allocated 
for the shift value when the CObject is destroyed. This destructor is very 
simple:
\begin{verbatim}
static void
_destructor(void* cobject, void *cdata)
{
  if (cdata)
    free(cdata);
}
\end{verbatim}
To use these functions, an extension module is build:
\begin{verbatim}
static PyMethodDef methods[] = {
  {"shift_function", (PyCFunction)py_shift_function, METH_VARARGS, ""},
  {NULL, NULL, 0, NULL}
};

void
initexample(void)
{
  Py_InitModule("example", methods);
}
\end{verbatim}
This extension can then be used in Python, for example:
\begin{verbatim}
>>> import example
>>> array = arange(12, shape=(4,3), type = Float64)
>>> fnc = example.shift_function(0.5)
>>> print geometric_transform(array, fnc)
[[ 0.      0.      0.    ]
 [ 0.      1.3625  2.7375]
 [ 0.      4.8125  6.1875]
 [ 0.      8.2625  9.6375]]
\end{verbatim}

C Callback functions for use with \module{nd\_image} functions must all be 
written according to this scheme. The next section lists the 
\module{nd\_image} functions that acccept a C callback function and gives 
the prototype of the callback function.

\subsection{Functions that support C callback functions}
The \module{nd\_image} functions that support C callback functions are 
described here. Obviously, the prototype of the function that is provided 
to these functions must match exactly that what they expect. Therefore we 
give here the prototypes of the callback functions. All these callback 
functions accept a void \var{callback_data} pointer that must be wrapped in 
a CObject using the Python \cfunction{PyCObject_FromVoidPtrAndDesc} 
function, which can also accept a pointer to a destructor function to free 
any memory allocated for \var{callback_data}. If \var{callback_data} is not 
needed, \cfunction{PyCObject_FromVoidPtr} may be used instead. The callback 
functions must return an integer error status that is equal to zero if 
something went wrong, or 1 otherwise. If an error occurs, you should 
normally set the python error status with an informative message before 
returning, otherwise, a default error message is set by the calling 
function.

The function \function{generic_filter} (see section 
\ref{sec:ndimage:genericfilters}) accepts a callback function with the 
following prototype:
\begin{cfuncdesc}{int}{FilterFunction}{double *buffer, int filter_size,
double *return_value, void *callback_data} The calling function iterates 
over the elements of the input and output arrays, calling the callback 
function at each element. The elements within the footprint of the filter 
at the current element are passed through the \var{buffer} parameter, and 
the number of elements within the footprint through \var{filter_size}. The 
calculated valued should be returned in the \var{return_value} argument.
\end{cfuncdesc}

The function \function{generic_filter1d} (see section 
\ref{sec:ndimage:genericfilters}) accepts a callback function with the 
following prototype: 
\begin{cfuncdesc}{int}{FilterFunction1D}{double *input_line, int 
input_length, double *output_line, int output_length, void *callback_data} 
The calling function iterates over the lines of the input and output 
arrays, calling the callback function at each line. The current line is 
extended according to the border conditions set by the calling function, 
and the result is copied into the array that is passed through the 
\var{input_line} array. The length of the input line (after extension) is 
passed through \var{input_length}. The callback function should apply the 
1D filter and store the result in the array passed through 
\var{output_line}. The length of the output line is passed through 
\var{output_length}.
\end{cfuncdesc}

The function \function{geometric_transform} (see section 
\ref{sec:ndimage:interpolation}) expects a function with the following 
prototype: 
\begin{cfuncdesc}{int}{MapCoordinates}{int *output_coordinates, 
double* input_coordinates, int output_rank, int input_rank, 
void *callback_data} The calling function iterates over the elements of the 
output array, calling the callback function at each element. The 
coordinates of the current output element are passed through 
\var{output_coordinates}. The callback function must return the coordinates 
at which the input must be interpolated in \var{input_coordinates}. The 
rank of the input and output arrays are given by \var{input_rank} and 
\var{output_rank} respectively.
\end{cfuncdesc}