def adjust_spot_positions(image, label_image, hp, debug=None):
    """Re-evaluate the spot positions based on the segmentation. 
        Parameters: 
        image: The original image (can be masked) that was sent to findspot3d
        label_image: the label image containing two labels 
        hp: the original hotpoints
        debug: set to true to write out an image debugimg.nii.gz with the stuff
        """
        
    struct2 = generate_binary_structure(3, 2)
    struct1 = generate_binary_structure(3, 1)
    peak_points =[] 

    if debug is None:
        temp_path = os.getenv("PYSBR_TEMP")
        if temp_path is not None:
            debug = os.path.join(temp_path, "debug-labels.nii.gz")
    
    if debug is not None:
            debimg = image.copy()

    nlabels = label_image.max()

    if nlabels!=len(hp):
        raise RuntimeError( 'number of labels and hotspots should be the same' )

    tins = []
    for n in range(nlabels):
        label = n+1
        area = binary_closing(label_image == label, struct2)
        thiniter = np.sum(area.reshape(-1)) / 1500 + 1
        csbr.thinning3d(area, thiniter)
        tins.append(area)
   
    for n in range(nlabels):
        label = n+1
        
        #avoid that a single pixel breaks the evaluation by running a closing 
        area = label_image == label
        
        #evaluate the boundary 
        dmask = binary_dilation(area, struct1)
        border = np.bitwise_xor(dmask, area)
        
        p = adjust_spot_position(image, border, image[tuple(hp[n])], tins[n], tins[(n + 1) % 2])
        peak_points.append(p)

        if debug is not None:
            debimg[border>0] = 196
            debimg[p] = 0
            nib.save(nib.Nifti1Image(debimg, global_affine), debug)

    peak_points = np.array( peak_points )
    return peak_points

def getBoundariesOfimage(image):
    """
    find edges by using erosion
    """
    if np.ndim(image) == 2:
        sElement = ndimage.generate_binary_structure(2, 1)
    else:
        sElement = ndimage.generate_binary_structure(3, 1)
    erode_im = scipy.ndimage.morphology.binary_erosion(image, sElement)
    b = image - erode_im
    return b

 def InDecPatch(self,which,amount):
     s = ndimage.generate_binary_structure(2,1) # taxi-cab struct
     if which == 0:
         ras = ndimage.binary_dilation(self.cl_array,s,iterations=amount,border_value=0)
     else:
         ras = ndimage.binary_erosion(self.cl_array,s,iterations=amount,border_value=0)
     return(ras)

def artifact_mask(imdata, airdata, distance, zscore=10.):
    """Computes a mask of artifacts found in the air region"""
    from statsmodels.robust.scale import mad

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    if np.sum((bg_img > 0).astype(np.uint8)) < 100:
        return np.zeros_like(airdata)

    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    bg_location = np.median(bg_img[bg_img > 0])
    bg_spread = mad(bg_img[bg_img > 0])
    bg_img[bg_img > 0] -= bg_location
    bg_img[bg_img > 0] /= bg_spread

    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > zscore] = 1
    qi1_img[distance < .10] = 0

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 1)
    qi1_img = nd.binary_opening(qi1_img, struc).astype(np.uint8)
    qi1_img[airdata <= 0] = 0

    return qi1_img

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 __init__(self, label_image=None, connectivity=1, data=None, **attr):

        super(RAG, self).__init__(data, **attr)
        if self.number_of_nodes() == 0:
            self.max_id = 0
        else:
            self.max_id = max(self.nodes_iter())

        if label_image is not None:
            fp = ndi.generate_binary_structure(label_image.ndim, connectivity)
            # In the next ``ndi.generic_filter`` function, the kwarg
            # ``output`` is used to provide a strided array with a single
            # 64-bit floating point number, to which the function repeatedly
            # writes. This is done because even if we don't care about the
            # output, without this, a float array of the same shape as the
            # input image will be created and that could be expensive in
            # memory consumption.
            ndi.generic_filter(
                label_image,
                function=_add_edge_filter,
                footprint=fp,
                mode='nearest',
                output=as_strided(np.empty((1,), dtype=np.float_),
                                  shape=label_image.shape,
                                  strides=((0,) * label_image.ndim)),
                extra_arguments=(self,))

 def f_returnInternalEdge(self,cl_array):
     # Internal edge: Count of neighboring non-zero cell       
     kernel = ndimage.generate_binary_structure(2, 1) # Make a kernel
     kernel[1, 1] = 0
     b = ndimage.convolve(cl_array, kernel, mode="constant")
     n_interior = b[cl_array != 0].sum() # Number of interiror edges
     return n_interior

