Python skimage.measure.label() Examples

def filter_cloudmask(cloudmask, threshold=1, connectivity=1):
    """Filter a given cloudmask for small cloud objects defined by their pixel
    number. 
    
    Parameters:
        cloudmask (ndarray): 2d binary cloud mask (optional with NaNs).
        threshold (int): minimum pixel number of objects remaining in cloudmask.
        connectivity (int):  Maximum number of orthogonal hops to consider
            a pixel/voxel as a neighbor (see :func:`skimage.measure.label`).
    
    Return:
        ndarray: filtered cloudmask without NaNs.
    """
    cloudmask[np.isnan(cloudmask)] = 0
    labels = measure.label(cloudmask, connectivity=connectivity)
    props = measure.regionprops(labels)
    area = [prop.area for prop in props]

    # Find objects < threshold pixle number, get their labels, set them to 0-clear.
    smallclouds = [t[0] for t in filter(lambda a: a[1] < threshold, enumerate(area, 1))]
    for label in smallclouds:
        cloudmask[labels == label] = 0

    return cloudmask 

def centroids_of_connected_components(bitmap, threshold=0.05, rescale=1.0):
  # TODO: don't do raw (binary) threshold; instead use P(y) as weighting for centroid
  #       e.g. https://arxiv.org/abs/1806.03413 sec 3.D
  #       update: this didn't help much :/ centroid weighted by intensities moved only up
  #       to a single pixel (guess centroids are already quite evenly dispersed)
  #       see https://gist.github.com/matpalm/20a3974ceb7f632f935285262fac4e98
  # TODO: hunt down the x/y swap between PIL and label db :/

  # threshold
  mask = bitmap > threshold
  bitmap = np.zeros_like(bitmap)
  bitmap[mask] = 1.0
  # calc connected components
  all_labels = measure.label(bitmap)
  # return centroids
  centroids = []
  for region in measure.regionprops(label_image=all_labels):
    cx, cy = map(lambda p: int(p*rescale), (region.centroid[0], region.centroid[1]))
    centroids.append((cx, cy))
  return centroids 

def segment_body (image, smooth=1, th=-300):
    blur = scipy.ndimage.filters.gaussian_filter(image, smooth, mode='constant')
    binary = np.array(blur < th, dtype=np.uint8)

    # body is a rough region covering human body
    body = np.zeros_like(binary)
    for i, sl in enumerate(binary):
        #H, W = sl.shape
        ll = measure.label(sl, background=1)   # connected components
        # biggest CC should be body
        pp = measure.regionprops(ll)
        boxes = [(x.area, x.bbox, x.filled_image) for x in pp if x.label != 0]  # label 0 is air
        boxes = sorted(boxes, key = lambda x: -x[0])
        if len(boxes) == 0:
            continue
        y0, x0, y1, x1 = boxes[0][1]
        body[i,y0:y1,x0:x1] = boxes[0][2]
        pass
    return body, None 

def filter_blank_slices_thick(img_vol, label_vol, weight_vol, threshold=50):
    """
    Function to filter blank slices from the volume using the label volume
    :param np.ndarray img_vol: orig image volume
    :param np.ndarray label_vol: label images (ground truth)
    :param np.ndarray weight_vol: weight corresponding to labels
    :param int threshold: threshold for number of pixels needed to keep slice (below = dropped)
    :return:
    """
    # Get indices of all slices with more than threshold labels/pixels
    select_slices = (np.sum(label_vol, axis=(0, 1)) > threshold)

    # Retain only slices with more than threshold labels/pixels
    img_vol = img_vol[:, :, select_slices, :]
    label_vol = label_vol[:, :, select_slices]
    weight_vol = weight_vol[:, :, select_slices]

    return img_vol, label_vol, weight_vol


# weight map generator 

