# Copyright (C) 2003-2005 Peter J. Verveer
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met: 
#
# 1. Redistributions of source code must retain the above copyright
#    notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above
#    copyright notice, this list of conditions and the following
#    disclaimer in the documentation and/or other materials provided
#    with the distribution.
#
# 3. The name of the author may not be used to endorse or promote
#    products derived from this software without specific prior
#    written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
# OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
# GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.      

import types
import math
import numarray
import _ni_support
import _nd_image

def correlate1d(input, weights, axis = -1, output = None,
            mode = "nearest", cval = 0.0, origin = 0, output_type = None):
    """Calculate a one-dimensional correlation along the given axis.

    The lines of the array along the given axis are correlated with
    the given weights. The weights parameter must be a one-dimensional
    sequence of numbers.  The type of the output can be forced by the
    type parameter. If equal to None, the output array type is equal
    to the input array type. An output array can optionally be
    provided. The origin parameter controls the placement of the
    filter. The mode parameter determines how the array borders are
    handled, where cval is the value when mode is equal to 'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    return _nd_image.correlate1d(input, weights, axis, output, mode,
                       cval, origin, _ni_support._type_to_num[output_type])


def correlate(input, weights, output = None, mode = "nearest", cval = 0.0,
              origin = 0, output_type = None):
    """Multi-dimensional correlation.

    The array is correlated with the given kernel. The type of the
    output can be forced by the type parameter. If equal to None, the
    output array type is equal to the input array type. An output
    array can optionally be provided.  The origin parameter controls
    the placement of the filter. The mode parameter determines how the
    array borders are handled, where cval is the value when mode is
    equal to 'constant'.
    """
    footprint = numarray.abs(weights) > 1e-15
    mode = _ni_support._extend_mode_to_code(mode)
    origin = _ni_support._normalize_sequence(origin, input)
    return _nd_image.correlate(input, weights, footprint, output, mode,
                       cval, origin, _ni_support._type_to_num[output_type])


def convolve1d(input, weights, axis = -1, output = None, mode = "nearest",
               cval = 0.0, origin = 0, output_type = None):
    """Calculate a one-dimensional convolution along the given axis.

    The lines of the array along the given axis are convoluted with
    the given weights. The weights parameter must be a one-dimensional
    sequence of numbers.  The type of the output can be forced by the
    type parameter. If equal to None, the output array type is equal
    to the input array type. An output array can optionally be
    provided. The origin parameter controls the placement of the
    filter. The mode parameter determines how the array borders are
    handled, where cval is the value when mode is equal to 'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    weights = weights[::-1]
    origin = -origin
    if not len(weights) & 1:
        origin -= 1
    return _nd_image.correlate1d(input, weights, axis, output, mode,
                      cval, origin, _ni_support._type_to_num[output_type])


