def clean_by_area(self, binary_image):
        image = binary_image.copy()
        image = ndi.binary_fill_holes(image)

        label_image = label(binary_image)
        initial_label = regionprops(label_image[0, :, :])[0].label

        for z in range(0, image.shape[0]):
            regions = regionprops(label_image[z, :, :])
            for region in regions:
                if region.label != initial_label:
                    for coords in region.coords:
                        image[z, coords[0], coords[1]] = 0

        for z in range(0, image.shape[0]):
            label_image = label(image[z, :, :], connectivity=1)
            regions = regionprops(label_image)
            if len(regions) > 1:
                max_area = np.max([r.area for r in regions])
                for region in regions:
                    if region.centroid[1] > 120 and region.area < max_area:
                        for coords in region.coords:
                            image[z, coords[0], coords[1]] = 0

        return image

def extract_cell_stats(img1_path, img2_path):

    # Function reads in the images and labels the cells. The features are
    # extracted from these labelled images.
    #
    # Inputs:   img1_path - path to previous image
    #           img2_path - path to current image
    #
    # Outputs:  out -   dict containing the relevant information
    #

    # TODO: be more accommodating with image types, RGB etc, tifffile warning
    # read image data
    img1 = skimage.io.imread(img1_path)
    img2 = skimage.io.imread(img2_path)

    # Image shape
    if img1.shape != img2.shape:
        warnings.warn('Caution: Comparing image frames of different sizes.')
    img_shape = img1.shape

    # Label pre-segmented images
    l_label, l_cell_total = label(img1, return_num=True)
    r_label, r_cell_total = label(img2, return_num=True)

    # Collect cell features if cell is of minimum size (not segmented debris)
    # TODO: clever way of setting this number
    l_cells = [cell for cell in regionprops(l_label) if cell['filled_area'] > 50]
    r_cells = [cell for cell in regionprops(r_label) if cell['filled_area'] > 50]

    # Output
    out = {'img1': l_cells, 'img2': r_cells, 'img_shape': img_shape}
    return out

def get_segmentation_features(im):
    dilwindow = [4, 4]
    imthr = np.where(im > np.mean(im), 0.0, 1.0)
    imdil = morphology.dilation(imthr, np.ones(dilwindow))
    labels = measure.label(imdil)
    labels = imthr * labels
    labels = labels.astype(int)
    regions = measure.regionprops(labels)
    numregions = len(regions)
    while len(regions) < 1:
        dilwindow[0] = dilwindow[0] - 1
        dilwindow[1] = dilwindow[1] - 1
        if dilwindow == [0, 0]:
            regions = None
            break
        imthr = np.where(im > np.mean(im), 0.0, 1.0)
        imdil = morphology.dilation(imthr, np.ones(dilwindow))
        labels = measure.label(imdil)
        labels = imthr * labels
        labels = labels.astype(int)
        regions = measure.regionprops(labels)
    regionmax = get_largest_region(regions, labels, imthr)

    if regionmax is None:
        return (np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan)
    eccentricity = regionmax.eccentricity
    convex_area = regionmax.convex_area
    convex_to_total_area = regionmax.convex_area / regionmax.area
    extent = regionmax.extent
    filled_area = regionmax.filled_area
    return (eccentricity, convex_area, convex_to_total_area, extent,
            filled_area, numregions)

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

        lnb = x.copy()
        lnb[0] = np.array([[1, 2, 2],
                           [1, 2, 2],
                           [2, 2, 2]])
        lnb[1] = np.array([[2, 2, 2],
                           [2, 1, 3],
                           [2, 2, 2]])
        lb = x.copy()
        lb[0] = np.array([[1,  BG, BG],
                          [1,  BG, BG],
                          [BG, BG, BG]])
        lb[1] = np.array([[BG, BG, BG],
                          [BG, 1,   2],
                          [BG, BG, BG]])

        assert_array_equal(label(x), lb)
        assert_array_equal(label(x, background=-1), lnb)

    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], 3)
        assert_array_equal(label(x, background=-1, return_num=True)[1], 4)

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

        lnb = x.copy()
        lnb[0] = np.array([[0, 1, 1],
                           [0, 1, 1],
                           [1, 1, 1]])
        lnb[1] = np.array([[1, 1, 1],
                           [1, 0, 2],
                           [1, 1, 1]])
        lb = x.copy()
        lb[0] = np.array([[0,  BG, BG],
                          [0,  BG, BG],
                          [BG, BG, BG]])
        lb[1] = np.array([[BG, BG, BG],
                          [BG, 0,   1],
                          [BG, BG, BG]])

        with expected_warnings(['`background`']):
            assert_array_equal(label(x), lnb)

        assert_array_equal(label(x, background=0), lb)

