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 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 detect_tc_in_step(self, nc, i, ecc_th=0.75):
     mask = self.ocean_mask.copy()
     uas = nc.variables["uas"][i].squeeze()
     vas = nc.variables["vas"][i].squeeze()
     wind_speed = numpy.sqrt(uas**2+vas**2)
     wind_mask = logical_and(self.ocean_mask, wind_speed > 20.)
     temp = nc.variables["ts"][i].squeeze()
     temp_mask = logical_and(self.ocean_mask, temp > 298.15)
     ps = nc.variables["ps"][i].squeeze()
     ps_mask = logical_and(self.ocean_mask, ps < 1005)
     mask = logical_or(wind_mask, logical_and(temp_mask, ps_mask))
     mask = remove_small_objects(mask, 20)
     lbl = label(mask)
     props_windspeed = regionprops(lbl, wind_speed)
     props_pressure = regionprops(lbl, ps)
     centroids = []
     for windspeed, pressure in zip(props_windspeed, props_pressure):
         max_wind_speed = windspeed["max_intensity"]
         min_pressure = pressure["min_intensity"]
         if windspeed["eccentricity"] > ecc_th or max_wind_speed<20.:
             lbl[lbl == windspeed["label"]]=0
         else:
             y, x = windspeed["centroid"]
             lon = float(self.idx_to_lon(x, y))
             lat = float(self.idx_to_lat(x, y))
             centroids.append([lon, lat, max_wind_speed, min_pressure])
     mask = lbl>0
     return mask, centroids

    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 filter_segments(labels, max_ecc, min_area, max_area, max_detect=None,
                    circ=None, intensity=None, **extra_args):
    """ filter_segments(labels, max_ecc=0.5, min_area=15, max_area=200) -> [Segment]
        Returns a list of Particles and masks out labels for
        particles not meeting acceptance criteria.
    """
    pts = []
    strengths = []
    centroid = 'Centroid' if intensity is None else 'WeightedCentroid'
    if skversion < version('0.10'):
        rprops = regionprops(labels, ['Area', 'Eccentricity', centroid], intensity)
    else:
        rprops = regionprops(labels, intensity)
    for rprop in rprops:
        area = rprop['area']
        if area < min_area or area > max_area:
            continue
        ecc = rprop['eccentricity']
        if ecc > max_ecc:
            continue
        x, y = rprop[centroid]
        if circ:
            co, ro = circ
            if (x - co[0])**2 + (y - co[1])**2 > ro**2:
                continue
        pts.append(Segment(x, y, rprop.label, ecc, area))
        if max_detect is not None:
            strengths.append(rprop['mean_intensity'])
    if max_detect is not None:
        pts = pts[np.argsort(-strengths)]
    return pts[:max_detect]

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 test_orientation():
    orientation = regionprops(SAMPLE, ['Orientation'])[0]['Orientation']
    # determined with MATLAB
    assert_almost_equal(orientation, 0.10446844651921)
    # test correct quadrant determination
    orientation2 = regionprops(SAMPLE.T, ['Orientation'])[0]['Orientation']
    assert_almost_equal(orientation2, math.pi / 2 - orientation)

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 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_bbox():
    bbox = regionprops(SAMPLE, ['BoundingBox'])[0]['BoundingBox']
    assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]))

    SAMPLE_mod = SAMPLE.copy()
    SAMPLE_mod[:, -1] = 0
    bbox = regionprops(SAMPLE_mod, ['BoundingBox'])[0]['BoundingBox']
    assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1))

def test_euler_number():
    en = regionprops(SAMPLE, ['EulerNumber'])[0]['EulerNumber']
    assert en == 0

    SAMPLE_mod = SAMPLE.copy()
    SAMPLE_mod[7, -3] = 0
    en = regionprops(SAMPLE_mod, ['EulerNumber'])[0]['EulerNumber']
    assert en == -1

def test_filled_area():
    area = regionprops(SAMPLE, ['FilledArea'])[0]['FilledArea']
    assert area == np.sum(SAMPLE)

    SAMPLE_mod = SAMPLE.copy()
    SAMPLE_mod[7, -3] = 0
    area = regionprops(SAMPLE_mod, ['FilledArea'])[0]['FilledArea']
    assert area == np.sum(SAMPLE)

