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 tissue_region_from_rgb(_img, _min_area=150, _g_th=None):
    """
    TISSUE_REGION_FROM_RGB detects the region(s) of the image containing the
    tissue. The original image is supposed to represent a haematoxylin-eosin
    -stained pathology slide.
    
    The main purpose of this function is to detect the parts of a large image
    which most probably contain tissue material, and to discard the background.
    
    Usage:
        tissue_mask = tissue_from_rgb(img, _min_area=150, _g_th=None)
        
    Args:
        img (numpy.ndarray): the original image in RGB color space
        _min_area (int, default: 150): any object with an area smaller than 
            the indicated value, will be discarded
        _g_th (int, default: None): the processing is done on the GREEN channel
            and all pixels below _g_th are considered candidates for "tissue
            pixels". If no value is given to _g_th, one is computed by K-Means
            clustering (K=2), and is returned.
        
    Returns:
        numpy.ndarray: a binary image containing the mask of the regions
            considered to represent tissue fragments
        int: threshold used for GREEN channel
    """
    
    if _g_th is None:
        # Apply vector quantization to remove the "white" background - work in the
        # green channel:
        vq = MiniBatchKMeans(n_clusters=2)
        _g_th = int(np.round(0.95 * np.max(vq.fit(_G(_img).reshape((-1,1)))
                                           .cluster_centers_.squeeze())))
    
    mask = _G(_img) < _g_th

    skm.binary_closing(mask, skm.disk(3), out=mask)
    mask = img_as_bool(mask)
    mask = skm.remove_small_objects(mask, min_size=_min_area, in_place=True)


    # Some hand-picked rules:
    # -at least 5% H and E
    # -at most 25% background
    # for a region to be considered tissue

    h, e, b = rgb2he2(_img)

    mask &= (h > np.percentile(h, 5)) | (e > np.percentile(e, 5))
    mask &= (b < np.percentile(b, 50))               # at most at 50% of "other components"

    mask = mh.close_holes(mask)

    return img_as_bool(mask), _g_th

