# Copyright (C) 2003 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 numarray
import _ni_support
import _nd_image
import filters

def iterate_structure(structure, iterations, origin = None):
    """Iterate a structure by dilating it with itself.

    If origin is None, only the iterated structure is returned. If
    not, a tuple of the iterated structure and the modified origin is
    returned.
    """

    structure = numarray.inputarray(structure)
    if iterations < 2:
        return structure.copy()
    ni = iterations - 1
    shape = [ii + ni * (ii - 1) for ii in structure.shape]
    pos = [ni * (structure.shape[ii] / 2) for ii in range(len(shape))]
    slc = [slice(pos[ii], pos[ii] + structure.shape[ii], None)
           for ii in range(len(shape))]
    out = numarray.zeros(shape, numarray.Bool)
    out[slc] = structure != 0
    out = binary_dilation(out, structure, iterations = ni)
    if origin == None:
        return out
    else:
        origin = _ni_support._normalize_sequence(origin, structure)
        origin = [iterations * o for o in origin]
        return out, origin

def generate_binary_structure(rank, connectivity):
    """Generate a binary structure for binary morphological operations.

    The inputs are the rank of the array to which the structure will
    be applied and the square of the connectivity of the structure.

    """
    if connectivity < 1:
        connectivity = 1
    if rank < 1:
        if connectivity < 1:
            return 0
        else:
            return 1
    output = numarray.zeros([3] * rank, numarray.Bool)
    output = numarray.abs(numarray.indices([3] * rank) - 1)
    output = numarray.add.reduce(output, 0)
    return output <= connectivity


def binary_erosion(input, structure = None, iterations = 1, mask = None,
                   output = None, border_value = 0, origin = 0):
    """Multi-dimensional binary erosion with the given structure.
    
    An output array can optionally be provided. The origin parameter
    controls the placement of the filter. If no structuring element is
    provided an element is generated with a squared connectivity equal
    to one. The border_value parameter gives the value of the array
    outside the border. The erosion operation is repeated iterations
    times. If iterations is less than 1, the erosion is repeated until
    the result does not change anymore. If a mask is given, only those
    elements with a true value at the corresponding mask element are
    modified at each iteration.

    """
    if structure == None:
        structure = generate_binary_structure(_ni_support._array_rank(input),
                                          1)
    origin = _ni_support._normalize_sequence(origin, input)
    inplace = isinstance(output, numarray.NumArray)
    if iterations == 1:
        result = _nd_image.binary_erosion(input, structure, mask, output,
                                          border_value, origin, 0, 0)
        if inplace:
            return None
        else:
            return result[0]
    else:
        tmp_in = numarray.zeros(input.shape, numarray.Bool)
        if inplace:
            tmp_out = output
        else:
            tmp_out = numarray.zeros(input.shape, numarray.Bool)
        if not iterations & 1:
            tmp_in, tmp_out = tmp_out, tmp_in
        changed = _nd_image.binary_erosion(input, structure, mask, tmp_out,
                                           border_value, origin, 0, 0)
        ii = 1
        while (ii < iterations) or (iterations < 1) and changed:
            tmp_in, tmp_out = tmp_out, tmp_in
            changed = _nd_image.binary_erosion(tmp_in, structure, mask,
                                               tmp_out, border_value,
                                               origin, 0, 0)
            ii += 1
        if not inplace:
            return tmp_out


