def label_image(image):
    
    ROI = np.zeros((470,400,3), dtype=np.uint8)
    for c in range(3):
        for i in range(50,520):
            for j in range(240,640):
                ROI[i-50,j-240,c] = image[i,j,c]

    
    gray_ROI = cv2.cvtColor(ROI,cv2.COLOR_BGR2GRAY)
    
    ROI_flou = cv2.medianBlur((ROI).astype('uint8'),3)
    
    Laser = Detecte_laser.Detect_laser(ROI_flou)
    
    open_laser = cv2.morphologyEx(Laser, cv2.MORPH_DILATE, disk(3))
    
    skel = skeletonize(open_laser > 0)
    
    tranche = Detecte_laser.tranche(skel,90,30)    
    
    ret, thresh = cv2.threshold(gray_ROI*tranche.astype('uint8'),0,1,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
    thresh01 = thresh<1.0
    
    open_thresh = cv2.morphologyEx(thresh01.astype('uint8'), cv2.MORPH_OPEN, disk(10))
    
    labelised = (label(open_thresh,8,0))+1
    
    return gray_ROI,labelised

def test_compare_8bit_vs_16bit():
    # filters applied on 8-bit image ore 16-bit image (having only real 8-bit
    # of dynamic) should be identical

    image8 = util.img_as_ubyte(data.camera())
    image16 = image8.astype(np.uint16)
    assert_equal(image8, image16)

    methods = [
        "autolevel",
        "bottomhat",
        "equalize",
        "gradient",
        "maximum",
        "mean",
        "subtract_mean",
        "median",
        "minimum",
        "modal",
        "enhance_contrast",
        "pop",
        "threshold",
        "tophat",
    ]

    for method in methods:
        func = getattr(rank, method)
        f8 = func(image8, disk(3))
        f16 = func(image16, disk(3))
        assert_equal(f8, f16)

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 returnProcessedImage(que,folder,img_flist):
	X = []
	for fname in img_flist:
		cur_img = imread(folder+'/'+fname , as_grey=True)
		cur_img = 1 - cur_img

		######## randomly add samples

		# random add contrast
		r_for_eq = random()
		cur_img = equalize_adapthist(cur_img,ntiles_x=8,ntiles_y=8,clip_limit=(r_for_eq+0.5)/3)

		
		#random morphological operation
		r_for_mf_1 = random()
		if 0.05 < r_for_mf_1 < 0.25: # small vessel
			selem1 = disk(0.5+r_for_mf_1)
			cur_img = dilation(cur_img,selem1)
			cur_img = erosion(cur_img,selem1)
		elif 0.25 < r_for_mf_1 < 0.5: # large vessel
			selem2 = disk(2.5+r_for_mf_1*3)
			cur_img = dilation(cur_img,selem2)
			cur_img = erosion(cur_img,selem2)
		elif 0.5 < r_for_mf_1 < 0.75: # exudate
			selem1 = disk(9.21)
			selem2 = disk(7.21)
			dilated1 = dilation(cur_img, selem1)
			dilated2 = dilation(cur_img, selem2)
			cur_img = np.subtract(dilated1, dilated2)
		
		cur_img = img_as_float(cur_img)
		X.append([cur_img.tolist()])
	# X = np.array(X , dtype = theano.config.floatX)
	que.put(X)
	return X

def find_grains(input_file, output_dir=None):
    """Return tuple of segmentaitons (grains, difficult_regions)."""
    name = fpath2name(input_file)
    name = "grains-" + name + ".png"
    if output_dir:
        name = os.path.join(output_dir, name)

    image = Image.from_file(input_file)
    intensity = mean_intensity_projection(image)
    image = remove_scalebar(intensity, np.mean(intensity))
    image = threshold_abs(image, 75)
    image = invert(image)
    image = fill_holes(image, min_size=500)
    image = erode_binary(image, selem=disk(4))
    image = remove_small_objects(image, min_size=500)
    image = dilate_binary(image, selem=disk(4))

    dist = distance(image)
    seeds = local_maxima(dist)
    seeds = dilate_binary(seeds)  # Merge spurious double peaks.
    seeds = connected_components(seeds, background=0)

    segmentation = watershed_with_seeds(dist, seeds=seeds, mask=image)
    # Need a copy to avoid over-writing original.
    initial_segmentation = np.copy(segmentation)

    # Remove spurious blobs.
    segmentation = remove_large_segments(segmentation, max_size=3000)
    segmentation = remove_small_segments(segmentation, min_size=500)
    props = skimage.measure.regionprops(segmentation)
    segmentation = remove_non_round(segmentation, props, 0.6)

    difficult = initial_segmentation - segmentation

    return segmentation, difficult

    def compute(self, src):
        image = img_as_ubyte(src)

        # denoise image
        denoised = denoise_tv_chambolle(image, weight=0.05)
        denoised_equalize= exposure.equalize_hist(denoised)

        # find continuous region (low gradient) --> markers
        markers = rank.gradient(denoised_equalize, disk(5)) < 10
        markers = ndi.label(markers)[0]

        # local gradient
        gradient = rank.gradient(denoised, disk(2))

        # labels
        labels = watershed(gradient, markers)

        # display results
        fig, axes = plt.subplots(2,3)
        axes[0, 0].imshow(image)#, cmap=plt.cm.spectral, interpolation='nearest')
        axes[0, 1].imshow(denoised, cmap=plt.cm.spectral, interpolation='nearest')
        axes[0, 2].imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
        axes[1, 0].imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest')
        axes[1, 1].imshow(labels, cmap=plt.cm.spectral, interpolation='nearest', alpha=.7)
        plt.show()

def find_grains(input_file, output_dir=None):
    """Return tuple of segmentaitons (grains, difficult_regions)."""
    name = fpath2name(input_file)
    name = "grains-" + name + ".png"
    if output_dir:
        name = os.path.join(output_dir, name)

    image = Image.from_file(input_file)
    intensity = mean_intensity_projection(image)

# Median filter seems more robust than Otsu.
#   image = threshold_otsu(intensity)
    image = threshold_median(intensity, scale=0.8)

    image = invert(image)
    image = erode_binary(image, selem=disk(2))
    image = dilate_binary(image, selem=disk(2))
    image = remove_small_objects(image, min_size=200)
    image = fill_holes(image, min_size=50)

    dist = distance(image)
    seeds = local_maxima(dist)
    seeds = dilate_binary(seeds)  # Merge spurious double peaks.
    seeds = connected_components(seeds, background=0)

    segmentation = watershed_with_seeds(dist, seeds=seeds, mask=image)

    # Remove spurious blobs.
    segmentation = remove_large_segments(segmentation, max_size=3000)
    segmentation = remove_small_segments(segmentation, min_size=100)

    return segmentation

def extract_region_opening(img, is_demo=False):
    """
    Extracts fingerprint region of image via mophological opening
    """

    after_median = skimage.filter.rank.median(img, skmorph.disk(9))
    after_erode = skmorph.erosion(after_median, skmorph.disk(11))
    after_dil = skmorph.dilation(after_erode, skmorph.disk(5))
    _, t_dil_img = cv2.threshold(after_dil, 240, 40, cv2.THRESH_BINARY)

    if is_demo:
        _, t_med_img = cv2.threshold(after_median, 240, 255, cv2.THRESH_BINARY)
        _, t_erd_img = cv2.threshold(after_erode, 240, 40, cv2.THRESH_BINARY)
        erd_gry = t_erd_img.astype(np.uint8) * 255
        rgb_erd = np.dstack((erd_gry, img, img))
        dil_gry = t_dil_img.astype(np.uint8) * 255
        rgb_dil = np.dstack((dil_gry, img, img))

        plt.subplot(2,2,1)
        plt.imshow(after_erode, cmap="gray", interpolation="nearest")

        plt.subplot(2,2,2)
        plt.imshow(rgb_erd, interpolation="nearest")

        plt.subplot(2,2,3)
        plt.imshow(after_dil, cmap="gray", interpolation="nearest")

        plt.subplot(2,2,4)
        plt.imshow(rgb_dil, interpolation="nearest")
        plt.show()

    return t_dil_img

    def split_object(self, labeled_image):
        """ split object when it's necessary
        """
        
        labeled_image = labeled_image.astype(np.uint16)

        labeled_mask = np.zeros_like(labeled_image, dtype=np.uint16)
        labeled_mask[labeled_image != 0] = 1

        #ift structuring element about center point. This only affects eccentric structuring elements (i.e. selem with even num===============================
       
        labeled_image = skr.median(labeled_image, skm.disk(4))
        labeled_mask = np.zeros_like(labeled_image, dtype=np.uint16)
        labeled_mask[labeled_image != 0] = 1
        distance = scipym.distance_transform_edt(labeled_image).astype(np.uint16)
        #=======================================================================
        # binary = np.zeros(np.shape(labeled_image))
        # binary[labeled_image > 0] = 1
        #=======================================================================
        distance = skr.mean(distance, skm.disk(15))
         
        l_max = skr.maximum(distance, skm.disk(5))
        #l_max = skf.peak_local_max(distance, indices=False,labels=labeled_image, footprint=np.ones((3,3)))
        l_max = l_max - distance <= 0
        
        l_max = skr.maximum(l_max.astype(np.uint8), skm.disk(6))
        
       
        
        marker = ndimage.label(l_max)[0]
        split_image = skm.watershed(-distance, marker)
        
        split_image[split_image[0,0] == split_image] = 0
        
        return split_image

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 watershed(image):
    hsv_image = color.rgb2hsv(image)

    low_res_image = rescale(hsv_image[:, :, 0], SCALE)
    local_mean = mean(low_res_image, disk(50))
    local_minimum_flat = np.argmin(local_mean)
    local_minimum = np.multiply(np.unravel_index(local_minimum_flat, low_res_image.shape), round(1 / SCALE))

    certain_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_bone_pixels[
    local_minimum[0] - INITIAL_WINDOW_SIZE/2:local_minimum[0]+INITIAL_WINDOW_SIZE/2,
    local_minimum[1] - INITIAL_WINDOW_SIZE/2:local_minimum[1]+INITIAL_WINDOW_SIZE/2
    ] = True

    certain_non_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_non_bone_pixels[0:BORDER_SIZE, :] = True
    certain_non_bone_pixels[-BORDER_SIZE:-1, :] = True
    certain_non_bone_pixels[:, 0:BORDER_SIZE] = True
    certain_non_bone_pixels[:, -BORDER_SIZE:-1] = True

    smoothed_hsv = median(hsv_image[:, :, 0], disk(50))
    threshold = MU * np.median(smoothed_hsv[certain_bone_pixels])

    possible_bones = np.zeros_like(hsv_image[:, :, 0])
    possible_bones[smoothed_hsv < threshold] = 1

    markers = np.zeros_like(possible_bones)
    markers[certain_bone_pixels] = 1
    markers[certain_non_bone_pixels] = 2

    labels = morphology.watershed(-possible_bones, markers)

    return labels

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 segm(img):
    den = rank.median(img, morph.disk(3))
    # continuous regions (low gradient)
    markers = rank.gradient(den, morph.disk(5)) < 10
    mrk, tot = ndimage.label(markers)
    grad = rank.gradient(den, morph.disk(2))
    labels = morph.watershed(grad, mrk)
    return seg_regions(img, labels, tot)

def watershed_segmentation(im):
    gray=skimage.color.rgb2gray(im)
    denoised = rank.median(gray, disk(2))
    # find continuous region (low gradient) --> markers
    markers = rank.gradient(denoised, disk(5)) < 10
    markers = nd.label(markers)[0]
    seg = watershed(gray,markers)
    return seg

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 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 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 smooth(self):
        # TODO: there is non nan in the ff img, or?
        mask = self.flatField == 0
        from skimage.filters.rank import median, mean
        from skimage.morphology import disk

        ff = mean(median(self.flatField, disk(5), mask=~mask),
                  disk(13), mask=~mask)

        return ff.astype(float) / ff.max(), mask

def EllipsoidFit(image,threshold = "otsu global",voxelscale=[1,1,1],
	outfile=None,minsize = 30,tol = 0.01,
	outputdata = True,output3D = False, outputhist = True,
	display3D=False, displayhist=False, hull = True,
	fittype = "ellipsoid",smooth = "median"):
	
	starttime = time.time()
	# If outfile is not provided, set automatically
	if outfile is None:
		if isinstance(image,str):
			outfile = image.split(".")[0]
		else:
			outfile = time.asctime(time.localtime(time.time()))
	# Read image into array
	if isinstance(image,str): 
		imgtiff = tiff.TIFFfile(image)
		image = imgtiff.asarray()

	ishape = image.shape #image shape
	if smooth == "median":
		for i in range(ishape[0]):
			image[i] = filters.rank.median(image[i],morphology.disk(1))
	elif smooth == "mean":
		for i in range(ishape[0]):
			image[i] = filters.rank.mean(image[i],morphology.disk(1))
	print image.shape
	# Rescale voxels to be cubic
	if len(set(voxelscale)) > 1:
		image = cubifyvoxels(image,voxelscale)
	print image.shape

	# Binarize images and label objects
	binimg, thresh = thresholding(image,threshold)
	lbimg,objnum = ndimage.label(binimg)
	print "Objects: ",objnum
	print "Fitting Ellipsoids:"
	# Fit ellipsoids
	ellipsoids = get_ellipsoids(lbimg,objnum,minsize=minsize,tol=tol,hull=hull,voxelscale=voxelscale)
	print " Done!"
	print "Calculating Ellipsoid Information:"
	# Calculate information on ellipsoids
	data = elldatacalc(ellipsoids)
	print "Done!"
	endtime = time.time()
	# Output various results of the fitting
	print "Graphing and Writing to files: "
	if outputdata:
		datafileprint(data,outfile)
		summaryfileprint(data,tol,minsize,outfile,threshold,hull,objnum,thresh,endtime-starttime,smooth)
	if display3D or output3D:
		graphrods(data,binimg,outfile,display3D,output3D,voxelscale=voxelscale)
	if displayhist or outputhist:
		if len(data)>1: graphdata(data,outfile,displayhist,outputhist)
	print "Done!"
	return data

    def test_compare_8bit_vs_16bit(self, method):
        # filters applied on 8-bit image ore 16-bit image (having only real 8-bit
        # of dynamic) should be identical
        image8 = util.img_as_ubyte(data.camera())[::2, ::2]
        image16 = image8.astype(np.uint16)
        assert_equal(image8, image16)

        func = getattr(rank, method)
        f8 = func(image8, disk(3))
        f16 = func(image16, disk(3))
        assert_equal(f8, f16)