def incorporate_cells(binary_image):    
    # invert input binary image
    inv_image = np.invert(binary_image)
    
    # matrix for binary_dilation
    struct = generate_binary_structure(2, 1)

    # do bunary dilation until the colony number even out
    plate_bin_dil = binary_dilation(inv_image, structure=struct)
    plate_dil_labels = label(plate_bin_dil)
    labels_number = len(np.unique(plate_dil_labels))  # initial number of colonies
    new_labels_number = labels_number - 1  # starting value
    cycle_number = 0  # starting value for dilation cycles
    while True:
        cycle_number += 1
        if cycle_number >= 30:
            break  # defence against infinite cycling
        else:
            if new_labels_number >= labels_number:
                break   # further dilation is useless (in theory)
            else:
                labels_number = new_labels_number
                plate_bin_dil = binary_dilation(plate_bin_dil, structure=struct)
                plate_dil_labels = label(plate_bin_dil)
                new_labels_number = len(np.unique(plate_dil_labels))
                
    return plate_bin_dil

def dice(img,y_true,y_pred):

    h, w = img.shape
    im_true = y_true.reshape(h, w)
    im_pred = y_pred.reshape(h, w)

    labels_true = measure.label(im_true)
    regions_true = regionprops(labels_true)

    labels_pred = measure.label(im_pred)
    regions_pred = regionprops(labels_pred)
    features = ['coords','area','dice']
    df = pd.DataFrame(columns=features)

    i=0
    for x_pred in regions_pred :
        centroid = (np.array(x_pred.centroid)).astype(int)
        if im_true[(centroid[0], centroid[1])] == 1:
            for x_true in regions_true:
               if centroid in x_true.coords:
                   A = np.zeros((img.shape[0], img.shape[1]))
                   B = np.zeros((img.shape[0], img.shape[1]))


                   A[x_pred.coords[:, 0], x_pred.coords[:, 1]] = 1
                   B[x_true.coords[:, 0], x_true.coords[:, 1]] = 1
                   intersect = float((sum(sum(B))))

                   D = intersect/(sum(sum(B))+ sum(sum(A)))
                   df.loc[i] = [x_pred.coords , x_pred.area, D]
                   break
        i+=1
    return df

    def start(self):
        """Segment the frame.

        The returned value is a labeled uint16 image.
        """
        background = np.bincount(self._frame.ravel()).argmax()  # Most common value.
        I_label = measure.label(self._frame, background=background)
        I_label += 1  # Background is labeled as -1, make it 0.
        I_bin = I_label > 0

        # Remove cells which are too small (leftovers).
        if self._a_min:
            I_label = mh.label(I_bin)[0]
            sizes = mh.labeled.labeled_size(I_label)
            too_small = np.where(sizes < self._a_min)
            I_cleanup = mh.labeled.remove_regions(I_label, too_small)
            I_bin = I_cleanup > 0

        # Fill holes.
        if self._fill:
            I_bin = ndimage.morphology.binary_fill_holes(I_bin)  # Small holes.
            # Bigger holes.
            labels = measure.label(I_bin)
            label_count = np.bincount(labels.ravel())
            background = np.argmax(label_count)
            I_bin[labels != background] = True

        I_label = mh.label(I_bin)[0].astype('uint16')
        return I_label

 def get_rough_detection(self, img, bigsize=40.0, smallsize=4.0, thresh = 0):
     diff = self.difference_of_gaussian(-img, bigsize, smallsize)
     diff[diff>thresh] = 1
     
     se = morphology.square(4)
     ero = morphology.erosion(diff, se)
     
     labimage = label(ero)
     #rec = morphology.reconstruction(ero, img, method='dilation').astype(np.dtype('uint8'))
     
     # connectivity=1 corresponds to 4-connectivity.
     morphology.remove_small_objects(labimage, min_size=600, connectivity=1, in_place=True)
     #res = np.zeros(img.shape)
     ero[labimage==0] = 0
     ero = 1 - ero
     labimage = label(ero)
     morphology.remove_small_objects(labimage, min_size=400, connectivity=1, in_place=True)
     ero[labimage==0] = 0
     res = 1 - ero
     res[res>0] = 255
     
     #temp = 255 - temp
     #temp = morphology.remove_small_objects(temp, min_size=400, connectivity=1, in_place=True)
     #res = 255 - temp
     
     return res

 def test_4_vs_8(self):
     x = np.zeros((2, 2, 2), int)
     x[0, 1, 1] = 1
     x[1, 0, 0] = 1
     label4 = x.copy()
     label4[1, 0, 0] = 2
     assert_array_equal(label(x, 4), label4)
     assert_array_equal(label(x, 8), x)

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

        with expected_warnings(['`background`']):
            assert_array_equal(label(x, return_num=True)[1], 4)

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

 def test_4_vs_8(self):
     x = np.zeros((2, 2, 2), int)
     x[0, 1, 1] = 1
     x[1, 0, 0] = 1
     label4 = x.copy()
     label4[1, 0, 0] = 2
     with expected_warnings(['`background`']):
         assert_array_equal(label(x, 4), label4)
         assert_array_equal(label(x, 8), x)

    def test_basic(self):
        assert_array_equal(label(self.x), self.labels)

        # Make sure data wasn't modified
        assert self.x[0, 2] == 3

        # Check that everything works if there is no background
        assert_array_equal(label(self.x, background=99), self.labels_nobg)
        # Check that everything works if background value != 0
        assert_array_equal(label(self.x, background=9), self.labels_bg_9)

