"""
`reconstruction` originally part of CellProfiler, code licensed under both GPL
and BSD licenses.

Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky

"""
import numpy as np

from skimage.filter.rank_order import rank_order


def reconstruction(image, mask, selem=None, offset=None):
    """Perform a morphological reconstruction of an image.

    Reconstruction requires a "seed" image and a "mask" image. Currently, this
    only implements reconstruction by dilation, such that the seed image is
    dilated until it is constrained by the mask. Thus, he "seed" and "mask"
    images will be the minimum and maximum possible values of the reconstructed
    image, respectively.

    Parameters
    ----------
    image : ndarray
        The seed image; a.k.a. marker image.
    mask : ndarray
        The maximum allowed value at each point.
    selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.

    Returns
    -------
    reconstructed : ndarray
       The result of morphological reconstruction.

    Notes
    -----
    The algorithm is taken from:
    Robinson, "Efficient morphological reconstruction: a downhill filter",
    Pattern Recognition Letters 25 (2004) 1759-1767.

    Applications for greyscale reconstruction are discussed in:

    [1] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis:
        Applications and Efficient Algorithms", IEEE Transactions on Image
        Processing (1993)

    [2] Soille, P., "Morphological Image Analysis: Principles and Applications",
        Chapter 6, 2nd edition (2003), ISBN 3540429883.

    Examples
    --------
    Uses for greyscale reconstruction are described in Vincent (1993). For
    example, let's try to extract the features of an image by subtracting a
    background image created by reconstruction.

    First, create an image where the "bumps" are the features that
    we want to extract:

    >>> import numpy as np
    >>> from scikits.image.morphology.grey import grey_reconstruction
    >>> y, x = np.mgrid[:20:0.5, :20:0.5]
    >>> bumps = np.sin(x) + np.sin(y)

    To create the background image, set the mask image to the original image,
    and the seed image to the original image with an intensity offset, `h`.

    >>> h = 0.3
    >>> seed = bumps - h
    >>> rec = grey_reconstruction(seed, bumps)

    The resulting reconstructed image looks exactly like the original image,
    but with the peaks of the bumps cut off. Subtracting this reconstructed
    image from the original image leaves just the peaks of the bumps

    >>> hdome = bumps - rec

    This operation is known as the h-dome of the image, which leaves features
    of height `h` in the subtracted image. The h-dome transform, and its
    inverse h-basin, are analogous to the white top-hat and black top-hat
    transforms, but don't require a structuring element.

    """
    assert tuple(image.shape) == tuple(mask.shape)
    assert np.all(image <= mask)
    try:
        from ._greyreconstruct import reconstruction_loop
    except ImportError:
        raise ImportError("_greyreconstruct extension not available.")

    if selem is None:
        selem = np.ones([3]*image.ndim, dtype=bool)
    else:
        selem = selem.copy()

    if offset == None:
        if not all([d % 2 == 1 for d in selem.shape]):
            ValueError("Footprint dimensions must all be odd")
        offset = np.array([d / 2 for d in selem.shape])
    # Cross out the center of the selem
    selem[[slice(d, d + 1) for d in offset]] = False

    # Make padding for edges of reconstructed image so we can ignore boundaries
    padding = (np.array(selem.shape) / 2).astype(int)
    dims = np.zeros(image.ndim + 1, dtype=int)
    dims[1:] = np.array(image.shape) + 2 * padding
    dims[0] = 2
    inside_slices = [slice(p, -p) for p in padding]
    # Set padded region to minimum image intensity and mask along first axis so
    # we can interleave image and mask pixels when sorting.
    values = np.ones(dims) * np.min(image)
    values[[0] + inside_slices] = image
    values[[1] + inside_slices] = mask

    # Create a list of strides across the array to get the neighbors within
    # a flattened array
    value_stride = np.array(values.strides[1:]) / values.dtype.itemsize
    image_stride = values.strides[0] / values.dtype.itemsize
    selem_mgrid = np.mgrid[[slice(-o, d - o)
                            for d, o in zip(selem.shape, offset)]]
    selem_offsets = selem_mgrid[:, selem].transpose()
    nb_strides = np.array([np.sum(value_stride * selem_offset)
                           for selem_offset in selem_offsets], np.int32)
    values = values.flatten()
    value_sort = np.lexsort([-values]).astype(np.int32)

    # Make a linked list of pixels sorted by value. -1 is the list terminator.
    prev = -np.ones(len(values), np.int32)
    next = -np.ones(len(values), np.int32)
    prev[value_sort[1:]] = value_sort[:-1]
    next[value_sort[:-1]] = value_sort[1:]

    # Create a rank-order value array so that the Cython inner-loop
    # can operate on a uniform data type
    rec_img, value_map = rank_order(values)
    current = value_sort[0]

    reconstruction_loop(rec_img, prev, next, nb_strides, current, image_stride)

    # Reshape reconstructed image to original image shape and remove padding.
    rec_img = value_map[rec_img[:image_stride]]
    rec_img.shape = np.array(image.shape) + 2 * padding
    return rec_img[inside_slices]