def erosion(image, selem=None, out=None, shift_x=False, shift_y=False):
    """Return greyscale morphological erosion of an image.

    Morphological erosion sets a pixel at (i,j) to the minimum over all pixels
    in the neighborhood centered at (i,j). Erosion shrinks bright regions and
    enlarges dark regions.

    Parameters
    ----------
    image : ndarray
        Image array.
    selem : ndarray, optional
        The neighborhood expressed as an array of 1's and 0's.
        If None, use cross-shaped structuring element (connectivity=1).
    out : ndarrays, optional
        The array to store the result of the morphology. If None is
        passed, a new array will be allocated.
    shift_x, shift_y : bool, optional
        shift structuring element about center point. This only affects
        eccentric structuring elements (i.e. selem with even numbered sides).

    Returns
    -------
    eroded : array, same shape as `image`
        The result of the morphological erosion.

    Notes
    -----
    For ``uint8`` (and ``uint16`` up to a certain bit-depth) data, the
    lower algorithm complexity makes the `skimage.filter.rank.minimum`
    function more efficient for larger images and structuring elements.

    Examples
    --------
    >>> # Erosion shrinks bright regions
    >>> import numpy as np
    >>> from skimage.morphology import square
    >>> bright_square = np.array([[0, 0, 0, 0, 0],
    ...                           [0, 1, 1, 1, 0],
    ...                           [0, 1, 1, 1, 0],
    ...                           [0, 1, 1, 1, 0],
    ...                           [0, 0, 0, 0, 0]], dtype=np.uint8)
    >>> erosion(bright_square, square(3))
    array([[0, 0, 0, 0, 0],
           [0, 0, 0, 0, 0],
           [0, 0, 1, 0, 0],
           [0, 0, 0, 0, 0],
           [0, 0, 0, 0, 0]], dtype=uint8)

    """
    selem = np.array(selem)
    selem = _shift_selem(selem, shift_x, shift_y)
    if out is None:
        out = np.empty_like(image)
    ndi.grey_erosion(image, footprint=selem, output=out)
    return out

def fan_structure(median=11):
    """function to analyse fan sub structure by setting non-fan pixels to 0 and histo-equalizing the remaining img with only fans. Also
    median-filtering is done to reduce noise."""
    xcoords=[388, 449,497]
    ycoords =[142, 254, 118]
    x2 = [403,590]
    y2 = [286,254]
    x3 = [403,459]
    y3 = [286,375]
    x4 = [1372,1420]
    y4 = [610,680]
    x5 = [1372,1467]
    y5 = [610,590]
    x6 = [1321,1422]
    y6 = [612,750]
    x7 = [1321,1439]
    y7 = [612,585]
    fig = plt.figure()
    ax = fig.add_subplot(111)
    data = get_data('ESP_011931_0945')
    data = nd.grey_erosion(data,footprint=np.ones((3,3)))
    data = nd.grey_erosion(data,footprint=np.ones((3,3)))
    data = nd.grey_dilation(data,footprint=np.ones((3,3)))
    data = nd.grey_dilation(data,footprint=np.ones((3,3)))
    threshold=0.045
    fans = data < threshold
    data = data*255/data.max()
    intfans = np.zeros(data.shape, dtype=np.uint16)
    intfans[fans] = np.round(data[fans])
    filtered = nd.median_filter(intfans,median)
    equ = hist_equal(filtered)
    im = ax.imshow(equ,cmap = cm.spectral,aspect='equal')
    ax.set_title('Fans within fans in Ithaca, filtered, opened and hist-equalized')
    ax.set_xlabel('0.5 m/pixel')
#    fig.savefig('Fans_within_fans.png')
#    cb =plt.colorbar(im,shrink=0.75)
#    cb.set_label('I/F')
    plt.plot(xcoords[:-1],ycoords[:-1],[xcoords[0],xcoords[2]],[ycoords[0],ycoords[2]],
            color='white',
            hold=True,
            scalex=False,scaley=False)
    plt.plot(x2,y2,color='white',hold=True,scalex=False,scaley=False)
    plt.plot(x3,y3,color='white',hold=True,scalex=False,scaley=False)
    plt.plot(x4,y4,color='white',hold=True,scalex=False,scaley=False)
    plt.plot(x5,y5,color='white',hold=True,scalex=False,scaley=False)

    plt.plot(x6,y6,color='white',hold=True,scalex=False,scaley=False)
    plt.plot(x7,y7,color='white',hold=True,scalex=False,scaley=False)

