def binaryClosing(binarydata, structure=None, iterations=1):
    
    result = np.zeros_like(binarydata)
    if structure is None:
        ndimage.binary_closing(binarydata, iterations=iterations, output=result)
    else:
        ndimage.binary_closing(binarydata, structure=structure, iterations=iterations, output=result)
    return result

def closing():
    image_list = get_one_imagefrom_mnist()
    image_array =np.asarray(image_list)
    image =image_array.reshape(28, 28)
    
    ndimage.binary_closing(image, structure=np.ones((2,2))).astype(int)
    plt.imshow(image, cmap=cm.binary)
    plt.show()

def refine_worm(image, initial_area, candidate_edges):
    # find strong worm edges (roughly equivalent to the edges found by find_initial_worm,
    # which are in candidate_edges): smooth the image, do canny edge-finding, and
    # then keep only those edges near candidate_edges
    smooth_image = restoration.denoise_tv_bregman(image, 140).astype(numpy.float32)
    smoothed, gradient, sobel = canny.prepare_canny(smooth_image, 8, initial_area)
    local_maxima = canny.canny_local_maxima(gradient, sobel)
    candidate_edge_region = ndimage.binary_dilation(candidate_edges, iterations=4)
    strong_edges = local_maxima & candidate_edge_region

    # Now threshold the image to find dark blobs as our initial worm region
    # First, find areas in the initial region unlikely to be worm pixels
    mean, std = mcd.robust_mean_std(smooth_image[initial_area][::4], 0.85)
    non_worm = (smooth_image > mean - std) & initial_area
    # now fit a smoothly varying polynomial to the non-worm pixels in the initial
    # region of interest, and subtract that from the actual image to generate
    # an image with a flat illumination field
    background = polyfit.fit_polynomial(smooth_image, mask=non_worm, degree=2)
    minus_bg = smooth_image - background
    # now recalculate a threshold from the background-subtracted pixels
    mean, std = mcd.robust_mean_std(minus_bg[initial_area][::4], 0.85)
    initial_worm = (minus_bg < mean - std) & initial_area
    # Add any pixels near the strong edges to our candidate worm position
    initial_worm |= ndimage.binary_dilation(strong_edges, iterations=3)
    initial_worm = mask.fill_small_radius_holes(initial_worm, 5)

    # Now grow/shrink the initial_worm region so that as many of the strong
    # edges from the canny filter are in contact with the region edges as possible.
    ac = active_contour.EdgeClaimingAdvection(initial_worm, strong_edges,
        max_region_mask=initial_area)
    stopper = active_contour.StoppingCondition(ac, max_iterations=100)
    while stopper.should_continue():
        ac.advect(iters=1)
        ac.smooth(iters=1, depth=2)
    worm_mask = mask.fill_small_radius_holes(ac.mask, 7)

    # Now, get edges from the image at a finer scale
    smoothed, gradient, sobel = canny.prepare_canny(smooth_image, 0.3, initial_area)
    local_maxima = canny.canny_local_maxima(gradient, sobel)
    strong_sum = strong_edges.sum()
    highp = 100 * (1 - 1.5*strong_sum/local_maxima.sum())
    lowp = max(100 * (1 - 3*strong_sum/local_maxima.sum()), 0)
    low_worm, high_worm = numpy.percentile(gradient[local_maxima], [lowp, highp])
    fine_edges = canny.canny_hysteresis(local_maxima, gradient, low_worm, high_worm)

    # Expand out the identified worm area to include any of these finer edges
    closed_edges = ndimage.binary_closing(fine_edges, structure=S)
    worm = ndimage.binary_propagation(worm_mask, mask=worm_mask|closed_edges, structure=S)
    worm = ndimage.binary_closing(worm, structure=S, iterations=2)
    worm = mask.fill_small_radius_holes(worm, 5)
    worm = ndimage.binary_opening(worm)
    worm = mask.get_largest_object(worm)
    # Last, smooth the shape a bit to reduce sharp corners, but not too much to
    # sand off the tail
    ac = active_contour.CurvatureMorphology(worm, max_region_mask=initial_area)
    ac.smooth(depth=2, iters=2)
    return strong_edges, ac.mask

