def compute_mask(self, params):
        """Creates the mask for the base image.
        Needs the base image, an instance of imageloaderparams
        and the clip area, which should be already defined
        by the load_base_image method.
        Creates the mask by improving the base mask created by the
        compute_base_mask method. Applies the mask closing, dilation and
        fill holes parameters.
        """
        self.compute_base_mask(params)

        mask = np.copy(self.base_mask)
        closing_matrix = np.ones((params.mask_closing, params.mask_closing))

        if params.mask_closing > 0:
            # removes small dark spots and then small white spots
            mask = img_as_float(morphology.closing(
                mask, closing_matrix))
            mask = 1 - \
                img_as_float(morphology.closing(
                    1 - mask, closing_matrix))

        for f in range(params.mask_dilation):
            mask = morphology.erosion(mask, np.ones((3, 3)))

        if params.mask_fill_holes:
            # mask is inverted
            mask = 1 - img_as_float(ndimage.binary_fill_holes(1.0 - mask))

        self.mask = mask

        self.overlay_mask_base_image()

def close_image(img, mask_length):
    # Morphological closing on greyscale/binary image
    
    img = img.astype(np.uint8)
    img = closing(img, rectangle(mask_length,1))
    
    return(img)

def getRegions():
    """Geocode address and retreive image centered
    around lat/long"""
    address = request.args.get('address')
    results = Geocoder.geocode(address)
    lat, lng = results[0].coordinates
    zip_code = results[0].postal_code

    map_url = 'https://maps.googleapis.com/maps/api/staticmap?center={0},{1}&size=640x640&zoom=19&sensor=false&maptype=roadmap&&style=visibility:simplified|gamma:0.1'
    request_url = map_url.format(lat, lng)
    req = urllib.urlopen(request_url)
    img = io.imread(req.geturl(),flatten=True)
    labels, numobjects = ndimage.label(img)
    image = filter.canny(img, sigma=3)
    thresh = threshold_otsu(image)
    bw = closing(image > thresh, square(3))

    # remove artifacts connected to image border
    cleared = bw.copy()
    clear_border(cleared)

    # label image regions
    label_image = label(cleared)
    borders = np.logical_xor(bw, cleared)
    label_image[borders] = -1
    image_label_overlay = label2rgb(label_image, image=image)

    fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
    ax.imshow(image_label_overlay)

def segment_out_cells(base):
    # TODO: try using OTSU for GFP thresholding

    sel_elem = disk(2)
    gfp_collector = np.sum(base, axis=0)
    gfp_clustering_markers = np.zeros(gfp_collector.shape, dtype=np.uint8)
    # random walker segment
    gfp_clustering_markers[gfp_collector > np.mean(gfp_collector) * 2] = 2
    gfp_clustering_markers[gfp_collector < np.mean(gfp_collector) * 0.20] = 1
    labels = random_walker(gfp_collector, gfp_clustering_markers, beta=10, mode='bf')
    # round up the labels and set the background to 0 from 1
    labels = closing(labels, sel_elem)
    labels -= 1
    # prepare distances for the watershed
    distance = ndi.distance_transform_edt(labels)
    local_maxi = peak_local_max(distance,
                                indices=False,  # we want the image mask, not peak position
                                min_distance=10,  # about half of a bud with our size
                                threshold_abs=10,  # allows to clear the noise
                                labels=labels)
    # we fuse the labels that are close together that escaped the min distance in local_maxi
    local_maxi = ndi.convolve(local_maxi, np.ones((5, 5)), mode='constant', cval=0.0)
    # finish the watershed
    expanded_maxi_markers = ndi.label(local_maxi, structure=np.ones((3, 3)))[0]
    segmented_cells_labels = watershed(-distance, expanded_maxi_markers, mask=labels)

    # log debugging data
    running_debug_frame.gfp_collector = gfp_collector
    running_debug_frame.gfp_clustering_markers = gfp_clustering_markers
    running_debug_frame.labels = labels
    running_debug_frame.segmented_cells_labels = segmented_cells_labels

    return gfp_collector, segmented_cells_labels

