def label_clusters(self):
        # by default only fully connected voxels and no diagonal connections
        # if needed structure must be passed, eg structure = np.ones((3,3,3))
        self.cluster_array_labelled, self.cluster_count = label(self.cluster_array_bool)
        # filter out clusters smaller than threshold
        if self.min_cluster_size is not None:
            # loop over clusters
            for i in range(self.cluster_count):
                n_cluster = i + 1
                cluster_ids = np.where(self.cluster_array_labelled == n_cluster)
                cluster_size = cluster_ids[0].size
                # if cluster below limit set it to np.nan in cluster_array
                if cluster_size < self.min_cluster_size:
                    self.cluster_array_labelled[cluster_ids] = 0
            self.cluster_array_labelled, self.cluster_count = label(self.cluster_array_labelled)

        if self.cluster_count == 0:
            raise NoClustersError('Exiting PearsonMerger: parameters too strict, no clusters detectable!')
        # overwrite zeros with nan for plotting
        zeros_ids = np.where(self.cluster_array_labelled == 0)
        # doesnt work on int array, thus convert to float type
        self.cluster_array_labelled = self.cluster_array_labelled.astype('float')
        self.cluster_array_labelled[zeros_ids] = np.nan
        self.cluster_img = nib.Nifti1Image(self.cluster_array_labelled,
                                           self.affine)

def formatTS(start, stop, r, useIntersects=True):
    aStarts = {}
    X = []
    Y = []
    coords = []
    sigNum = []
    strcArr = [[1, 1, 1], [1, 1, 1], [1, 1, 1]]
    for sig in xrange(start, stop):
        if sig % 50 == 0:
            print sig
        path = np.array(sigs[sig])
        imgO = imread(getFilePath(sig), as_grey=True)
        xs = [n for n in xrange()]
        plt.subplot(1, 2, 1)
        plt.imshow(imgO)
        plt.subplot(1, 2, 2)
        plt.imshow(guassian_filter(imgO, 1))
        plt.show()
        thresh = 0.9
        img = gaussian_filter(imgO, 1) < thresh
        imgX, imgY = np.nonzero(img)
        imgW = imgX.max() - imgX.min()
        imgH = imgY.max() - imgY.min()
        img = skeletonize(img)
        img = img[imgX.min() : imgX.max(), imgY.min() : imgY.max()]
        skltnSegs, nSkltnSegs = label(img, structure=strcArr)
        # print nSkltnSegs
        # plt.imshow(skltnSegs)
        # plt.show()
        imgO = imgO[imgX.min() : imgX.max(), imgY.min() : imgY.max()]
        sumArr = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
        summed = convolve(img, sumArr, mode="constant", cval=0)
        corners = (summed == 2) & (img == 1)
        startx = path[0][0] * img.shape[1]
        starty = path[0][1] * img.shape[0]
        aStarts[sig] = [startx, starty]
        if useIntersects:
            intersects = (summed >= 4) & (img == 1)
            labeled, nLabels = label(intersects, structure=strcArr)
            itrscts = []
            for l in xrange(1, nLabels + 1):
                intersect = np.array((labeled == l), dtype=int)
                posX, posY = np.nonzero(intersect)
                xC, yC = np.array([[np.sum(posX) / posX.size], [np.sum(posY) / posY.size]])
                itrscts.append([xC[0], yC[0]])
            itrscts = np.array(itrscts)
        corners = np.transpose(np.array(np.nonzero(corners)))
        if useIntersects:
            try:
                corners = np.vstack((itrscts, corners))
            except:
                print "something went wrong at", sig
        corners[:, [0, 1]] = corners[:, [1, 0]]
        for i, corner in enumerate(corners):
            x, y = corner[0], corner[1]
            x = getFeatures(imgO, x, y, r, imgH, imgW, skltnSegs)
            X.append(x)
            sigNum.append(sig)
            coords.append([x, y])
    return np.array(X), np.array(sigNum), np.array(coords)