def get_uv_mask(vertices_vis, triangles, uv_coords, h, w, resolution):
    triangles = triangles.T
    vertices_vis = vertices_vis.astype(np.float32)
    uv_mask = render_texture(uv_coords.T, vertices_vis[np.newaxis, :], triangles, resolution, resolution, 1)
    uv_mask = np.squeeze(uv_mask > 0)
    uv_mask = ndimage.binary_closing(uv_mask)
    uv_mask = ndimage.binary_erosion(uv_mask, structure = np.ones((4,4)))  
    uv_mask = ndimage.binary_closing(uv_mask)
    uv_mask = ndimage.binary_erosion(uv_mask, structure = np.ones((4,4)))  
    uv_mask = ndimage.binary_erosion(uv_mask, structure = np.ones((4,4)))  
    uv_mask = ndimage.binary_erosion(uv_mask, structure = np.ones((4,4)))  
    uv_mask = uv_mask.astype(np.float32)

    return np.squeeze(uv_mask)

def plot_mask(mask, plot_axis=None, color='#ff0000', closing_iteration=None, **kwargs):
    '''
    plot mask (ROI) borders by using pyplot.contour function. all the 0s and Nans in the input mask will be considered
    as background, and non-zero, non-nan pixel will be considered in ROI.
    '''
    if not check_binary_2d_array(mask):
        raise(ValueError, 'input mask should be a 2d binary numpy.ndarray with dtype as integer and contains '
                          'only 0s and 1s.')

    if not plot_axis:
        f = plt.figure()
        plot_axis = f.add_subplot(111)

    if closing_iteration is not None:
        ploting_mask = ni.binary_closing(mask, iterations=closing_iteration).astype(np.uint8)
    else:
        ploting_mask = mask

    currfig = plot_axis.contourf(ploting_mask, levels=[0.5, 1], colors=color, **kwargs)

    # put y axis in decreasing order
    y_lim = list(plot_axis.get_ylim())
    y_lim.sort()
    plot_axis.set_ylim(y_lim[::-1])

    plot_axis.set_aspect('equal')

    return currfig

def main():
	
	usage="Look for holes in a 3d density map.\n\tfindholesinmap.py map.mrc\n**********scipy is required!*********"
	parser = EMArgumentParser(usage=usage,version=EMANVERSION)
	parser.add_argument("--thr", type=float,help="Threshold for the isosurface", default=1)
	parser.add_argument("--closeiter", type=int,help="Number of iterations for the closing operation", default=10)
	parser.add_argument("--filter_res", type=float,help="Resolution for the final filter", default=10)
	parser.add_argument("--output", type=str,help="output file name", default=None)
	(options, args) = parser.parse_args()
	logid=E2init(sys.argv)
	
	e=EMData(args[0])
	img=e.numpy()
	if options.output==None:
		options.output=args[0][:-4]+"_holes.hdf"
	
	apix=e["apix_x"]
	img_open=ndimage.binary_closing(img>options.thr,iterations=options.closeiter)
	m=img.copy()
	m[m<0]=0
	m/=np.max(m)
	hole=img_open-m
	a=from_numpy(hole)

	
	a["apix_x"]=apix
	a["apix_y"]=apix
	a["apix_z"]=apix
	a.process_inplace("filter.lowpass.gauss",{"cutoff_freq":1./options.filter_res})

	a.write_image(options.output)
	
	E2end(logid)

  def findCom(self,data): 
    data = data.astype(int)
    if self.update_counter >= 5:
      self.update_counter = 0

      ##########################################################################
      ## Update the background image, adding a new image and removing the oldest.
      ##########################################################################
      self.background_list.insert(0,data)
      self.background_list.pop()
      background = np.zeros((480, 640, 3), dtype=int)
      for b in self.background_list:
        background += b
      self.background = background/len(self.background_list)
    
    ############################################################################
    ## Detect foreground by looking at difference from mean.
    ############################################################################
    foreground = np.sum(np.abs(np.subtract(self.background,data)),axis=2)
    falseImage = foreground
    ## clean foreground image
    falseImage[falseImage > 100] = 255
    falseImage[falseImage < 101] = 0
    falseImage = ndimage.binary_opening(falseImage)
    falseImage = ndimage.binary_closing(falseImage)
    com = ndimage.measurements.center_of_mass(falseImage)
    self.update_counter += 1
    return com

