def watershed(base_image, seed_image=None, threshold_distance=80):
    """ execute watershed with chosen seeds
    """
    from scipy import ndimage as ndi
    
    from skimage.morphology import watershed
    from skimage.feature import peak_local_max
    from skimage.morphology import label
    import matplotlib.pyplot as plt
    
    
    distance = ndi.distance_transform_edt(base_image)
    #    imgplot = plt.imshow(distance)
    fig = plt.figure()  # a new figure window
    ax = fig.add_subplot(1, 1, 1)  # specify (nrows, ncols, axnum)
    #    ax.imshow(distance>threshold_distance, cmap='Greys')
    
    thresh = distance > threshold_distance
    ax.imshow(thresh, cmap='Greys') 
    
    
    
    #local_maxi = peak_local_max(distance, labels=jac, footprint=np.ones((100, 100)), indices=False)
    if seed_image is None: 
        markers = label(thresh)
    else:
        markers = label(seed_image)

    #    imgplot = plt.imshow(markers)
     
      
    watersh = watershed(-distance, markers, mask=base_image) 
    #    plt.imshow(watersh, cmap=plt.cm.viridis, interpolation='nearest')
    
    return watersh

    def test_return_num(self):
        x = np.array([[1, 0, 6],
                      [0, 0, 6],
                      [5, 5, 5]])

        assert_array_equal(label(x, return_num=True)[1], 4)
        assert_array_equal(label(x, background=0, return_num=True)[1], 3)

def segmentationize(imageSe):
    """
    Divides coherent forms of an image in smaller groups of type integer.
    """
    
    # create an matrix of distances to the next sourrounding area
    distance = ndimage.distance_transform_edt(imageSe, sampling=3)
    erosed = ndimage.binary_erosion(imageSe, iterations=8).astype(imageSe.dtype)
    distanceE = ndimage.distance_transform_edt(erosed, sampling=3)
    distance += (2 * distanceE)
    labels, num = label(imageSe, background=0, return_num='True')
    sizes_image = ndimage.sum(imageSe, labels, range(num))
    sizes_image = np.sort(sizes_image, axis=None)
    pos = int(0.4 * num)
    areal = int(sizes_image[pos] ** 0.5)
    if areal <= 10:
        areal = 10
    elif (areal % 2) != 0:
        areal += 1
    footer = circarea(areal) # draw circle area
    
    # find the positions of the maxima from the distances
    local_maxi = peak_local_max(distance, indices=False, footprint=footer, labels=imageSe)
    markers = label(local_maxi)
    
    # watershed algorithm starts at the maxima and returns labels of particles
    simplefilter("ignore", FutureWarning)   # avoid warning in watershed method
    labels_ws = watershed(-distance, markers, mask=imageSe)
    simplefilter("default", FutureWarning)
    
    return labels, labels_ws, local_maxi

    def clear_body(self, body, minimum_object_size_px=2400):
        """ Vycisti obraz od stolu a ostatnich veci okolo tela a zanecha pouze a jen telo """
        body = scipy.ndimage.filters.gaussian_filter(copy.copy(body).astype(float), sigma=[15, 0, 0]) > 0.7

        # fallowing lines are here to supress warning "Only one label was provided to `remove_small_objects`. "
        blabeled = morphology.label(body)
        if np.max(blabeled) > 1:
            body = morphology.remove_small_objects(morphology.label(blabeled), minimum_object_size_px)
        del blabeled
        
        body[0] = False
        body[-1] = False
        
        body = scipy.ndimage.filters.gaussian_filter(copy.copy(body.astype(float)), sigma=[1, 3, 3]) > 0.2
    
        bodylabel = label(body)
    
        n_of_pixels = [np.count_nonzero(bodylabel==i) for i in range(len(np.unique(bodylabel)))]
        labelsort = np.argsort(n_of_pixels)[::-1]
    
        newbody = np.zeros(body.shape)
        newbody[bodylabel==labelsort[0]] = body[(bodylabel==labelsort[0])]
        newbody[bodylabel==labelsort[1]] = body[(bodylabel==labelsort[1])]
    
        return newbody.astype(bool)

 def test_4_vs_8(self):
     x = np.array([[0, 1],
                   [1, 0]], dtype=int)
     assert_array_equal(label(x, 4),
                        [[0, 1],
                         [2, 3]])
     assert_array_equal(label(x, 8),
                        [[0, 1],
                         [1, 0]])