def GetCC(match1, match2, mask1, mask2, score) : 
	match1_arr = np.array(match1).astype(int)
	mask1[match1_arr[:, 0].flatten(), match1_arr[:, 1].flatten()] = 1
	label1, nb_label1 = measure.label(mask1, connectivity=2, return_num=True)
	dict_match = dict(zip(match1, match2))
	dict_score = dict(zip(match1, score))
	
	CC = []
	CC_score = np.zeros(nb_label1)
	for i in range(1, nb_label1 + 1) : 
		CC.append({})
		posx, posy = np.where(label1 == i)
		tmp_match1 = [(posx[j], posy[j]) for j in range(len(posx))]
		tmp_match2 = [dict_match[item] for item in tmp_match1]
		CC[i-1] = dict(zip(tmp_match1, tmp_match2))
		CC[i-1]['mask2shape'] = mask2.shape
		CC_score[i -1 ] = sum(dict_score[item] for item in tmp_match1)
	return CC, CC_score 

def extend_nonforest(config, buffered_array, full_array):

   """ 
   After buffering forest, get changes that are spatially connected 
   to previous change or non-forest
   """

   # convert full array into connected labels
   ones_array = np.copy(full_array[0,:,:])
   ones_array[ones_array > 0] = 1
   ones_array[ones_array != 1] = 0

   label_image = label(ones_array, neighbors=8)
   labels = np.unique(label_image)

   for _label in labels:
	indices = np.where(label_image == _label)
        if ones_array[indices][0] > 0:
	    connected_parts = buffered_array[0,:,:][indices].sum()
	    if connected_parts > 0:
	        for i in range(full_array.shape[0]):
	            buffered_array[i,:,:][indices] = full_array[i,:,:][indices]

   return buffered_array 

def debug_img(img, bitmap, logistic_output):
  # create a debug image with three columns; 1) original RGB. 2) black/white
  # bitmap of labels 3) black/white bitmap of predictions (with centroids coloured
  # red.
  h, w, _channels = bitmap.shape
  canvas = Image.new('RGB', (w*3, h), (50, 50, 50))
  # original input image on left
  img = zero_centered_array_to_pil_image(img)
  img = img.resize((w, h))
  canvas.paste(img, (0, 0))
  # label bitmap in center
  canvas.paste(bitmap_to_pil_image(bitmap), (w, 0))
  # logistic output on right
  canvas.paste(bitmap_to_pil_image(logistic_output), (w*2, 0))
  # draw red dots on right hand side image corresponding to
  # final thresholded prediction
  draw = ImageDraw.Draw(canvas)
  for y, x in centroids_of_connected_components(logistic_output):
    draw.rectangle((w*2+x,y,w*2+x,y), fill='red')
  # finally draw blue lines between the three to delimit boundaries
  draw.line([w,0,w,h], fill='blue')
  draw.line([2*w,0,2*w,h], fill='blue')
  draw.line([3*w,0,3*w,h], fill='blue')
  # done
  return canvas 

def __init__(self, threshold=0.4, iou_range=(0.5, 1.0), ignore_index=-1, min_instance_size=None,
                 use_last_target=False):
        """
        :param threshold: probability value at which the input is going to be thresholded
        :param iou_range: compute ROC curve for the the range of IoU values: range(min,max,0.05)
        :param ignore_index: label to be ignored during computation
        :param min_instance_size: minimum size of the predicted instances to be considered
        :param use_last_target: if True use the last target channel to compute AP
        """
        self.threshold = threshold
        # always have well defined ignore_index
        if ignore_index is None:
            ignore_index = -1
        self.iou_range = iou_range
        self.ignore_index = ignore_index
        self.min_instance_size = min_instance_size
        self.use_last_target = use_last_target 

def _find_overlapping_target(self, predicted_label, predicted, target, min_iou):
        """
        Return ground truth label which overlaps by at least 'min_iou' with a given input label 'p_label'
        or None if such ground truth label does not exist.
        """
        mask_predicted = predicted == predicted_label
        overlapping_labels = target[mask_predicted]
        labels, counts = np.unique(overlapping_labels, return_counts=True)
        # retrieve the biggest overlapping label
        target_label_ind = np.argmax(counts)
        target_label = labels[target_label_ind]
        # return target label if IoU greater than 'min_iou'; since we're starting from 0.5 IoU there might be
        # only one target label that fulfill this criterion
        mask_target = target == target_label
        # return target_label if IoU > min_iou
        if self._iou(mask_predicted, mask_target) > min_iou:
            return target_label
        return None 