def xy_map_to_np_image(xy_map,m_per_pixel,dilation_count=0,padding=50):
    ''' returns binary numpy image. (255 for occupied
        pixels, 0 for unoccupied)
        2d array
    '''
    min_x = np.min(xy_map[0,:])
    max_x = np.max(xy_map[0,:])
    min_y = np.min(xy_map[1,:])
    max_y = np.max(xy_map[1,:])
    br = np.matrix([min_x,min_y]).T

    n_x = int(round((max_x-min_x)/m_per_pixel)) + padding
    n_y = int(round((max_y-min_y)/m_per_pixel)) + padding
    img = np.zeros((n_x+padding,n_y+padding),dtype='int')
    occupied_pixels = np.matrix([n_x,n_y]).T - np.round((xy_map-br)/m_per_pixel).astype('int')
    
    if dilation_count == 0:
        img[(occupied_pixels[0,:],occupied_pixels[1,:])] = 255
    else:
        img[(occupied_pixels[0,:],occupied_pixels[1,:])] = 1
        connect_structure = np.empty((3,3),dtype='int')
        connect_structure[:,:] = 1
        img = ni.binary_closing(img,connect_structure,iterations=dilation_count)
        img = ni.binary_dilation(img,connect_structure,iterations=1)
        img = img*255

    return img,n_x,n_y,br

def get_difference_spots(pix):
    bpix = pix > 20
    bpix = ndimage.binary_opening(bpix)
    bpix = ndimage.binary_closing(bpix)
    labels, n = ndimage.measurements.label(bpix)
    clicks = ndimage.measurements.center_of_mass(pix, labels, range(1, n+1))
    return clicks