#    plt.close(fig)
    plt.show()

def rolling_ball_filter(data, ball_radius, spacing=None, top=False, **kwargs):
    """Rolling ball filter implemented with morphology operations

    This implenetation is very similar to that in ImageJ and uses a top hat transform
    with a ball shaped structuring element
    https://en.wikipedia.org/wiki/Top-hat_transform

    Parameters
    ----------
    data : ndarray type uint8
        image data (assumed to be on a regular grid)
    ball_radius : float
        the radius of the ball to roll
    spacing : int or sequence
        the spacing of the image data
    top : bool
        whether to roll the ball on the top or bottom of the data
    kwargs : key word arguments
        these are passed to the ndimage morphological operations

    Returns
    -------
    data_nb : ndarray
        data with background subtracted as uint8
    bg : ndarray
        background that was subtracted from the data
    """
    data = data.astype(np.int16)
    ndim = data.ndim
    if spacing is None:
        spacing = _normalize_sequence(1, ndim)
    else:
        spacing = _normalize_sequence(spacing, ndim)

    radius = np.asarray(_normalize_sequence(ball_radius, ndim))
    mesh = np.array(np.meshgrid(*[np.arange(-r, r + s, s) for r, s in zip(radius, spacing)], indexing="ij"))
    structure = 2 * np.sqrt(2 - ((mesh / radius.reshape(-1, *((1,) * ndim)))**2).sum(0))
    structure[~np.isfinite(structure)] = 0
    if not top:
        # ndi.white_tophat(y, structure=structure, output=background)
        background = ndi.grey_erosion(data, structure=structure, **kwargs)
        background = ndi.grey_dilation(background, structure=structure, **kwargs)
    else:
        # ndi.black_tophat(y, structure=structure, output=background)
        background = ndi.grey_dilation(data, structure=structure, **kwargs)
        background = ndi.grey_erosion(background, structure=structure, **kwargs)

    data_corr = data - background
    data_corr[data_corr<0] = 0

    return data_corr.astype(np.uint8), background.astype(np.uint8)

def artifact_mask(imdata, airdata):
    """Computes a mask of artifacts found in the air region"""

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    hist, bin_edges = np.histogram(bg_img[bg_img > 0], bins=128)
    bg_threshold = np.mean(bin_edges[np.argmax(hist)])


    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > bg_threshold] = bg_img[bg_img > bg_threshold]

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 2)

    # Perform an a grayscale erosion operation.
    qi1_img = nd.grey_erosion(qi1_img, structure=struc).astype(np.float32)
    # Binarize and binary dilation
    qi1_img[qi1_img > 0.] = 1
    qi1_img[qi1_img < 1.] = 0
    qi1_img = nd.binary_dilation(qi1_img, structure=struc).astype(np.uint8)
    return qi1_img

 def generateData(self):
     input = self.getInput(0).getData()/255
     fill = ndimage.grey_erosion(input, size = (3,3))
     
     output = input - fill
     self.getOutput(0).setData(output*255)
     self.getOutput(1).setData(input*255)        

def multi_label_edge_detection(data):
    f = nd.generate_binary_structure(len(data.shape), 1)
    bound = (nd.grey_erosion(data,footprint=f) != nd.grey_dilation(data,footprint=f)) - \
            (nd.binary_dilation(data.astype(np.bool)) - data.astype(np.bool))     # the unwanted thick bounds
    data=bound.astype(data.dtype)
    
    return data