def _filter_instances(self, input):
        """
        Filters instances smaller than 'min_instance_size' by overriding them with 'ignore_index'
        :param input: input instance segmentation
        :return: tuple: (instance segmentation with small instances filtered, set of unique labels without the 'ignore_index')
        """
        if self.min_instance_size is not None:
            labels, counts = np.unique(input, return_counts=True)
            for label, count in zip(labels, counts):
                if count < self.min_instance_size:
                    mask = input == label
                    input[mask] = self.ignore_index

        labels = set(np.unique(input))
        labels.discard(self.ignore_index)

        return input, labels 

def plot_regions(self, fill=True, bgimage=None, alpha=0.5):
        import pylab as pl
        ax = pl.gca()
        assert isinstance(ax, pl.Axes)

        colors = i12.JET_12

        self._plot_background(bgimage)

        for label in self.regions:
            color = colors[i12.LABELS.index(label)] / 255.

            for region in self.regions[label]:
                t = region['top']
                l = self.facade_left + region['left']
                b = region['bottom']
                r = self.facade_left + region['right']
                patch = pl.Rectangle((l, t), r - l, b - t, color=color, fill=fill, alpha=alpha)
                ax.add_patch(patch) 

def DynamicWatershedAlias(p_img, lamb, p_thresh = 0.5):
    """
    Applies our dynamic watershed to 2D prob/dist image.
    """
    b_img = (p_img > p_thresh) + 0
    Probs_inv = PrepareProb(p_img)

    Hrecons = HreconstructionErosion(Probs_inv, lamb)
    markers_Probs_inv = find_maxima(Hrecons, mask = b_img)
    markers_Probs_inv = label(markers_Probs_inv)
    ws_labels = watershed(Hrecons, markers_Probs_inv, mask=b_img)
    arrange_label = ArrangeLabel(ws_labels)
    wsl = generate_wsl(arrange_label)
    arrange_label[wsl > 0] = 0
    
    return arrange_label

#### 

def computeEvaluationMask(maskDIR, resolution, level):
    """Computes the evaluation mask.
    
    Args:
        maskDIR:    the directory of the ground truth mask
        resolution: Pixel resolution of the image at level 0
        level:      The level at which the evaluation mask is made
        
    Returns:
        evaluation_mask
    """
    slide = openslide.open_slide(maskDIR)
    dims = slide.level_dimensions[level]
    pixelarray = np.zeros(dims[0]*dims[1], dtype='uint')
    pixelarray = np.array(slide.read_region((0,0), level, dims))
    distance = nd.distance_transform_edt(255 - pixelarray[:,:,0])
    Threshold = 75/(resolution * pow(2, level) * 2) # 75µm is the equivalent size of 5 tumor cells
    binary = distance < Threshold
    filled_image = nd.morphology.binary_fill_holes(binary)
    evaluation_mask = measure.label(filled_image, connectivity = 2) 
    return evaluation_mask 

