def imfill(image):
    """
    Replicates the imfill function available within MATLAB.
    Based on the example provided in http://scikit-image.org/docs/dev/auto_examples/plot_holes_and_peaks.html#example-plot-holes-and-peaks-py.

    :param image):
        A 2D numpy array containing "holes". Darker pixels surrounded by brighter pixels.

    :return:
        A 2D numpy array of type Float with the "holes" from image filled.

    """

    seed = image.copy()

    # Define seed points and the start points for the erosion process.
    seed[1:-1, 1:-1] = image.max()

    # Define the mask; Probably unneeded.
    mask = image

    # Fill the holes
    filled = morphology.reconstruction(seed, mask, method='erosion')

    return filled

def bin_to_color(RGB_image, bin_image, feat_vec, marge=None, do_label=True):
    bin_image_copy = bin_image.copy()
    res = np.zeros(shape=(bin_image.shape[0], bin_image.shape[1], 3))
    if marge is not None:
        seed = np.zeros_like(bin_image_copy)
        seed[marge:-marge, marge:-marge] = 1
        mask = bin_image_copy.copy()
        mask[ mask > 0 ] = 1
        mask[marge:-marge, marge:-marge] = 1
        reconstructed = reconstruction(seed, mask, 'dilation')
        bin_image_copy[reconstructed == 0] = 0
    if do_label:
        bin_image_copy = label(bin_image_copy)

    if len(np.unique(bin_image_copy)) != 2:
        if len(np.unique(bin_image_copy)) == 1:
            if 0 in bin_image_copy:
                print "Return blank matrix."
                return bin_image_copy
            else:
                print "Error, must give a bin image."
    RegProp = regionprops(bin_image_copy)
    for i in range(len(RegProp)):
        mini_reg = RegProp[i]
        bin_image_copy[mini_reg.coords] = feat_vec[i]
    return bin_image_copy

def detectOpticDisc(image):
    kernel = octagon(10, 10)
    thresh = threshold_otsu(image[:,:,1])
    binary = image > thresh
    print binary.dtype
    luminance = convertToHLS(image)[:,:,2]
    t = threshold_otsu(luminance)
    t = erosion(luminance, kernel)
    
    
    labels = segmentation.slic(image[:,:,1], n_segments = 3)
    out = color.label2rgb(labels, image[:,:,1], kind='avg')
    skio.imshow(out)
    
    x, y = computeCentroid(t)
    print x, y
    rows, cols, _ = image.shape
    p1 = closing(image[:,:,1],kernel)
    p2 = opening(p1, kernel)
    p3 = reconstruction(p2, p1, 'dilation')
    p3 = p3.astype(np.uint8)
    #g = dilation(p3, kernel)-erosion(p3, kernel)
    #g = rank.gradient(p3, disk(5))
    g = cv2.morphologyEx(p3, cv2.MORPH_GRADIENT, kernel)
    #markers = rank.gradient(p3, disk(5)) < 10
    markers = drawCircle(rows, cols, x, y, 85)
    #markers = ndimage.label(markers)[0]
    #skio.imshow(markers)
    g = g.astype(np.uint8)
    #g = cv2.cvtColor(g, cv2.COLOR_GRAY2RGB)
    w = watershed(g, markers)
    print np.max(w), np.min(w)
    w = w.astype(np.uint8)
    #skio.imshow(w)
    return w