def edge_matrix(labels, connectivity=1):
    """Generate a COO matrix containing the coordinates of edge pixels.

    Parameters
    ----------
    labels : array of int
        An array of labeled pixels (or voxels).
    connectivity : int in {1, ..., labels.ndim}
        The square connectivity for considering neighborhood.

    Returns
    -------
    edges : sparse.coo_matrix
        A COO matrix where (i, j) indicate neighboring labels and the
        corresponding data element is the linear index of the edge pixel
        in the labels array.
    """
    conn = ndi.generate_binary_structure(labels.ndim, connectivity)
    eroded = ndi.grey_erosion(labels, footprint=conn).ravel()
    dilated = ndi.grey_dilation(labels, footprint=conn).ravel()
    labels = labels.ravel()
    boundaries0 = np.flatnonzero(eroded != labels)
    boundaries1 = np.flatnonzero(dilated != labels)
    labels_small = np.concatenate((eroded[boundaries0], labels[boundaries1]))
    labels_large = np.concatenate((labels[boundaries0], dilated[boundaries1]))
    n = np.max(labels_large) + 1
    data = np.concatenate((boundaries0, boundaries1))
    sparse_graph = sparse.coo_matrix((data, (labels_small, labels_large)),
                                     dtype=np.int_, shape=(n, n))
    return sparse_graph

def findendsjunctions(img, disp=None):
    if disp is None:
        disp = 0

    # Create a look up table to find junctions.
    # lut = ndimage.grey_erosion(junction(img), size=(3,3))
    junctions = ndimage.grey_erosion(img, size=(3, 3))
    # Row and column coordinates of junction points in the image.
    rjcj = np.where(junctions)

    ends = ndimage.grey_erosion(img, size=(3, 3))
    # Row and column coordinates of end points in the image.
    rece = np.where(ends)

    # Display image with the junctions and endings marked.
    if disp:
        plt.imshow(img)

def apply(array, **kwargs):
    """
    Apply a set of standard filter to array data: 
    
    Call: apply(array-data, <list of key=value arguments>)

    The list of key-value define the filtering to be done and should be given in
    the order to be process. Possible key-value are:
    
      * smooth:  gaussian filtering, value is the sigma parameter (scalar or tuple)
      * uniform: uniform  filtering (2)
      * max:     maximum  filtering (1)
      * min:     minimum  filtering (1)
      * median:  median   filtering (1)
      
      * dilate: grey dilatation (1)
      * erode:  grey erosion    (1)
      * close:  grey closing    (1)
      * open:   grey opening    (1)
      
      * linear_map: call linear_map(), value is the tuple (min,max)   (3)
      * normalize:  call normalize(),  value is the method            (3)
      * adaptive:   call adaptive(),   value is the sigma             (3)
      * adaptive_:  call adaptive(),   with uniform kernel            (3)
          
    The filtering is done using standard scipy.ndimage functions.
    
    (1) The value given (to the key) is the width of the the filter: 
        the distance from the center pixel (the size of the filter is thus 2*value+1)
        The neighborhood is an (approximated) boolean circle (up to discretization)
    (2) Same as (*) but the neighborhood is a complete square
    (3) See doc of respective function
    """
    for key in kwargs:
        value = kwargs[key]
        if key not in ('smooth','uniform'):
            fp = _kernel.distance(array.ndim*(2*value+1,))<=value  # circular filter
            
        if   key=='smooth' : array = _nd.gaussian_filter(array, sigma=value)
        elif key=='uniform': array = _nd.uniform_filter( array, size=2*value+1)
        elif key=='max'    : array = _nd.maximum_filter( array, footprint=fp)
        elif key=='min'    : array = _nd.minimum_filter( array, footprint=fp)
        elif key=='median' : array = _nd.median_filter(  array, footprint=fp)

        elif key=='dilate' : array = _nd.grey_dilation(  array, footprint=fp)
        elif key=='erode'  : array = _nd.grey_erosion(   array, footprint=fp)
        elif key=='open'   : array = _nd.grey_opening(   array, footprint=fp)
        elif key=='close'  : array = _nd.grey_closing(   array, footprint=fp)
        
        elif key=='linear_map': array = linear_map(array, min=value[0], max=value[1])
        elif key=='normalize' : array = normalize( array, method = value)
        elif key=='adaptive'  : array = adaptive(  array, sigma  = value, kernel='gaussian')
        elif key=='adaptive_' : array = adaptive(  array, sigma  = value, kernel='uniform')
        else: 
            print '\033[031mUnrecognized filter :', key
            
    return array