def split_exclusions(image, labels, exclusions, dilation=0, connectivity=1,
    standard_seeds=False):
    """Ensure that no segment in 'labels' overlaps more than one exclusion."""
    labels = labels.copy()
    cur_label = labels.max()
    dilated_exclusions = exclusions.copy()
    foot = generate_binary_structure(exclusions.ndim, connectivity)
    for i in range(dilation):
        dilated_exclusions = grey_dilation(exclusions, footprint=foot)
    hashed = labels * (exclusions.max() + 1) + exclusions
    hashed[exclusions == 0] = 0
    violations = bincount(hashed.ravel()) > 1
    violations[0] = False
    if sum(violations) != 0:
        offending_labels = labels[violations[hashed]]
        mask = zeros(labels.shape, dtype=bool)
        for offlabel in offending_labels:
            mask += labels == offlabel
        if standard_seeds:
            seeds = label(mask * (image == 0))[0]
        else:
            seeds = label(mask * dilated_exclusions)[0]
        seeds[seeds > 0] += cur_label
        labels[mask] = watershed(image, seeds, connectivity, mask)[mask]
    return labels

def segs_to_raveler(sps, bodies, **kwargs):
    import morpho
    min_sp_size = kwargs.get('min_sp_size', 16)
    sps_out = []
    sps_per_plane = []
    sp_to_segment = []
    segment_to_body = [array([[0,0]])]
    total_nsegs = 0
    for i, (sp_map, body_map) in enumerate(zip(sps, bodies)):
        sp_map, nsps = label(
            morpho.remove_small_connected_components(sp_map, min_sp_size, True)
        )
        segment_map, nsegs = label(body_map)
        segment_map += total_nsegs
        segment_map *= sp_map.astype(bool)
        total_nsegs += nsegs
        sps_out.append(sp_map[newaxis,...])
        sps_per_plane.append(nsps)
        valid = (sp_map != 0) + (segment_map == 0)
        sp_to_segment.append(unique(
            zip(it.repeat(i), sp_map[valid], segment_map[valid])))
        valid = segment_map != 0
        logging.debug('plane %i done'%i)
        segment_to_body.append(unique(
                                zip(segment_map[valid], body_map[valid])))
    logging.info('total superpixels before: ' + str(len(unique(sps))) +
                    'total superpixels after: ' + str(sum(sps_per_plane)))
    sps_out = concatenate(sps_out, axis=0)
    sp_to_segment = concatenate(sp_to_segment, axis=0)
    segment_to_body = concatenate(segment_to_body, axis=0)
    return sps_out, sp_to_segment, segment_to_body

def ucm_to_raveler(ucm, sp_threshold=0.0, body_threshold=0.1, **kwargs):
    """Return Raveler map from a UCM.
    
    Parameters
    ----------
    ucm : numpy ndarray, shape (M, N, P)
        An ultrametric contour map. This is a map of scored segment boundaries
        such that if A, B, and C are segments, then 
        score(A, B) = score(B, C) >= score(A, C), for some permutation of
        A, B, and C.
        A hierarchical agglomeration process produces a UCM.
    sp_threshold : float, optional (default: 0.0)
        The value for which to threshold the UCM to obtain the superpixels.
    body_threshold : float, optional (default: 0.1)
        The value for which to threshold the UCM to obtain the segments/bodies.
        The condition `body_threshold >= sp_threshold` should hold in order
        to obtain sensible results.
    **kwargs : dict, optional
        Keyword arguments to be passed through to `segs_to_raveler`.

    Returns
    -------
    superpixels : numpy ndarray, shape (M, N, P)
        The superpixel map. Non-zero superpixels are unique to each plane.
        That is, `np.unique(superpixels[i])` and `np.unique(superpixels[j])` 
        have only 0 as their intersection.
    sp_to_segment : numpy ndarray, shape (Q, 3)
        The superpixel to segment map. Segments are unique to each plane. The
        first number on each line is the plane number.
    segment_to_body : numpy ndarray, shape (R, 2)
        The segment to body map.
    """
    sps = label(ucm < sp_threshold)[0]
    bodies = label(ucm <= body_threshold)[0]
    return segs_to_raveler(sps, bodies, **kwargs)