def dilation(data):
    """
    Use dilation to define the background. Not working too well...
    """
    image = gaussian_filter(data, 1)

    h = 1
    seed = np.copy(image) - h
    mask = image
    dilated = reconstruction(seed, mask, method='dilation')
    dilated = data - dilated

    fileIO.writeFITS(dilated, 'dilation.fits', int=False)

    #plot the image
    fig = plt.figure(figsize=(18, 10))
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)

    ax1.set_title('Data')
    ax2.set_title('Background')

    im1 = ax1.imshow(np.log10(data), origin='lower', vmin=2., vmax=3.5, interpolation='none')
    im2 = ax2.imshow(dilated, origin='lower', interpolation='none')
    c1 = plt.colorbar(im1, ax=ax1, orientation='horizontal')
    c2 = plt.colorbar(im2, ax=ax2, orientation='horizontal')
    c1.set_label('$\log_{10}$(Counts [ADU])')
    c2.set_label('Dilation')
    plt.savefig('dilation.png')
    plt.close()

    return dilated.ravel()

 def FilterMaxima(self, Minimum=None):
     if Minimum == None:
         Minimum = self.ManipulatedData.min()
     Seed = self.ManipulatedData - Minimum
     Mask = self.ManipulatedData
     self.Dilated = reconstruction(Seed, Mask, method="dilation")
     self.ManipulatedData = self.ManipulatedData - self.Dilated

 def h_minima(self, img, h):
     img_shift = img.copy()
     img_shift[img_shift >= 255 - h] = 255-h
     img_shift = img_shift + h
     rec = morphology.reconstruction(img_shift, img, method='erosion').astype(np.dtype('uint8'))
     diff = rec - img
     return diff

 def morpho_rec2(self, img, size=10):
     # internal gradient of the cells: 
     se = morphology.diamond(size)
     dil = morphology.dilation(img, se)
     rec = morphology.reconstruction(dil, img, method='erosion').astype(np.dtype('uint8'))
             
     return rec

def get_stomata(max_proj_image, min_obj_size=200, max_obj_size=1000):
    """Performs image segmentation from a max_proj_image.
     Disposes of objects in range min_obj_size to
    max_obj_size

    :param max_proj_image: the maximum projection image
    :type max_proj_image: numpy.ndarray, uint16
    :param min_obj_size: minimum size of object to keep
    :type min_obj_size: int
    :param max_obj_size: maximum size of object to keep
    :type max_obj_size: int
    :returns: list of [ [coordinates of kept objects - list of slice objects],
                        binary object image - numpy.ndarray,
                        labelled object image - numpy.ndarray
                     ]

    """

    # pore_margin = 10
    # max_obj_size = 1000
    # min_obj_size = 200
    # for prop, value in segment_options:
    #     if prop == 'pore_margin':
    #         pore_margin = value
    #     if prop == 'max_obj_size':
    #         max_obj_size = value
    #     if prop == 'min_obj_size':
    #         min_obj_size = value
    #
    # print(pore_margin)
    # print(max_obj_size)
    # print(min_obj_size)

    #rescale_min = 50
    #rescale_max= 100
    #rescaled = exposure.rescale_intensity(max_proj_image, in_range=(rescale_min,rescale_max))
    rescaled = max_proj_image
    seed = np.copy(rescaled)
    seed[1:-1, 1:-1] = rescaled.max()
    #mask = rescaled
    #if gamma != None:
    #    rescaled = exposure.adjust_gamma(max_proj_image, gamma)
    #filled = reconstruction(seed, mask, method='erosion')
    closed = dilation(rescaled)
    seed = np.copy(closed)
    seed[1:-1, 1:-1] = closed.max()
    mask = closed


    filled = reconstruction(seed, mask, method='erosion')
    label_objects, nb_labels = ndimage.label(filled)
    sizes = np.bincount(label_objects.ravel())
    mask_sizes = sizes
    mask_sizes = (sizes > min_obj_size) & (sizes < max_obj_size)
    #mask_sizes = (sizes > 200) & (sizes < 1000)
    mask_sizes[0] = 0
    big_objs = mask_sizes[label_objects]
    stomata, _ = ndimage.label(big_objs)
    obj_slices = ndimage.find_objects(stomata)
    return [obj_slices, big_objs, stomata]

