Python scipy.ndimage 模块，label() 实例源码

def loop_connected_components(mask):
    """
    @type mask: np.ndarray
    @rtype (np.ndarray, np.ndarray, np.ndarray)
    """

    c = np.array([])
    init = np.array([])
    fin = np.array([])

    if mask.sum() > 0:
        labeled = sp_image.label(mask)[0]
        components = sp_image.measurements.find_objects(labeled)
        c_fin = [(s[0].stop - s[0].start, s[0].stop - 1) for s in components]
        if len(c_fin) > 1 and mask[0] and mask[-1]:
            c_fin[0] = c_fin[0][0] + c_fin[-1][0], c_fin[0][1]
            c_fin = c_fin[0:-1]

        c, fin = zip(*c_fin)
        c = np.array(c, dtype=int)
        fin = np.array(fin, dtype=int)
        init = (fin - c + 1) % mask.shape[0]
    return c, init, fin

def remove_artifacts(self, image):
        """
        Remove the connected components that are not within the parameters
        Operates in place
        :param image: sudoku's thresholded image w/o grid
        :return: None
        """
        labeled, features = label(image, structure=CROSS)
        lbls = np.arange(1, features + 1)
        areas = extract_feature(image, labeled, lbls, np.sum,
                                np.uint32, 0)
        sides = extract_feature(image, labeled, lbls, min_side,
                                np.float32, 0, True)
        diags = extract_feature(image, labeled, lbls, diagonal,
                                np.float32, 0, True)

        for index in lbls:
            area = areas[index - 1] / 255
            side = sides[index - 1]
            diag = diags[index - 1]
            if side < 5 or side > 20 \
                    or diag < 15 or diag > 25 \
                    or area < 40:
                image[labeled == index] = 0
        return None

def remove_artifacts(self, image):
        """
        Remove the connected components that are not within the parameters
        Operates in place
        :param image: sudoku's thresholded image w/o grid
        :return: None
        """
        labeled, features = label(image, structure=CROSS)
        lbls = np.arange(1, features + 1)
        areas = extract_feature(image, labeled, lbls, np.sum,
                                np.uint32, 0)
        sides = extract_feature(image, labeled, lbls, min_side,
                                np.float32, 0, True)
        diags = extract_feature(image, labeled, lbls, diagonal,
                                np.float32, 0, True)

        for index in lbls:
            area = areas[index - 1] / 255
            side = sides[index - 1]
            diag = diags[index - 1]
            if side < 5 or side > 20 \
                    or diag < 15 or diag > 25 \
                    or area < 40:
                image[labeled == index] = 0
        return None

def adjust_prediction(self, prediction):

        new_prediction = prediction

        labeled, n_objects = ndimage.label(prediction > 0)

        max_volume = 0

        volumes = {}

        for object_n in range(1, n_objects + 1):
            volume = np.sum(prediction[labeled == object_n])
            if volume > max_volume:
                max_volume = volume
            volumes.update({object_n: volume})

        for object_n, volume in volumes.iteritems():
            if volume < self.threshold * max_volume:
                new_prediction[labeled == object_n] = 0
        return new_prediction

def write_file(filename, resolution, **kwargs):
        h5file = h5py.File(filename, 'w')
        config = {'hdf5_file': os.path.basename(filename)}
        channels = ['image', 'label', 'mask']
        default_datasets = {
            'image': 'volumes/raw',
            'label': 'volumes/labels/neuron_ids',
            'mask': 'volumes/labels/mask',
        }
        for channel in channels:
            data = kwargs.get('{}_data'.format(channel), None)
            dataset_name = kwargs.get('{}_dataset'.format(channel), default_datasets[channel])
            if data is not None:
                dataset = h5file.create_dataset(dataset_name, data=data, dtype=data.dtype)
                dataset.attrs['resolution'] = resolution
                config['{}_dataset'.format(channel)] = dataset_name

        h5file.close()

        return config

def find_windows_from_heatmap(image):
    hot_windows = []
    # Threshold the heatmap
    thres = 0
    image[image <= thres] = 0
    # Set labels
    labels = ndi.label(image)
    # iterate through labels and find windows
    for car_number in range(1, labels[1]+1):
        # Find pixels with each car_number label value
        nonzero = (labels[0] == car_number).nonzero()
        # Identify x and y values of those pixels
        nonzeroy = np.array(nonzero[0])
        nonzerox = np.array(nonzero[1])
        # Define a bounding box based on min/max x and y
        bbox = ((np.min(nonzerox), np.min(nonzeroy)), (np.max(nonzerox), np.max(nonzeroy)))
        hot_windows.append(bbox)
    return hot_windows

def calc_detection_rate(prediction,label):
    '''
    Calculate the detection True Positives, False Negatives and False Positives.
    TP occurs whenever the prediction region intersects the ground truth region.
    FP occurs 
    '''
    label_temp=np.copy(label)   
    TP=FN=FP=0.0
    # pattern for ndimage neighbouring pixels
    s = [[1,1,1],
        [1,1,1],
        [1,1,1]]
    labeled_prediction, num_features_prediction = ndimage.label(prediction, structure=s)
    labeled_label, num_features_label = ndimage.label(label_temp, structure=s)
    for i in range(1,num_features_prediction+1):
        intersection=np.sum(labeled_prediction[label==1]==i)
        if intersection>0:
            TP+=1
        else:
            FP+=1
    for i in range(1,num_features_label+1):
        intersection=np.sum(labeled_label[prediction==1]==i)
        if intersection==0:
            FN=+1
    return TP,FN,FP