def getPara(predict, true, threshold, resolution, windowsize):
    (TP, FP, TN, FN, class_lable) = perf_measure(true, predict, threshold)
    if((TP + FN) == 0):
        TPR = 0
    else:
        TPR = np.float(TP) / (TP + FN)

    class_lable = class_lable.astype(bool).reshape(250,  130)
    true = true.astype(bool).reshape((250,  130))

    num = 2
    x = np.arange( -num , num+1, 1)
    xx, yy  = np.meshgrid( x, x )
    struc = (xx * xx + yy * yy)<= num * num
    class_lable = binary_dilation(class_lable, struc)
    class_lable = binary_erosion(class_lable, struc)

    # predict2 = remove_small_objects(class_lable, windowsize * resolution, in_place=False)
    predict2 = remove_small_objects(class_lable, windowsize, in_place=False)
    labeled_array1, num_features1 = label(predict2)
    labeled_array2, num_features2 = label(true)
    FP_num = num_features1 - num_features2
    if FP_num < 0:
        FP_num = 0
    return TPR, FP_num

def __distinct_binary_object_correspondences(reference, result, connectivity=1):
    """
    Determines all distinct (where connectivity is defined by the connectivity parameter
    passed to scipy's `generate_binary_structure`) binary objects in both of the input
    parameters and returns a 1to1 mapping from the labelled objects in reference to the
    corresponding (whereas a one-voxel overlap suffices for correspondence) objects in
    result.

    All stems from the problem, that the relationship is non-surjective many-to-many.

    @return (labelmap1, labelmap2, n_lables1, n_labels2, labelmapping2to1)
    """
    result = np.atleast_1d(result.astype(np.bool))
    reference = np.atleast_1d(reference.astype(np.bool))

    # binary structure
    footprint = generate_binary_structure(result.ndim, connectivity)

    # label distinct binary objects
    labelmap1, n_obj_result = label(result, footprint)
    labelmap2, n_obj_reference = label(reference, footprint)

    # find all overlaps from labelmap2 to labelmap1; collect one-to-one relationships and store all one-two-many for later processing
    slicers = find_objects(labelmap2)  # get windows of labelled objects
    mapping = dict()  # mappings from labels in labelmap2 to corresponding object labels in labelmap1
    used_labels = set()  # set to collect all already used labels from labelmap2
    one_to_many = list()  # list to collect all one-to-many mappings
    for l1id, slicer in enumerate(slicers):  # iterate over object in labelmap2 and their windows
        l1id += 1  # labelled objects have ids sarting from 1
        bobj = (l1id) == labelmap2[slicer]  # find binary object corresponding to the label1 id in the segmentation
        l2ids = np.unique(labelmap1[slicer][
                                 bobj])  # extract all unique object identifiers at the corresponding positions in the reference (i.e. the mapping)
        l2ids = l2ids[0 != l2ids]  # remove background identifiers (=0)
        if 1 == len(
                l2ids):  # one-to-one mapping: if target label not already used, add to final list of object-to-object mappings and mark target label as used
            l2id = l2ids[0]
            if not l2id in used_labels:
                mapping[l1id] = l2id
                used_labels.add(l2id)
        elif 1 < len(l2ids):  # one-to-many mapping: store relationship for later processing
            one_to_many.append((l1id, set(l2ids)))

    # process one-to-many mappings, always choosing the one with the least labelmap2 correspondences first
    while True:
        one_to_many = [(l1id, l2ids - used_labels) for l1id, l2ids in
                       one_to_many]  # remove already used ids from all sets
        one_to_many = [x for x in one_to_many if x[1]]  # remove empty sets
        one_to_many = sorted(one_to_many, key=lambda x: len(x[1]))  # sort by set length
        if 0 == len(one_to_many):
            break
        l2id = one_to_many[0][1].pop()  # select an arbitrary target label id from the shortest set
        mapping[one_to_many[0][0]] = l2id  # add to one-to-one mappings
        used_labels.add(l2id)  # mark target label as used
        one_to_many = one_to_many[1:]  # delete the processed set from all sets

    return labelmap1, labelmap2, n_obj_result, n_obj_reference, mapping

def label(image,**kw):
    """Redefine the scipy.ndimage.measurements.label function to
    work with a wider range of data types.  The default function
    is inconsistent about the data types it accepts on different
    platforms."""
    try: return measurements.label(image,**kw)
    except: pass
    types = ["int32","uint32","int64","unit64","int16","uint16"]
    for t in types:
        try: return measurements.label(array(image,dtype=t),**kw) 
        except: pass
    # let it raise the same exception as before
    return measurements.label(image,**kw)

