def segment_cells(frame, mask=None):
    """
    Compute the initial segmentation based on ridge detection + watershed.
    This works reasonably well, but is not robust enough to use by itself.
    """
    
    blurred = filters.gaussian_filter(frame, 2)
    ridges = enhance_ridges(frame)
    
    # threshold ridge image
    thresh = filters.threshold_otsu(ridges)
    thresh_factor = 0.6
    prominent_ridges = ridges > thresh_factor*thresh
    prominent_ridges = morphology.remove_small_objects(prominent_ridges, min_size=256)
    prominent_ridges = morphology.binary_closing(prominent_ridges)
    prominent_ridges = morphology.binary_dilation(prominent_ridges)
    
    # skeletonize
    ridge_skeleton = morphology.medial_axis(prominent_ridges)
    ridge_skeleton = morphology.binary_dilation(ridge_skeleton)
    ridge_skeleton *= mask
    ridge_skeleton -= mask
    
    # label
    cell_label_im = measure.label(ridge_skeleton)
    
    # morphological closing to fill in the cracks
    for cell_num in range(1, cell_label_im.max()+1):
        cell_mask = cell_label_im==cell_num
        cell_mask = morphology.binary_closing(cell_mask, disk(3))
        cell_label_im[cell_mask] = cell_num
    
    return cell_label_im 

def findROI(footprint):
    
    dims = footprint.shape
    footprint = footprint.reshape(dims[0],dims[1],-1)
    ROI = zeros(footprint.shape)
    
    mxs = []
    for i in range(np.size(footprint[0,0,:])):
        mxs.append(footprint.reshape(-1,1).max())
    
    for i in range(np.size(footprint[0,0,:])):
        img = footprint[:,:,i].reshape(dims[0],dims[1])
        mx = mxs[i]
        thresh = img>0.4*mx
        thresh2 = binary_dilation(binary_dilation(thresh))
        lbls,marks = label(thresh2)
        rgs = regionprops(lbls)
        if np.size(rgs)>0:
            szs = []
            for prop in rgs:
                szs.append(prop.area)
            ind = np.argmax(szs)
            if rgs[ind].area>100:
                region = lbls==ind+1
                out = zeros([dims[0],dims[1]])
                out[region] = 1
                ROI[:,:,i] = out
                
    ROI = ROI.reshape(dims)
    
    return ROI

def shapesPlot(shapes,inds,fig,ax):
    
    from skimage.measure import label,regionprops
    from skimage import feature
    from skimage.morphology import binary_dilation
    from skimage.segmentation import find_boundaries
    import pylab as plt
    import numpy as np
    
    #fig = plt.figure()
    #ax = fig.add_subplot(111)
    sz = np.int32(shapes.shape)
    
    
    for i in inds:
        img = shapes[i,:,:]
        mx = img[:].max()
        test = img>0.4*mx
        test2 = binary_dilation(binary_dilation(test))
        lbls = label(test2)
        rgs = regionprops(lbls)
        if np.size(rgs)>0:
            szs = []
            for prop in rgs:
                szs.append(prop.area)
            ind = np.argmax(szs)
            if rgs[ind].area>100:
                pt = rgs[ind].centroid
                region = lbls==ind+1
                edges = find_boundaries(region)
                eln = edges.nonzero()
                ax.scatter(eln[1],eln[0],marker='.',color='r',linewidths=0.01)
                ax.text(pt[1]-4,pt[0]+4,'%i' % i,fontsize=14,color='k')
    
    return fig,ax

def parasites(image, cells, voronoi):
    img = Functions.global_otsu(image)
    cells = Functions.global_otsu(cells)
    s_elem = Functions.fig(Functions.fig_size)

    # Remove cells

    for i in range(Functions.iterations):
        cells = binary_dilation(cells, s_elem)
    return_image = Functions.subtraction(img, cells)

    # Remove stuff from cells

    for i in range(Functions.iterations-1):
        return_image = binary_erosion(return_image)
    return_image = binary_opening(return_image)
    for i in range(Functions.iterations - 1):
        return_image = binary_dilation(return_image)

    # Remove bigger objects

    removal_image = return_image.copy()

    for i in range(Functions.iterations + 5):
        removal_image = binary_erosion(removal_image)
    removal_image = binary_opening(removal_image)
    for i in range(Functions.iterations + 10):
        removal_image = binary_dilation(removal_image)

    return_image = Functions.subtraction(return_image, removal_image)

    # Remove voronoi lines for better quality
    return Functions.subtraction(return_image, voronoi)