def edge_matrix(labels, connectivity=1):
    """Generate a COO matrix containing the coordinates of edge pixels.

    Parameters
    ----------
    labels : array of int
        An array of labeled pixels (or voxels).
    connectivity : int in {1, ..., labels.ndim}
        The square connectivity for considering neighborhood.

    Returns
    -------
    edges : sparse.coo_matrix
        A COO matrix where (i, j) indicate neighboring labels and the
        corresponding data element is the linear index of the edge pixel
        in the labels array.
    """
    conn = ndi.generate_binary_structure(labels.ndim, connectivity)
    eroded = ndi.grey_erosion(labels, footprint=conn).ravel()
    dilated = ndi.grey_dilation(labels, footprint=conn).ravel()
    labels = labels.ravel()
    boundaries0 = np.flatnonzero(eroded != labels)
    boundaries1 = np.flatnonzero(dilated != labels)
    labels_small = np.concatenate((eroded[boundaries0], labels[boundaries1]))
    labels_large = np.concatenate((labels[boundaries0], dilated[boundaries1]))
    n = np.max(labels_large) + 1
    data = np.concatenate((boundaries0, boundaries1))
    sparse_graph = sparse.coo_matrix((data, (labels_small, labels_large)),
                                     dtype=np.int_, shape=(n, n))
    return sparse_graph

def objextract(Fg):

    s = nd.generate_binary_structure(2,2)
    labeled_array, num_features = nd.measurements.label(Fg, structure=s)
    coor = []
    cnt = []

    if num_features == 0:
       idx = []
    else:    
        lth = 200   # label pixel number less than lth will be removed
        Lth = 6500
        for i in range(1,num_features+1):
            coor.append(np.where(labeled_array==i))
            cnt.append(len(np.where(labeled_array==i)[1]))

        cnt = array(cnt)
        idx = arange(num_features)
        idx = idx[(cnt<Lth)&(cnt>lth)]
      
        if len(idx)==0:
            idx = []
        elif len(idx)>1:              
            #idx = [idx[cnt[idx].argmax()]]
            idx = sorted(range(len(cnt)),key=lambda x:cnt[x])[::-1][0:2]

    return idx,labeled_array,coor,cnt

def find_local_max(img, d_rad, threshold=1e-15):
    """
    This is effectively a replacement for pkfnd in the matlab/IDL code.

    The output of this function is meant to be feed into :py:func:`~subpixel_centroid`

    The magic of numpy means this should work for any dimension data.

    :param img: an ndarray representing the data to find the local maxes
    :param d_rad: the radius of the dilation, the smallest possible spacing between local maximum
    :param threshold: optional, voxels < threshold are ignored.

    :rtype: (d,N) array of the local maximums.
    """
    d_rad = int(d_rad)
    img = np.array(np.squeeze(img))       # knock out singleton dimensions
    img[img < threshold] = -np.inf        # mask out pixels below threshold
    dim = img.ndim                        # get the dimension of data

    # make structuring element
    s = ndimage.generate_binary_structure(dim, 1)
    # scale it up to the desired size
    d_struct = ndimage.iterate_structure(s, int(d_rad))
    dilated_img = ndimage.grey_dilation(img,
                                        footprint=d_struct,
                                        cval=0,
                                        mode='constant')   # do the dilation

    # find the locations that are the local maximum
    # TODO clean this up
    local_max = np.where(np.exp(img - dilated_img) > (1 - 1e-15))
    # the extra [::-1] is because matplotlib and ndimage disagree an xy vs yx
    return np.vstack(local_max[::-1])

def thresholding(img, thresh, size=9):
    """
    Segment using a thresholding algorithm
    
    Input:
     - img  ndarray : Image array (ndim=2)
     - thresh float : Threshold value for pixels selectino
     - size     int : Minimum size a group of pixels must have
    
    Output:
     - regions : Binary array for each segmented region
    
    ---
    """

    logging.debug("Threshold: %.2f", thresh)
    logging.debug("Objects min size: %d", size)

    # Take the binary image thresholded
    img_bin = img > thresh

    # And use (MO) binary opening (erosion + dilation) for cleaning spurious Trues
    strct = ndi.generate_binary_structure(2, 2)
    img_bin = ndi.binary_opening(img_bin, strct)

    # Label each group/region (value==True) of pixels
    regions, nlbl = ndi.label(img_bin)
    for i in xrange(1, nlbl + 1):
        inds = np.where(regions == i)
        if inds[0].size < size:
            regions[inds] = 0

    logging.debug("Threshold labels: %s", np.unique(regions))

    return regions.astype(np.bool)