def translate_back(outputs,threshold=0.7):
    """Translate back. Thresholds on class 0, then assigns
    the maximum class to each region."""
    labels,n = measurements.label(outputs[:,0]<threshold)
    mask = tile(labels.reshape(-1,1),(1,outputs.shape[1]))
    maxima = measurements.maximum_position(outputs,mask,arange(1,amax(mask)+1))
    return [c for (r,c) in maxima]

def manual_split(probs, seg, body, seeds, connectivity=1, boundary_seeds=None):
    """Manually split a body from a segmentation using seeded watershed.

    Input:
        - probs: the probability of boundary in the volume given.
        - seg: the current segmentation.
        - body: the label to be split.
        - seeds: the seeds for the splitting (should be just two labels).
        [-connectivity: the connectivity to use for watershed.]
        [-boundary_seeds: if not None, these locations become inf in probs.]
    Value:
        - the segmentation with the selected body split.
    """
    struct = generate_binary_structure(seg.ndim, connectivity)
    body_pixels = seg == body
    bbox = find_objects(body_pixels)[0]
    body_pixels = body_pixels[bbox]
    body_boundary = binary_dilation(body_pixels, struct) - body_pixels
    non_body_pixels = True - body_pixels - body_boundary
    probs = probs.copy()[bbox]
    probs[non_body_pixels] = probs.min()-1
    if boundary_seeds is not None:
        probs[boundary_seeds[bbox]] = probs.max()+1
    probs[body_boundary] = probs.max()+1
    seeds = label(seeds.astype(bool)[bbox], struct)[0]
    outer_seed = seeds.max()+1 # should be 3
    seeds[non_body_pixels] = outer_seed
    seg_new = watershed(probs, seeds, 
        dams=(seg==0).any(), connectivity=connectivity, show_progress=True)
    seg = seg.copy()
    new_seeds = unique(seeds)[:-1]
    for new_seed, new_label in zip(new_seeds, [0, body, seg.max()+1]):
        seg[bbox][seg_new == new_seed] = new_label
    return seg

 def set_array(self, arr, **kwds):
     kwds['structure'] = kwds.get('structure', np.ones((3,3,3)))
     self._structure = kwds['structure'] # keep record of this
     self._arr = np.asarray(arr)
     labels, num = label(self._arr, **kwds)
     self._labels = labels
     self._nlabels = num

def blob_define(array, thresh, min_area=None, max_area=None, minmax_area=None):
    array[array >= thresh] = 0  # T threshold maskout
    array[np.isnan(array)] = 0  # set ocean nans to 0

    labels, numL = label(array)

    u, inv = np.unique(labels, return_inverse=True)
    n = np.bincount(inv)

    goodinds = u[u!=0]

    if min_area != None:
        goodinds = u[(n>=min_area) & (u!=0)]
        badinds = u[n<min_area]

        for b in badinds:
            pos = np.where(labels==b)
            labels[pos]=0

    if max_area != None:
        goodinds = u[(n<=max_area)  & (u!=0)]
        badinds = u[n>max_area]

    if minmax_area != None:
        goodinds = u[(n <= minmax_area[1]) & (u != 0) & (n>=minmax_area[0])]
        badinds = u[(n > minmax_area[1]) | (n < minmax_area[0])]

        for b in badinds:
            pos = np.where(labels==b)
            labels[pos]=0

    return labels, goodinds

    def snippets_except_labels(self, exclude_labels=None):
        """
        Returns list of snippets which have labels not in the 'exclude_labels'
        list
        TODO - perhaps this should instead return the entire part?
        if exclude_labels is None, then return *all* unlabelled sections
        """

        if isinstance(exclude_labels, basestring):
            exclude_labels = [exclude_labels]

        # create vector which is one whenever one of the exclude labels occurs
        is_label = np.zeros(self.series.shape[0])

        for annotation in self.annotations:
            if annotation['Label'] in exclude_labels:
                start_point = self.time_to_position(annotation['LabelStartTime_Seconds'])
                end_point = self.time_to_position(annotation['LabelEndTime_Seconds'])
                is_label[start_point:end_point] += 1

        # now extract the snippets from which there are no labels
        labs, nlabels = measurements.label(is_label==0)

        snippets = []
        for lab in range(1, nlabels+1):
            snippet = Wave(self.series[labs==lab], self.sample_rate)
            snippets.append(snippet)

        return snippets

