\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
arbritrary dimensionality. The packages currently includes functions for linear
and non-linear filtering, binary morphology, B-spline interpolation, and object
measurements.

\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 arbritrary 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='nearest', 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='nearest', 
    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='nearest', 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 results in a shift in the opposite directions.}
\end{funcdesc}

\begin{funcdesc}{convolve}{input, weights, output=None, mode='nearest', 
    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: it 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='nearest', 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}. A filter 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.
\end{funcdesc}

\begin{funcdesc}{gaussian_filter}{input, sigma, order=0, output=None, 
    mode='nearest', 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}{boxcar_filter1d}{input, size, axis=-1, output=None,
    mode='nearest', cval=0.0, origin=0, output_type=None} The
  \function{boxcar_filter1d} function calculates a one-dimensional boxcar
  filter of the given \var{size} along the given axis.
\end{funcdesc}

\begin{funcdesc}{boxcar_filter}{input, size, output=None, mode='nearest', 
    cval=0.0, origin=0, output_type=None} The \function{boxcar_filter}
  implements a multi-dimensional boxcar filter.  The sizes of the boxcar 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 boxcar 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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., 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, mask=None, 
    output=None, mode='nearest', 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.,
  percentile=-20 equals 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, mask=None, output=None, 
    mode='nearest', 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='nearest', 
    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='nearest', 
    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, output=None, mode='nearest', 
    cval=0.0, output_type=None, derivative2=None}
  
  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, output_type)}
  
  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} and \var{output_type} have the usual meaning.
\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='nearest', 
    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, output=None, mode='nearest', 
    cval=0.0, output_type=None} The function \function{gaussian_laplace}
  calculates the Laplace using \function{gaussian_filter} to calculate the
  second derivatives.
\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, output=None, 
    mode='nearest', cval=0.0, output_type=None, derivative=None} 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, output_type)}
  
  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} and \var{output_type} have the usual meaning.
\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, output=None, 
    mode='nearest', 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.
\end{funcdesc}

\subsection{Fourier space 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_boxcar}{input, size, n=-1, axis=-1, output=None}
  The \function{fourier_boxcar} function multiplies the input array with the
  multi-dimensional Fourier transform of a boxcar 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 \var{output_type} 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 problems 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} The \function{geometric_transform} function applies an
  arbritrary 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.
\end{funcdesc}

\begin{funcdesc}{map_indices}{input, indices, output_type=None, output=None, 
    order=3, mode='constant', cval=0.0, prefilter=True} The function
  \function{map_indices} applies an arbritrary coordinate transformation using
  the given array of indices. The parameter \var{indices} is used to find for
  each point in the output the corresponding indices in the input. The values
  of \var{indices} along the first axis are the coordinates in the input array
  at which the output value is found. (See also the numarray \function{indices}
  function.) Since the indices may be non-integer coordinates, the value of the
  input at these indices is determined by spline interpolation of the requested
  order.

\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{Morphology}
\label{sec:ndimage:morphology}

\subsection{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 Euclidian 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
  arbritrary 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 of a dilation with a given
structure $n$ times is equivalent to an erosions or a dilation with a structure
that is $n - 1$ times dilated with itself. For a low number of iterations this
is more efficient, and for this reason 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 of 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 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_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}

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

Grey-scale morphology operations are the equivalents of binary morphology
operations that operate on arrays with arbritrary 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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='nearest', 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{Object measurements}
This section describes some basic functions for performing measurements on
objects whose shape is defined by some mask. For instance, assume we have some
two-dimensional array and a mask that represents the objects in the image:
\begin{verbatim}
[[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 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 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}

Given an array of labeled objects the \function{find_objects} functions 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.

\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}{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}

%% Local Variables:
%% mode: LaTeX
%% mode: auto-fill
%% fill-column: 79
%% indent-tabs-mode: nil
%% ispell-dictionary: "american"
%% reftex-fref-is-default: nil
%% TeX-auto-save: t
%% TeX-command-default: "pdfeLaTeX"
%% TeX-master: "numarray"
%% TeX-parse-self: t
%% End:
% LocalWords:  latexonly numarray ndimage Verveer cval gaussian Unser vol pp
% LocalWords:  prefilter prefiltered american reftex fref pdfeLaTeX Prewitt fft
% LocalWords:  Sobel prewitt sobel tophat laplace fourier