def calc_detection_rate(prediction,label):
    '''
    Calculate the detection True Positives, False Negatives and False Positives.
    TP occurs whenever the prediction region intersects the ground truth region.
    FP occurs 
    '''
    label_temp=np.copy(label)   
    TP=FN=FP=0.0
    # pattern for ndimage neighbouring pixels
    s = [[1,1,1],
        [1,1,1],
        [1,1,1]]
    labeled_prediction, num_features_prediction = ndimage.label(prediction, structure=s)
    labeled_label, num_features_label = ndimage.label(label_temp, structure=s)
    for i in range(1,num_features_prediction+1):
        intersection=np.sum(labeled_prediction[label==1]==i)
        if intersection>0:
            TP+=1
        else:
            FP+=1
    for i in range(1,num_features_label+1):
        intersection=np.sum(labeled_label[prediction==1]==i)
        if intersection==0:
            FN=+1
    return TP,FN,FP

def calc_detection_rate(prediction,label):
    '''
    Calculate the detection True Positives, False Negatives and False Positives.
    TP occurs whenever the prediction region intersects the ground truth region.
    FP occurs 
    '''
    label_temp=np.copy(label)   
    TP=FN=FP=0.0
    # pattern for ndimage neighbouring pixels
    s = [[1,1,1],
        [1,1,1],
        [1,1,1]]
    labeled_prediction, num_features_prediction = ndimage.label(prediction, structure=s)
    labeled_label, num_features_label = ndimage.label(label_temp, structure=s)
    for i in range(1,num_features_prediction+1):
        intersection=np.sum(labeled_prediction[label==1]==i)
        if intersection>0:
            TP+=1
        else:
            FP+=1
    for i in range(1,num_features_label+1):
        intersection=np.sum(labeled_label[prediction==1]==i)
        if intersection==0:
            FN=+1
    return TP,FN,FP

def calc_detection_rate(prediction,label):
    '''
    Calculate the detection True Positives, False Negatives and False Positives.
    TP occurs whenever the prediction region intersects the ground truth region.
    FP occurs 
    '''
    label_temp=np.copy(label)   
    TP=FN=FP=0.0
    # pattern for ndimage neighbouring pixels
    s = [[1,1,1],
        [1,1,1],
        [1,1,1]]
    labeled_prediction, num_features_prediction = ndimage.label(prediction, structure=s)
    labeled_label, num_features_label = ndimage.label(label_temp, structure=s)
    for i in range(1,num_features_prediction+1):
        intersection=np.sum(labeled_prediction[label==1]==i)
        if intersection>0:
            TP+=1
        else:
            FP+=1
    for i in range(1,num_features_label+1):
        intersection=np.sum(labeled_label[prediction==1]==i)
        if intersection==0:
            FN=+1
    return TP,FN,FP

def calc_detection_rate(prediction,label):
    '''
    Calculate the detection True Positives, False Negatives and False Positives.
    TP occurs whenever the prediction region intersects the ground truth region.
    FP occurs 
    '''
    label_temp=np.copy(label)   
    TP=FN=FP=0.0
    # pattern for ndimage neighbouring pixels
    s = [[1,1,1],
        [1,1,1],
        [1,1,1]]
    labeled_prediction, num_features_prediction = ndimage.label(prediction, structure=s)
    labeled_label, num_features_label = ndimage.label(label_temp, structure=s)
    for i in range(1,num_features_prediction+1):
        intersection=np.sum(labeled_prediction[label==1]==i)
        if intersection>0:
            TP+=1
        else:
            FP+=1
    for i in range(1,num_features_label+1):
        intersection=np.sum(labeled_label[prediction==1]==i)
        if intersection==0:
            FN=+1
    return TP,FN,FP

def get_centers_of_mass_from_blobs(segmentation_layer, iterations=3):
    """
    Determine the centers of each object in an image

    ::param segmentation_layer: NxM ndarray image mask of all target objects
    ::param iterations: threshold for removal of small non-target objects
    ::return centers_of_mass: a np ndarray of x,y coordinates for the center of each target object

    """

    segmentation_layer = ndimage.binary_opening(segmentation_layer, iterations=iterations)  # remove small objects
    labels, label_number = ndimage.measurements.label(segmentation_layer)  # label remaining blobs

    centers_of_mass = np.zeros((label_number, 2))

    for i in range(label_number):
        idx = np.where(labels == i+1)  # calculate the center of mass for each blob
        centers_of_mass[i, 1] = np.mean(idx[1].astype(float))  # must be float
        centers_of_mass[i, 0] = np.mean(idx[0].astype(float))  # must be float

    return centers_of_mass

def remove_small_blobs(centers_of_mass, segmentation_layer):
    """
    removes non-overlapping pixel-islands and cell centres (centers_of_mass)
    :param segmentation_layer: NxM ndarray image mask of all target objects
    :param centers_of_mass: a np ndarray of x,y coordinates for the center of each target object
    :return updated_labels:
    """
    labels, label_number = ndimage.label(segmentation_layer)  # label all pixel islands
    updated_labels = np.zeros_like(labels)
    centers_of_mass_to_keep = np.zeros(len(centers_of_mass))  # holder

    for i in range(label_number):
        idx = np.where(labels == i+1)
        for j, c in enumerate(centers_of_mass.astype(int)):
            if labels[c[0], c[1]] == i+1:  # if the center_of_mass is in the blob
                updated_labels[idx] = 1  # add the blob
                centers_of_mass_to_keep[j] = 1  # add the center_of_mass

    centers_of_mass_idx = np.where(centers_of_mass_to_keep == 1)
    updated_centers_of_mass = centers_of_mass[centers_of_mass_idx]

    return updated_labels, updated_centers_of_mass