def character_seg_erosion(grey_scale_image, max_w_h_ratio=0.85):
    bin_img = grey_scale_image > 0
    labels, num_labels = label(binary_erosion(bin_img > 0))
    for span, mask in _create_spans_and_masks(labels, num_labels):
        char_img = grey_scale_image[:, span[0]:span[1]].copy()
        char_img[mask == False] = 0
        yield char_img

def clean_cc_mask(mask):
    """
    Cleans a segmentation of the corpus callosum so no random pixels are included.

    Parameters
    ----------
    mask : ndarray
        Binary mask of the coarse segmentation.

    Returns
    -------
    new_cc_mask : ndarray
        Binary mask of the cleaned segmentation.
    """

    from scipy.ndimage.measurements import label

    new_cc_mask = np.zeros(mask.shape)

    # Flood fill algorithm to find contiguous regions.
    labels, numL = label(mask)

    volumes = [len(labels[np.where(labels == l_idx+1)]) for l_idx in np.arange(numL)]
    biggest_vol = np.arange(numL)[np.where(volumes == np.max(volumes))] + 1
    new_cc_mask[np.where(labels == biggest_vol)] = 1

    return new_cc_mask

def calc_modes(N2, bottom_depth, z_bins):
    """Wave velocity and structure of first three modes"""

    dz = np.mean(np.diff(z_bins))

    # Truncate N2 to appropriate length based on depth and dz
    Nz = (bottom_depth/dz).astype(int)
    N2 = N2[:Nz]

    # Find indices of start and end of finite values
    finite_vals = nan_or_masked(N2) == 0
    labels = label(finite_vals)[0]
    main_data = np.where(labels == mode(labels[finite_vals]))[1]
    start_ind, end_ind = main_data[0], main_data[-1]

    # Fill in NaN values with start or end values
    N2[:start_ind] = N2[start_ind]
    N2[end_ind + 1:] = N2[end_ind]

    # Preallocate arrays for horizontal and vertical structure
    hori = np.full((len(z_bins) - 1, 3), np.nan)
    vert = hori.copy()

    hori[:len(N2), :], vert[:len(N2), :], c, _ = vertModes(N2, dz, 3)

    return hori, vert, c[:3]

def _filter_small_slopes(hgt, dx, min_slope=0):
    """Masks out slopes with NaN until the slope if all valid points is at 
    least min_slope (in degrees).
    """

    min_slope = np.deg2rad(min_slope)
    slope = np.arctan(-np.gradient(hgt, dx))  # beware the minus sign
    # slope at the end always OK
    slope[-1] = min_slope

    # Find the locs where it doesn't work and expand till we got everything
    slope_mask = np.where(slope >= min_slope, slope, np.NaN)
    r, nr = label(~np.isfinite(slope_mask))
    for objs in find_objects(r):
        obj = objs[0]
        i = 0
        while True:
            i += 1
            i0 = objs[0].start-i
            if i0 < 0:
                break
            ngap =  obj.stop - i0 - 1
            nhgt = hgt[[i0, obj.stop]]
            current_slope = np.arctan(-np.gradient(nhgt, ngap * dx))
            if i0 <= 0 or current_slope[0] >= min_slope:
                break
        slope_mask[i0:obj.stop] = np.NaN
    out = hgt.copy()
    out[~np.isfinite(slope_mask)] = np.NaN
    return out