def multi_label_edge_detection(data):
    """Detect the edge in the image with multi-labels."""
    f = nd.generate_binary_structure(len(data.shape), 1)
    # the unwanted thick bounds
    bound = (nd.grey_erosion(data, footprint=f) != nd.grey_dilation(data, footprint=f)) - (
        nd.binary_dilation(data.astype(np.bool)) - data.astype(np.bool)
    )
    data = bound.astype(data.dtype)
    return data

    def outline_segments(self, mask_background=False):
        """
        Outline the labeled segments.

        The "outlines" represent the pixels *just inside* the segments,
        leaving the background pixels unmodified.

        Parameters
        ----------
        mask_background : bool, optional
            Set to `True` to mask the background pixels (labels = 0) in
            the returned image.  This is useful for overplotting the
            segment outlines on an image.  The default is `False`.

        Returns
        -------
        boundaries : 2D `~numpy.ndarray` or `~numpy.ma.MaskedArray`
            An image with the same shape of the segmentation image
            containing only the outlines of the labeled segments.  The
            pixel values in the outlines correspond to the labels in the
            segmentation image.  If ``mask_background`` is `True`, then
            a `~numpy.ma.MaskedArray` is returned.

        Examples
        --------
        >>> from photutils import SegmentationImage
        >>> segm = SegmentationImage([[0, 0, 0, 0, 0, 0],
        ...                           [0, 2, 2, 2, 2, 0],
        ...                           [0, 2, 2, 2, 2, 0],
        ...                           [0, 2, 2, 2, 2, 0],
        ...                           [0, 2, 2, 2, 2, 0],
        ...                           [0, 0, 0, 0, 0, 0]])
        >>> segm.outline_segments()
        array([[0, 0, 0, 0, 0, 0],
               [0, 2, 2, 2, 2, 0],
               [0, 2, 0, 0, 2, 0],
               [0, 2, 0, 0, 2, 0],
               [0, 2, 2, 2, 2, 0],
               [0, 0, 0, 0, 0, 0]])
        """

        from scipy.ndimage import grey_erosion, grey_dilation

        # mode='constant' ensures outline is included on the image borders
        selem = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
        eroded = grey_erosion(self.data, footprint=selem, mode='constant',
                              cval=0.)
        dilated = grey_dilation(self.data, footprint=selem, mode='constant',
                                cval=0.)

        outlines = ((dilated != eroded) & (self.data != 0)).astype(int)
        outlines *= self.data

        if mask_background:
            outlines = np.ma.masked_where(outlines == 0, outlines)

        return outlines

 def gradient(self):
     Type=self.getEdgeType()
     Size=self.image_temp.shape
     EdgeImage=np.zeros(Size)
     if (Type=='S'):
         EdgeImage=(ndimage.grey_dilation(self.image_temp,footprint=np.ones([3,3]))-ndimage.grey_erosion(self.image_temp,footprint=np.ones([3,3])))/2
     if (Type=='I'):
         EdgeImage=(self.image_temp-ndimage.grey_erosion(self.image_temp,footprint=np.ones([3,3])))/2
     if (Type=='E'):
         EdgeImage=(ndimage.grey_dilation(self.image_temp,footprint=np.ones([3,3]))-self.image_temp)/2
     self.showView4(EdgeImage)

def calcSupNp(preprob, posp, lungs, imscan, pat, midx, psp, dictSubP, dimtabx):
    '''calculate the number of reticulation and HC in subpleural'''
    print 'number of subpleural for : ', pat, psp
    imgngray = np.copy(lungs)
    np.putmask(imgngray, imgngray == 1, 0)
    np.putmask(imgngray, imgngray > 0, 1)

# subErosion=  in mm
#avgPixelSpacing=0.734 in mm/ pixel
    subErosionPixel = int(round(2 * subErosion / avgPixelSpacing))
    erosion = ndimage.grey_erosion(
        imgngray, size=(
            subErosionPixel, subErosionPixel))

    np.putmask(erosion, erosion > 0, 1)
    mask_inv = np.bitwise_not(erosion)
    subpleurmask = np.bitwise_and(imgngray, mask_inv)

    ill = 0
    for ll in posp:

        xpat = ll[0]
        ypat = ll[1]

        proba = preprob[ill]
        prec, mprobai = maxproba(proba)
        classlabel = fidclass(prec, classif)

        if xpat >= midx:
            pospr = 1
            pospl = 0
        else:
            pospr = 0
            pospl = 1
        if classlabel == pat and mprobai > thrprobaUIP:
            tabpatch = np.zeros((dimtabx, dimtabx), np.uint8)
            tabpatch[ypat:ypat + dimpavy, xpat:xpat + dimpavx] = 1
            tabsubpl = np.bitwise_and(subpleurmask, tabpatch)

            area = tabsubpl.sum()