def remove_small_objects(ar, min_size=64, connectivity=1, in_place=False):
    # Should use `issubdtype` for bool below, but there's a bug in numpy 1.7
    if not (ar.dtype == bool or np.issubdtype(ar.dtype, np.integer)):
        raise TypeError("Only bool or integer image types are supported. "
                        "Got %s." % ar.dtype)

    if in_place:
        out = ar
    else:
        out = ar.copy()

    if min_size == 0:  # shortcut for efficiency
        return out

    if out.dtype == bool:
        selem = nd.generate_binary_structure(ar.ndim, connectivity)
        ccs = np.zeros_like(ar, dtype=np.int32)
        nd.label(ar, selem, output=ccs)
    else:
        ccs = out

    try:
        component_sizes = np.bincount(ccs.ravel())
    except ValueError:
        raise ValueError("Negative value labels are not supported. Try "
                         "relabeling the input with `scipy.ndimage.label` or "
                         "`skimage.morphology.label`.")

    too_small = component_sizes < min_size
    too_small_mask = too_small[ccs]
    out[too_small_mask] = 0

    return out

def artifact_mask(imdata, airdata, distance):
    """Computes a mask of artifacts found in the air region"""
    import nibabel as nb

    if not np.issubdtype(airdata.dtype, np.integer):
        airdata[airdata < .95] = 0
        airdata[airdata > 0.] = 1

    bg_img = imdata * airdata
    # Find the background threshold (the most frequently occurring value
    # excluding 0)
    # CHANGED - to the 75 percentile
    bg_threshold = np.percentile(bg_img[airdata > 0], 75)

    # Apply this threshold to the background voxels to identify voxels
    # contributing artifacts.
    qi1_img = np.zeros_like(bg_img)
    qi1_img[bg_img > bg_threshold] = 1
    qi1_img[distance < .10] = 0

    # Create a structural element to be used in an opening operation.
    struc = nd.generate_binary_structure(3, 1)
    qi1_img = nd.binary_opening(qi1_img, struc).astype(np.uint8)
    qi1_img[airdata <= 0] = 0

    return qi1_img

def compute_sparsity(im):
    l_x = len(im)
    X, Y = np.ogrid[:l_x, :l_x]
    mask = ((X - l_x/2)**2 + (Y - l_x/2)**2 <= (l_x/2)**2)
    grad1 = ndimage.morphological_gradient(im, footprint=np.ones((3, 3)))
    grad2 = ndimage.morphological_gradient(im, footprint=ndimage.generate_binary_structure(2, 1))
    return (grad1[mask] > 0).mean(), (grad2[mask] > 0).mean() 

def ClusterizeImage(image,thresh=None,connectivity=3):
    if thresh is None:
        thresh = 0
    image[np.where(image<=thresh)] = 0
    s = generate_binary_structure(3,connectivity)
    larray, nf = label(image,s)
    return larray

def test_labeling():
    "Test cluster labeling"
    shape = flat_shape = (4, 20)
    pmap = np.empty(shape, np.float_)
    out = np.empty(shape, np.uint32)
    bin_buff = np.empty(shape, np.bool_)
    int_buff = np.empty(shape, np.uint32)
    struct = ndimage.generate_binary_structure(2, 1)
    struct[::2] = False
    conn = np.array([(0, 1), (0, 3), (1, 2), (2, 3)], np.uint32)
    criteria = None

    # some clusters
    pmap[:] = [[ 3, 3, 0, 0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 0],
               [ 0, 1, 0, 0, 0, 0, 8, 0, 0, 4, 4, 4, 0, 0, 0, 0, 0, 0, 4, 0],
               [ 0, 3, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 4, 4],
               [ 0, 0, 0, 0, 0, 0, 8, 0, 0, 0, 0, 0, 4, 4, 4, 0, 0, 0, 0, 0]]
    cids = _label_clusters(pmap, out, bin_buff, int_buff, 2, 0, struct, False,
                           flat_shape, conn, criteria)
    assert_equal(len(cids), 6)

    # some other clusters
    pmap[:] = [[ 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0],
               [ 0, 4, 0, 0, 0, 0, 0, 4, 0, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0],
               [ 0, 0, 4, 4, 0, 4, 4, 0, 4, 0, 0, 0, 4, 4, 1, 0, 4, 4, 0, 0],
               [ 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 0]]
    cids = _label_clusters(pmap, out, bin_buff, int_buff, 2, 0, struct, False,
                           flat_shape, conn, criteria)
    assert_equal(len(cids), 6)