def calculate_distance(centers_of_mass, image):
    """
    takes the centers of each blob, and an image to be segmented. Divides the image according to the center of masses
    by a random walk

    :param image: a binarised image to be segmented
    :param centers_of_mass: the centres that will define the maxima of the watershed segmentation
    :return segmentation_labels: a labelled image/segmentation, where each index belongs do a different center of mass

    """
    # random walk segmentation of 2D image-mask array
    distance = ndimage.distance_transform_edt(np.abs(image-1))
    local_maxi = np.zeros_like(image)

    for c in centers_of_mass:
        local_maxi[int(c[0]), int(c[1])] = 1

    markers = ndimage.label(local_maxi)[0]
    segmentation_labels = segmentation.random_walker(distance, markers, beta=60)

    return segmentation_labels

def blob_labels(centers_of_mass, blob_image):
    """
    label nuclei with segmentation - so labels are in the same order as the outer layer

    :param list centers_of_mass: centers of target blobs/cells
    :param np.array blob_image: image of the target cells, or nuclei of cells
    :return segmented_blobs: a labelled image where each index is a different cell
    :return distance: image where each pixel's value is related to how far away from a blob center it is
    """

    image = np.abs(blob_image-1)
    distance = ndimage.distance_transform_edt(np.abs(image-1))
    local_maxi = np.zeros_like(image)
    for c in centers_of_mass:
        local_maxi[int(c[0]), int(c[1])] = 1
    markers = ndimage.label(local_maxi)[0]
    segmented_blobs = segmentation.random_walker(distance, markers, beta=20)

    return segmented_blobs, distance

def getAlignImg(t,label = None):#!!!notice, only take uint8 type for the imrotate function!!!
    f = lambda x:np.asarray([float(a) for a in x]);
    o = f(t.ImageOrientationPatient);
    o1 = o[:3];
    o2 = o[3:];
    oh = np.cross(o1,o2);
    or1 = np.asarray([0.6,0.6,-0.2]);
    o2new = np.cross(oh,or1);
    theta = np.arccos(np.dot(o2,o2new)/np.sqrt(np.sum(o2**2)*np.sum(o2new**2)))*180/3.1416;
    theta = theta * np.sign(np.dot(oh,np.cross(o2,o2new)));
    im_max = np.percentile(t.pixel_array.flatten(),99);
    res = imrotate(np.array(np.clip(np.array(t.pixel_array,dtype=np.float)/im_max*256,0,255),dtype=np.uint8),theta);
    if label is None:
        return res;
    else:
        lab = imrotate(label,theta);
        return res,lab

def getAlignImg(t,label = None):#!!!notice, only take uint8 type for the imrotate function!!!
    print 'CALLING GET_ALIGN'
    f = lambda x:np.asarray([float(a) for a in x]);
    o = f(t.ImageOrientationPatient);
    o1 = o[:3];
    o2 = o[3:];
    oh = np.cross(o1,o2);
    or1 = np.asarray([0.6,0.6,-0.2]);
    o2new = np.cross(oh,or1);
    theta = np.arccos(np.dot(o2,o2new)/np.sqrt(np.sum(o2**2)*np.sum(o2new**2)))*180/3.1416;
    theta = theta * np.sign(np.dot(oh,np.cross(o2,o2new)));
    im_max = np.percentile(t.pixel_array.flatten(),99);
    res = imrotate(np.array(np.clip(np.array(t.pixel_array,dtype=np.float)/im_max*256,0,255),
        dtype=np.uint8),theta);
    if label is None:
        return res;
    else:
        lab = imrotate(label,theta);
        return res,lab

def check(n):
    a = [(n>>i) & 1 for i in range(8)]
    a.insert(4, 0) # make the 3x3 unit
    # if up, down, left, right all are 1, you cannot make a hole
    if a[1] & a[3] & a[5] & a[7]:return False
    #if sum(a)==1: return False
    a = np.array(a).reshape((3,3))
    # segments
    n = label(a, strc)[1]
    # if sum is 0, it is a isolate point, you cannot remove it.
    # if number of segments > 2, you cannot split them.
    return n<2
    return a.sum()>1 and n<2
    if a.sum()==1 or n>2: return 2
    if a.sum()>1 and n<2: return 1
    return 0

def run(self, ips, snap, img, para = None):
        self.ips.lut[:] = self.buflut
        ndimg.gaussian_filter(snap, para['sigma'], output=img)
        mark = img<para['thr'] if para['ud'] else img>para['thr']

        markers, n = ndimg.label(mark, np.ones((3,3)), output=np.uint16)
        if not para['ud']:img[:] = 255-img
        mark = watershed(img, markers, line=True, conn=para['con']+1)
        mark = np.multiply((mark==0), 255, dtype=np.uint8)
        if para['type'] == 'white line':
            img[:] = mark
        if para['type'] == 'gray line':
            np.minimum(snap, mark, out=img)
        if para['type'] == 'white line on ori':
            #img //=2
            np.maximum(snap, mark, out=img)