def bin_analyser(RGB_image, bin_image, list_feature, marge=None, pandas_table=False, do_label=True):
    bin_image_copy = bin_image.copy()

    p = 0
    for feat in list_feature:
        p += feat.size

    if marge is not None and marge != 0:
        seed = np.zeros_like(bin_image_copy)
        seed[marge:-marge, marge:-marge] = 1
        mask = bin_image_copy.copy()
        mask[ mask > 0 ] = 1
        mask[marge:-marge, marge:-marge] = 1
        reconstructed = reconstruction(seed, mask, 'dilation')
        bin_image_copy[reconstructed == 0] = 0
    if do_label:
        bin_image_copy = label(bin_image_copy)

    if len(np.unique(bin_image_copy)) != 2:
        if len(np.unique(bin_image_copy)) == 1:
            if 0 in bin_image_copy:
                print "Return blank matrix. Change this shit"
                white_npy = np.zeros(shape=(1, p))
                if not pandas_table:
                    return white_npy
                else:
                    names = GetNames(list_feature) 
                    return pd.DataFrame(white_npy, columns=names)
            else:
                print "Error, must give a bin image."
    
    GrowRegion_N = NeededGrownRegion(list_feature)
    img = {0: bin_image_copy}
    RegionProp = {0: regionprops(bin_image_copy)}
    for val in GrowRegion_N:
        if val != 0:
            img[val] = GrowRegion(bin_image_copy, val)
            RegionProp[val] = regionprops(img[val])
    
    n = len(RegionProp[0])

    TABLE = np.zeros(shape=(n,p))
    for i in range(n):
        offset_ALL = 0
        for j, feat in enumerate(list_feature):
            tmp_regionprop = RegionProp[feat._return_n_extension()][i]
            off_tmp = feat.size      
            TABLE[i, (j + offset_ALL):(j + offset_ALL + off_tmp)] = feat._apply_region(tmp_regionprop ,RGB_image)

            offset_ALL += feat.size - 1

    if pandas_table:
        names = GetNames(list_feature)
        return pd.DataFrame(TABLE, columns=names)
    else:
        return TABLE

def shadow_morphology(image):
    image_shape = image.shape
    _image = np.zeros((image_shape[0], image_shape[1]))
    _image = image
    seed = np.copy(_image)
    seed[1:-1, 1:-1] = image.max()
    mask = _image

    filled = reconstruction(seed, mask, method='erosion')
    result = filled - _image
    return result

 def homogenize(self, img):
     img_sub = img.astype(np.dtype('float')) - self.settings.homogenize_settings['h']
     img_sub[img_sub<0] = 0
     img_sub = img_sub.astype(img.dtype)
     
     # seed and then mask
     print 'img: ', np.min(img), np.max(img)
     print 'img_sub: ', np.min(img_sub), np.max(img_sub)
     temp = morphology.reconstruction(img_sub, img)
     res = img - temp.astype(img.dtype)
     return res

def temblob(image,ind):

	"""
	Laplacian of gaussian blob detection for TEM images
	"""

	org = image[4:-256,4:-4]


	with warnings.catch_warnings():
		warnings.simplefilter("ignore")
		warnings.warn("user", UserWarning)

		igray = img_as_ubyte(rgb2gray(org))
		iinv = np.invert(igray)


	igaus = img_as_float(iinv)
	igaus = gaussian_filter(igaus, 1)
	h = 0.5
	sd = igaus - h
	msk = igaus
	dilat = reconstruction(sd, msk, method='dilation')
	hdome = igaus - dilat

	if ind == 'AgNP':
		kwargs = {}
		kwargs['threshold'] = 0.01
		kwargs['overlap'] = 0.4
		kwargs['min_sigma'] = 25
		kwargs['max_sigma'] = 50
		kwargs['num_sigma'] = 25
		calib = 500/(969-26)
		
	elif ind == 'AuNP':
		kwargs = {}
		kwargs['threshold'] = 0.01
		kwargs['overlap'] = 0.4
		kwargs['min_sigma'] = 18
		kwargs['max_sigma'] = 23
		kwargs['num_sigma'] = 5
		calib = 200/(777-23)
		
	else:
		warnmsg='Unable to identify keyword: {:}'.format(ind)
		warnings.warn(warnmsg, UserWarning)

	blobs = blob_log(hdome, **kwargs)
	diam = 2*sqrt(2)*blobs[:,-1]
	npdiam = [ind]
	npdiam.extend(calib*diam)

	return(npdiam)

    def local_maxima(self, img, h=1):
        img_sub = img.astype(np.dtype('float')) - h
        img_sub[img_sub<0] = 0
        img_sub = img_sub.astype(img.dtype)
        
        # seed and then mask
        temp = morphology.reconstruction(img_sub, img)
        res = img - temp.astype(img.dtype)
        res[res>0] = 255        
        res = res.astype(np.dtype('uint8'))

        return res    