def convolve(input, weights, output = None, mode = "nearest", cval = 0.0,
             origin = 0, output_type = None):
    """Multi-dimensional convolution.

    The array is convolved with the given kernel. The type of the
    output can be forced by the type parameter. If equal to None, the
    output array type is equal to the input array type. An output
    array can optionally be provided.  The origin parameter controls
    the placement of the filter. The mode parameter determines how the
    array borders are handled, where cval is the value when mode is
    equal to 'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    origin = _ni_support._normalize_sequence(origin, input)
    weights = numarray.asarray(weights)
    weights = weights[tuple([slice(None, None, -1)] * weights.rank)]
    for ii in range(len(origin)):
        origin[ii] = -origin[ii]
        if not weights.shape[ii] & 1:
            origin[ii] -= 1
    footprint = numarray.abs(weights) > 1e-15
    return _nd_image.correlate(input, weights, footprint, output, mode,
                    cval, origin, _ni_support._type_to_num[output_type])


def gaussian_filter1d(input, sigma, axis = -1, order = 0, output = None,
                      mode = "nearest", cval = 0.0, output_type = None):
    """One-dimensional Gaussian filter.

    The standard-deviation of the Gaussian filter is given by
    sigma. 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 data type of the output array
    can be specified. The data type of the output is equal to the
    input type if the output_type argument is equal to
    None. Optionally an output array can be provided.
    """
    if output_type == None:
        if isinstance(input, numarray.NumArray):
            output_type = input.type()
        else:
            output_type = numarray.Float64
    if isinstance(output_type, numarray.ComplexType):
        raise RuntimeError, "complex array types not supported"
    sd = float(sigma)
    # make the length of the filter equal to 4 times the standard
    # deviations:
    lw = int(4.0 * sd + 0.5)
    weights = [0.0] * (2 * lw + 1)
    weights[lw] = 1.0

    sum = 1.0
    sd = sd * sd
    # calculate the kernel:
    for ii in range(1, lw + 1):
        tmp = math.exp(-0.5 * float(ii * ii) / sd)
        weights[lw + ii] = tmp
        weights[lw - ii] = tmp
        sum += 2.0 * tmp
    for ii in range(2 * lw + 1):
        weights[ii] /= sum
    # implement first, second and third order derivatives:
    if order == 1 : # first derivative
        weights[lw] = 0.0
        for ii in range(1, lw + 1):
            x = float(ii)
            tmp = -x / sd * weights[lw + ii]
            weights[lw + ii] = tmp
            weights[lw - ii] = -tmp
    elif order == 2: # second derivative
        weights[lw] = -1.0
        for ii in range(1, lw + 1):
            x = float(ii)
            tmp = (2.0 * x * x / sd - 1.0) * weights[lw + ii]
            weights[lw + ii] = tmp
            weights[lw - ii] = tmp
    elif order == 3: # third derivative
        weights[lw] = 0.0
        sd2 = sd * sd
        for ii in range(1, lw + 1):
            x = float(ii)
            tmp = (6.0 - 4.0 * x * x / sd) * x * weights[lw + ii]
            weights[lw + ii] = tmp
            weights[lw - ii] = -tmp
    return correlate1d(input, weights, axis, origin = 0, mode = mode,
                   cval = cval, output_type = output_type, output = output)


def gaussian_filter(input, sigma, order = 0, output = None,
                    mode = "nearest", cval = 0.0, output_type = None):
    """Multi-dimensional Gaussian filter.

    The standard-deviations of the Gaussian filter are given for each
    axis as a sequence, or as a single number, in which case it is
    equal for all axes. The order of the filter along each axis is
    given as a sequence of integers, or as a single number. 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 data type of the output array can be
    specified. The data type of the output is equal to the input type
    if the output_type argument is equal to None. Optionally an output
    array can be provided.

    Note: The multi-dimensional filter is implemented as a sequence of
    one-dimensional convolution filters. The intermediate arrays are
    stored in the same data type as the output. Therefore, for output
    types with a limited precision, the results may be imprecise
    because intermediate results may be stored with insufficient
    precision.
    """
    sigma = _ni_support._normalize_sequence(sigma, input)
    order = _ni_support._normalize_sequence(order, input)
    inplace = isinstance(output, numarray.NumArray)
    if output_type == None:
        if isinstance(input, numarray.NumArray):
            output_type = input.type()
        else:
            output_type = numarray.Float64
    if isinstance(output_type, numarray.ComplexType):
        raise RuntimeError, "complex array types not supported"
    done = False
    # loop over the dimensions
    for axis in range(len(sigma)):
        sd = float(sigma[axis])
        # small values of the standard deviation are skipped:
        if sd < 1e-15:
            continue
        if done:
            # if one or more axes are already processed, work in-place:
            gaussian_filter1d(output, sd, axis, order[axis], output,
                              mode, cval, output_type)
        else:
            if inplace:
                gaussian_filter1d(input, sd, axis, order[axis], output,
                                  mode, cval, output_type)
            else:
                output = gaussian_filter1d(input, sd, axis, order[axis],
                                           None, mode, cval, output_type)
        done = True
    if not done: # if no filter was applied, make a copy:
        if inplace:
            output[...] = input.astype(output_type)[...]
        else:
            output = input.copy().astype(output_type)
    if not inplace:
        return output
        

def boxcar_filter1d(input, size, axis = -1, output = None,
                    mode = "nearest", cval = 0.0, origin = 0,
                    output_type = None):
    """Calculates a one-dimensional boxcar filter along the given axis.

    The lines of the array along the given axis are filtered by a
    boxcar filter. The type of the output can be forced by the type
    parameter. If equal to None, the output array type is equal to the
    input array type. An output array can optionally be provided. The
    origin parameter controls the placement of the filter. The mode
    parameter determines how the array borders are handled, where
    cval is the value when mode is equal to 'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    return _nd_image.boxcar_filter1d(input, size, axis, output, mode,
                      cval, origin, _ni_support._type_to_num[output_type])