def binary_dilation(input, structure = None, iterations = 1, mask = None,
                    output = None, border_value = 0, origin = 0):
    """Multi-dimensional binary dilation with the given structure.
    
    An output array can optionally be provided. The origin parameter
    controls the placement of the filter. If no structuring element is
    provided an element is generated with a squared connectivity equal
    to one. The erosion operation is repeated iterations times.  If
    iterations is less than 1, the dilation is repeated until the
    result does not change anymore.  If a mask is given, only those
    elements with a true value at the corresponding mask element are
    modified at each iteration.

    """

    if structure == None:
        structure = generate_binary_structure(_ni_support._array_rank(input),
                                              1)
    origin = _ni_support._normalize_sequence(origin, input)
    structure = numarray.inputarray(structure)
    structure = structure[tuple([slice(None, None, -1)] *
                                len(structure.shape))]
    for ii in range(len(origin)):
        origin[ii] = -origin[ii]
        if not structure.shape[ii] & 1:
            origin[ii] -= 1
    # implemented by eroding the background with the reflected structure:
    inplace = isinstance(output, numarray.NumArray)
    if iterations == 1:
        result = _nd_image.binary_erosion(input, structure, mask, output,
                                          border_value, origin, 1, 1)
        if inplace:
            return None
        else:
            return result[0]
    else:
        tmp_in = numarray.zeros(input.shape, numarray.Bool)
        if inplace:
            tmp_out = output
        else:
            tmp_out = numarray.zeros(input.shape, numarray.Bool)
        if not iterations & 1:
            tmp_in, tmp_out = tmp_out, tmp_in
        changed = _nd_image.binary_erosion(input, structure, mask,
                                           tmp_out, border_value,
                                           origin, 1, 1)
        ii = 1
        while (ii < iterations) or (iterations < 1) and changed:
            tmp_in, tmp_out = tmp_out, tmp_in
            changed = _nd_image.binary_erosion(tmp_in, structure, mask,
                                               tmp_out, border_value,
                                               origin, 1, 1)
            ii += 1
        if not inplace:
            return tmp_out


def binary_opening(input, structure = None, iterations = 1, output = None,
                   origin = 0):
    """Multi-dimensional binary opening with the given structure.

    An output array can optionally be provided. The origin parameter
    controls the placement of the filter. If no structuring element is
    provided an element is generated with a squared connectivity equal
    to one. The iterations parameter gives the number of times the
    erosions and then the dilations are done.

    """
    if structure == None:
        structure = generate_binary_structure(_ni_support._array_rank(input),
                                              1)
    tmp = binary_erosion(input, structure, iterations, None, None, 0, origin)
    return binary_dilation(tmp, structure, iterations, None, output, 0, origin)


def binary_closing(input, structure = None, iterations = 1, output = None,
                   origin = 0):
    """Multi-dimensional binary closing with the given structure.

    An output array can optionally be provided. The origin parameter
    controls the placement of the filter. If no structuring element is
    provided an element is generated with a squared connectivity equal
    to one. The iterations parameter gives the number of times the
    erosions and then the dilations are done.

    """
    if structure == None:
        structure = generate_binary_structure(_ni_support._array_rank(input),
                                              1)
    tmp = binary_dilation(input, structure, iterations, None, None, 0, origin)
    return binary_erosion(tmp, structure, iterations, None, output, 0, origin)


def binary_hit_or_miss(input, structure1 = None, structure2 = None,
                       output = None, origin1 = 0, origin2 = None):
    """Multi-dimensional binary hit-or-miss transform with given structures.
    
    An output array can optionally be provided. The origin parameters
    controls the placement of the structuring elements. If the first
    structuring element is not given one is generated with a squared
    connectivity equal to one. If the second structuring element is
    not provide, it set equal to the inverse of the first structuring
    element. If the origin for the second structure is equal to None
    it is set equal to the origin of the first.

    """
    if structure1 == None:
        structure1 = generate_binary_structure(_ni_support._array_rank(input),
                                               1)
    if structure2 == None:
        structure2 = numarray.logical_not(structure1)

    origin1 = _ni_support._normalize_sequence(origin1, input)
    if origin2 == None:
        origin2 = origin1
    else:
        origin2 = _ni_support._normalize_sequence(origin2, input)

    tmp1, changed = _nd_image.binary_erosion(input, structure1, None, None,
                                             0, origin1, 0, 0)
    inplace = isinstance(output, numarray.NumArray)
    result= _nd_image.binary_erosion(input, structure2, None, output,
                                     0, origin2, 1, 0)
    if inplace:
        numarray.logical_and(tmp1, output, output)
    else:
        return numarray.logical_and(tmp1, result[0])


def grey_erosion(input,  size = None, footprint = None, structure = None,
                 output = None, mode = "nearest", cval = 0.0, origin = 0):
    """Calculate a grey values erosion.

    ither a size or a footprint, or the structure 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 filters._minimum_or_maximum_filter(True, input, size, footprint,
                                              structure, output, mode, cval,
                                              origin)


def grey_dilation(input,  size = None, footprint = None, structure = None,
                  output = None, mode = "nearest", cval = 0.0, origin = 0):
    """Calculate a grey values dilation.
    
    Either a size or a footprint, or the structure 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 structure != None:
        structure = numarray.inputarray(structure)
        structure = structure[tuple([slice(None, None, -1)] *
                                    len(structure.shape))]
    if footprint != None:
        footprint = numarray.inputarray(footprint)
        footprint = footprint[tuple([slice(None, None, -1)] *
                                    len(footprint.shape))]
    origin = _ni_support._normalize_sequence(origin, input)
    for ii in range(len(origin)):
        origin[ii] = -origin[ii]
        if footprint != None:
            sz = footprint.shape[ii]
        else:
            sz = size[ii]
        if not sz & 1:
            origin[ii] -= 1
    return filters._minimum_or_maximum_filter(False, input, size, footprint,
                                              structure, output, mode, cval,
                                              origin)