def closing3D(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_closing(data[i, :, :], selem)
        elif sliceId == 2:
            for i in range(data.shape[2]):
                data[:, :, i] = skimor.binary_closing(data[:, :, i], selem)
    else:
        data = scindimor.binary_closing(data, selem)
    return data

def label_particles_edge(im, sigma=2, closing_size=0, **extra_args):
    """ Segment image using Canny edge-finding filter.

        parameters
        ----------
        im : image in which to find particles
        sigma : size of the Canny filter
        closing_size : size of the closing filter

        returns
        -------
        labels : an image array of uniquely labeled segments
    """
    from skimage.morphology import square, binary_closing, skeletonize
    if skimage_version < StrictVersion('0.11'):
        from skimage.filter import canny
    else:
        from skimage.filters import canny
    edges = canny(im, sigma=sigma)
    if closing_size > 0:
        edges = binary_closing(edges, square(closing_size))
    edges = skeletonize(edges)
    labels = sklabel(edges)
    print "found {} segments".format(labels.max())
    # in ma.array mask, False is True, and vice versa
    labels = np.ma.array(labels, mask=edges == 0)
    return labels

def needles2tips(needles, mri, number_of_slices=3):
    """
    This function accepts a list of volumes (each volumes containing a needle) and returns a single volume containing just the tips.
    Here there are a plenty of attempts of removing untrusted data.
    """
    tips = np.zeros_like(mri).astype(np.int32)
    #print(tips.shape)
    for needle in needles:
        needle = needle.astype(np.int32)
        #print("MIN %f, MAX %f" % (needle.min(), needle.max()))
        if np.sum(needle) < (np.shape(needle)[0] * np.shape(needle)[1] * np.shape(needle)[2]):
            #print("Valid needle")
            needle = binary_closing(needle, selem=np.ones((3,3,3)))
            needle[needle!=0]=1
            #print(" after closing: MIN %f, MAX %f " % (needle.min(), needle.max()))
            for z in range(np.shape(mri)[0]-1, 0, -1):
                if 200 > np.sum(needle[z,:,:]) > 0.5 and z-number_of_slices-1 >= 0:
                    #print(" valid slice %d" % z)
                    tmp = deepcopy(needle)
                    tmp[:z-number_of_slices-1,:,:] = 0
                    tips[tmp!=0] = 1
                    del(tmp)
                    break
    
    tips[tips!=0]=1
        
    return tips.astype(np.int32)

def removeChessboard(img):

    # Get the major lines in the image
    edges, dilatedEdges, (h, theta, d) = findLines(img)

    # Create image with ones to fill inn lines
    lines = np.ones(img.shape[:2])

    # Add lines to image as zeroes
    for _, angle, dist in zip(*hough_line_peaks(h, theta, d)):
        y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
        y1 = (dist - img.shape[1] * np.cos(angle)) / np.sin(angle)
        x, y = line(int(y1), 0, int(y0), img.shape[1] - 1)
        x = np.clip(x, 0, img.shape[0] - 1)
        y = np.clip(y, 0, img.shape[1] - 1)
        lines[x, y] = 0

    # Remove border edges from image with all edges
    w = 4
    edges = np.pad(edges[w:img.shape[0] - w, w:img.shape[1] - w], w, mode='constant')

    # Erode the lines bigger, such that they cover the original lines
    lines = binary_erosion(lines, square(13))

    # Remove major lines and close shape paths
    removedChessboard = binary_closing(edges * lines, square(8))

    return removedChessboard

    def get_compactness(self, labels):
        nlabels = labels.max() + 1
        # nlabels = len(np.unique(labels)) - 1
        eccs = np.zeros(nlabels)

        for lab in range(nlabels):
            obj = labels == lab
            if labels.ndim == 2:
                strel = np.ones((3, 3), dtype=np.bool)
                obj_c = skimor.binary_closing(obj, strel)
                if obj_c.sum() >= obj.sum():
                    obj = obj_c
            else:
                strel = np.ones((3, 3, 3), dtype=np.bool)
                obj_c = tools.closing3D(obj, strel)
                if obj_c.sum() >= obj.sum():
                    obj = obj_c

            if labels.ndim == 2:
                ecc = skimea.regionprops(obj)[0]['eccentricity']
            else:
                ecc = tools.get_zunics_compatness(obj)
                # if np.isnan(ecc):
                #     ecc = 0
            eccs[lab] = ecc

        return eccs

def image_filter(img):
  img2 = img.copy();
  img2[img2 < 30] = 100;
  img2 = exp.smooth_image(img2, sigma = 1.0);
  #plt.figure(6); plt.clf();
  #plt.imshow(img2);

  # threshold image and take zero smaller components..
  
  th = img2 < 92;

  th2 = morph.binary_closing(th, morph.diamond(1))
  
  label = meas.label(th2, background=0)
  #plt.imshow(mask)
  
  bs = meas.regionprops(label+1);
  area = np.array([prop.area for prop in bs]);
  if len(area) > 0:
    mask = np.logical_and(label > -1, label != np.argsort(area)[-1]);
  
    img2[mask] = 100;
  
  img2[:2,:] = 100; img2[-2:,:] = 100;
  img2[:,:2] = 100; img2[:,-2:] = 100;
  
  #plt.figure(6); plt.clf();
  #plt.subplot(1,2,1);
  #plt.imshow(img2, vmin = 84, vmax = 92, cmap = plt.cm.gray)
  #plt.subplot(1,2,2);
  #plt.imshow(img2);
  return img2;

def segment_roi(roi):
    # step 1. phase congruency (edge detection)
    Mm = phasecong_Mm(roi)
    # step 2. hysteresis thresholding (of edges)
    B = hysthresh(Mm,HT_T1,HT_T2)
    # step 3. trim pixels off border
    B[B[:,1]==0,0]=0
    B[B[:,-2]==0,-1]=0
    B[0,B[1,:]==0]=0
    B[-1,B[-2,:]==0]=0
    # step 4. threshold to find dark areas
    dark = dark_threshold(roi, DARK_THRESHOLD_ADJUSTMENT)
    # step 5. add dark areas back to blob
    B = B | dark
    # step 6. binary closing
    B = binary_closing(B,SE3)
    # step 7. binary dilation
    B = binary_dilation(B,SE2)
    # step 8. thinning
    B = bwmorph_thin(B,3)
    # step 9. fill holes
    B = binary_fill_holes(B)
    # step 10. remove blobs smaller than BLOB_MIN
    B = remove_small_objects(B,BLOB_MIN,connectivity=2)
    # done.
    return B

def closing(gray_img, kernel=None):
    """Wrapper for scikit-image closing functions. Opening can remove small dark spots (i.e. pepper).

    Inputs:
    gray_img = input image (grayscale or binary)
    kernel   = optional neighborhood, expressed as an array of 1s and 0s. If None, use cross-shaped structuring element.

    :param gray_img: ndarray
    :param kernel = ndarray
    :return filtered_img: ndarray
    """

    params.device += 1

    # Make sure the image is binary/grayscale
    if len(np.shape(gray_img)) != 2:
        fatal_error("Input image must be grayscale or binary")

    # If image is binary use the faster method
    if len(np.unique(gray_img)) == 2:
        bool_img = morphology.binary_closing(image=gray_img, selem=kernel)
        filtered_img = np.copy(bool_img.astype(np.uint8) * 255)
    # Otherwise use method appropriate for grayscale images
    else:
        filtered_img = morphology.closing(gray_img, kernel)

    if params.debug == 'print':
        print_image(filtered_img, os.path.join(params.debug_outdir, str(params.device) + '_opening' + '.png'))
    elif params.debug == 'plot':
        plot_image(filtered_img, cmap='gray')

    return filtered_img

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 refine_aseg(aseg, ball_size=4):
    """
    First step to reconcile ANTs' and FreeSurfer's brain masks.

    Here, the ``aseg.mgz`` mask from FreeSurfer is refined in two
    steps, using binary morphological operations:

      1. With a binary closing operation the sulci are included
         into the mask. This results in a smoother brain mask
         that does not exclude deep, wide sulci.

      2. Fill any holes (typically, there could be a hole next to
         the pineal gland and the corpora quadrigemina if the great
         cerebral brain is segmented out).


    """
    # Read aseg data
    bmask = aseg.copy()
    bmask[bmask > 0] = 1
    bmask = bmask.astype(np.uint8)

    # Morphological operations
    selem = sim.ball(ball_size)
    newmask = sim.binary_closing(bmask, selem)
    newmask = binary_fill_holes(newmask.astype(np.uint8), selem).astype(np.uint8)

    return newmask.astype(np.uint8)

def get_segmented_lungs(im, plot=False):
    # Step 1: Convert into a binary image.
    binary = im < -400
    # Step 2: Remove the blobs connected to the border of the image.
    cleared = clear_border(binary)
    # Step 3: Label the image.
    label_image = label(cleared)
    # Step 4: Keep the labels with 2 largest areas.
    areas = [r.area for r in regionprops(label_image)]
    areas.sort()
    if len(areas) > 2:
        for region in regionprops(label_image):
            if region.area < areas[-2]:
                for coordinates in region.coords:
                       label_image[coordinates[0], coordinates[1]] = 0
    binary = label_image > 0
    # Step 5: Erosion operation with a disk of radius 2. This operation is seperate the lung nodules attached to the blood vessels.
    selem = disk(2)
    binary = binary_erosion(binary, selem)
    # Step 6: Closure operation with a disk of radius 10. This operation is    to keep nodules attached to the lung wall.
    selem = disk(10) # CHANGE BACK TO 10
    binary = binary_closing(binary, selem)
    # Step 7: Fill in the small holes inside the binary mask of lungs.
    edges = roberts(binary)
    binary = ndi.binary_fill_holes(edges)
    # Step 8: Superimpose the binary mask on the input image.
    get_high_vals = binary == 0
    im[get_high_vals] = -2000
    return im, binary

def get_segmented_lungs(im):

    binary = im < -320
    cleared = clear_border(binary) 
    cleared=morph(cleared,5)
    label_image = label(cleared)
  
    areas = [r.area for r in regionprops(label_image)]
    areas.sort()
    if len(areas) > 2:
        for region in regionprops(label_image):
            if region.area < areas[-2]:
                for coordinates in region.coords:
                       label_image[coordinates[0], coordinates[1]] = 0
    binary = label_image > 0  
    selem = disk(2)
    binary = binary_erosion(binary, selem)
 
    selem = disk(10)
    binary = binary_closing(binary, selem)
    edges = roberts(binary)
    binary = ndi.binary_fill_holes(edges)
 
    get_high_vals = binary == 0
    im[get_high_vals] = 0
  
    binary = morphology.dilation(binary,np.ones([5,5]))
    return binary

def fill_gaps(image, closing_radius=0, min_hole_size=0, median_radius=0.6):
    """
    This function closes small gaps between and within objects and smooths edges. It is a
    'finishing' step before skeletonization, and improves the quality of the skeleton by removing
    gaps and minimizing bumps. It also enables removing close, parallel objects such as appear under
    the microscope as a single, long, clear object with sharp, parallel edges. These spurrious objects would
    otherwise pass earlier filters but are, in fact, spurrious. The function itself is a wrapper for
    `skimage.morphology.binary_closing`, `skimage.morphology.remove_small_holes`, and `skimage.filters.median`
    on a binary image.

    Parameters
    ----------
    image : ndarray
        Binary image of candidate objects
    closing_radius : ndarray
        Binary structure to perform binary closing. Defaults to 0 (skips).
    min_hole_size : int
        Holes with areas smaller than this (in pixels) are removed. Defaults to 0 (skips).
    median_radius : ndarray
        Binary structure to use for a median filter. Defaults at 0.6, giving square connectivity of 1
        (manhattan = 1). 0 to skip.

    Returns
    -------
    ndarray : Binary image of candidate objects.

    """
    closing_structure = _disk(closing_radius)
    median_structure = _disk(median_radius)

    out = morphology.binary_closing(image, closing_structure)
    out = morphology.remove_small_holes(out, min_size=min_hole_size)
    out = filters.median(out, selem=median_structure)

    return(out)

def get_rois(im_array, markers):
    """Return the regions of interest."""
    salem = disk(3)
    im = get_thresholded_image(im_array)
    im = remove_small_objects(im, min_size=50)
    im = binary_closing(im, salem)
    im = watershed(-im_array, markers, mask=im)
    return im

def outlineSmoothing(image, radius=1):
    """Smooth outlines of foreground object using morphological operations."""

    struct = median(disk(radius), disk(1))
    label_image = binary_closing(image, struct)
    label_image = binary_opening(label_image, struct)

    return label_image.astype(image.dtype)

def mask_numpy_array(channel1, channel2, thresh_o = 20):

    channel2_Otsu = threshold_otsu(channel2)
    channel2_thresh = channel2 > ((channel2_Otsu*thresh_o)/100)
    channel2_thresh = binary_closing(channel2_thresh, square(5))
    channel1_channel2_masked = ma.masked_array(channel1, mask=~channel2_thresh)
    channel1_channel2 = channel1_channel2_masked.filled(0)
    return channel1_channel2, channel2_thresh

def rotate_blob(blob, theta):
    """rotate a blob and smooth out rotation artifacts"""
    blob = rotate(blob,-1*theta,resize=True)
    blob = binary_closing(blob,SE3)
    blob = binary_dilation(blob,SE2)
    # note that H Sosik's version does one iteration
    # of thinning but 3 is closer to area-preserving
    blob = bwmorph_thin(blob,3)
    return blob

def ProcessImage(im, targetDim = 250, doDenoiseOpening = True):

	#Resize to specified pixels max edge size
	scaling = 1.
	if im.shape[0] > im.shape[1]:
		if im.shape[0] != targetDim:
			scaling = float(targetDim) / im.shape[0]
			im = misc.imresize(im, (targetDim, int(round(im.shape[1] * scaling))))
	else:
		if im.shape[1] != targetDim:
			scaling = float(targetDim) / im.shape[1]
			im = misc.imresize(im, (int(round(im.shape[0] * scaling)), targetDim))
	#print "scaling", scaling

	greyim = 0.2126 * im[:,:,0] + 0.7152 * im[:,:,1] + 0.0722 * im[:,:,2]

	#Highlight number plate
	imnorm = np.array(greyim, dtype=np.uint8)
	se = np.ones((3, 30), dtype=np.uint8)
	opim = morph.opening(imnorm, se)
	diff = greyim - opim + 128.

	misc.imsave("diff.png", diff)

	#Binarize image
	vals = diff.copy()
	vals = vals.reshape((vals.size))

	meanVal = vals.mean()
	stdVal = vals.std()
	threshold = meanVal + stdVal

	#print "Threshold", threshold

	binIm = diff > threshold
	misc.imsave("threshold.png", binIm)
	#print vals.shape
	#plt.plot(vals)
	#plt.show()

	#Denoise
	diamond = morph.diamond(2)
	if doDenoiseOpening:
		currentIm = morph.binary_opening(binIm, diamond)
	else:
		currentIm = binIm
	denoiseIm2 = morph.binary_closing(currentIm, np.ones((3, 13)))

	#print "currentIm", currentIm.min(), currentIm.max(), currentIm.mean()
	#print "denoiseIm2", denoiseIm2.min(), denoiseIm2.max(), currentIm.mean()
	#misc.imsave("denoised1.png", currentIm * 255)
	#misc.imsave("denoised2.png", denoiseIm2 * 255)

	#Number candidate regions
	#print "Numbering regions"
	numberedRegions, maxRegionNum = morph.label(denoiseIm2, 4, 0, return_num = True)
	return numberedRegions, scaling