def Mask_ROI_cl(im,disk_size,thresh=None,black_spots=None,with_morph=False):
    l=np.array([0,0])
    if not isinstance(im,l.__class__):
        numpy_array=np.array(im)
    else:
        numpy_array=im 
    if len(numpy_array.shape)==3:
        numpy_array=numpy_array[:,:,0:3].mean(axis=2)

    selem = disk(disk_size)
    closed = closing(numpy_array, selem)
    if thresh is None:
        thresh = threshold_otsu(closed)
    binary = closed > thresh
    if binary.dtype=='bool':
        binary=binary+0
    if black_spots is not None:    
        binary2 = closed > black_spots
        binary2 = binary2 + 0
        binary = binary - binary2 
    else:
        binary -=1
    binary=binary * -1
    if with_morph:
        return(binary,closed)
    else:
        return(binary)

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 test_large_radius():
    ''' Compare execution time against scikit: single closing case
    Here, our implementation does not take advantage of smaller radius results
    so ours is slower than scikit, but it uses significantly less memory.
    '''
    base_dir = '/home/omar/data/DATA_NeoBrainS12/'
    neo_subject = '30wCoronal/example2/'

    # Read subject files
    t2CurrentSubjectName  = base_dir + 'trainingDataNeoBrainS12/'+neo_subject+'T2_1-1.nii.gz'
    t2CurrentSubject_data = nib.load(t2CurrentSubjectName).get_data()
    affineT2CS            = nib.load(t2CurrentSubjectName).get_affine()
    zoomsT2CS             = nib.load(t2CurrentSubjectName).get_header().get_zooms()[:3]
    # Step 1.4 - Resampling for isotropic voxels

    n_zooms = (zoomsT2CS[0],zoomsT2CS[0],zoomsT2CS[0])
    t2CurrentSubject_data,affineT2CS = reslice(t2CurrentSubject_data,affineT2CS,zoomsT2CS,n_zooms)

    S = t2CurrentSubject_data.astype(np.float64)

    ###########compare times#########
    # in-house
    radius = 15
    start = time.time()
    d = isotropic_dilation(S, radius)
    c = isotropic_erosion(d, radius)
    end = time.time()
    print('Elapsed (in-home): %f'%(end-start,))

    # scikit
    start = time.time()
    expected = closing(S, ball(radius))
    end = time.time()
    print('Elapsed (scikit): %f'%(end-start,))

def threshold_image(image, threshold=0):
	"""
	This function takes out any values in an image's RGB matrix that are
	below the threshold value.

	Inputs:
	- image: a matrix describing an image with only one channel represented.
	- threshold: a value, between 0 and 1, for which if an image matrix's
				 value is below, will be set to 0, and if above, will be 
				 set to 1.

				 If the threshold is set to 0, then an Otsu thresholding will
				 be returned.

	Outputs:
	- thresholded_image: a matrix representation of the thresholded image.
						 this is essentially a black and white image.
	- thresh: the threshold value

	To screen: the black-and-white image representation.
	- 
	"""
	if threshold == 0:
		thresh = threshold_otsu(image)

	if threshold != 0:
		thresh = threshold

	thresholded_image = closing(image > thresh, square(3), out=None)
	imshow(thresholded_image)

	return thresholded_image, thresh

def main():
    for file_path in glob.glob("/home/lucas/Downloads/Lucas/GSK 10uM/*.JPG"):

        img = data.imread(file_path, as_grey=True)

        img = transform.resize(img, [600, 600])
        img_color = transform.resize(data.imread(file_path), [600, 600])

        img[img >img.mean()-0.1] = 0

        # io.imshow(img)
        # io.show()
        #
        edges = canny(img)
        bordas_fechadas = closing(img > 0.1, square(15)) # fechando gaps
        fill_cells = ndi.binary_fill_holes(bordas_fechadas)
        # io.imshow(fill_cells)
        # io.show()
        img_label = label(fill_cells, background=0)
        n= 0
        for  x in regionprops(img_label):
            if x.area < 2000 and x.area > 300:
                n +=1
                print x.area
                minr, minc, maxr, maxc = x.bbox
                try:
                    out_path_name = file_path.split("/")[-1].rstrip(".JPG")
                    io.imsave("out/cell_{}_pic_{}_area_{}.png".format(n, out_path_name, str(round(x.area))),img_color[minr-3: maxr+3, minc-3: maxc+3])
                    #io.show()
                except:
                    pass

 def morphOps( imgIn, sizeSE, sizeCC ):
     imgOut = imgIn.astype(bool) #boolean image
     imgOut = ~imgOut #img negative
     imgOut = morphology.remove_small_objects( imgOut, sizeCC ) #cclargest
     SE = morphology.selem.disk( sizeSE ) #structuring element
     imgOut = morphology.closing(imgOut, SE)
     return imgOut

    def get_valley_image(self, image):
        valley_img = np.zeros_like(image)
        for z in range(0, image.shape[0]):
            valley_img[z, :, :] = closing(image[z, :, :], disk(5))
        valley_img -= image

        return valley_img