def gen_data(xsize, ysize, nstars=3, starradius=10, brightness=2000):
    # 1) lots of stars, big
    # 2) lots of tiny stars
    # 3) few stars, big
    # 4) few stars, tiny

    footprint = ndimage.generate_binary_structure(2,1)

    ret = numpy.zeros((xsize, ysize))
    for star in xrange(nstars):
        xcenter = random.randint(0, xsize-1)
        ycenter = random.randint(0, ysize-1)
        for x in xrange(xcenter-1, xcenter+2):
            for y in xrange(ycenter-1, ycenter+2):
                if x >= 0 and y >= 0 and x < xsize and y < ysize:
                    ret[x,y] = brightness / 3
        ret[xcenter, ycenter] = brightness
    for i in xrange(starradius):
        ret = ndimage.grey_dilation(ret, footprint=footprint)

    # add some cosmic rays (single points)
    for i in xrange(30):
        xcenter = random.randint(0, xsize-1)
        ycenter = random.randint(0, ysize-1)
        ret[xcenter, ycenter] = brightness

    return ret

def get_largest_two_component(img, prt = False, threshold = None):
    s = ndimage.generate_binary_structure(3,2) # iterate structure
    labeled_array, numpatches = ndimage.label(img,s) # labeling
    sizes = ndimage.sum(img,labeled_array,range(1,numpatches+1))
    sizes_list = [sizes[i] for i in range(len(sizes))]
    sizes_list.sort()
    #if(prt):
    #    print('component size', sizes_list, flush = True)
    if(len(sizes) == 1):
        return img
    else:
        if(threshold):
            out_img = np.zeros_like(img)
            for temp_size in sizes_list:
                if(temp_size > threshold):
                    temp_lab = np.where(sizes == temp_size)[0] + 1
                    temp_cmp = labeled_array == temp_lab
                    out_img = (out_img + temp_cmp) > 0
            return out_img
        else:
            max_size1 = sizes_list[-1]
            max_size2 = sizes_list[-2]
            max_label1 = np.where(sizes == max_size1)[0] + 1
            max_label2 = np.where(sizes == max_size2)[0] + 1
            component1 = labeled_array == max_label1
            component2 = labeled_array == max_label2
            #if(prt):
            #    print(max_size2, max_size1, max_size2/max_size1, flush = True)
            if(max_size2*10 > max_size1):
                component1 = (component1 + component2) > 0

            return component1

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 art_qi2(img, airmask, ncoils=12, erodemask=True):
    """
    Calculates **qi2**, the distance between the distribution
    of noise voxel (non-artifact background voxels) intensities, and a
    centered Chi distribution.

    :param numpy.ndarray img: input data
    :param numpy.ndarray airmask: input air mask without artifacts

    """
    from matplotlib import rc
    import seaborn as sn
    import matplotlib.pyplot as plt
    rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
    # rc('text', usetex=True)

    if erodemask:
        struc = nd.generate_binary_structure(3, 2)
        # Perform an opening operation on the background data.
        airmask = nd.binary_erosion(airmask, structure=struc).astype(np.uint8)

    # Artifact-free air region
    data = img[airmask > 0]
    data = data[data < np.percentile(data, 99.5)]
    maxvalue = int(data.max())
    nbins = maxvalue if maxvalue < 100 else 100

    # Estimate data pdf
    hist, bin_edges = np.histogram(data, density=True, bins=nbins)
    bin_centers = [np.mean(bin_edges[i:i+1]) for i in range(len(bin_edges)-1)]
    max_pos = np.argmax(hist)

    # Fit central chi distribution
    param = chi.fit(data, 2*ncoils, loc=bin_centers[max_pos])
    pdf_fitted = chi.pdf(bin_centers, *param[:-2], loc=param[-2], scale=param[-1])

    # Write out figure of the fitting
    out_file = op.abspath('background_fit.png')
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    sn.distplot(data, bins=nbins, norm_hist=True, kde=False, ax=ax1)
    #_, bins, _ = ax1.hist(data, nbins, normed=True, color='gray', linewidth=0)
    ax1.plot(bin_centers, pdf_fitted, 'k--', linewidth=1.2)
    fig.suptitle('Noise distribution on the air mask, and fitted chi distribution')
    ax1.set_xlabel('Intensity')
    ax1.set_ylabel('Frequency')
    fig.savefig(out_file, format='png', dpi=300)
    plt.close()

    # Find t2 (intensity at half width, right side)
    ihw = 0.5 * hist[max_pos]
    t2idx = 0
    for i in range(max_pos + 1, len(bin_centers)):
        if hist[i] < ihw:
            t2idx = i
            break

    # Compute goodness-of-fit (gof)
    return (float(np.abs(hist[t2idx:] - pdf_fitted[t2idx:]).sum() /
                  len(pdf_fitted[t2idx:])), out_file)