def skeleton(seg):
    skel, dist = skmorph.medial_axis(seg, return_distance=True)
    node, edge, leaf = (spim.label(g, np.ones((3, 3), bool))[0] for g in skel2graph(skel))

    trim_edge = (edge != 0) & ~(skmorph.binary_dilation(node != 0, np.ones((3, 3), bool)) != 0)
    trim_edge = spim.label(trim_edge, np.ones((3, 3), bool))[0]

    leaf_edge_vals = skmorph.binary_dilation(leaf != 0, np.ones((3, 3), bool)) != 0
    leaf_edge_vals = np.unique(trim_edge[leaf_edge_vals])
    leaf_edge_vals = leaf_edge_vals[leaf_edge_vals > 0]
    leaf_edge = leaf != 0

    trim_edge = ndshm.fromndarray(trim_edge)
    leaf_edge = ndshm.fromndarray(leaf_edge)
    Parallel()(delayed(set_msk)(leaf_edge, trim_edge, l) for l in leaf_edge_vals)
    trim_edge = np.copy(trim_edge)
    leaf_edge = np.copy(leaf_edge)

    leaf_edge[(skmorph.binary_dilation(leaf_edge, np.ones((3, 3), bool)) != 0) & (edge != 0)] = True
    leaf_edge = spim.label(leaf_edge, np.ones((3, 3), bool))[0]

    leaf_edge_node = skmorph.binary_dilation(leaf_edge != 0, np.ones((3, 3), bool)) != 0
    leaf_edge_node = ((node != 0) & leaf_edge_node) | leaf_edge
    leaf_edge_node = spim.label(leaf_edge_node, np.ones((3, 3), bool))[0]

    cand_node = leaf_edge_node * (node != 0)
    cand_node = cand_node.nonzero()
    cand_node = np.transpose((leaf_edge_node[cand_node],) + cand_node + (2 * dist[cand_node],))

    cand_leaf = leaf_edge_node * (leaf != 0)
    cand_leaf = cand_leaf.nonzero()
    cand_leaf = np.transpose((leaf_edge_node[cand_leaf],) + cand_leaf)

    if len(cand_node) > 0 and len(cand_leaf) > 0:
        cand_leaf = ndshm.fromndarray(cand_leaf)
        cand_node = ndshm.fromndarray(cand_node)
        pruned = Parallel()(delayed(prune_leaves)(cand_leaf, cand_node, j) for j in np.unique(cand_node[:, 0]))
        cand_leaf = np.copy(cand_leaf)
        cand_node = np.copy(cand_node)

        pruned_ind = []
        for p in pruned:
            pruned_ind.extend(p)
        pruned_ind = tuple(np.transpose(pruned_ind))

        pruned = ~skel

        pruned = ndshm.fromndarray(pruned)
        leaf_edge = ndshm.fromndarray(leaf_edge)
        Parallel()(delayed(set_msk)(pruned, leaf_edge, l) for l in np.unique(leaf_edge[pruned_ind]))
        pruned = np.copy(pruned)
        leaf_edge = np.copy(leaf_edge)

        pruned = ~pruned
    else:
        pruned = skel

    return pruned

def db_eval_boundary(foreground_mask,gt_mask,bound_th=0.008):
	"""
	Compute mean,recall and decay from per-frame evaluation.
	Calculates precision/recall for boundaries between foreground_mask and
	gt_mask using morphological operators to speed it up.

	Arguments:
		foreground_mask (ndarray): binary segmentation image.
		gt_mask         (ndarray): binary annotated image.

	Returns:
		F (float): boundaries F-measure
		P (float): boundaries precision
		R (float): boundaries recall
	"""
	assert np.atleast_3d(foreground_mask).shape[2] == 1

	bound_pix = bound_th if bound_th >= 1 else \
			np.ceil(bound_th*np.linalg.norm(foreground_mask.shape))

	# Get the pixel boundaries of both masks
	fg_boundary = seg2bmap(foreground_mask);
	gt_boundary = seg2bmap(gt_mask);

	from skimage.morphology import binary_dilation,disk

	fg_dil = binary_dilation(fg_boundary,disk(bound_pix))
	gt_dil = binary_dilation(gt_boundary,disk(bound_pix))

	# Get the intersection
	gt_match = gt_boundary * fg_dil
	fg_match = fg_boundary * gt_dil

	# Area of the intersection
	n_fg     = np.sum(fg_boundary)
	n_gt     = np.sum(gt_boundary)

	#% Compute precision and recall
	if n_fg == 0 and  n_gt > 0:
		precision = 1
		recall = 0
	elif n_fg > 0 and n_gt == 0:
		precision = 0
		recall = 1
	elif n_fg == 0  and n_gt == 0:
		precision = 1
		recall = 1
	else:
		precision = np.sum(fg_match)/float(n_fg)
		recall    = np.sum(gt_match)/float(n_gt)

	# Compute F measure
	if precision + recall == 0:
		F = 0
	else:
		F = 2*precision*recall/(precision+recall);

	return F