def boxcar_filter(input, size = 3, output = None, mode = "nearest",
                  cval = 0.0, origin = 0, output_type = None):
    """Multi-dimensional boxcar filter.

    The sizes of the boxcar filter are given for each axis as a
    sequence, or as a single number, in which case the size is equal
    for all axes. The data type of the output array can be
    specified. The data type of the output is equal to the input type
    if the output_type argument is equal to None. Optionally an output
    array can be provided. The origin parameter controls the placement
    of the filter. The mode parameter determines how the array borders
    are handled, where cval is the value when mode is equal to
    'constant'.

    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 limited precision, the results may be imprecise because
    intermediate results may be stored with insufficient precision.
    """
    inplace = isinstance(output, numarray.NumArray)
    if output_type == None:
        if isinstance(input, numarray.NumArray):
            output_type = input.type()
        else:
            output_type = numarray.Float64
    if isinstance(output_type, numarray.ComplexType):
        raise RuntimeError, "complex array types not supported"
    origin = _ni_support._normalize_sequence(origin, input)
    size = _ni_support._normalize_sequence(size, input)
    done = False
    # loop over the dimensions
    for axis in range(len(size)):
        fs = int(size[axis])
        # zero filter sizes are skipped:
        if fs <= 1:
            continue
        if done:
            # if one or more axes are already processed, work in-place:
            boxcar_filter1d(output, fs, axis, origin = origin[axis],
                            mode = mode, cval = cval,
                            output_type = None, output = output)
        else:
            if inplace:
                boxcar_filter1d(input, fs, axis, origin = origin[axis],
                            mode = mode, cval = cval,
                            output_type = output_type, output = output)
            else:
                output = boxcar_filter1d(input, fs, axis,
                                origin = origin[axis], mode = mode,
                                cval = cval, output_type = output_type)
        done = True
    if not done: # if no filter was applied, make a copy:
        if inplace:
            output[...] = numarray.array(input).astype(output_type)[...]
        else:
            output = numarray.array(input).copy().astype(output_type)
    if not inplace:
        return output
    

def minimum_filter1d(input, size, axis = -1, output = None,
                     mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a minimum filter along the given axis.

    The lines of the array along the given axis are processed by a
    one-dimensional minimum filter of given size. An output array can
    optionally be provided. The origin parameter controls the placement
    of the filter. The mode parameter determines how the array borders
    are handled, where cval is the value when mode is equal to
    'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    return _nd_image.minimum_maximum_filter1d(input, size, axis, output,
                                              mode, cval, origin, 1)


def maximum_filter1d(input, size, axis = -1, output = None,
                     mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a maximum filter along the given axis.

    The lines of the array along the given axis are processed by a
    one-dimensional maximum filter of given size. An output array can
    optionally be provided. The origin parameter controls the placement
    of the filter. The mode parameter determines how the array borders
    are handled, where cval is the value when mode is equal to
    'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    return _nd_image.minimum_maximum_filter1d(input, size, axis, output,
                                              mode, cval, origin, 0)