def bg_sub_signPreserveNorm(imageFile, path, extent, extension):
    vec = signPreserveNorm(imageFile, path, extent, extension)
    image = np.reshape(vec, (20,20), order="F")

    image = gaussian_filter(image, 1)
    seed = np.copy(image)
    seed[1:-1, 1:-1] = image.min()
    mask = image

    dilated = reconstruction(seed, mask, method='dilation')
    
    return np.ravel(image - dilated, order="F")

def opening_by_reconstruction(img, se):
    diffImg = True
    orig_img = img.copy()
    last_rec = img.copy()
    while diffImg:
        er_img = morphology.erosion(img, se)
        img = morphology.reconstruction(er_img, orig_img)
        if np.array_equal(last_rec, img):
            diffImg = False
        else:
            last_rec = img.copy()
    return last_rec

def change_color(table, bin, color_vec, indice, res):
    x_cent, y_cent = table[indice, 3:5]
    X, Y = int(x_cent), int(y_cent)
    only_cell = np.zeros_like(bin)
    if bin[X, Y] != 1:
        X, Y = find_closest_cc(X, Y, bin)


    only_cell[X, Y] = 1
    only_cell = reconstruction(only_cell, bin)
    x, y = np.where(only_cell == 1)
    res[x, y] = color_vec

def HreconstructionErosion(prob_img, h):

    def making_top_mask(x, lamb=h):
	   return min(255, x + lamb)

    f = np.vectorize(making_top_mask)
    shift_prob_img = f(prob_img)

    seed = shift_prob_img
    mask = prob_img
    recons = reconstruction(
        seed, mask, method='erosion').astype(np.dtype('ubyte'))
    return recons

def ifcb_segment(img):
    Y = img_as_float(img)
    # step 1. local variance
    Yv = rescale_intensity(generic_filter(Y, np.var, footprint=disk(3)))
    # step 2. threshold local variance, aggressively
    Ye = Yv > (threshold_otsu(Yv) / 2.)
    # step 3. dark areas
    Yt = Y < threshold_otsu(Y)
    thin_blob = Ye | Yt
    # step 4. morphological reconstruction
    seed = np.copy(thin_blob)
    seed[1:-1,1:-1] = 1
    four=np.disk(1).astype(np.bool)
    return reconstruction(seed,thin_blob,method='erosion',selem=four)

def analyze_image():  
    
    global sin_par, dipole, multiple, contours
    
    # Change to the gray image
    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        
    # Adaptive Thresholding
    block_size = 201
    img_binary = threshold_adaptive(img_gray, block_size, offset=7)    
    img_binary = img_binary.astype(dtype='uint8')*255
    
    # Remove the salt and pepper noise
    img_binary = cv2.medianBlur(img_binary,7)    
    
    # Fill the hole
    img_binary = 255-img_binary # Change the white/black for filling the hole     
    seed = np.copy(img_binary)
    seed[1:-1, 1:-1] = img_binary.max()
    mask = img_binary
    filled = reconstruction(seed, mask, method='erosion')
    filled = filled.astype(dtype = 'uint8')

    # Find the contours
    # For terminal use    
#    image, contours, hierarchy = cv2.findContours(filled,
#                                 cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)         
    # For Anaconda use    
    hierarchy, contours, hierarchy = cv2.findContours(filled,
                                  cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)        

    
    # Calculate the area (= len) of each contour
    area = np.zeros(len(contours))
    area_contour = np.zeros(len(contours))

    for kk in range(0,len(contours)):
        area[kk] = len(contours[kk])
        area_contour[kk] = cv2.contourArea(contours[kk])
  
    single = np.where((area_contour>200)&(area_contour<=624))[0]
    dipole = np.where((area_contour>624)&(area_contour<=936))[0]
    multiple = np.where(area_contour>936)[0]    
    
    # Analyze the single particle     
    sin_par = analyze_single(single,filled)
    
    # Draw the contours
    refreshFigure_contours(img,sin_par,dipole,multiple,contours)    

def remove_bg(args):
    """
    Remove background from reconstructions.
    """
    data, args, ind_start, ind_end = args
    
    for m in range(ind_end-ind_start):
        img = data[m, :, :]
        
        # first remove background.
        seed = np.copy(img)
        seed[1:-1, 1:-1] = img.min()
        img -= reconstruction(seed, img, method='dilation')
        
        data[m, :, :] = img
    return ind_start, ind_end, data