def dilating3D(data, selem=skimor.disk(3), slicewise=False, sliceId=0):
    if slicewise:
        if sliceId == 0:
            for i in range(data.shape[0]):
                data[i, :, :] = skimor.binary_dilation(data[i, :, :], selem)
        elif sliceId == 2:
            for i in range(data.shape[2]):
                data[:, :, i] = skimor.binary_dilation(data[:, :, i], selem)
    else:
        data = scindimor.binary_dilation(data, selem)
    return data

def process_cell(img):

    # la binariza en caso de que sea escala de grises
    if not img.dtype == 'bool':
        img = img > 0  # Binarizar

    # Calcular máscaras para limpiar lineas largas verticales
    h_k = 0.8
    sum0 = np.sum(img, 0)  # Aplastar la matriz a una fila con las sumas de los valores de cada columna.
    thr0 = sum0 < h_k * img.shape[0]
    thr0 = thr0.reshape(len(thr0), 1) # Convertirlo a vector de una dimensión

    # Calcular máscaras para limpiar lineas largas horizontales
    w_k = 0.5
    sum1 = np.sum(img, 1)
    thr1 = sum1 < w_k * img.shape[1]
    thr1 = thr1.reshape(len(thr1), 1)

    mask = thr0.transpose() * thr1 # Generar máscara final para la celda
    mask_lines = mask.copy()

    elem = morphology.square(5)
    mask = morphology.binary_erosion(mask, elem) # Eliminar ruido

    img1 = np.bitwise_and(mask, img) # Imagen filtrada

    # segmentación del bloque de números
    kerw = 5  # Kernel width
    thr_k = 0.8

    # Calcular mascara para marcar inicio y fin de región con dígitos horizontalmente
    sum0 = np.sum(img1, 0)
    sum0 = signal.medfilt(sum0, kerw)
    thr0 = sum0 > thr_k * np.median(sum0)
    thr0 = np.bitwise_and(thr0.cumsum() > 0, np.flipud(np.flipud(thr0).cumsum() > 0))
    thr0 = thr0.reshape(len(thr0), 1)

    # Calcular mascara para marcar inicio y fin de región con dígitos verticalmente
    sum1 = np.sum(img1, 1)
    sum1 = signal.medfilt(sum1, kerw)
    thr1 = sum1 > thr_k * np.median(sum1)
    thr1 = np.bitwise_and(thr1.cumsum() > 0, np.flipud(np.flipud(thr1).cumsum() > 0))
    thr1 = thr1.reshape(len(thr1), 1)

    # Mascara final para inicio y fin de caracteres (bounding box of digit region)
    mask = thr0.transpose() * thr1
    mask = morphology.binary_dilation(mask, morphology.square(2))


    img = np.bitwise_and(mask_lines.astype(img.dtype), img)  # Aplicar máscara para quitar lineas
    img = morphology.binary_dilation(img, morphology.disk(1)) # Dilatación para unir números quebrados por la máscara anterior
    img = morphology.binary_erosion(img, morphology.disk(1)) # Volver a la fomorma 'original' con los bordes unidos

    return np.bitwise_and(mask, img)

def double_dilation(binary, selem):
    '''Returns the result of two sequential binary dilations'''

    for i in (1,2):
        binary = binary_dilation(binary, selem)

    return binary

    def compute_fluor_baseline(self, mask, fluor, margin):
        """mask and fluor are the global images
           NOTE: mask is 0 (black) at cells and 1 (white) outside
        """

        x0, y0, x1, y1 = self.box
        wid, hei = mask.shape
        x0 = max(x0 - margin, 0)
        y0 = max(y0 - margin, 0)
        x1 = min(x1 + margin, wid - 1)
        y1 = min(y1 + margin, hei - 1)
        mask_box = mask[x0:x1, y0:y1]

        count = 0

        inverted_mask_box = 1 - mask_box

        while count < 5:
            inverted_mask_box = morphology.binary_dilation(inverted_mask_box)
            count += 1

        mask_box = 1 - inverted_mask_box

        fluor_box = fluor[x0:x1, y0:y1]
        self.stats["Baseline"] = np.median(mask_box[mask_box > 0] * fluor_box[mask_box > 0])