def convert_image_to_text(path):
    img = cv2.imread(path, 0)
    new_mask = np.zeros(img.shape, dtype=np.uint8)
    im2, contours, hierarchy = cv2.findContours(img.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    print('Total contours: {}'.format(len(contours)))
    small_area = 0
    for c in contours:
        area = cv2.contourArea(c)
        # print(area)
        if area < 100:
            small_area += 1
            continue
        cv2.drawContours(new_mask, [c], 0, (255, 255, 255), -1)

    # show_resized_image(new_mask)
    arr = measure.label(new_mask, connectivity=1)
    str1 = rle(arr)
    # print(str1)
    # tmp_mask = np.zeros(img.shape, dtype=np.uint8)
    # tmp_mask[arr == 1] = 255
    # show_resized_image(tmp_mask)
    # print(arr.min(), arr.max())
    print('Small contours: {}'.format(small_area))
    return str1 

def extractCornersFromSegmentation(segmentation, cornerTypeRange=[0, 13]):
    from skimage import measure
    orientationPoints = []
    for heatmapIndex in xrange(cornerTypeRange[0], cornerTypeRange[1]):
        heatmap = segmentation == heatmapIndex
        #heatmap = cv2.dilate(cv2.erode(heatmap, kernel), kernel)
        components = measure.label(heatmap, background=0)
        points = []
        for componentIndex in xrange(components.min()+1, components.max() + 1):
            ys, xs = (components == componentIndex).nonzero()
            points.append((xs.mean(), ys.mean()))
            continue
        orientationPoints.append(points)
        continue
    return orientationPoints

#Extract corners from heatmaps 

def LoadGTBatch(path, normalize=True):
    """
    Generic function to be used in multi-threaded batch read for GT.
    """
    image = ni.load(path)
    img = image.get_data()
    img = measure.label(img)
    wsl = generate_wsl(img[:,:,0])
    img[ img > 0 ] = 1
    wsl[ wsl > 0 ] = 1
    img[:,:,0] = img[:,:,0] - wsl
    if len(img.shape) == 3:
        img = img[:, :, 0].transpose()
    else:
        img = img.transpose()
    return img 

def Process(self, img_lbl_Mwgt, f, crop_n):
        """
        Apply necessary transformation to input image and label.
        """

        if self.UNet:
            img_lbl_Mwgt = self.Unet_cut(*img_lbl_Mwgt)

        img_lbl_Mwgt = f._apply_(*img_lbl_Mwgt)  # change _apply_

        if self.crop != 1:
            img_lbl_Mwgt = self.DivideImage(crop_n, *img_lbl_Mwgt)

        if self.random_crop:
            img_lbl_Mwgt = self.CropImgLbl(*img_lbl_Mwgt)

        if self.UNet:
            CheckNumberForUnet(img_lbl_Mwgt[0].shape[0])
            img_lbl_Mwgt = self.ReduceDimUNet(*img_lbl_Mwgt)    

        if self.Mean:
            img_lbl_Mwgt = self.SubtractMean(*img_lbl_Mwgt)

        return img_lbl_Mwgt 

def LoadGT(self, path, normalize=True):
        """
        Loads one single GT image.
        """
        image = ni.load(path)
        img = image.get_data()
        img = measure.label(img)
        wsl = generate_wsl(img[:,:,0])
        img[ img > 0 ] = 1
        wsl[ wsl > 0 ] = 1
        img[:,:,0] = img[:,:,0] - wsl
        if len(img.shape) == 3:
            img = img[:, :, 0].transpose()
        else:
            img = img.transpose()
        return img 

def DynamicWatershedAlias(p_img, lamb, p_thresh = 0.5):
    """
    Applies our dynamic watershed to 2D prob/dist image.
    """
    b_img = (p_img > p_thresh) + 0
    Probs_inv = PrepareProb(p_img)


    Hrecons = HreconstructionErosion(Probs_inv, lamb)
    markers_Probs_inv = find_maxima(Hrecons, mask = b_img)
    markers_Probs_inv = label(markers_Probs_inv)
    ws_labels = watershed(Hrecons, markers_Probs_inv, mask=b_img)
    arrange_label = ArrangeLabel(ws_labels)
    wsl = generate_wsl(arrange_label)
    arrange_label[wsl > 0] = 0
    

    return arrange_label 

def keep_largest_connected_components(mask):
    '''
    Keeps only the largest connected components of each label for a segmentation mask.
    '''

    out_img = np.zeros(mask.shape, dtype=np.uint8)

    for struc_id in [1, 2, 3]:

        binary_img = mask == struc_id
        blobs = measure.label(binary_img, connectivity=1)

        props = measure.regionprops(blobs)

        if not props:
            continue

        area = [ele.area for ele in props]
        largest_blob_ind = np.argmax(area)
        largest_blob_label = props[largest_blob_ind].label

        out_img[blobs == largest_blob_label] = struc_id

    return out_img 

def segment_body (image, smooth=1, th=-300):
    blur = scipy.ndimage.filters.gaussian_filter(image, smooth, mode='constant')
    binary = np.array(blur < th, dtype=np.uint8)

    # body is a rough region covering human body
    body = np.zeros_like(binary)
    for i, sl in enumerate(binary):
        #H, W = sl.shape
        ll = measure.label(sl, background=1)   # connected components
        # biggest CC should be body
        pp = measure.regionprops(ll)
        boxes = [(x.area, x.bbox, x.filled_image) for x in pp if x.label != 0]  # label 0 is air
        boxes = sorted(boxes, key = lambda x: -x[0])
        if len(boxes) == 0:
            continue
        y0, x0, y1, x1 = boxes[0][1]
        body[i,y0:y1,x0:x1] = boxes[0][2]
        pass
    return body, None 

def label_mask(pred, main_threshold=0.3, seed_threshold=0.7, w_pixel_t=20, pixel_t=100):
    av_pred = pred / 255.
    av_pred = av_pred[..., 0] * (1 - av_pred[..., 2])
    av_pred = 1 * (av_pred > seed_threshold)
    av_pred = av_pred.astype(np.uint8)

    y_pred = measure.label(av_pred, neighbors=8, background=0)
    props = measure.regionprops(y_pred)
    for i in range(len(props)):
        if props[i].area < w_pixel_t:
            y_pred[y_pred == i + 1] = 0
    y_pred = measure.label(y_pred, neighbors=8, background=0)

    nucl_msk = (255 - pred[..., 0])
    nucl_msk = nucl_msk.astype('uint8')
    y_pred = watershed(nucl_msk, y_pred, mask=(pred[..., 0] > main_threshold * 255), watershed_line=True)

    props = measure.regionprops(y_pred)

    for i in range(len(props)):
        if props[i].area < pixel_t:
            y_pred[y_pred == i + 1] = 0
    y_pred = measure.label(y_pred, neighbors=8, background=0)
    return y_pred 

def get_thick_slices(img_data, slice_thickness=3):
    """
    Function to extract thick slices from the image 
    (feed slice_thickness preceeding and suceeding slices to network, 
    label only middle one)
    :param np.ndarray img_data: 3D MRI image read in with nibabel 
    :param int slice_thickness: number of slices to stack on top and below slice of interest (default=3) 
    :return: 
    """
    h, w, d = img_data.shape
    img_data_pad = np.expand_dims(np.pad(img_data, ((0, 0), (0, 0), (slice_thickness, slice_thickness)), mode='edge'),
                                  axis=3)
    img_data_thick = np.ndarray((h, w, d, 0), dtype=np.uint8)
    
    for slice_idx in range(2 * slice_thickness + 1):
        img_data_thick = np.append(img_data_thick, img_data_pad[:, :, slice_idx:d + slice_idx, :], axis=3)

    return img_data_thick 

def BW_img(input, thresholding):
    if input.max() > thresholding:
        binary = input > thresholding
    else:
        binary = input > input.max() / 2.0

    label_image = label(binary)
    regions = regionprops(label_image)
    area_list = []
    for region in regions:
        area_list.append(region.area)
    if area_list:
        idx_max = np.argmax(area_list)
        binary[label_image != idx_max + 1] = 0
    return scipy.ndimage.binary_fill_holes(np.asarray(binary).astype(int)) 

def elastic_transform(image, label, alpha, sigma, random_state=None):
    """Elastic deformation of images as described in [Simard2003]_.
    .. [Simard2003] Simard, Steinkraus and Platt, "Best Practices for
       Convolutional Neural Networks applied to Visual Document Analysis", in
       Proc. of the International Conference on Document Analysis and
       Recognition, 2003.
    """
    seed = random.random()
    if seed > 0.5:
        assert len(image.shape) == 3

        if random_state is None:
            random_state = np.random.RandomState(None)

        shape = image.shape[0:2]
        dx = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma, mode="constant", cval=0) * alpha
        dy = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma, mode="constant", cval=0) * alpha

        x, y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]), indexing='ij')
        indices = np.reshape(x + dx, (-1, 1)), np.reshape(y + dy, (-1, 1))
        transformed_image = np.zeros(image.shape)
        transformed_label = np.zeros(image.shape)
        for i in range(image.shape[-1]):
            transformed_image[:, :, i] = map_coordinates(image[:, :, i], indices, order=1).reshape(shape)
            if label is not None:
                transformed_label[:, :, i] = map_coordinates(label[:, :, i], indices, order=1).reshape(shape)
            else:
                transformed_label = None
        transformed_image = transformed_image.astype(np.uint8)
        if label is not None:
            transformed_label = transformed_label.astype(np.uint8)
        return transformed_image, transformed_label
    else:
        return image, label 