def run(self, ips, snap, img, para = None):
        #denoised = rank.median(img, disk(para['sigma']))
        #gradient = rank.gradient(denoised, disk(para['gdt']))
        ndimg.gaussian_filter(snap, para['sigma'], output=img)

        markers, n = ndimg.label(ips.get_msk(), np.ones((3,3)), output=np.uint16)
        if not para['ud']:img[:] = 255-img
        mark = watershed(img, markers, line=True, conn=para['con']+1)
        mark = np.multiply((mark==0), 255, dtype=np.uint8)

        if para['type'] == 'white line':
            img[:] = mark
        if para['type'] == 'gray line':
            np.minimum(snap, mark, out=img)
        if para['type'] == 'white line on ori':
            np.maximum(snap, mark, out=img)

def run(self, ips, imgs, para = None):
        k, unit = ips.unit
        strc = generate_binary_structure(3, 1 if para['con']=='4-connect' else 2)

        lab, n = label(imgs==0 if para['inv'] else imgs, strc, output=np.uint16)
        idx = (np.ones(n+1)*(0 if para['inv'] else para['front'])).astype(np.uint8)
        ls = regionprops(lab)

        for i in ls:
            if para['vol'] == 0: break
            if para['vol']>0:
                if i.area*k**3 < para['vol']: idx[i.label] = para['back']
            if para['vol']<0:
                if i.area*k**3 >= -para['vol']: idx[i.label] = para['back']

        for i in ls:
            if para['dia'] == 0: break
            d = norm(np.array(i.bbox[:3]) - np.array(i.bbox[3:]))
            if para['dia']>0:
                if d*k < para['dia']: idx[i.label] = para['back']
            if para['dia']<0:
                if d*k >= -para['dia']: idx[i.label] = para['back']

        idx[0] = para['front'] if para['inv'] else 0
        imgs[:] = idx[lab]

def __init__(self,path):
        blur_radius = 1.0
        threshold = 50

        img = misc.imread(path)
        self.origin = img
        # smooth the image (to remove small objects)
        imgf = ndimage.gaussian_filter(img, blur_radius)
        threshold = 50
        # find connected components
        self.img, self.count = ndimage.label(imgf > threshold)
        self.labels = self.calculate_labels()


    # Return
    # ----------
    # labels : dictionary, key = label, value = list of bounding in order y1, y2, x1, x2