def _properties2d(image, dim):
    """
    Compute shape property of the input 2D image. Accounts for partial volume information.
    :param image: 2D input image in uint8 or float (weighted for partial volume) that has a single object.
    :param dim: [px, py]: Physical dimension of the image (in mm). X,Y respectively correspond to AP,RL.
    :return:
    """
    upscale = 5  # upscale factor for resampling the input image (for better precision)
    pad = 3  # padding used for cropping
    # Check if slice is empty
    if not image.any():
        logging.debug('The slice is empty.')
        return None
    # Normalize between 0 and 1 (also check if slice is empty)
    image_norm = (image - image.min()) / (image.max() - image.min())
    # Convert to float64
    image_norm = image_norm.astype(np.float64)
    # Binarize image using threshold at 0. Necessary input for measure.regionprops
    image_bin = np.array(image_norm > 0.5, dtype='uint8')
    # Get all closed binary regions from the image (normally there is only one)
    regions = measure.regionprops(image_bin, intensity_image=image_norm)
    # Check number of regions
    if len(regions) > 1:
        logging.debug('There is more than one object on this slice.')
        return None
    region = regions[0]
    # Get bounding box of the object
    minx, miny, maxx, maxy = region.bbox
    # Use those bounding box coordinates to crop the image (for faster processing)
    image_crop = image_norm[np.clip(minx-pad, 0, image_bin.shape[0]): np.clip(maxx+pad, 0, image_bin.shape[0]),
                 np.clip(miny-pad, 0, image_bin.shape[1]): np.clip(maxy+pad, 0, image_bin.shape[1])]
    # Oversample image to reach sufficient precision when computing shape metrics on the binary mask
    image_crop_r = transform.pyramid_expand(image_crop, upscale=upscale, sigma=None, order=1)
    # Binarize image using threshold at 0. Necessary input for measure.regionprops
    image_crop_r_bin = np.array(image_crop_r > 0.5, dtype='uint8')
    # Get all closed binary regions from the image (normally there is only one)
    regions = measure.regionprops(image_crop_r_bin, intensity_image=image_crop_r)
    region = regions[0]
    # Compute area with weighted segmentation and adjust area with physical pixel size
    area = np.sum(image_crop_r) * dim[0] * dim[1] / upscale ** 2
    # Compute ellipse orientation, rotated by 90deg because image axis are inverted, modulo pi, in deg, and between [0, 90]
    orientation = _fix_orientation(region.orientation)
    # Find RL and AP diameter based on major/minor axes and cord orientation=
    [diameter_AP, diameter_RL] = \
        _find_AP_and_RL_diameter(region.major_axis_length, region.minor_axis_length, orientation,
                                 [i / upscale for i in dim])
    # TODO: compute major_axis_length/minor_axis_length by summing weighted voxels along axis
    # Fill up dictionary
    properties = {'area': area,
                  'diameter_AP': diameter_AP,
                  'diameter_RL': diameter_RL,
                  'centroid': region.centroid,
                  'eccentricity': region.eccentricity,
                  'orientation': orientation,
                  'solidity': region.solidity  # convexity measure
    }

    return properties

def test_get_boundaries_of_image_3d():
    # Test if equivalent diameter of the maximum intensity project of edges of the object is same
    # as the input sphere, measure.regionprops, 3D perimeter parameter not implemented in skimage
    radius = 4
    binary = morphology.ball(radius)
    boundary = radius_skeleton.get_boundaries_of_image(binary)
    maxip = np.amax(boundary, 0)
    nose.tools.assert_almost_equal(measure.regionprops(binary)[0].equivalent_diameter,
                                   measure.regionprops(maxip)[0].equivalent_diameter, places=1)