def smooth(image):
  filename_split = os.path.splitext(image)
  filename_zero, fileext = filename_split
  basename = os.path.basename(filename_zero)
  im = np.array(Image.open(image))
  with rasterio.open(image) as r:
    im = r.read()
    p = r.profile
  im = im.squeeze()
  selem = disk(1)
  print("image min and max: ", im.min(),im.max())
  dilated = skimage.morphology.dilation(im, selem)
  print("dilated image min and max: ", dilated.min(),dilated.max())
  eroded = skimage.morphology.erosion(im, selem)
  print("eroded image min and max: ", eroded.min(),eroded.max())
  opened = opening(im, selem)
  print("opened image min and max: ", opened.min(),opened.max())
  closed = closing(im, selem)
  print("closed image min and max: ", closed.min(),closed.max())
  #im[im==1]=0
  #im[im==2]=1
  median = cv2.medianBlur(im,9)
  average = cv2.blur(im,(9,9))
  #gaussian = cv2.GaussianBlur(im,(9,9),0)
  gaussian = cv2.GaussianBlur(dilated,(9,9),0)
  #bilateral = cv2.bilateralFilter(im,9,75,75)
  bilateral = cv2.bilateralFilter(gaussian,9,75,75)
  with rasterio.open(outPath+basename+fileext, 'w', **p) as dst:
      dst.write(bilateral, 1)
  color_outPath = outPath+'color/'
  if not os.path.exists(color_outPath):
        os.mkdir(color_outPath) 
  colored_image = color_outPath+basename+'.png'
  os.system("gdaldem color-relief", bilateral, colorfile, colored_image)
  return im, dilated, eroded, opened, closed, median, average, gaussian, bilateral

def upsample_smooth(image):
  filename_split = os.path.splitext(image)
  filename_zero, fileext = filename_split
  basename = os.path.basename(filename_zero)
  upsampleRes = int(upsampleRes)
  upsampled_image = outPath+basename+'_cubicSpline.png'
  os.system("gdalwarp -tr", upsample_res, upsample_res," -r cubicspline ", image, upsampled_image)
  im = np.array(Image.open(upsampled_image))
  with rasterio.open(image) as r:
    im = r.read()
    p = r.profile
  im = im.squeeze()
  selem = disk(1)
  print("image min and max: ", im.min(),im.max())
  dilated = skimage.morphology.dilation(im, selem)
  print("dilated image min and max: ", dilated.min(),dilated.max())
  eroded = skimage.morphology.erosion(im, selem)
  print("eroded image min and max: ", eroded.min(),eroded.max())
  opened = opening(im, selem)
  print("opened image min and max: ", opened.min(),opened.max())
  closed = closing(im, selem)
  print("closed image min and max: ", closed.min(),closed.max())
  dilated = Image.fromarray(dilated)
  dilated.save(outPath+basename+'.png')
  with rasterio.open(outPath+basename+fileext, 'w', **p) as dst:
      dst.write(dilated, 1)
  color_outPath = outPath+'color/'
  if not os.path.exists(color_outPath):
        os.mkdir(color_outPath) 
  colored_image = color_outPath+basename+'.png'
  os.system("gdaldem color-relief", dilated, colorfile, colored_image)
  return im, dilated, eroded, opened, closed

def scikit_example_plot_label():
    image = data.coins()[50:-50, 50:-50]
    
    # apply threshold
    thresh = threshold_otsu(image)
    bw = closing(image > thresh, square(3))
    
    # remove artifacts connected to image border
    cleared = bw.copy()
    clear_border(cleared)
    
    # label image regions
    label_image = label(cleared)
    borders = np.logical_xor(bw, cleared)
    label_image[borders] = -1
    
    fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
    ax.imshow(label_image, cmap='jet')
    
    for region in regionprops(label_image, ['Area', 'BoundingBox']):
    
        # skip small images
        if region['Area'] < 100:
            continue
    
        # draw rectangle around segmented coins
        minr, minc, maxr, maxc = region['BoundingBox']
        rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
                                  fill=False, edgecolor='red', linewidth=2)
        ax.add_patch(rect)
    
    plt.show()