def calculate_masked_stats():
    plate_no = "59798"
    parsed = get_plate_files(plate_no)
    for w in ['w2']:
        files = filter(lambda f: f.wave == w[1], parsed)
        # accum = np.zeros((2160, 2160), dtype=np.uint32)
        # files = filter(lambda x: 's1' not in x and 's7' not in x, all_files)
        nof = len(files)
        for i, frame in enumerate(files[0:5], 1):
            LogHelper.logText(frame.fullpath)
            img = imread(frame.fullpath)
            t = filters.threshold_yen(img)
            b1 = img > t
            b2 = binary_erosion(b1, square(2))
            b3 = binary_dilation(b2, square(10))
            b4 = binary_closing(b3, square(3))
            imm = np.ma.masked_where(b4, img)
            mn, mx = np.percentile(imm, (1, 99))
            LogHelper.logText(
                '%3d of %d, %4d-%4d-%4d-%5d, %.0f-%.0f'
                % (i, nof, imm.min(), mn, mx, imm.max(), imm.mean(), imm.std())
            )
            im2 = imm.filled(int(imm.mean()))
            out_name = "{0}\\{5}-{1}{2}-{3}-{4}.tif".format(ROOT_DIR, frame.row, frame.column, frame.site, LogHelper.init_ts, frame.experiment)
            imsave(out_name, im2)

def draw_gray_tree(frame):
    """
    use a grayscale copy of the frame to draw a quadtree on the original frame
    """
    tree = trees.tree_edges(grayscale(frame))
    tree = morphology.binary_dilation(tree)
    return color_mask(frame, np.logical_not(tree))

def draw_tree(frame):
    """
    draw the edges of a quadtree on the frame
    """
    tree = trees.tree_edges(frame)
    tree = morphology.binary_dilation(tree)
    return color_mask(frame, np.logical_not(tree))

def neg_tree(frame):
    """
    draw a tree in negative
    """
    tree = trees.tree_edges(frame)
    tree = morphology.binary_dilation(tree)
    return color_mask(frame, tree)

def estimate_rotation(img):
    assert(img.dtype == 'bool')

    # elimina bloques rellenos para acelerar la deteccion de lineas
    elem = morphology.square(2)
    aux = morphology.binary_dilation(img, elem) - morphology.binary_erosion(img, elem)

    # Detección de lineas usando transformada de Hough probabilística
    thres = 50
    minlen = 0.1 * min(aux.shape)
    maxgap = 0.01 * minlen
    lines = transform.probabilistic_hough(aux, threshold=thres, line_length=minlen, line_gap=maxgap)

    # me aseguro que el primer punto de cada línea sea el más próximo al origen
    for lin in lines:
        (x0,y0), (x1,y1) = lin
        if x1*x1+y1*y1 < x0*x0+y0*y0:
            (x0, x1) = (x1, x0)
            (y0, y1) = (y1, y0)

    # orientación dominante
    angle_half_range = np.math.pi / 4
    nbins = int(2 * angle_half_range * (180./np.math.pi) / 0.2)

    orient = []
    for lin in lines:
        (x0,y0), (x1,y1) = lin
        orient.append(np.math.atan2(y1-y0, x1-x0))

    (h, binval) = np.histogram(orient, range=(-angle_half_range, angle_half_range), bins=nbins)
    alpha = binval[h.argmax()] * (180./ np.math.pi)
    return alpha + 0.5 * (binval[1] - binval[0]) * (180./ np.math.pi)