def analyseClusters(binary, newlabels):
    """
    Calculates the sizes and porosities of the clusters.
    """
    
    # dilate particles to find cluster
    dilated = ndimage.binary_dilation(binary, iterations=_DILATIONFACTOR_TO_FIND_CLUSTER)
    labels, num = label(dilated, background=0, return_num=True)
    pxArea = (_CONVERSIONFACTOR_FOR_PIXEL) ** 2
    outputImage = labels.copy()
    clusterAreas = np.zeros(num)
    porosities = np.zeros(num)
    circumference = np.zeros(num)
    fcirc = np.zeros(num)
    particlesPerCluster = np.zeros(num)
    illegalIndex = []
    
    for i in range(num):
        cluster = labels == i
        cluster = ndimage.binary_fill_holes(cluster)
        helper = np.zeros_like(newlabels)
        helper[cluster] = newlabels[cluster]
        newLabel, particleNum = label(helper, background=0, return_num=True)
        particlesPerCluster[i] = particleNum
        particleArea = float(np.sum(binary[cluster].astype(np.int)))
        
        # cluster area and porosity
        outputImage[cluster] = i
        helper = ndimage.binary_erosion(cluster, iterations=_DILATIONFACTOR_TO_FIND_CLUSTER-3, border_value=1)        
        helper = ndimage.binary_erosion(helper, iterations=3, border_value=0)
        fl = float(np.sum(helper[cluster].astype(np.int)))
        clusterAreas[i] = fl * pxArea
        porosity = (fl - particleArea)/ fl
        porosity = porosity if porosity >= 0 else 0.0  # porosity can not be less than 0
        porosities[i] = porosity
        
        # circumference
        new = np.zeros((helper.shape[0],helper.shape[1],3), dtype=np.uint8)
        new[:,:,1] = helper
        gray = cv2.cvtColor(new, cv2.COLOR_RGB2GRAY)
        gray[gray > 0] = 255
        blur = cv2.GaussianBlur(gray,(5,5),0)
        gray = cv2.Canny(blur, 10, 200)
        contours, hierarchy = cv2.findContours(gray, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
        arclength = 0
        for con in contours:
            arclength += cv2.arcLength(con,True)
        circumference[i] = arclength * _CONVERSIONFACTOR_FOR_PIXEL
        fcirc[i] = (4. * np.pi * fl) / arclength**2
        
        if fcirc[i] > 1.0:  # fcirc can not be greater than 1
            illegalIndex.append(i)
    
    fcirc = np.delete(fcirc, illegalIndex)
    clusterData = {'areas':clusterAreas,'circ':circumference,'ppc':particlesPerCluster,'fcirc':fcirc,'porosities':porosities}
    return outputImage, clusterData, num

 def test_4_vs_8(self):
     x = np.array([[0, 1],
                   [1, 0]], dtype=int)
     with catch_warnings():
         assert_array_equal(label(x, 4),
                            [[0, 1],
                             [2, 3]])
         assert_array_equal(label(x, 8),
                            [[0, 1],
                             [1, 0]])

def label_img(img):
    # Labelling the nests is done using connected components
    img = create_bin(img)

    labeled_img = label(input=img, connectivity=2, background=0)
    #min size holes ina nest
    rem_holes = remove_small_holes(labeled_img, min_size=100, connectivity=2)
    #min size of a nest
    labeled_img1 = remove_small_objects(rem_holes, min_size=70, connectivity=2)
    labeled = label(labeled_img1, connectivity=2, background=0)

    print labeled
    return labeled

    def test_background(self):
        x = np.array([[1, 0, 0],
                      [1, 1, 5],
                      [0, 0, 0]])

        assert_array_equal(label(x), [[0, 1, 1],
                                      [0, 0, 2],
                                      [3, 3, 3]])

        assert_array_equal(label(x, background=0),
                           [[0, -1, -1],
                            [0,  0,  1],
                            [-1, -1, -1]])

    def run(self, workspace):
        cell_object = workspace.object_set.get_objects(self.primary_objects.value)
        cell_labeled = cell_object.get_segmented()
        cell_image = cell_object.get_parent_image()

        cell_image = (cell_image * 1000).astype(np.uint16)
        #        object_count = cell_labeled.max()
        maxi = skr.maximum(cell_image.astype(np.uint8), skm.disk(10))
        local_max = maxi - cell_image < 10

        local_max_labelize, object_count = scipy.ndimage.label(local_max, np.ones((3, 3), bool))
        histo_local_max, not_use = np.histogram(local_max_labelize, range(object_count + 2))
        old = local_max_labelize.copy()

        # filter in intensity mean
        # =======================================================================
        #
        # regionprops_result = skmes.regionprops(local_max_labelize, intensity_image=cell_image)
        #
        # for region in regionprops_result:
        #     if region["mean_intensity"]
        # =======================================================================

        # filter on size
        for i in range(object_count + 1):
            value = histo_local_max[i]
            if value > self.range_size.max or value < self.range_size.min:
                local_max_labelize[local_max_labelize == i] = 0

        # split granule for each cell
        cell_labeled = skm.label(cell_labeled)
        cell_count = np.max(cell_labeled)

        for cell_object_value in range(1, cell_count):
            cell_object_mask = cell_labeled == cell_object_value
            granule_in_cell = np.logical_and(cell_object_mask, local_max_labelize)
            granule_in_cell = skm.label(granule_in_cell)
            # ===================================================================
            # plt.imshow(granule_in_cell + cell_object_mask)
            # plt.show()
            # ===================================================================
        #
        # get the filename
        #
        measurements = workspace.measurements
        file_name_feature = self.source_file_name_feature
        filename = measurements.get_current_measurement("Image", file_name_feature)
        print "filename = ", filename

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 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()

    def separate_segments(self):
        """
        Perform image segmentation on the "segment image", and remove any
        segments that aren't the right part of the track.
        """

        # binary image
        binary_segment_image = (
            self.end_segment_image > self.options.low_threshold_kev)
        # segmentation: labeled regions, 8-connectivity
        labels = morph.label(binary_segment_image, connectivity=2)
        x1 = self.end_coordinates[0] - self.end_segment_offsets[0]
        y1 = self.end_coordinates[1] - self.end_segment_offsets[1]
        x2 = self.start_coordinates[0] - self.end_segment_offsets[0]
        y2 = self.start_coordinates[1] - self.end_segment_offsets[1]
        chosen_label = labels[x1, y1]
        if labels[x2, y2] != chosen_label:
            # this happens with 4-connectivity. need to use 8-connectivity
            raise RuntimeError('What the heck happened?')
        binary_again = (labels == chosen_label)
        # dilate this region, in order to capture information below threshold
        #  (it won't include the other regions, because there must be a gap
        #   between)
        pix_to_keep = morph.binary_dilation(binary_again)
        self.end_segment_image[np.logical_not(pix_to_keep)] = 0

 def FilterChangeByArea(self,area):
     size = np.round(area / (0.4**2)) #devide with square of image resolution
     
     
     
     BinArray = gfh.ReadData(self.chng)
     
     BigChng = np.zeros(np.shape(BinArray))
     
     LabelArray = morphology.label(np.squeeze(BinArray))        
     
     N_chng = np.max(LabelArray)
     if N_chng > 0:
         k = 1
         for i in range(N_chng):
             indx,indy = np.where(LabelArray == i)    
             #larger than x m² and smaller than half the image
             if (len(indx) >= size) & (len(indx) < (self.Ndata/2)): 
                 BigChng[indx,indy] = k
                 k += 1     
         
         print str(k)
         ROI = gfh.ReadData(self.ROI)
         xROI,yROI,_ = np.where(ROI == 0)
         BigChng[xROI,yROI] = -9999
         
         gfh.WriteToFile(self.chng,BigChng)
         
         gfh.TIFF2PNG(self.chng,self.png_path + 'chng_AreaSelect' +
         self.png_fn_tag,'Fit')
     else:
         print 'no changes found, test thresholdvalues!'
     return 0

def process_img (timestamp,img, point):
    img_gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    thr = cv2.adaptiveThreshold(img_gray,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,\
            cv2.THRESH_BINARY,11,2)
    
    thr = 255 - thr
    kernel = np.ones((1,7),np.uint8)
    dilated_img = cv2.dilate(thr, kernel,iterations = 1)
    L = morphology.label(dilated_img)
    #cv2.circle(img,point,5,(255,0,255),3)
    #cv2.imshow('img',img)
    
    cnt = 1
    for region in regionprops(L):
        minr, minc, maxr, maxc = region.bbox
        
        if maxr - minr > 200:
            continue

        if point[1]>=minr and point[1]<=maxr and point[0]>=minc and point[0]<=maxc:
            subimg = img[minr:maxr,minc:maxc]
            
            imgName = 'sub/' + timestamp + '_' + str(cnt) + '.png'
            cnt = cnt + 1
            cv2.imwrite(imgName,subimg)

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""" 
     
    #  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
    from skimage.io import imsave
    markers, n_label = label(small_cc_cleaned_mask, connectivity=1, background=0, return_num=True)
    
    #for each cc, defines a region    
    region_prop = regionprops(markers, (base_image*255).astype(np.uint8))
    
    #for each region, do something
    
    base_path = '/media/sf_RemiCura/PROJETS/belleepoque/extract_data_from_old_paris_map/jacoubet/results/annotations/'
    for region in region_prop: 
        #print(region.bbox, region.area)
        imsave(base_path+str(region.bbox)+'.png', region.intensity_image) 
        
    return region_prop

def get_cells(image):
    '''
    Get cellls from the polygon.
    '''
    new_image=np.ones([3,image.shape[0],image.shape[1]],dtype=float)
    # apply threshold
    thresh = threshold_otsu(image)
    bw=image

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

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

    #extract the regions and get a polygon per region
    polygons=[]
    for i,region in enumerate(regionprops(label_image)):
        # skip small images
        if region.area < 100:
            continue
        #polygons.append(matplotlib.path.Path(region.coords))
        print (region.coords.shape)
        a=np.zeros(region.coords.shape)
        a[:,0]=region.coords[:,1]
        a[:,1]=region.coords[:,0]
        polygons.append(a)   
    return polygons

def large_contiguous_regions(thresholded_image, min_area):
	"""
	NOTE: THIS FUNCTION IS SPECIALLY WRITTEN FOR TRACHEAL SECTIONS.

	This function gets the largest contiguous region of an image, based on 
	its labels.
	"""

	labelled_image = label(thresholded_image)

	regions = np.zeros((np.shape(thresholded_image)[0], np.shape(thresholded_image)[1]))
	maxarea = 0

	for region in np.unique(labelled_image):
		image_area = np.shape(labelled_image)[0] ** 2
		
		region_matrix = labelled_image == region
		region_in_thresholded = np.logical_and(region_matrix, thresholded_image)
		region_in_thresholded_area = np.sum(region_in_thresholded)

		if region_in_thresholded_area > image_area * min_area:

			regions = regions + region_in_thresholded

	return regions		

def remove_none_edge_intersecting(img, edge=0, width=1):

    mask = np.zeros(img.shape,dtype=int)
    out = np.zeros(img.shape,dtype=int)
    # print '--->', img.sum()

    if edge == 0:
        mask[:,0:0+width] = 1
    elif edge == 1:
        mask[:,-1-width:-1] = 1
    elif edge == 2:
        mask[0:width,:] = 1
    elif edge == 3:
        mask[-1-width:-1,:] = 1
    else:
        raise ValueError('Edge is duff')

    s = label(img.astype(int))
    s_set = np.unique(s * mask)
    if s_set.sum() > 0:
        for v in s_set:
            q = (s == v)
            if np.all(img[q]):
                out[s == v] = 1

    return out

def human_afterloop(output_directory, pre_time, fle_name, buffer_directory):
    start2 = time()

    d_c = import_edited(buffer_directory)
    rebw = load(open(buffer_directory+'-'+'DO_NOT_TOUCH_ME.dmp','rb'))
    seg_dc = (label(d_c,neighbors=4)+1)*d_c
    if np.max(seg_dc)<4:
        return 'FAILED: mask for %s looks unsegmented' % fle_name

    colormap = repaint_culsters(int(np.max(seg_dc)))

    segs = len(set(seg_dc.flatten().tolist()))-1

    # shows the result before saving the clustering and printing to the user the number of the images
    plt.subplot(1,2,1)
    plt.title(fle_name)
    plt.imshow(rebw, cmap='gray', interpolation='nearest')

    plt.subplot(1,2,2)
    plt.title('Segmentation - clusters: %s'%str(segs))
    plt.imshow(mark_boundaries(rebw, d_c))
    plt.imshow(seg_dc, cmap=colormap, interpolation='nearest', alpha=0.3)

    plt.show()

    plt.imshow(mark_boundaries(rebw, d_c))
    plt.imshow(seg_dc, cmap=colormap, interpolation='nearest', alpha=0.3)

    plt.savefig(path.join(output_directory, fle_name+'_%s_clusters.png'%str(segs)), dpi=500, bbox_inches='tight', pad_inches=0.0)

    return fle_name+'\t clusters: %s,\t total time : %s'%(segs, "{0:.2f}".format(time()-start2+pre_time))