def _minimum_or_maximum_filter(do_minimum, input,  size = None,
                               footprint = None, structure = None,
                               output = None, mode = "nearest",
                               cval = 0.0, origin = 0):
    """Low level code for minimum or maximum filter"""
    if structure == None:
        if footprint == None:
            if size == None:
               raise RuntimeError, "no structure provided" 
            separable= True
        else:
            footprint = numarray.array(footprint, numarray.Bool)
            if numarray.alltrue(numarray.ravel(footprint)):
                size = footprint.shape
                footprint = None
                separable = True
            else:
                separable = False
    else:
        separable = False
        if footprint == None:
            footprint = numarray.ones(structure.shape, numarray.Bool)
    inplace = isinstance(output, numarray.NumArray)
    origin = _ni_support._normalize_sequence(origin, input)
    if separable:
        size = _ni_support._normalize_sequence(size, input)
        done = False
        # loop over the dimensions
        for axis in range(len(size)):
            fs = int(size[axis])
            # small values of the standard deviation are skipped:
            if fs <= 1:
                continue
            if done:
                # if one or more axes are already processed, 
                # work in-place:
                if do_minimum:
                    minimum_filter1d(output, fs, axis,
                                     origin = origin[axis], mode = mode,
                                     cval = cval, output = output)
                else:
                     maximum_filter1d(output, fs, axis,
                                      origin = origin[axis], mode = mode,
                                      cval = cval, output = output) 
            else:
                if inplace:
                    if do_minimum:
                        minimum_filter1d(input, fs, axis,
                                    origin = origin[axis], mode = mode,
                                    cval = cval, output = output)
                    else:
                        maximum_filter1d(input, fs, axis,
                                    origin = origin[axis], mode = mode,
                                    cval = cval, output = output)
                else:
                    if do_minimum:
                        output = minimum_filter1d(input, fs, axis,
                                                origin = origin[axis],
                                                mode = mode, cval = cval)
                    else:
                        output = maximum_filter1d(input, fs, axis,
                                                origin = origin[axis],
                                                mode = mode, cval = cval)
            done = True
        if not done: # if no filter was applied, make a copy:
            if inplace:
                output[...] = input.astype(output_type)[...]
            else:
                output = input.copy().astype(output_type)
    else:
        mode = _ni_support._extend_mode_to_code(mode)
        if do_minimum:
            return _nd_image.minimum_maximum_filter(input, footprint,
                                                    structure, output,
                                                    mode, cval, origin, 1)
        else:
            return _nd_image.minimum_maximum_filter(input, footprint,
                                                    structure, output,
                                                    mode, cval, origin, 0)
    if not inplace:
        return output
    