def eval_area(pred_ror, true_ror):
    '''
        面積に対する評価をし、[true_size, pred_size, score]の配列を返す

        それぞれのiごとで、
            T     : 岩の大きさの真値
            P     : pred_rorの大きさの新地
            TnorP : 岩領域で検出できていないピクセル数
            TandP : 検出できているピクセル数
            PnorT : 後検出のピクセル数
    '''
    detect = 0
    total_pre = 0
    total_recall = 0

    pred = sk.label(pred_ror, return_num = False, background=None)
    true = sk.label(true_ror, return_num = False, background=0)

    for i in range(0, np.max(true_ror)+1):

        TandP = np.count_nonzero(pred[true == i]) # pred_ror[true_ror == i]でzeroじゃなかった数 / iの領域のサイズ
        
        if TandP !=0: # もし検出できていれば

            ## Get P
            non = np.nonzero(pred[true == i]) 
            p = np.unique(pred[true == i][non]) ## 被っている領域のpredの番号
            P = 0 # Initialization
            for i2 in p:
                P += (pred == i2).sum()

            ## Get others
            T = (true == i).sum()
            TnorP = (true == i).sum() - np.count_nonzero(pred[true == i])
            PnorT = P - TandP

            ## Get score
            pre = 1. * TandP / P
            recall = 1. * TandP / T
            
            ## renew total score
            total_pre += pre
            total_recall += recall
            detect += 1
            
            ## Draw
            plt.scatter(pre, recall, color = 'b')
            # print T,P,TandP,TnorP,PnorT

    pre_ave    = 1. * total_pre   / detect
    recall_ave = 1. * total_recall/ detect

    return pre_ave, recall_ave