def label_lesion(self):
        printv('\nLabel connected regions of the masked image...', self.verbose, 'normal')
        im = Image(self.fname_mask)
        im_2save = im.copy()
        im_2save.data = label(im.data, connectivity=2)
        im_2save.save(self.fname_label)

        self.measure_pd['label'] = [l for l in np.unique(im_2save.data) if l]
        printv('Lesion count = ' + str(len(self.measure_pd['label'])), self.verbose, 'info') 

def get_largest_fillhole(binary):
    label_image = label(binary)
    regions = regionprops(label_image)
    area_list = []
    for region in regions:
        area_list.append(region.area)
    if area_list:
        idx_max = np.argmax(area_list)
        binary[label_image != idx_max + 1] = 0
    return scipy.ndimage.binary_fill_holes(np.asarray(binary).astype(int)) 

def projectSegmentation(self, pointSegmentation, min_x, max_x, min_y, max_y):
        if max_x - min_x == 1 and max_y - min_y == 1:
            segments, counts = np.unique(pointSegmentation[:, 2], return_counts=True)
            segmentList = zip(segments.tolist(), counts.tolist())
            segmentList = [segment for segment in segmentList if segment[0] not in [0, 2]]
            label = 0
            if 2 in segments:
                label = 2
                pass
            if len(segmentList) > 0:
                segment = max(segmentList, key=lambda x: x[1])
                if segment[1] > 0:
                    label = segment[0]
                    pass
                pass
            self.segmentation[min_y][min_x] = label
        elif max_x - min_x >= max_y - min_y:
            middle_x = int((max_x + min_x + 1) / 2)
            mask_1 = pointSegmentation[:, 1] < middle_x
            self.projectSegmentation(pointSegmentation[mask_1], min_x, middle_x, min_y, max_y)
            mask_2 = pointSegmentation[:, 1] >= middle_x
            self.projectSegmentation(pointSegmentation[mask_2], middle_x, max_x, min_y, max_y)
        else:
            middle_y = int((max_y + min_y + 1) / 2)
            mask_1 = pointSegmentation[:, 0] < middle_y
            self.projectSegmentation(pointSegmentation[mask_1], min_x, max_x, min_y, middle_y)
            mask_2 = pointSegmentation[:, 0] >= middle_y
            self.projectSegmentation(pointSegmentation[mask_2], min_x, max_x, middle_y, max_y)
            pass
        return 