def bin_analyser(RGB_image, bin_image, list_feature, marge=None, pandas_table=False, do_label=True):
    bin_image_copy = bin_image.copy()

    p = 0
    for feat in list_feature:
        p += feat.size

    if marge is not None and marge != 0:
        seed = np.zeros_like(bin_image_copy)
        seed[marge:-marge, marge:-marge] = 1
        mask = bin_image_copy.copy()
        mask[ mask > 0 ] = 1
        mask[marge:-marge, marge:-marge] = 1
        reconstructed = reconstruction(seed, mask, 'dilation')
        bin_image_copy[reconstructed == 0] = 0
    if do_label:
        bin_image_copy = label(bin_image_copy)

    if len(np.unique(bin_image_copy)) != 2:
        if len(np.unique(bin_image_copy)) == 1:
            if 0 in bin_image_copy:
                print "Return blank matrix. Change this shit"
                white_npy = np.zeros(shape=(1, p))
                if not pandas_table:
                    return white_npy
                else:
                    names = GetNames(list_feature) 
                    return pd.DataFrame(white_npy, columns=names)
            else:
                print "Error, must give a bin image."
    
    GrowRegion_N = NeededGrownRegion(list_feature)
    img = {0: bin_image_copy}
    RegionProp = {0: regionprops(bin_image_copy)}
    for val in GrowRegion_N:
        if val != 0:
            img[val] = GrowRegion(bin_image_copy, val)
            RegionProp[val] = regionprops(img[val])
    
    n = len(RegionProp[0])

    TABLE = np.zeros(shape=(n,p))
    for i in range(n):
        offset_ALL = 0
        for j, feat in enumerate(list_feature):
            tmp_regionprop = RegionProp[feat._return_n_extension()][i]
            off_tmp = feat.size      
            TABLE[i, (j + offset_ALL):(j + offset_ALL + off_tmp)] = feat._apply_region(tmp_regionprop ,RGB_image)

            offset_ALL += feat.size - 1

    if pandas_table:
        names = GetNames(list_feature)
        return pd.DataFrame(TABLE, columns=names)
    else:
        return TABLE

def rag_solidity(labels, connectivity=2):

    graph = RAG()

    # The footprint is constructed in such a way that the first
    # element in the array being passed to _add_edge_filter is
    # the central value.
    fp = ndi.generate_binary_structure(labels.ndim, connectivity)
    for d in range(fp.ndim):
        fp = fp.swapaxes(0, d)
        fp[0, ...] = 0
        fp = fp.swapaxes(0, d)

    # For example
    # if labels.ndim = 2 and connectivity = 1
    # fp = [[0,0,0],
    #       [0,1,1],
    #       [0,1,0]]
    #
    # if labels.ndim = 2 and connectivity = 2
    # fp = [[0,0,0],
    #       [0,1,1],
    #       [0,1,1]]

    ndi.generic_filter(
        labels,
        function=_add_edge_filter,
        footprint=fp,
        mode='nearest',
        output=np.zeros(labels.shape, dtype=np.uint8),
        extra_arguments=(graph,))

    regions = regionprops(labels)
    regions = {r.label: r for r in regionprops(labels)}

    graph.remove_node(0)

    for n in graph:
        region = regions[n]
        graph.node[n].update({'labels': [n],
                              'solidity': region['solidity'],
                              'mask': labels == region.label})

    for x, y, d in graph.edges_iter(data=True):
        new_mask = np.logical_or(graph.node[x]['mask'], graph.node[y]['mask'])
        new_solidity = 1. * new_mask.sum() / convex_hull_image(new_mask).sum()
        org_solidity = np.mean([graph.node[x]['solidity'],
                                graph.node[y]['solidity']])
        d['weight'] = org_solidity / new_solidity

    return graph

def extractSegmentPropertiesRGB(segments,img):
	"""Function to extract 6 properties from each segment of a segmented image.

    Returns a numpy array with 6 columns and a row per segment.
    [:,1] : Mean red (r) intensity of segment
    [:,2] : Mean green (g) intensity of segment
    [:,3] : Mean blue (b) intensity of segment
    [:,4] : (b-r)/(b+r)
    [:,5] : (b-g)/(b+g)
    [:,6] : (g-r)/(2b - g -r)

    Parameters
    ----------
    segments : ndarray, shape(M, N)
        The segmented image
    img : ndarray, shape(M, N, 3)
        The RGB image that has been segmented.


    Returns
    -------
    properties : ndarray, shape(6, n)
        A numpy array.
    """

	# Extract feature properties
	# Mean RGB intensities
	props = measure.regionprops(segments,img[:,:,0], cache=False)
	r=[prop.mean_intensity for prop in props]
	props = measure.regionprops(segments,img[:,:,1], cache=False)
	g=[prop.mean_intensity for prop in props]
	props = measure.regionprops(segments,img[:,:,2], cache=False)
	b=[prop.mean_intensity for prop in props]

	# RGB Ratios
	denom = np.add(b,r)
	denom[denom==0] = 0.05
	Ratio1 = np.divide(np.subtract(b,r) , denom)

	denom = np.add(b,g)
	denom[denom==0] = 0.05
	Ratio2 = np.divide(np.subtract(b,g) , denom)

	denom = np.subtract(np.subtract(np.multiply(b,2),g),r)
	denom[denom==0] = 0.05
	Ratio3 = np.divide(np.subtract(g,r) , denom)

	# Stack properties to array
	properties = np.column_stack((r,g,b,Ratio1,Ratio2,Ratio3))

	return properties