def minimum_filter(input,  size = None, footprint = None, output = None,
                   mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a multi-dimensional minimum filter.
    
    Either a size or a footprint with the filter must be
    provided. An output array can optionally be provided. The origin
    parameter controls the placement of the filter. The mode parameter
    determines how the array borders are handled, where cval is the
    value when mode is equal to 'constant'.
    """
    return _minimum_or_maximum_filter(True, input, size, footprint, None,
                                origin = origin, mode = mode, cval = cval)


def maximum_filter(input,  size = None, footprint = None, output = None,
                   mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a multi-dimensional maximum filter.
    
    Either a size or a footprint with the filter must be
    provided. An output array can optionally be provided. The origin
    parameter controls the placement of the filter. The mode parameter
    determines how the array borders are handled, where cval is the
    value when mode is equal to 'constant'.
    """
    return _minimum_or_maximum_filter(False, input, size, footprint, None,
                                origin = origin, mode = mode, cval = cval)


def rank_filter(input, rank, size = None, footprint = None, output = None,
                mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a multi-dimensional rank filter.
    
    The rank parameter may be less then zero, i.e., rank = -1
    indicates the larges element. Either a size or a footprint with
    the filter must be provided. An output array can optionally be
    provided. The origin parameter controls the placement of the
    filter. The mode parameter determines how the array borders are
    handled, where cval is the value when mode is equal to 'constant'.
    """
    if rank == 0:
        return minimum_filter(input, size, footprint, origin = origin,
                              mode = mode, cval = cval, output = output)
    if rank == -1:
        return maximum_filter(input, size, footprint, origin = origin,
                              mode = mode, cval = cval, output = output)
    if ((footprint == None and size == None) or
        (footprint != None and size != None)):
        msg = "filter dimensions or filter shape must be given"
        raise RuntimeError, msg
    if footprint == None:
        size = _ni_support._normalize_sequence(size, input)
        footprint = numarray.ones(size, numarray.Bool)
    mode = _ni_support._extend_mode_to_code(mode)
    origin = _ni_support._normalize_sequence(origin, input)
    return _nd_image.rank_filter(input, rank, 0.0, footprint, output,
                                 mode, cval, origin, 0)


def percentile_filter(input, percentile, size = None, footprint = None,
                output = None, mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a multi-dimensional percentile filter.
    
    The percentile parameter may be less then zero, i.e., percentile =
    -20 equals percentile = 80. Either a size or a footprint with the
    filter must be provided. An output array can optionally be
    provided. The origin parameter controls the placement of the
    filter. The mode parameter determines how the array borders are
    handled, where cval is the value when mode is equal to 'constant'.
    """
    if percentile == 0.0:
        return minimum_filter(input, size, footprint, origin = origin,
                              mode = mode, cval = cval,output = output)
    if percentile == 100.0:
        return maximum_filter(input, size, footprint, origin= origin,
                              mode = mode, cval = cval, output = output)
    if ((footprint == None and size == None) or
        (footprint != None and size != None)):
        msg = "filter dimensions or filter shape must be given"
        raise RuntimeError, msg
    if footprint == None:
        size = _ni_support._normalize_sequence(size, input)
        footprint = numarray.ones(size, numarray.Bool)
    mode = _ni_support._extend_mode_to_code(mode)
    origin = _ni_support._normalize_sequence(origin, input)
    return _nd_image.rank_filter(input, 0, percentile, footprint, output,
                                 mode, cval, origin, 1)


def median_filter(input, size = None, footprint = None, output = None,
                  mode = "nearest", cval = 0.0, origin = 0):
    """Calculates a multi-dimensional median filter.
    
    Either a size or a footprint with the filter must be provided. An
    output array can optionally be provided. The origin parameter
    controls the placement of the filter. The mode parameter
    determines how the array borders are handled, where cval is the
    value when mode is equal to 'constant'.
    """
    if ((footprint == None and size == None) or
        (footprint != None and size != None)):
        msg = "filter dimensions or filter shape must be given"
        raise RuntimeError, msg
    if footprint == None:
        size = _ni_support._normalize_sequence(size, input)
        footprint = numarray.ones(size, numarray.Bool)
    mode = _ni_support._extend_mode_to_code(mode)
    origin = _ni_support._normalize_sequence(origin, input)
    return _nd_image.rank_filter(input, 0, 50.0, footprint, output, mode,
                                 cval, origin, 1)


def prewitt(input, axis = -1, output = None, mode = "nearest", cval = 0.0,
            output_type = None):
    """Calculate a Prewitt filter.

    The type of the output can be forced by the type parameter. If
    equal to None, the output array type is equal to the input array
    type. An output array can optionally be provided. The mode
    parameter determines how the array borders are handled, where cval
    is the value when mode is equal to 'constant'.
    """
    inplace = isinstance(output, numarray.NumArray)
    input = numarray.asarray(input)
    if axis < 0:
        axis += input.rank
    if axis < 0 or axis >= input.rank:
        raise RuntimeError, 'invalid axis specification'
    if inplace:
        correlate1d(input, [1, 0, -1], axis, output, mode, cval, 0,
                    output_type)
    else:
        output = correlate1d(input, [1, 0, -1], axis, None, mode, cval, 0,
                             output_type)
    for ii in range(input.rank):
        if ii != axis:
            correlate1d(output, [1, 1, 1], ii, output, mode, cval, 0,
                        output_type)
    if not inplace:
        return output


def sobel(input, axis= - 1, output = None, mode = "nearest", cval = 0.0,
          output_type = None):
    """Calculate a Prewitt filter.

    The type of the output can be forced by the type parameter. If
    equal to None, the output array type is equal to the input array
    type. An output array can optionally be provided. The mode
    parameter determines how the array borders are handled, where cval
    is the value when mode is equal to 'constant'.
    """
    inplace = isinstance(output, numarray.NumArray)
    input = numarray.asarray(input)
    if axis < 0:
        axis += input.rank
    if axis < 0 or axis >= input.rank:
        raise RuntimeError, 'invalid axis specification'
    if inplace:
        correlate1d(input, [1, 0, -1], axis, output, mode, cval, 0,
                    output_type)
    else:
        output = correlate1d(input, [1, 0, -1], axis, None, mode, cval, 0,
                             output_type)
    for ii in range(input.rank):
        if ii != axis:
            correlate1d(output, [1, 2, 1], ii, output, mode, cval, 0,
                        output_type)
    if not inplace:
        return output


def generic_laplace(input, derivative2, output = None,
                    mode = "nearest", cval = 0.0, output_type = None):
    """Calculate a multidimensional laplace filter using the provided
    second derivative function.

    The type of the output can be forced by the type parameter. If
    equal to None, the output array type is equal to the input array
    type. An output array can optionally be provided. The mode
    parameter determines how the array borders are handled, where cval
    is the value when mode is equal to 'constant'. The derivative2
    parameter must be a callable with the following signature:

    derivative2(input, axis, output, mode, cval, output_type)
    """
    inplace = isinstance(output, numarray.NumArray)
    if inplace:
        derivative2(input, 0, output, mode, cval, output_type)
    else:
        output = derivative2(input, 0, None, mode, cval, output_type)
    for ii in range(1, input.rank):
        tmp = derivative2(input, ii, None, mode, cval, output_type)
        output += tmp
    if not inplace:
        return output
        

def laplace(input, output = None, mode = "nearest", cval = 0.0,
            output_type = None):
    """Calculate a multidimensional laplace filter.

    The type of the output can be forced by the type parameter. If
    equal to None, the output array type is equal to the input array
    type. An output array can optionally be provided. The mode
    parameter determines how the array borders are handled, where cval
    is the value when mode is equal to 'constant'.
    """
    def derivative2(input, axis, output, mode, cval, output_type):
        return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0,
                           output_type)
    return generic_laplace(input, derivative2, output, mode, cval,
                           output_type)


def gaussian_laplace(input, sigma, output = None, mode = "nearest",
                     cval = 0.0, output_type = None):
    """Calculate a multidimensional laplace filter using gaussian
    second derivatives.

    The standard-deviations of the Gaussian filter are given for each
    axis as a sequence, or as a single number, in which case it is
    equal for all axes.  The type of the output can be forced by the
    type parameter. If equal to None, the output array type is equal
    to the input array type. An output array can optionally be
    provided. The mode parameter determines how the array borders are
    handled, where cval is the value when mode is equal to 'constant'.
    """
    sigma = _ni_support._normalize_sequence(sigma, input)
    def derivative2(input, axis, output, mode, cval, output_type):
        input = numarray.asarray(input)
        order = [0] * input.rank
        order[axis] = 2
        return gaussian_filter(input, sigma, order, output, mode, cval,
                               output_type)
    return generic_laplace(input, derivative2, output, mode, cval,
                           output_type)


def generic_gradient_magnitude(input, derivative, output = None,
                               mode = "nearest", cval = 0.0,
                               output_type = None):
    """Calculate a gradient magnitude using the provdide function for
    the gradient.

    The type of the output can be forced by the type parameter. If
    equal to None, the output array type is equal to the input array
    type. An output array can optionally be provided. The mode
    parameter determines how the array borders are handled, where cval
    is the value when mode is equal to 'constant'.  The derivative
    parameter must be a callable with the following signature:
                    
    derivative(input, axis, output, mode, cval, output_type)
    """
    inplace = isinstance(output, numarray.NumArray)
    if inplace:
        derivative(input, 0, output, mode, cval, output_type)
    else:
        output = derivative(input, 0, None, mode, cval, output_type)
    numarray.multiply(output, output, output)
    for ii in range(1, input.rank):
        tmp = derivative(input, ii, None, mode, cval, output_type)
        numarray.multiply(tmp, tmp, tmp)
        output += tmp
    numarray.sqrt(output, output)    
    if not inplace:
        return output


def gaussian_gradient_magnitude(input, sigma, output = None,
                                mode = "nearest", cval = 0.0,
                                output_type = None):
    """Calculate a multidimensional gradient magnitude using gaussian
    derivatives.

    The standard-deviations of the Gaussian filter are given for each
    axis as a sequence, or as a single number, in which case it is
    equal for all axes.  The type of the output can be forced by the
    type parameter. If equal to None, the output array type is equal
    to the input array type. An output array can optionally be
    provided. The mode parameter determines how the array borders are
    handled, where cval is the value when mode is equal to 'constant'.
    """
    sigma = _ni_support._normalize_sequence(sigma, input)
    def derivative(input, axis, output, mode, cval, output_type):
        input = numarray.asarray(input)
        order = [0] * input.rank
        order[axis] = 1
        return gaussian_filter(input, sigma, order, output, mode, cval,
                               output_type)
    return generic_gradient_magnitude(input, derivative, output, mode,
                                      cval, output_type)
                                      
def generic_filter1d(input, function, filter_size, axis = -1,
                     output = None, mode = "nearest", cval = 0.0,
                     origin = 0, output_type = None):
    """Calculate a one-dimensional filter along the given axis.

    The function iterates over the lines of the array, calling the given
    function at each line. The arguments of the line are the input line,
    the output line. The input and output lines are 1D double arrays. The 
    input line is extended appropiately according the filter size and 
    origin. The output line must be modified in-place with the result. The 
    type of the output can be forced by the type parameter. If equal to 
    None, the output array type is equal to the input array type. An
    output array can optionally be provided. The origin parameter controls 
    the placement of the filter. The mode parameter determines how the 
    array borders are handled, where cval is the value when mode is equal 
    to 'constant'.
    """
    mode = _ni_support._extend_mode_to_code(mode)
    return _nd_image.generic_filter1d(input, function, filter_size, axis,
         output, mode, cval, origin, _ni_support._type_to_num[output_type])

def generic_filter(input, function, size = None, footprint = None,
                   output = None, mode = "nearest", cval = 0.0, 
                   origin = 0, output_type = None):
    """Calculates a multi-dimensional filter using the given function.
    
    At each element the provided function is called. The input values
    within the filter footprint at that element are passed to the function
    as a 1D array of Float64 values.
    current element are passed as a tuple.
       
    Either a size or a footprint with the filter must be provided. The 
    type of the output can be forced by the type parameter. If equal to 
    None, the output array type is equal to the input array type. An
    output array can optionally be provided. The origin parameter
    controls the placement of the filter. The mode parameter
    determines how the array borders are handled, where cval is the
    value when mode is equal to 'constant'.
    """
    if ((footprint == None and size == None) or
        (footprint != None and size != None)):
        msg = "filter dimensions or filter shape must be given"
        raise RuntimeError, msg
    if footprint == None:
        size = _ni_support._normalize_sequence(size, input)
        footprint = numarray.ones(size, numarray.Bool)
    mode = _ni_support._extend_mode_to_code(mode)
    origin = _ni_support._normalize_sequence(origin, input)
    return _nd_image.generic_filter(input, function, footprint, output,
                mode, cval, origin, _ni_support._type_to_num[output_type])


                                      