def fix_elevations(z, riv_i, riv_j, ch_depth, sea_level, slope, dx, max_rand, SLRR):

    test_elev = z - sea_level
    max_cell_h = slope * dx
    riv_prof = test_elev[riv_i, riv_j]
    test_elev[riv_i, riv_j] += 2*ch_depth

    # set new subaerial cells to marsh elevation
    test_elev[test_elev == 0] = max_cell_h

    # make mask for depressions
    ocean_mask = test_elev < max_cell_h
    labeled_ponds, ocean = measurements.label(ocean_mask)

    # # fill in underwater spots that are below SL (?)
    # below_SL = [z <= sea_level]
    # underwater_cells, big_ocean = measurements.label(below_SL)
    # underwater_cells[underwater_cells == big_ocean] = 0
    # test_elev[underwater_cells > 0] = max_cell_h + SLRR + (np.random.rand() * max_rand)

    # create an ocean and shoreline mask
    ocean_and_shore = np.copy(labeled_ponds)

    # create an ocean mask
    # ocean_cells = np.copy(ocean_and_shore)
    # ocean_and_shore[test_elev > 0] = 0

    # create mask for pond cells and fix them
    area = measurements.sum(ocean_mask, labeled_ponds, index=np.arange(labeled_ponds.max() + 1))
    areaPonds = area[labeled_ponds]
    labeled_ponds[areaPonds == areaPonds.max()] = 0

    #finish creating ocean and shoreline mask
    ocean_and_shore[areaPonds != areaPonds.max()] = 0

    # something here to get rid of ocean cells
    test_elev[labeled_ponds > 0] = max_cell_h + SLRR + (np.random.rand() * max_rand)

    # raise cells close to sea level above it
    test_elev[(test_elev >= max_cell_h) & (test_elev <= (max_cell_h + SLRR))] = \
        (max_cell_h + SLRR + (np.random.rand() * max_rand))

    riv_buffer = np.zeros_like(test_elev)
    riv_buffer[riv_i, riv_j] = 1
    riv_buffer[riv_i[1:]-1, riv_j[1:]] = 1

    for i in xrange(1, test_elev.shape[0]-1):
        for j in xrange(test_elev.shape[1]):
            if (not ocean_and_shore[i, j]
                and not ocean_and_shore[i-1, j]
                and not ocean_and_shore[i+1, j]
                and not riv_buffer[i, j]):
                if test_elev[i+1, j] >= test_elev[i, j]:
                    test_elev[i, j] = test_elev[i+1, j] + (np.random.rand() * slope)
    
    test_elev[riv_i, riv_j] = riv_prof

    z = test_elev + sea_level

    return z

def mod_regions(composite, mod_id, algorithm):
  # paths
  template = composite.templates.get(name='region') # REGION TEMPLATE

  # get region img set that has the region template
  region_img_set = composite.gons.filter(channel__name='-regionimg', template__name='region')

  # channel
  region_channel, region_channel_created = composite.channels.get_or_create(name='-regions')

  # iterate
  for t in range(composite.series.ts):
    print(t)
    region_img = region_img_set.filter(t=t)
    if region_img.count()==0:
      region_img = region_img_set.get(t=t-1)
    else:
      region_img = region_img_set.get(t=t)

    # for each image, determine unique values of labelled array
    # make gon with label array and save

    region_gon = composite.gons.create(experiment=composite.experiment, series=composite.series, channel=region_channel, template=template)
    region_gon.set_origin(0, 0, 0, t)
    region_gon.set_extent(composite.series.rs, composite.series.cs, 1)

    # modify image
    region_array = region_img.load()
    region_array = region_array[:,:,0]
    region_array[region_array>0] = 1
    region_array, n = label(region_array)

    region_gon.array = region_array.copy()
    region_gon.save_array(composite.experiment.composite_path, template)
    region_gon.save()

def rand_by_threshold(ucm, gt, npoints=None):
    """Compute Rand and Adjusted Rand indices for each threshold of a UCM

    Parameters
    ----------
    ucm : np.ndarray, arbitrary shape
        An Ultrametric Contour Map of region boundaries having specific
        values. Higher values indicate higher boundary probabilities.
    gt : np.ndarray, int type, same shape as ucm
        The ground truth segmentation.
    npoints : int, optional
        If provided, only compute values at npoints thresholds, rather than
        all thresholds. Useful when ucm has an extremely large number of
        unique values.

    Returns
    -------
    ris : np.ndarray of float, shape (3, len(np.unique(ucm))) or (3, npoints)
        The rand indices of the segmentation induced by thresholding and
        labeling `ucm` at different values. The 3 rows of `ris` are the values
        used for thresholding, the corresponding Rand Index at that threshold,
        and the corresponding Adjusted Rand Index at that threshold.
    """
    ts = np.unique(ucm)[1:]
    if npoints is None:
        npoints = len(ts)
    if len(ts) > 2 * npoints:
        ts = ts[np.arange(1, len(ts), len(ts) / npoints)]
    result = np.zeros((2, len(ts)))
    for i, t in enumerate(ts):
        seg = label(ucm < t)[0]
        result[0, i] = rand_index(seg, gt)
        result[1, i] = adj_rand_index(seg, gt)
    return np.concatenate((ts[np.newaxis, :], result), axis=0)