def extractCornersFromHeatmaps(heatmaps, heatmapThreshold=0.5, numPixelsThreshold=5, returnRanges=True):
    from skimage import measure
    heatmaps = (heatmaps > heatmapThreshold).astype(np.float32)
    orientationPoints = []
    #kernel = np.ones((3, 3), np.float32)
    for heatmapIndex in xrange(0, heatmaps.shape[-1]):
        heatmap = heatmaps[:, :, heatmapIndex]
        #heatmap = cv2.dilate(cv2.erode(heatmap, kernel), kernel)
        components = measure.label(heatmap, background=0)
        points = []
        for componentIndex in xrange(components.min() + 1, components.max() + 1):
            ys, xs = (components == componentIndex).nonzero()
            if ys.shape[0] <= numPixelsThreshold:
                continue
            #print(heatmapIndex, xs.shape, ys.shape, componentIndex)
            if returnRanges:
                points.append(((xs.mean(), ys.mean()), (xs.min(), ys.min()), (xs.max(), ys.max())))
            else:
                points.append((xs.mean(), ys.mean()))
                pass
            continue
        orientationPoints.append(points)
        continue
    return orientationPoints

#Extract corners from heatmaps 

def __init__(self, ignore_index=None):
        """
        :param ignore_index: id of the label to be ignored from IoU computation
        """
        self.ignore_index = ignore_index