def largest_region(imData):

    belowMeanFilter = np.where(imData > np.mean(imData), 0., 1.0)
    dialated = morphology.dilation(belowMeanFilter, np.ones((3, 3)))
    regionLabels = (belowMeanFilter * measure.label(dialated)).astype(int)

    # calculate common region properties for each region within the segmentation
    regions = measure.regionprops(regionLabels)
    areas = [(None
              if sum(belowMeanFilter[regionLabels == region.label]) * 1.0 / region.area < 0.50
              else region.filled_area)
             for region in regions]

    if len(areas) > 0:

        regionMax = regions[np.argmax(areas)]

        # trim image to the max region
        regionMaxImg = trim_image(
            np.minimum(
                imData*np.where(regionLabels == regionMax.label, 1, 255),
                255))

        # rotate
        angle = intertial_axis(regionMaxImg)[2]
        rotatedRegionMaxImg = ndimage.rotate(regionMaxImg, np.degrees(angle))
        rotatedRegionMaxImg = trim_image(trim_image(rotatedRegionMaxImg, 0), 255)

    else:
        regionMax = None
        rotatedRegionMaxImg = None
        angle = 0

    return regionMax, rotatedRegionMaxImg, angle, regionLabels, regions, areas, belowMeanFilter, dialated

def get_largest_cc(u,v):
    """
    Return mask with largest connected component in u,v

    """

    if not skimage_available:
        print('*** skimage is not available. get_larget_cc() will not work. ***')
        return np.ones_like(u).astype('bool')
    
    fxx = np.array([[1,-2.0,1.0]])
    fxy = np.array([[-0.25,0,0.25],[0.0,0,0],[0.25,0,-0.25]])
    fyy = fxx.T

    u_ = u.astype('float32')
    v_ = v.astype('float32')
    uxx = cv2.filter2D(u_,-1,fxx)
    uxy = cv2.filter2D(u_,-1,fxy)
    uyy = cv2.filter2D(u_,-1,fyy)

    vxx = cv2.filter2D(v_,-1,fxx)
    vxy = cv2.filter2D(v_,-1,fxy)
    vyy = cv2.filter2D(v_,-1,fyy)

    THRESH=0.1
    ue = np.logical_or(np.logical_or(np.abs(uxx)>THRESH, np.abs(uxy)>THRESH),np.abs(uyy)>THRESH)
    ve = np.logical_or(np.logical_or(np.abs(vxx)>THRESH, np.abs(vxy)>THRESH),np.abs(vyy)>THRESH)
    edg = np.logical_or(ue,ve)
    
    L = measure.label(edg.astype('int32'),neighbors=4)
    
    sums = np.bincount(L.ravel())
    biggest_cc = L==np.argmax(sums)
    return biggest_cc

    def get_maxima(self, src):
        '''
        入力された画像を領域分割し、各領域の最大値を算出する。
        src: 1ch-img
        dst: 領域の最大値のピクセルにのみその値が格納された画像。
        '''

        img = copy.deepcopy(src)
        img[img!=0] = 255

        # 各領域にラベルをつける
        labels, num = sk.label(img, return_num = True) 

        seed_img = np.zeros_like(src)
        
        # 各領域の最大値を求める
        for i in range(1,num+1):

            # iの領域だけ残す
            img = copy.deepcopy(src) # 初期に戻す
            img[labels!=i] = 0 # これで残った領域の最大値求める
            
            # 最大値を求める,1行にしたときの値が出てくるからこんなになってる
            y = np.argmax(img)/len(img)
            x = np.argmax(img)%len(img)

            if img[y,x] != 0: # 中に空いた穴じゃなければ種にする
                seed_img[y,x] = src[y,x]
        
        return seed_img

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 class_and_replace(org_label, classes):
    """
    CLASSIFIES AND REPLACES VALUES IN NUMPY ARRAY
    
    org_label = numpy array with original labels
    classes = array of same length as len of org_label with new labels
    
    classified = numpy array with new classified values
    classified_ID = merged regions with a unique ID
    
    """
    # Classify and replace
    keys = np.unique(org_label)
    values = classes #0 for non veg, 1 for veg
    dictionary = dict(zip(keys, values))
    classified = replace(org_label, dictionary)
    
    # Label merged regions with unique ID
    labeld = label(classified)
    number_of_labels = np.unique(labeld)
    newvals = np.asarray(list(range(1,(len(number_of_labels) + 1))))
    keys_ID = number_of_labels
    values_ID = newvals
    dictionary_ID = dict(zip(keys_ID, values_ID))
    classified_ID = replace(labeld, dictionary_ID)
    
    del(labeld)
    
    return classified, classified_ID