def remove_abnormal_samples(seg, sigma=0.9):
    '''removes abnormal samples based on sig. deviation from mean centroid,
       then removes based on mean size'''
    # Get centroids of samples
    centroids = []
    for idx, lab in enumerate(measure.regionprops(scipy.ndimage.label(seg)[0])):
        centroids.append(lab.centroid)    
    row_vals = [x[0] for x in centroids]
    # Get relevant stats of row values
    mean = np.mean(row_vals)
    std = np.std(row_vals)
    hi = mean + sigma*std
    lo = mean - sigma*std
    # Eliminate sig. deviation from mean
    for idx, lab in enumerate(measure.regionprops(scipy.ndimage.label(seg)[0])):
        centroid = lab.centroid
        if centroid[0] < lo or centroid[0] > hi:
            seg = floodfill_fast(
                    seg, 
                    x=int(centroid[0]), 
                    y=int(centroid[1]), 
                    value=0, 
                    border_color=0, 
                    dtype=np.uint8
                    )[0]
              
    # Get sizes of samples
    areas = []
    for idx, lab in enumerate(measure.regionprops(scipy.ndimage.label(seg)[0])):
        areas.append(lab.filled_area) 
    mean = np.mean(areas)
    std = np.std(areas)
    hi = mean + 3*sigma*std
    lo = mean - 3*sigma*std  
    # Eliminate sig. deviation from mean
    for idx, lab in enumerate(measure.regionprops(scipy.ndimage.label(seg)[0])):
        area = lab.filled_area
        centroid = lab.centroid
        if area < lo or area > hi:
            seg = floodfill_fast(
                    seg, 
                    x=int(centroid[0]), 
                    y=int(centroid[1]), 
                    value=0, 
                    border_color=0, 
                    dtype=np.uint8
                    )[0]      
    return seg

def computeITCList(evaluation_mask, resolution, level):
    """Compute the list of labels containing Isolated Tumor Cells (ITC)
    
    Description:
        A region is considered ITC if its longest diameter is below 200µm.
        As we expanded the annotations by 75µm, the major axis of the object 
        should be less than 275µm to be considered as ITC (Each pixel is 
        0.243µm*0.243µm in level 0). Therefore the major axis of the object 
        in level 5 should be less than 275/(2^5*0.243) = 35.36 pixels.
        
    Args:
        evaluation_mask:    The evaluation mask
        resolution:         Pixel resolution of the image at level 0
        level:              The level at which the evaluation mask was made
        
    Returns:
        Isolated_Tumor_Cells: list of labels containing Isolated Tumor Cells
    """
    max_label = np.amax(evaluation_mask)    
    properties = measure.regionprops(evaluation_mask)
    Isolated_Tumor_Cells = [] 
    threshold = 275/(resolution * pow(2, level))
    for i in range(0, max_label):
        if properties[i].major_axis_length < threshold:
            Isolated_Tumor_Cells.append(i+1)
    return Isolated_Tumor_Cells

def getImgData(img):
    labelPos = np.array((img < 0.9), dtype=int)
    # labelNeg = np.array((img >= 0.9),dtype=int)
    props = [
        "Area",
        "Centroid",
        "WeightedCentroid",
        "WeightedMoments",
        "WeightedHuMoments",
        "HuMoments",
        "EulerNumber",
        "Eccentricity",
        "EquivDiameter",
        "Extent",
        "MeanIntensity",
        "MinIntensity",
        "MaxIntensity",
        "Orientation",
        "Solidity",
    ]
    # props = ['Centroid','WeightedCentroid']
    dataPos = regionprops(labelPos, properties=props, intensity_image=img)[0]
    del dataPos["Label"]
    # dataNeg = regionprops(labelNeg,properties=props,intensity_image=img)[0]
    # del dataNeg['Label']
    return dataPos  # ,dataNeg