def detect_objects_heatmap(heatmap):
    data = 256 * heatmap
    data_max = filters.maximum_filter(data, 3)
    maxima = (data == data_max)
    data_min = filters.minimum_filter(data, 3)
    diff = ((data_max - data_min) > 0.3)
    maxima[diff == 0] = 0
    labeled, num_objects = ndimage.label(maxima)
    slices = ndimage.find_objects(labeled)
    objects = np.zeros((num_objects, 2), dtype=np.int32)
    pidx = 0
    for (dy, dx) in slices:
        pos = [(dy.start + dy.stop - 1) // 2, (dx.start + dx.stop - 1) // 2]
        if heatmap[pos[0], pos[1]] > config.CENTER_TR:
            objects[pidx, :] = pos
            pidx += 1
    return objects[:pidx]

def findMaximaOnFG(self, param):
        self.defineFG(param)
        #self.smooth_corr()
        self.coorExtract = [0, 0]
        xmin, ymin = self.coorExtract


        img=self.candi
        img [self.FG ==0] =0

        im = img_as_float(img)
        image_max = ndimage.maximum_filter(im, size=10, mode='constant')
        coordinates = peak_local_max(im, min_distance=10)

        tep=np.zeros(self.candi.shape)
        for i,ide in enumerate(coordinates):
            tep[ide[0],ide[1]] = self.candi[ide[0],ide[1]]
        lbl = ndimage.label(tep)[0]
        centerc = np.round(ndimage.measurements.center_of_mass(tep, lbl, range(1,np.max(lbl)+1)))
        if centerc.size > 0:
            self.centersX = centerc[:,0].astype(int)
            self.centersY = centerc[:,1].astype(int)
        self.nComponents = len(self.centersX)

def extract_digits(self, image):
        """
        Extract digits from a binary image representing a sudoku
        :param image: binary image/sudoku
        :return: array of digits and their probabilities
        """
        prob = np.zeros(4, dtype=np.float32)
        digits = np.zeros((4, 9, 9), dtype=object)
        for i in range(4):
            labeled, features = label(image, structure=CROSS)
            objs = find_objects(labeled)
            for obj in objs:
                roi = image[obj]
                # center of bounding box
                cy = (obj[0].stop + obj[0].start) / 2
                cx = (obj[1].stop + obj[1].start) / 2
                dists = cdist([[cy, cx]], CENTROIDS, 'euclidean')
                pos = np.argmin(dists)
                cy, cx = pos % 9, pos / 9
                # 28x28 image, center relative to sudoku
                prediction = self.classifier.classify(morph(roi))
                if digits[i, cy, cx] is 0:
                    # Newly found digit
                    digits[i, cy, cx] = prediction
                    prob[i] += prediction[0, 0]
                elif prediction[0, 0] > digits[i, cy, cx][0, 0]:
                    # Overlapping! (noise), choose the most probable prediction
                    prob[i] -= digits[i, cy, cx][0, 0]
                    digits[i, cy, cx] = prediction
                    prob[i] += prediction[0, 0]
            image = np.rot90(image)
        logging.info(prob)
        return digits[np.argmax(prob)]

def extract_digits(self, image):
        """
        Extract digits from a binary image representing a sudoku
        :param image: binary image/sudoku
        :return: array of digits and their probabilities
        """
        prob = np.zeros(4, dtype=np.float32)
        digits = np.zeros((4, 9, 9), dtype=object)
        for i in range(4):
            labeled, features = label(image, structure=CROSS)
            objs = find_objects(labeled)
            for obj in objs:
                roi = image[obj]
                # center of bounding box
                cy = (obj[0].stop + obj[0].start) / 2
                cx = (obj[1].stop + obj[1].start) / 2
                dists = cdist([[cy, cx]], CENTROIDS, 'euclidean')
                pos = np.argmin(dists)
                cy, cx = pos % 9, pos / 9
                # 28x28 image, center relative to sudoku
                prediction = self.classifier.classify(morph(roi))
                if digits[i, cy, cx] is 0:
                    # Newly found digit
                    digits[i, cy, cx] = prediction
                    prob[i] += prediction[0, 0]
                elif prediction[0, 0] > digits[i, cy, cx][0, 0]:
                    # Overlapping! (noise), choose the most probable prediction
                    prob[i] -= digits[i, cy, cx][0, 0]
                    digits[i, cy, cx] = prediction
                    prob[i] += prediction[0, 0]
            image = np.rot90(image)
        logging.info(prob)
        return digits[np.argmax(prob)]

def filter_isolated_cells(array, struct):

    filtered_array = np.copy(array)
    id_regions, num_ids = ndimage.label(filtered_array, structure=struct)
    id_sizes = np.array(ndimage.sum(array, id_regions, range(num_ids + 1)))
    area_mask = (id_sizes == 1)
    filtered_array[area_mask[id_regions]] = 0

    return filtered_array

#Decide if given spectrum shows bird sounds or noise only

def hysteresis_thresholding(probs: np.array, candidates: np.array, low_threshold: float, high_threshold: float):
    low_mask = candidates & (probs > low_threshold)
    # Connected components extraction
    label_components, count = label(low_mask, np.ones((3, 3)))
    # Keep components with high threshold elements
    good_labels = np.unique(label_components[low_mask & (probs > high_threshold)])
    label_masks = np.zeros((count + 1,), bool)
    label_masks[good_labels] = 1
    return label_masks[label_components]

def has_uniform_seed_margin(self, seed_margin=20.0):
        """Test if a subvolume has a margin of uniform label around its seed.

        Parameters
        ----------
        seed_margin : float, optional
            The minimum acceptable margin of uniform target label around the seed
            voxel (in nm, default 20.0).

        Returns
        -------
        bool
            True if the rectangular margin around the seed position is uniform.
        """
        margin = np.ceil(np.reciprocal(np.array(CONFIG.volume.resolution),
                                       dtype=np.float64) * seed_margin).astype(np.int64)

        mask_target = self.label_mask
        # If data is unlabeled, can not test so always succeed.
        if mask_target is None:
            return True
        # Seed location in the mask accounting for offset of label from image.
        ctr = self.seed - (np.asarray(self.image.shape) - np.asarray(mask_target.shape)) // 2
        seed_fov = (ctr - margin, ctr + margin + 1)
        seed_region = mask_target[seed_fov[0][0]:seed_fov[1][0],
                                  seed_fov[0][1]:seed_fov[1][1],
                                  seed_fov[0][2]:seed_fov[1][2]]
        return np.all(seed_region)

def __next__(self):
        subv = six.next(self.subvolume_generator)

        label_im, _ = ndimage.label(subv.label_mask)
        label_axis_margin = (np.array(subv.image.shape) - np.array(subv.label_mask.shape)) // 2
        seed_label = label_im[tuple(subv.seed - label_axis_margin)]

        subv.label_mask = label_im == seed_label

        return subv

def get_subvolume(self, bounds):
        subvol_shape = bounds.stop - bounds.start
        label_shape = subvol_shape - 2 * bounds.label_margin
        parent_bounds = SubvolumeBounds(self.world_coord_to_parent(bounds.start),
                                        self.world_coord_to_parent(bounds.stop),
                                        label_margin=self.world_coord_to_parent(bounds.label_margin))
        subvol = self.parent.get_subvolume(parent_bounds)
        subvol.image = subvol.image.reshape(
                [subvol_shape[0], self.scale[0],
                 subvol_shape[1], self.scale[1],
                 subvol_shape[2], self.scale[2]]).mean(5).mean(3).mean(1)

        if subvol.label_mask is not None:
            # Downsample body mask by considering blocks where the majority
            # of voxels are in the body to be in the body. Alternatives are:
            # - Conjunction (tends to introduce false splits)
            # - Disjunction (tends to overdilate and merge)
            # - Mode label (computationally expensive)
            if CONFIG.volume.label_downsampling == 'conjunction':
                subvol.label_mask = subvol.label_mask.reshape(
                        [label_shape[0], self.scale[0],
                         label_shape[1], self.scale[1],
                         label_shape[2], self.scale[2]]).all(5).all(3).all(1)
            else:
                subvol.label_mask = subvol.label_mask.reshape(
                        [label_shape[0], self.scale[0],
                         label_shape[1], self.scale[1],
                         label_shape[2], self.scale[2]]).mean(5).mean(3).mean(1) > 0.5

        # Note that this is not a coordinate xform to parent in the typical
        # sense, just a rescaling of the coordinate in the subvolume-local
        # coordinates. Hence no similar call in VolumeView.get_subvolume.
        subvol.seed = self.parent_coord_to_world(subvol.seed)

        return subvol

def get_largest_component(self, closing_shape=None):
        mask, bounds = self._get_bounded_mask(closing_shape)

        label_im, num_labels = ndimage.label(mask)
        label_sizes = ndimage.sum(mask, label_im, range(num_labels + 1))
        label_im[(label_sizes < label_sizes.max())[label_im]] = 0
        label_im = np.minimum(label_im, 1)

        if label_im[tuple(self.seed - bounds[0])] == 0:
            logging.warning('Seed voxel ({}) is not in connected component.'.format(np.array_str(self.seed)))

        return label_im, bounds

def get_seeded_component(self, closing_shape=None):
        mask, bounds = self._get_bounded_mask(closing_shape)

        label_im, _ = ndimage.label(mask)
        seed_label = label_im[tuple(self.seed - bounds[0])]
        if seed_label == 0:
            raise ValueError('Seed voxel (%s) is not in body.', np.array_str(self.seed))
        label_im[label_im != seed_label] = 0
        label_im[label_im == seed_label] = 1

        return label_im, bounds

def compute_centroids(image):
    """Find the centroids of all nonzero connected blobs in `image`.

    Parameters
    ----------
    image : ndarray
        The input image.

    Returns
    -------
    label_image : ndarray of int
        The input image, with each connected region containing a different
        integer label.

    Examples
    --------
    >>> image = np.array([[1, 0, 1, 0, 0, 1, 1],
    ...                   [1, 0, 0, 1, 0, 0, 0]])
    >>> labels, centroids = compute_centroids(image)
    >>> print(labels)
    [[1 0 2 0 0 3 3]
     [1 0 0 2 0 0 0]]
    >>> centroids
    array([[ 0.5,  0. ],
           [ 0.5,  2.5],
           [ 0. ,  5.5]])
    """
    connectivity = np.ones((3,) * image.ndim)
    labeled_image = ndi.label(image, connectivity)[0]
    nz = np.nonzero(labeled_image)
    nzpix = labeled_image[nz]
    sizes = np.bincount(nzpix)
    coords = np.transpose(nz)
    grouping = np.argsort(nzpix)
    sums = np.add.reduceat(coords[grouping], np.cumsum(sizes)[:-1])
    means = sums / sizes[1:, np.newaxis]
    return labeled_image, means

def mesh_sizes(skeleton):
    """Compute the area in pixels of the spaces between skeleton branches.

    This only makes sense for 2D images.

    Parameters
    ----------
    skeleton : array, shape (M, N)
        An image containing a single-pixel-wide closed skeleton.

    Returns
    -------
    sizes : array of int, shape (P,)
        The sizes of all spaces delineated by the skeleton *not* touching
        the borders.

    Examples
    --------
    >>> image = np.array([[0, 0, 1, 0, 0],
    ...                   [0, 0, 1, 1, 1],
    ...                   [0, 0, 1, 0, 0],
    ...                   [0, 1, 0, 1, 0]])
    >>> print(mesh_sizes(image))
    []
    >>> from skan.nputil import pad
    >>> image2 = pad(image, 1)  # make sure mesh not touching border
    >>> print(mesh_sizes(image2))  # sizes in row order of first pixel in space
    [7 2 3 1]
    """
    spaces = ~skeleton.astype(bool)
    labeled = ndi.label(spaces)[0]
    touching_border = np.unique(np.concatenate((labeled[0], labeled[-1],
                                                labeled[:, 0],
                                                labeled[:, -1])))
    sizes = np.bincount(labeled.flat)
    sizes[touching_border] = 0
    sizes = sizes[sizes != 0]
    return sizes

def image_summary(skeleton, *, spacing=1):
    """Compute some summary statistics for an image.

    Parameters
    ----------
    skeleton : array, shape (M, N)
        The input image.

    Other Parameters
    ----------------
    spacing : float or array of float, shape (`skeleton.ndim`,)
        The resolution along each axis of `skeleton`.

    Returns
    -------
    stats : pandas.DataFrame
        Selected statistics about the image.
    """
    stats = pd.DataFrame()
    stats['scale'] = [spacing]
    g, coords, degimg = csr.skeleton_to_csgraph(skeleton, spacing=spacing)
    degrees = np.diff(g.indptr)
    num_junctions = np.sum(degrees > 2)
    stats['number of junctions'] = num_junctions
    pixel_area = (spacing ** skeleton.ndim if np.isscalar(spacing) else
                  np.prod(spacing))
    stats['area'] = np.prod(skeleton.shape) * pixel_area
    stats['junctions per unit area'] = (stats['number of junctions'] /
                                        stats['area'])
    sizes = mesh_sizes(skeleton)
    stats['average mesh area'] = np.mean(sizes)
    stats['median mesh area'] = np.median(sizes)
    stats['mesh area standard deviation'] = np.std(sizes)

    structure = np.ones((3,) * skeleton.ndim)
    stats['number of disjoint skeletons'] = ndi.label(skeleton, structure)[1]

    return stats

def find_peaks_minmax(z, separation=5., threshold=10.):
    """
    Method to locate the positive peaks in an image by comparing maximum
    and minimum filtered images.
    Parameters
    ----------
    z : ndarray
        Matrix of image intensities.
    separation : float
        Expected distance between peaks.
    threshold : float
        Minimum difference between maximum and minimum filtered images.
    Returns
    -------
    peaks: array with dimensions (npeaks, 2) that contains the x, y coordinates
           for each peak found in the image.
    """
    data_max = ndi.filters.maximum_filter(z, separation)
    maxima = (z == data_max)
    data_min = ndi.filters.minimum_filter(z, separation)
    diff = ((data_max - data_min) > threshold)
    maxima[diff == 0] = 0
    labeled, num_objects = ndi.label(maxima)
    peaks = np.array(
        ndi.center_of_mass(z, labeled, range(1, num_objects + 1)))

    return clean_peaks(peaks)

def candidates_to_image(cands,radius):
    image_names = []
    for subset in xrange(0,10):
        subset_names = glob.glob("../data/subset{0}/*.mhd".format(subset))
        names = [os.path.split(x)[1].replace('.mhd','') for x in subset_names]
        image_names.append(names)
    previous_candidate = ""
    images = []
    image = []
    origin = []
    spacing = []
    number = 0
    for candidate in tqdm(cands.values):
        if candidate[0] != previous_candidate:
            number = 0
            previous_candidate = candidate[0]
            for image_subset in xrange(0,10):
                if candidate[0] in image_names[image_subset]:
                    image,origin,spacing = load_itk_image("../data/subset{0}/{1}.mhd".format(image_subset,candidate[0]))
                    break
        coords = world_2_voxel([candidate[3],candidate[2],candidate[1]],origin,spacing)
        im = image_part_from_candidate(image,coords,radius)
        #images.append(im)
        if candidate[4]:
            label = "true"
        else:
            label = "false"
        scipy.misc.imsave('../data/samples/{0}_{1}_{2}.jpg'.format(candidate[0],number,label), im)
        number += 1
    return images

def average_boxes(hot_windows, old_boxes, 
                  image_shape):
    hot_boxes = initialize_center_box(hot_windows)
    new_boxes = add_center_box(hot_boxes, old_boxes)
    # print('new_boxes', new_boxes)
    filtered_boxes = []
    for new_box in new_boxes:
        # Draw boxes only if the confidence level is above 1
        if new_box[-1] > 2:
            filtered_boxes.append(new_box)
    new_windows = []
    # convert center-width-height to lefttop-rightbottom format
    for filtered_box in filtered_boxes:
        new_center, new_width, new_height,new_move, new_prob = filtered_box
        new_windows.append(((int(new_center[0]-new_width), int(new_center[1]-new_height)), 
            (int(new_center[0]+new_width), int(new_center[1]+new_height))))
    # Create a heatmap
    heatmap = create_heatmap(new_windows, image_shape)
    # Check if there is any overlap of windows
    if np.unique(heatmap)[-1] >= 2:
        labels = ndi.label(heatmap)[0]
        heatmap_2 = np.zeros_like(heatmap)
        heatmap_2[heatmap>=2] = 1
        labels_2 = ndi.label(heatmap_2)
        array_2 = np.argwhere(labels_2[0])
        for car_number in range(1, labels_2[1]+1):
            # Find pixels with each car_number label value
            nonzero = (labels_2[0] == car_number).nonzero()
            # Identify x and y values of those pixels
            num = labels[nonzero[0][0], nonzero[1][0]]
            labels[labels == num] = 0
        heatmap = labels + heatmap_2
        new_windows = find_windows_from_heatmap(heatmap)

    # return the boxes with high confidence and new set of probability array
    return new_windows, new_boxes

def _run_interface(self, runtime):
        in_file = nb.load(self.inputs.in_file)
        data = in_file.get_data()
        mask = np.zeros_like(data, dtype=np.uint8)
        mask[data <= 0] = 1

        # Pad one pixel to control behavior on borders of binary_opening
        mask = np.pad(mask, pad_width=(1,), mode='constant', constant_values=1)

        # Remove noise
        struc = nd.generate_binary_structure(3, 2)
        mask = nd.binary_opening(mask, structure=struc).astype(
            np.uint8)

        # Remove small objects
        label_im, nb_labels = nd.label(mask)
        if nb_labels > 2:
            sizes = nd.sum(mask, label_im, list(range(nb_labels + 1)))
            ordered = list(reversed(sorted(zip(sizes, list(range(nb_labels + 1))))))
            for _, label in ordered[2:]:
                mask[label_im == label] = 0

        # Un-pad
        mask = mask[1:-1, 1:-1, 1:-1]

        # If mask is small, clean-up
        if mask.sum() < 500:
            mask = np.zeros_like(mask, dtype=np.uint8)

        out_img = in_file.__class__(mask, in_file.affine, in_file.header)
        out_img.header.set_data_dtype(np.uint8)

        out_file = fname_presuffix(self.inputs.in_file,
                                   suffix='_rotmask', newpath='.')
        out_img.to_filename(out_file)
        self._results['out_file'] = out_file
        return runtime

def makeMask(self):
        #Narrow down to pixels which measured something every frame
        counts = np.sum(self.Zs > 0, 2)
        self.Mask = (counts == self.Zs.shape[2])
        #Find biggest connected component out of remaining pixels
        ILabel, NLabels = ndimage.label(self.Mask)
        idx = np.argmax(ndimage.sum(self.Mask, ILabel, range(NLabels+1)))
        self.Mask = (ILabel == idx)
        plt.imshow(self.Mask)
        plt.show()

    #Actually sets up all of the vertex, face, and adjacency structures in the 
    #PolyMeshObject

def binary_volume_opening(vol, minvol):
    lb_vol, num_objs = label(vol)
    lbs = np.arange(1, num_objs + 1)
    v = labeled_comprehension(lb_vol > 0, lb_vol, lbs, np.sum, np.int, 0)
    ix = np.isin(lb_vol, lbs[v >= minvol])
    newvol = np.zeros(vol.shape)
    newvol[ix] = vol[ix]
    return newvol

def autofind_pois(self, neighborhood_size=1, min_threshold=10000, max_threshold=1e6):
        """Automatically search the xy scan image for POIs.

        @param neighborhood_size: size in microns.  Only the brightest POI per neighborhood will be found.

        @param min_threshold: POIs must have c/s above this threshold.

        @param max_threshold: POIs must have c/s below this threshold.
        """

        # Calculate the neighborhood size in pixels from the image range and resolution
        x_range_microns = np.max(self.roi_map_data[:, :, 0]) - np.min(self.roi_map_data[:, :, 0])
        y_range_microns = np.max(self.roi_map_data[:, :, 1]) - np.min(self.roi_map_data[:, :, 1])
        y_pixels = len(self.roi_map_data)
        x_pixels = len(self.roi_map_data[1, :])

        pixels_per_micron = np.max([x_pixels, y_pixels]) / np.max([x_range_microns, y_range_microns])
        # The neighborhood in pixels is nbhd_size * pixels_per_um, but it must be 1 or greater
        neighborhood_pix = int(np.max([math.ceil(pixels_per_micron * neighborhood_size), 1]))

        data = self.roi_map_data[:, :, 3]

        data_max = filters.maximum_filter(data, neighborhood_pix)
        maxima = (data == data_max)
        data_min = filters.minimum_filter(data, 3 * neighborhood_pix)
        diff = ((data_max - data_min) > min_threshold)
        maxima[diff is False] = 0

        labeled, num_objects = ndimage.label(maxima)
        xy = np.array(ndimage.center_of_mass(data, labeled, range(1, num_objects + 1)))

        for count, pix_pos in enumerate(xy):
            poi_pos = self.roi_map_data[pix_pos[0], pix_pos[1], :][0:3]
            this_poi_key = self.add_poi(position=poi_pos, emit_change=False)
            self.rename_poi(poikey=this_poi_key, name='spot' + str(count), emit_change=False)

        # Now that all the POIs are created, emit the signal for other things (ie gui) to update
        self.signal_poi_updated.emit()

def candidates_to_image(cands,radius):
    image_names = []
    for subset in xrange(0,10):
        subset_names = glob.glob("../data/subset{0}/*.mhd".format(subset))
        names = [os.path.split(x)[1].replace('.mhd','') for x in subset_names]
        image_names.append(names)
    previous_candidate = ""
    images = []
    image = []
    origin = []
    spacing = []
    number = 0
    for candidate in tqdm(cands.values):
        if candidate[0] != previous_candidate:
            number = 0
            previous_candidate = candidate[0]
            for image_subset in xrange(0,10):
                if candidate[0] in image_names[image_subset]:
                    image,origin,spacing = load_itk_image("../data/subset{0}/{1}.mhd".format(image_subset,candidate[0]))
                    break
        coords = world_2_voxel([candidate[3],candidate[2],candidate[1]],origin,spacing)
        im = image_part_from_candidate(image,coords,radius)
        #images.append(im)
        if candidate[4]:
            label = "true"
        else:
            label = "false"
        scipy.misc.imsave('../data/samples/{0}_{1}_{2}.jpg'.format(candidate[0],number,label), im)
        number += 1
    return images

def count(n):
    a = [(n>>i) & 1 for i in range(8)]
    if sum(a)<=1:return False
    if a[1] & a[3] & a[5] & a[7]:return False
    a = np.array([[a[0],a[1],a[2]],
                  [a[7],  0 ,a[3]],
                  [a[6],a[5],a[4]]])
    n = label(a, strc)[1]
    return n<2

def run(self, ips, snap, img, para = None):
        if para == None: para = self.para
        ips.lut = self.lut
        msk = img>para['thr']
        con = 1 if para['con']=='4-Connect' else 2
        strc = generate_binary_structure(2, con)
        msk = label(msk, strc, output = np.int16)

        IPy.show_img([msk[0]], ips.title+'-label')

def run(self, ips, snap, img, para = None):
        intenimg = WindowsManager.get(para['inten']).ips.img
        strc = ndimage.generate_binary_structure(2, 1 if para['con']=='4-connect' else 2)
        buf, n = ndimage.label(snap, strc, output=np.uint16)
        index = range(1, n+1)
        idx = (np.ones(n+1)*para['front']).astype(np.uint8)
        msk = np.ones(n, dtype=np.bool)

        if para['mean']>0: msk *= ndimage.mean(intenimg, buf, index)>=para['mean']
        if para['mean']<0: msk *= ndimage.mean(intenimg, buf, index)<-para['mean']
        if para['max']>0: msk *= ndimage.maximum(intenimg, buf, index)>=para['max']
        if para['max']<0: msk *= ndimage.maximum(intenimg, buf, index)<-para['max']
        if para['min']>0: msk *= ndimage.minimum(intenimg, buf, index)>=para['min']
        if para['min']<0: msk *= ndimage.minimum(intenimg, buf, index)<-para['min']
        if para['sum']>0: msk *= ndimage.sum(intenimg, buf, index)>=para['sum']
        if para['sum']<0: msk *= ndimage.sum(intenimg, buf, index)<-para['sum']
        if para['std']>0: msk *= ndimage.standard_deviation(intenimg, buf, index)>=para['std']
        if para['std']<0: msk *= ndimage.standard_deviation(intenimg, buf, index)<-para['std']


        xy = ndimage.center_of_mass(intenimg, buf, index)
        xy = np.array(xy).round(2).T

        idx[1:][~msk] = para['back']
        idx[0] = 0
        img[:] = idx[buf]

        WindowsManager.get(para['inten']).ips.mark = RGMark((xy.T, msk))
        WindowsManager.get(para['inten']).ips.update = True