def grey_opening(input, size = None, footprint = None, structure = None,
                 output = None, mode = "nearest", cval = 0.0, origin = 0):
    """Multi-dimensional grey valued opening.

    Either a size or a footprint, or the structure 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'.

    """
    tmp = grey_erosion(input, size, footprint, structure, None, mode, cval,
                       origin)
    return grey_dilation(tmp, size, footprint, structure, output, mode, cval,
                         origin)


def grey_closing(input, size = None, footprint = None, structure = None,
                 output = None, mode = "nearest", cval = 0.0, origin = 0):
    """Multi-dimensional grey valued closing.

    Either a size or a footprint, or the structure 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'.

    """
    tmp = grey_dilation(input, size, footprint, structure, None, mode, cval,
                        origin)
    return grey_erosion(tmp, size, footprint, structure, output, mode, cval,
                        origin)


def morphological_gradient(input, size = None, footprint = None,
                           structure = None, output = None, mode = "nearest",
                           cval = 0.0, origin = 0):
    """Multi-dimensional morphological gradient.

    Either a size or a footprint, or the structure 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'.

    """
    tmp = grey_dilation(input, size, footprint, structure, None, mode, cval,
                        origin)
    if isinstance(output, numarray.NumArray):
        grey_erosion(input, size, footprint, structure, output, mode, cval,
                     origin)
        return numarray.subtract(tmp, output, output)
    else:
        return (tmp - grey_erosion(input, size, footprint, structure, None,
                                    mode, cval, origin))
    

def morphological_laplace(input, size = None, footprint = None,
                          structure = None, output = None, mode = "nearest",
                          cval = 0.0, origin = 0):
    """Multi-dimensional morphological laplace.

    Either a size or a footprint, or the structure 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'.

    """
    tmp1 = grey_dilation(input, size, footprint, structure, None, mode, cval,
                         origin)
    if isinstance(output, numarray.NumArray):
        grey_erosion(input, size, footprint, structure, output, mode, cval,
                     origin)
        numarray.add(tmp1, output, output)
        del tmp1
        numarray.subtract(output, input, output)
        return numarray.subtract(output, input, output)
    else:
        tmp2 = grey_erosion(input, size, footprint, structure, None, mode,
                            cval, origin)
        numarray.add(tmp1, tmp2, tmp2)
        del tmp1
        numarray.subtract(tmp2, input, tmp2)
        numarray.subtract(tmp2, input, tmp2)
        return tmp2


def white_tophat(input, size = None, footprint = None,
                 structure = None, output = None, mode = "nearest",
                 cval = 0.0, origin = 0):
    """Multi-dimensional white tophat filter.

    Either a size or a footprint, or the structure 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'.

    """
    tmp = grey_erosion(input, size, footprint, structure, None, mode, cval,
                       origin)
    if isinstance(output, numarray.NumArray):
        grey_dilation(tmp, size, footprint, structure, output, mode, cval,
                      origin)
        del tmp
        return numarray.subtract(input, output, output)
    else:
        tmp = grey_dilation(tmp, size, footprint, structure, None, mode,
                            cval, origin)
        return input - tmp


def black_tophat(input, size = None, footprint = None,
                 structure = None, output = None, mode = "nearest",
                 cval = 0.0, origin = 0):
    """Multi-dimensional black tophat filter.

    Either a size or a footprint, or the structure 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'.

    """
    tmp = grey_dilation(input, size, footprint, structure, None, mode, cval,
                        origin)
    if isinstance(output, numarray.NumArray):
        grey_erosion(tmp, size, footprint, structure, output, mode, cval,
                      origin)
        del tmp
        return numarray.subtract(output, input, output)
    else:
        tmp = grey_erosion(tmp, size, footprint, structure, None, mode,
                           cval, origin)
        return tmp - input