def test_accuracy():
    ''' Verify that our implementation returns exactly the same as scikit
    '''
    base_dir = '/home/omar/data/DATA_NeoBrainS12/'
    neo_subject = '30wCoronal/example2/'

    # Read subject files
    t2CurrentSubjectName  = base_dir + 'trainingDataNeoBrainS12/'+neo_subject+'T2_1-1.nii.gz'
    t2CurrentSubject_data = nib.load(t2CurrentSubjectName).get_data()
    affineT2CS            = nib.load(t2CurrentSubjectName).get_affine()
    zoomsT2CS             = nib.load(t2CurrentSubjectName).get_header().get_zooms()[:3]

    n_zooms = (zoomsT2CS[0],zoomsT2CS[0],zoomsT2CS[0])
    t2CurrentSubject_data,affineT2CS = reslice(t2CurrentSubject_data,affineT2CS,zoomsT2CS,n_zooms)

    S = t2CurrentSubject_data.astype(np.float64)

    max_radius = 4
    D = SequencialSphereDilation(S)
    for r in range(1, 1+max_radius):
        D.expand(S)
        expected = dilation(S, ball(r))
        actual = D.get_current_dilation()
        assert_array_equal(expected, actual)
        expected = closing(S, ball(r))
        actual = D.get_current_closing()
        assert_array_equal(expected, actual)

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 = erosion(lines, square(13))

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

    return removedChessboard

def Image_ws_tranche(image):
    
    laser = Detect_laser(image)
    laser_tranche = tranche_image(laser,60)
    
    image_g = skimage.color.rgb2gray(image)
    image_g = image_g * laser_tranche
    
    image_med = rank2.median((image_g*255).astype('uint8'),disk(8))
    
    image_clahe = exposure.equalize_adapthist(image_med, clip_limit=0.03)
    image_clahe_stretch = exposure.rescale_intensity(image_clahe, out_range=(0, 256))

    image_grad = rank2.gradient(image_clahe_stretch,disk(3))
    
    image_grad_mark = image_grad<20
    image_grad_forws = rank2.gradient(image_clahe_stretch,disk(1))
    
    image_grad_mark_closed = closing(image_grad_mark,disk(1))
    
    Labelised = (skimage.measure.label(image_grad_mark_closed,8,0))+1
    Watersheded  = watershed(image_grad_forws,Labelised)
    
    cooc = coocurence_liste(Watersheded,laser,3)
    
    x,y = compte_occurences(cooc)
    
    return x,y

def im_proc(im):
    """Apply series of morphological procedures on image."""
    th = threshold_otsu(im)
    im_bin = im > th
    return(ndi.binary_fill_holes(
                morphology.closing(
                im_bin,np.ones((3,3)))))

    def plot_preprocessed_image(self):
        """
        plots pre-processed image. The plotted image is the same as obtained at the end
        of the get_text_candidates method.
        """
        image = restoration.denoise_tv_chambolle(self.image, weight=0.1)
        thresh = threshold_otsu(image)
        bw = closing(image > thresh, square(2))
        cleared = bw.copy()

        label_image = measure.label(cleared)
        borders = np.logical_xor(bw, cleared)

        label_image[borders] = -1
        image_label_overlay = label2rgb(label_image, image=image)

        fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(12, 12))
        ax.imshow(image_label_overlay)

        for region in regionprops(label_image):
            if region.area < 10:
                continue

            minr, minc, maxr, maxc = region.bbox
            rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
                                      fill=False, edgecolor='red', linewidth=2)
            ax.add_patch(rect)

        plt.show()

def roofRegion(edge):
    """Estimate region based on edges of roofRegion
    """
    # apply threshold
    thresh = threshold_otsu(image)
    bw = closing(image > thresh, square(3))

    # remove artifacts connected to image border
    cleared = bw.copy()
    clear_border(cleared)

    # label image regions
    label_image = label(cleared)
    borders = np.logical_xor(bw, cleared)
    label_image[borders] = -1
    image_label_overlay = label2rgb(label_image, image=image)

    fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
    ax.imshow(image_label_overlay)

    for region in regionprops(label_image):

        # skip small images
        if region.area < 100:
            continue

        # draw rectangle around segmented coins
        minr, minc, maxr, maxc = region.bbox
        rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
                                  fill=False, edgecolor='red', linewidth=2)
        ax.add_patch(rect)

    plt.show()