#                    check if area above threshold
            targ = float(area) / pxy

            if targ > thrpatch:
                dictSubP[pat]['all'] = (
                    dictSubP[pat]['all'][0] + pospl,
                    dictSubP[pat]['all'][1] + pospr)
                dictSubP[pat][psp] = (
                    dictSubP[pat][psp][0] + pospl,
                    dictSubP[pat][psp][1] + pospr)

        ill += 1
    return dictSubP

def PreProcess(im):
   im=NMode(im)
   im1 = ndimage.grey_erosion(im, size=(10,10))
   scipy.misc.imsave("eroded.jpg",im1)
   im1= Image.open("eroded.jpg")
 
   im=ImageOps.equalize(im,0)
   im=ImageChops.difference(im1, im)
   #print ("image height %d and width %d\n"%(imh,imw))
 
 
   im=GBinarization(im)#binarize the image
   return im

 def fit_background(self, image, subscriber=0):
     poldegree = int(self.paramPoldegree.value)
     swidth = int(self.paramSwidth.value)
     sheight = int(self.paramSheight.value)
     threshold = int(self.paramThreshold.value)
     mediansize = int(self.paramMediansize.value)
     medianruns = int(self.paramMedianruns.value)
     darksize = int(self.paramDarksize.value)
     darkruns = int(self.paramDarkruns.value)
     brightsize = int(self.paramBrightsize.value)
     brightruns = int(self.paramBrightruns.value)
     dopreview = self.paramDopreview.value
     data = image.data
     # Median:
     for run in xrange(medianruns):
         data = ndimage.median_filter(data, size=mediansize)
     # Suspend dark spots:
     for run in xrange(darkruns):
         data = 255 - ndimage.grey_erosion(255 - data, size=darksize)
     # Suspend features:
     for run in xrange(brightruns):
         data = ndimage.grey_erosion(data, size=brightsize)
     # Fit background:
     if not dopreview:
         data = self.fit(data, poldegree, swidth, sheight, threshold)
     longname = "FitBackground"
     result = DataContainer.FieldContainer(
         data,
         copy.deepcopy(image.unit),
         copy.deepcopy(image.error),
         copy.deepcopy(image.mask),
         copy.deepcopy(image.dimensions),
         longname,
         image.shortname,
         copy.deepcopy(image.attributes),
         False,
     )
     result.seal()
     return result

def erode(im):
    """`scipy.ndimage.grey_erosion` with size set to 3 on each axis.

    Parameters
    ----------
    im : np.ndarray, arbitrary type and shape.
        The input image.

    Returns
    -------
    out : np.ndarray, same type and shape as `im`
        The eroded image.
    """
    return nd.grey_erosion(im, size=[3] * im.ndim)

def calib(extra=False):
    global bins,refe,mall,dall,nomal
    #calibration plate
    sele=(array(mall)<20)*(array(dall)>700)
    from scipy import ndimage
    sele2=ndimage.grey_erosion(sele,2)
    ia,ib=where(sele2>0.5)
    nomal=mean([zzall[ia[i]][ib[i]] for i in range(100)],0)
    if extra:
        #testing individual lines:
        ps=array([sum(ia<i) for i in range(30,41)])
        koral=[mean([zzall[i+30][b] for b in ib[ps[i]:ps[i+1]]],0) for i in range(10)]
        return nomal,koral
    return nomal