def detect_vortices(cloud, radius=70, showplots=False):
    """
    Detects whether there are vortex-like features within a given radius
    of the peak density in the TOF image of an expanded BEC
    """
    OD = cloud.get_OD()    
    peak_coord = cloud.results['peak coordinates']
    center_region = ROI(center=peak_coord,
                        size=(1.5 * radius, 1.5 * radius)).slices
    smooth_cloud = ndi.median_filter(OD[center_region], size=4)
    minOD = smooth_cloud.min()
    maxOD = smooth_cloud.max()
    cloud_median = ndi.median_filter(smooth_cloud, size=10)
    belowthresh = where(smooth_cloud < cloud_median * 0.75, 1, 0)
    opened = ndi.binary_opening(belowthresh, iterations=1)
    closed = ndi.binary_closing(opened, iterations=1)
    vort_found = ndi.label(closed)[1]
    cloud.results['vort_found'] = vort_found
    if showplots == True:
        fig = plt.figure(1999)
        fig.add_subplot(221, xticks=[], yticks=[])   
        plt.imshow(smooth_cloud, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(222, xticks=[], yticks=[]) 
        plt.imshow(cloud_median, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(223, xticks=[], yticks=[]) 
        plt.imshow(closed, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
        fig.add_subplot(224, xticks=[], yticks=[]) 
        plt.imshow(belowthresh, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
    return vort_found

def detect_current(cloud, showplots=False):
    """
    Detects whether there is a vortex-like signature of persistent
    current in the center of a TOF image of an expanded ring BEC
    """
    OD = cloud.get_OD()    
    peak_coord = cloud.results['peak coordinates']
    center_region = ROI(center=peak_coord, size=(40, 40)).slices
    cloud_center = ndi.median_filter(OD[center_region], size=2)
    minOD = cloud_center.min()
    maxOD = cloud_center.max()
    cloud_median = ndi.median_filter(cloud_center, size=10)
    belowthresh = where(cloud_center < cloud_median * 0.75, 1, 0)
    opened = ndi.binary_opening(belowthresh, iterations=1)
    closed = ndi.binary_closing(opened, iterations=3)
    current_found = ndi.label(closed)[1]
    cloud.results['current_found'] = current_found
    if showplots == True:
        fig = plt.figure(1999)
        fig.add_subplot(221, xticks=[], yticks=[])   
        plt.imshow(cloud_center, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(222, xticks=[], yticks=[]) 
        plt.imshow(cloud_median, interpolation='nearest', vmin=minOD,
                                   vmax=maxOD)
        fig.add_subplot(223, xticks=[], yticks=[]) 
        plt.imshow(closed, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
        fig.add_subplot(224, xticks=[], yticks=[]) 
        plt.imshow(belowthresh, interpolation='nearest',
                   cmap=plt.cm.get_cmap('binary'))
    return current_found, asum(closed)     

def findNeuron(bgImage, threshold, xNeuron,  yNeuron):
    """find bright object in small roi image."""
    mask = np.where(bgImage > threshold[0], 1, 0)
    mask = ndimage.binary_opening(mask,structure = np.ones((2,2)))
    mask = ndimage.binary_closing(mask)
        # --- Individually label all connected regions and get their center of mass
    label_im, nb_labels = ndimage.label(mask)
    centroids = ndimage.measurements.center_of_mass(bgImage, label_im, xrange(1,nb_labels+1))
    # --- select brightest object by default (mean brightness)
    meanBrightness = ndimage.measurements.mean(bgImage, label_im, xrange(1,nb_labels+1))
#        # --- Calculate the distance of each new cms to the old neuron position
#        # --- and select the new neuron position to be the object closest to
#        # --- the old location
#    dist = []
#    for coords in centroids:
#        dist.append((coords[0]-yNeuron)**2 + (coords[1]-xNeuron)**2)
#    if len(dist)==0:
#        yNewNeuron,xNewNeuron = yNeuron,  xNeuron
#    else:
#        loc = np.argmin(dist)
#        yNewNeuron,xNewNeuron = centroids[loc]
    if nb_labels >1:
        loc = np.argmax(meanBrightness)
        yNewNeuron,xNewNeuron = centroids[loc]
    else:
        yNewNeuron,xNewNeuron = yNeuron,  xNeuron
        loc = -1  
    neuronObject = np.where(label_im == loc+1,0,1)
    neuronArea = np.sum(neuronObject)
        # --- Get average of the neuron fluoresence --- 
    tmp_neuron = np.ma.masked_array(bgImage, neuronObject)
    newNeuronAverage = np.ma.average(tmp_neuron[tmp_neuron>threshold[1]])
        
    return yNewNeuron,xNewNeuron, newNeuronAverage,neuronArea,  neuronObject

def fetch_icbm152_brain_gm_mask(data_dir=None, threshold=0.2, resume=True,
                                verbose=1):
    """Downloads ICBM152 template first, then loads 'gm' mask image.

    .. versionadded:: 0.2.5

    Parameters
    ----------
    data_dir: str, optional
        Path of the data directory. Used to force storage in a specified
        location. Defaults to None.

    threshold: float, optional
        The parameter which amounts to include the values in the mask image.
        The values lies above than this threshold will be included. Defaults
        to 0.2 (one fifth) of values.

    resume: bool, optional
        If True, try resuming partially downloaded data. Defaults to True.

    verbose: int, optional
        verbosity level (0 means no message).

    Returns
    -------
    gm_mask_img: Nifti image
        Corresponding to brain grey matter from ICBM152 template.

    Notes
    -----
    This function relies on ICBM152 templates where we particularly pick
    grey matter template and threshold the template at .2 to take one fifth
    of the values. Then, do a bit post processing such as binary closing
    operation to more compact mask image.

    Note: It is advised to check the mask image with your own data processing.

    See Also
    --------
    nilearn.datasets.fetch_icbm152_2009: for details regarding the ICBM152
        template.

    nilearn.datasets.load_mni152_template: for details about version of MNI152
        template and related.

    """
    # Fetching ICBM152 grey matter mask image
    icbm = fetch_icbm152_2009(data_dir=data_dir, resume=resume, verbose=verbose)
    gm = icbm['gm']
    gm_img = check_niimg(gm)
    gm_data = niimg._safe_get_data(gm_img)

    # getting one fifth of the values
    gm_mask = (gm_data > threshold)

    gm_mask = ndimage.binary_closing(gm_mask, iterations=2)
    gm_mask_img = new_img_like(gm_img, gm_mask)
    return gm_mask_img

def morph_sequence(pix, *param):
    for oc, wd, ht in param:
        logi(" Performing Morph : ", oc, wd, ht)
        structure = np.ones((ht, wd))
        if oc == "c":
            pix = binary_closing(pix, structure)
        elif oc == "o":
            pix = binary_opening(pix, structure)
    return pix

def masked_slic(img, mask, n_segments, compactness):
    labels = slic(img, n_segments=n_segments, compactness=compactness)
    labels += 1
    n_labels = len(np.unique(labels))
    try:
        mask = ndi.binary_closing(mask, structure=np.ones((3, 3)), iterations=1)
    except IndexError, e:
        rospy.logerr(e)
        return

def adaptive_segment(args):
    """
    Applies an adaptive threshold to reconstructed data.

    Also known as local or dynamic thresholding
    where the threshold value is the weighted mean
    for the local neighborhood of a pixel subtracted
    by constant. Alternatively the threshold can be
    determined dynamically by a given function using
    the 'generic' method.

    Parameters
    ----------
    data : ndarray, float32
        3-D reconstructed data with dimensions:
        [slices, pixels, pixels]

    block_size : scalar, int
        Uneven size of pixel neighborhood which is
        used to calculate the threshold value
        (e.g. 3, 5, 7, ..., 21, ...).

    offset : scalar, float
         Constant subtracted from weighted mean of
         neighborhood to calculate the local threshold
         value. Default offset is 0.

    Returns
    -------
    output : ndarray
        Thresholded data.

    References
    ----------
    - `http://scikit-image.org/docs/dev/auto_examples/plot_threshold_adaptive.html \
    <http://scikit-image.org/docs/dev/auto_examples/plot_threshold_adaptive.html>`_
    """
    # Arguments passed by multi-processing wrapper
    ind, dshape, inputs = args

    # Function inputs
    data = mp.tonumpyarray(mp.shared_arr, dshape)  # shared-array
    block_size, offset = inputs

    for m in ind:
        img = data[m, :, :]

        # Perform scikit adaptive thresholding.
        img = threshold_adaptive(img, block_size=block_size, offset=offset)

        # Remove small white regions
        img = ndimage.binary_opening(img)

        # Remove small black holes
        img = ndimage.binary_closing(img)

        data[m, :, :] = img

def breakup_region(component):
    distance = ndi.distance_transform_edt(component)
    skel = skeletonize(component)
    skeldist = distance*skel
    local_maxi = peak_local_max(skeldist, indices=False, footprint=disk(10))
    local_maxi=ndi.binary_closing(local_maxi,structure = disk(4),iterations = 2)
    markers = ndi.label(local_maxi)[0]
    labels = watershed(-distance, markers, mask=component)
    return(labels)

def main():
	i, h = load(sys.argv[1])
	
	i = i.copy()
	i = binary_closing(i, iterations=1)
	i = morphology2d(binary_closing, i, iterations=4)
	i = fill2d(i)
	

	save(i, sys.argv[1], h)

def masked_slic(img, mask, n_segments, compactness):
    labels = slic(img, n_segments=n_segments, compactness=compactness)
    labels += 1
    n_labels = len(np.unique(labels))
    mask = ndi.binary_closing(mask, structure=np.ones((3, 3)), iterations=1)
    labels[mask == 0] = 0  # set bg_label
    if len(np.unique(labels)) < n_labels - 2:
        sys.stderr.write('WARNING: number of label differs after masking.'
                         ' Maybe this is not good for RAG construction.\n')
    return labels

    def removeGrid(self,cs,removeGrid):
        """
        Detect the grid of the phantom and remove it from the image
        """
        shift = int(1./self.pixDim(cs)+.5)
        
        # try to find a threshold on pixelvalue to define a value representing the grid
        maskval = 0.75*cs.pixeldataIn.mean()
        # make a mask of grid-like values
        mask = (cs.pixeldataIn < maskval)

        # hole closing of the mask
        mask = scind.binary_closing(mask,structure=np.ones((5,5)))
        mask = scind.binary_opening(mask,structure=np.ones((5,5)))
        mask = scind.binary_dilation(mask)
        # new since 20150211
        mask = scind.binary_dilation(mask)

        # fill the gridlines with the median values of the pixels around it
        medimage = np.roll(cs.pixeldataIn,shift,axis=0).astype(float)
        dest = cs.pixeldataIn+mask*(medimage-cs.pixeldataIn)
        # repeat to remove propagated mask # new since 20150211
        medimage = np.roll(dest,shift,axis=0).astype(float)
        dest = cs.pixeldataIn+mask*(medimage-cs.pixeldataIn)
        medimage = None
        cs.gridimage = mask.astype(float)
        mask = None
        
        # find gridobject
        gridobject         = scind.binary_fill_holes(cs.gridimage)
        label_im,nb_labels = scind.label(gridobject)
        sizes = scind.sum(gridobject, label_im, range(nb_labels + 1))
        gridobject = None
        
        #Clean up small connect components:
        mask_size = sizes < max(sizes) #(100/self.pixDim())**2
        remove_pixel = mask_size[label_im]
        label_im[remove_pixel] = 0

        # Now reassign labels with np.searchsorted:
        labels = np.unique(label_im)
        label_im = np.searchsorted(labels, label_im)

        medimage = np.roll(dest,shift,axis=0).astype(float)
        dest += cs.gridimage*(medimage-dest)
        medimage = None
        cs.gridimage *= label_im
        if -1>0: # remove everything outside grid
            wid = dest.shape[0]
            hei = dest.shape[1]
            mv = np.mean(dest[wid/4:3*wid/4,hei/4:3*hei/4])
            dest = label_im*(dest-mv)+mv

        if removeGrid:
            cs.pixeldataIn = dest