def split_image_into_sudoku_pieces_adaptive_global(image, otsu_local=False, apply_gaussian=False):
    L = image.shape[0]
    d = int(np.ceil(L / 9))
    dd = d // 5
    output = []
    if apply_gaussian:
        image = gaussian_filter(image, sigma=1.0)
    if not otsu_local:
        image = to_binary_adaptive(image)
    for k in range(9):
        this_row = []
        start_row_i = max([k * d - dd, 0])
        stop_row_i = min([(k + 1) * d + dd, L])
        for kk in range(9):
            start_col_i = max([kk * d - dd, 0])
            stop_col_i = min([(kk + 1) * d + dd, L])
            i = image[start_row_i:stop_row_i, start_col_i:stop_col_i].copy()
            if otsu_local:
                i = to_binary_otsu(i)
            i = binary_opening(i)
            i = to_binary_otsu(i)
            if apply_gaussian:
                i = to_binary_otsu(binary_dilation(i))
            this_row.append(i)
        output.append(this_row)
    return output, image

def detect_edges(image_array):
    """ Detect edges in a given image
    Takes a numpy.array representing an image,
    apply filters and edge detection and return a numpy.array

    Parameters
    ----------
    image_array : ndarray (2D)
        Image data to be processed. Detect edges on this 2D array representing the image

    Returns
    -------
    edges : ndarray (2D)
        Edges of an image.
    """
    #Transform image into grayscale
    img = rgb2gray(image_array)
    #Remove some noise from the image
    img = denoise_tv_chambolle(img, weight=0.55)
    #Apply canny
    edges = filter.canny(img, sigma=3.2)
    #Clear the borders
    clear_border(edges, 15)
    #Dilate edges to make them more visible and connected
    edges = binary_dilation(edges, selem=diamond(3))
    return edges

def binarize_canny(pic_source, sensitivity = 5.):

    ht = 5. + ((10 - sensitivity)/5.)*20.

#    print ht

    edges = canny_filter(pic_source, sigma = 3, high_threshold = ht, low_threshold = 2.)

    selem_morph = np.array([0,1,0,1,1,1,0,1,0], dtype=bool).reshape((3,3))

    for i in (1,2):
        edges = binary_dilation(edges, selem_morph)

#    misc.imsave('/home/varnivey/Data/Biophys/Burnazyan/Experiments/fluor_calc/test/edges.jpg', edges)

#    binary = ndimage.binary_fill_holes(edges)

    labels = measure_label(edges)

    labelcount = np.bincount(labels.ravel())

    bg = np.argmax(labelcount)

    edges[labels != bg] = 255

    selem_med = np.ones((3,3), dtype = bool)

    binary = median_filter(edges, selem_med)

    for i in (1,2,3):
        binary = binary_erosion(edges, selem_morph)

    return edges

    def mask(self, dims=None, binary=True, outline=False):
        """
        Construct a mask from a source, either locally or within a larger image.

        Parameters
        ----------
        dims : list or tuple, optional, default = None
            Dimensions of large image in which to draw mask. If none, will restrict
            to the bounding box of the region.

        binary : boolean, optional, deafult = True
            Whether to incoporate values or only show a binary mask

        outline : boolean, optional, deafult = False
            Whether to only show outlines (derived using binary dilation)
        """
        coords = self.coordinates

        if dims is None:
            extent = self.bbox[len(self.center):] - self.bbox[0:len(self.center)] + 1
            m = zeros(extent)
            coords = (coords - self.bbox[0:len(self.center)])
        else:
            m = zeros(dims)

        if hasattr(self, 'values') and self.values is not None and binary is False:
            m[coords.T.tolist()] = self.values
        else:
            m[coords.T.tolist()] = 1

        if outline:
            from skimage.morphology import binary_dilation
            m = binary_dilation(m, ones((3, 3))) - m

        return m

def single_out_annotation(base_image, small_cc_image):
    """ extracting individual annotations :
    starting from potential annotation + noise, we remove the noise and
     consolidate annotation area, then return the coordinates of center of
     potential annotations"""
    import numpy as np

    # remove small stuff
    filtered_small_cc, removed_small_cc_small = remove_small_ccomponents(
        small_cc_image, size_closing=5, hist_thres=120)
    # plot_image(removed_small_cc_small)

    # dilate
    from skimage.morphology import binary_dilation, disk
    dilation_radius = 10
    small_cc_cleaned_mask = binary_dilation(filtered_small_cc, disk(dilation_radius))
    # plot_image(small_cc_cleaned_mask)

    # label connected compoenents
    from skimage.morphology import label
    from skimage.measure import regionprops

    markers, n_label = label(small_cc_cleaned_mask, connectivity=1, background=0, return_num=True)

    # for each cc, defines a region
    image_for_region = (base_image*255).astype(np.uint8)
    region_prop = regionprops(markers, image_for_region)

    # for each region, do something

    return region_prop