def fast_rag(labels, connectivity=1):
    """Build a data-free region adjacency graph quickly.

    Parameters
    ----------
    labels : array of int
        Image pre-segmentation or segmentation
    connectivity : int in {1, ..., labels.ndim}, optional
        Use square connectivity equal to `connectivity`. See
        `scipy.ndimage.generate_binary_structure` for more.

    Returns
    -------
    g : networkx Graph
        A graph where nodes represent regions in `labels` and edges
        indicate adjacency.

    Examples
    --------
    >>> labels = np.array([1, 1, 5, 5], dtype=np.int_)
    >>> fast_rag(labels).edges()
    [(1, 5)]
    >>> labels = np.array([[1, 1, 1, 2, 2],
    ...                    [1, 1, 1, 2, 2],
    ...                    [3, 3, 4, 4, 4],
    ...                    [3, 3, 4, 4, 4]], dtype=np.int_)
    >>> sorted(fast_rag(labels).edges())
    [(1, 2), (1, 3), (1, 4), (2, 4), (3, 4)]
    """
    conn = ndi.generate_binary_structure(labels.ndim, connectivity)
    eroded = ndi.grey_erosion(labels, footprint=conn)
    dilated = ndi.grey_dilation(labels, footprint=conn)
    boundaries0 = (eroded != labels)
    boundaries1 = (dilated != labels)
    labels_small = np.concatenate((eroded[boundaries0], labels[boundaries1]))
    labels_large = np.concatenate((labels[boundaries0], dilated[boundaries1]))
    n = np.max(labels_large) + 1
    # use a dummy broadcast array as data for RAG
    data = np.broadcast_to(np.ones((1,), dtype=np.int_),
                           labels_small.shape)
    sparse_graph = sparse.coo_matrix((data, (labels_small, labels_large)),
                                     dtype=np.int_, shape=(n, n)).tocsr()
    rag = nx.from_scipy_sparse_matrix(sparse_graph, edge_attribute='count')
    return rag

def morphop(im, operation='open', radius='5'):
    """Perform a morphological operation with spherical structuring element.

    Parameters
    ----------
    im : array, shape (M, N[, P])
        2D or 3D grayscale image.
    operation : string, optional
        The operation to perform. Choices are 'opening', 'closing',
        'erosion', and 'dilation'. Imperative verbs also work, e.g.
        'dilate'.
    radius : int, optional
        The radius of the structuring element (disk or ball) used.

    Returns
    -------
    imout : array, shape (M, N[, P])
        The transformed image.

    Raises
    ------
    ValueError : if the image is not 2D or 3D.
    """
    if im.ndim == 2:
        selem = skmorph.disk(radius)
    elif im.ndim == 3:
        selem = skmorph.ball(radius)
    else:
        raise ValueError("Image input to 'morphop' should be 2D or 3D"
                         ", got %iD" % im.ndim)
    if operation.startswith('open'):
        imout = nd.grey_opening(im, footprint=selem)
    elif operation.startswith('clos'):
        imout = nd.grey_closing(im, footprint=selem)
    elif operation.startswith('dila'):
        imout = nd.grey_dilation(im, footprint=selem)
    elif operation.startswith('ero'):
        imout = nd.grey_erosion(im, footprint=selem)
    return imout

def remove_merged_boundaries(labels, connectivity=1):
    """Remove boundaries in a label field when they separate the same region.

    By convention, the boundary label is 0, and labels are positive.

    Parameters
    ----------
    labels : array of int
        The label field to be processed.
    connectivity : int in {1, ..., labels.ndim}, optional
        The morphological connectivity for considering neighboring voxels.

    Returns
    -------
    labels_out : array of int
        The same label field, with unnecessary boundaries removed.

    Examples
    --------
    >>> labels = np.array([[1, 0, 1], [0, 1, 0], [2, 0, 3]], np.int)
    >>> remove_merged_boundaries(labels)
    array([[1, 1, 1],
           [0, 1, 0],
           [2, 0, 3]])
    """
    boundary = 0
    labels_out = labels.copy()
    is_boundary = (labels == boundary)
    labels_complement = labels.copy()
    labels_complement[is_boundary] = labels.max() + 1
    se = nd.generate_binary_structure(labels.ndim, connectivity)
    smaller_labels = nd.grey_erosion(labels_complement, footprint=se)
    bigger_labels = nd.grey_dilation(labels, footprint=se)
    merged = is_boundary & (smaller_labels == bigger_labels)
    labels_out[merged] = smaller_labels[merged]
    return labels_out