Python skimage.morphology.skeletonize() Examples

def test():
    from skimage.morphology import skeletonize
    import numpy as np
    from skimage import data
    import matplotlib.pyplot as plt

    img = data.horse()
    ske = skeletonize(~img).astype(np.uint16)
    graph = build_sknw(ske)
    plt.imshow(img, cmap='gray')
    for (s,e) in graph.edges():
        ps = graph[s][e]['pts']
        plt.plot(ps[:,1], ps[:,0], 'green')

    nodes = graph.nodes()
    ps = np.array([nodes[i]['o'] for i in nodes])
    plt.plot(ps[:,1], ps[:,0], 'r.')
    plt.title('Build Graph')
    plt.show() 

def metroize(img, mgrid=32, threshold=0.5):
    threshold = translate_threshold(img, threshold=threshold)
    img = skeletonize(img > threshold)
    img = skeletonize(rebin(img, [mgrid, mgrid]) > 0)
    return img 

def process_single_image(filename, image_format, scale_metadata_path,
                         threshold_radius, smooth_radius,
                         brightness_offset, crop_radius, smooth_method):
    image = imageio.imread(filename, format=image_format)
    scale = _get_scale(image, scale_metadata_path)
    if crop_radius > 0:
        c = crop_radius
        image = image[c:-c, c:-c]
    pixel_threshold_radius = int(np.ceil(threshold_radius / scale))

    pixel_smoothing_radius = smooth_radius * pixel_threshold_radius
    thresholded = pre.threshold(image, sigma=pixel_smoothing_radius,
                                radius=pixel_threshold_radius,
                                offset=brightness_offset,
                                smooth_method=smooth_method)
    quality = shape_index(image, sigma=pixel_smoothing_radius,
                          mode='reflect')
    skeleton = morphology.skeletonize(thresholded) * quality
    framedata = csr.summarise(skeleton, spacing=scale)
    framedata['squiggle'] = np.log2(framedata['branch-distance'] /
                                    framedata['euclidean-distance'])
    framedata['scale'] = scale
    framedata.rename(columns={'mean pixel value': 'mean shape index'},
                     inplace=True)
    framedata['filename'] = filename
    return image, thresholded, skeleton, framedata 

def test_skeleton(test_thresholded):
    skeleton = morphology.skeletonize(test_thresholded)
    return skeleton 

def thin(guess_in):
    """
    Use scikit-image skeletonize to 'thin' the guess raster.

    Parameters
    ----------
    guess_in : path-like or 2D array
        Output from threshold().

    Returns
    -------
    guess_skel : numpy array
        Thinned version.
    affine : Affine
        Only if path-like supplied.
    """

    if isinstance(guess_in, (str, Path)):
        guess_rd = rasterio.open(guess_in)
        guess_arr = guess_rd.read(1)
        affine = guess_rd.transform

        guess_skel = skeletonize(guess_arr)
        guess_skel = guess_skel.astype("int32")

        return guess_skel, affine

    elif isinstance(guess_in, np.ndarray):
        guess_skel = skeletonize(guess_in)
        guess_skel = guess_skel.astype("int32")

        return guess_skel

    else:
        raise ValueError 

def compute_precision_recall(pred, gt):
    pred_skel = skeletonize(pred)
    pred_dil = dilation(pred_skel, square(5))
    gt_skel = skeletonize(gt)
    gt_dil = dilation(gt_skel, square(5))
    return compute_metrics([pred_skel], [gt_skel], [pred_dil], [gt_dil]) 

def scale(self, ratio, interpolation='LINEAR'):
        h, w = self.data.shape[:2]
        tw = int(ratio * w)
        th = int(ratio * h)

        # solve the missed edges
        if ratio > 1:
            im = cv2.resize(self.data, dsize=(tw, th), interpolation=cv2.INTER_LINEAR_EXACT)
            im[im > 0.2] = 1
            im = skeletonize(im)
        else:
            im = cv2.resize(self.data, dsize=(tw, th), interpolation=cv2.INTER_LINEAR_EXACT)
            im[im > 0.4] = 1
            im = skeletonize(im)
        self.data = im.copy() 

def place_routers_on_skeleton(d, cmethod):
    wireless = np.where(d["graph"] == Cell.Wireless, 1, 0)
    # perform skeletonization
    skeleton = skeletonize(wireless)
    med_axis = medial_axis(wireless)

    skel = skeleton
    # skel = med_axis
    # get all skeleton positions
    pos = []
    for i in range(skel.shape[0]):
        for j in range(skel.shape[1]):
            if skel[i][j]:
                pos.append((i, j))

    budget = d['budget']
    shuffle(pos)

    max_num_routers = min([int(d['budget'] / d['price_router']), len(pos)])
    print("Num of routers constrained by:")
    print(" budget:   %d" % int(int(d['budget'] / d['price_router'])))
    print(" skeleton: %d" % len(pos))

    for i in tqdm(range(max_num_routers), desc="Placing Routers"):
        new_router = pos[i]
        a, b = new_router

        # check if remaining budget is enough
        d["graph"][a][b] = Cell.Router
        d, ret, cost = _add_cabel(d, new_router, budget)
        budget -= cost

        if not ret:
            break

    return d 

def place_routers_on_skeleton_iterative(d, cmethod):
    budget = d['budget']
    R = d['radius']
    max_num_routers = int(d['budget'] / d['price_router'])
    coverage = np.where(d["graph"] == Cell.Wireless, 1, 0).astype(np.bool)

    pbar = tqdm(range(max_num_routers), desc="Placing Routers")
    while budget > 0:
        # perform skeletonization
        # skeleton = skeletonize(coverage)
        skeleton = medial_axis(coverage)
        # get all skeleton positions
        pos = np.argwhere(skeleton > 0).tolist()
        # escape if no positions left
        if not len(pos):
            break
        # get a random position
        shuffle(pos)
        a, b = pos[0]
        # place router
        d["graph"][a][b] = Cell.Router
        d, ret, cost = _add_cabel(d, (a, b), budget)
        if not ret:
            print("No budget available!")
            break
        budget -= cost
        # refresh wireless map by removing new coverage
        m = wireless_access(a, b, R, d['graph']).astype(np.bool)
        coverage[(a - R):(a + R + 1), (b - R):(b + R + 1)] &= ~m
        pbar.update()
    pbar.close()

    return d 

def skeletonize(mask):
    """Reduces binary objects to 1 pixel wide representations (skeleton)

    Inputs:
    mask       = Binary image data

    Returns:
    skeleton   = skeleton image

    :param mask: numpy.ndarray
    :return skeleton: numpy.ndarray
    """

    # Convert mask to boolean image, rather than 0 and 255 for skimage to use it
    skeleton = skmorph.skeletonize(mask.astype(bool))

    skeleton = skeleton.astype(np.uint8) * 255

    # Auto-increment device
    params.device += 1

    if params.debug == 'print':
        print_image(skeleton, os.path.join(params.debug_outdir, str(params.device) + '_skeleton.png'))
    elif params.debug == 'plot':
        plot_image(skeleton, cmap='gray')

    return skeleton 

def breakup_region(component):
    distance = ndi.distance_transform_edt(component)
    skel = skeletonize(component)
    skeldist = distance*skel
    local_maxi = peak_local_max(skeldist, indices=False, footprint=disk(10))
    local_maxi=ndi.binary_closing(local_maxi,structure = disk(4),iterations = 2)
    markers = ndi.label(local_maxi)[0]
    labels = watershed(-distance, markers, mask=component)
    return(labels) 

def skeleton_transform(label, relabel=True):
    resolution = (1.0, 1.0)
    alpha = 1.0
    beta = 0.8

    if relabel == True: # run connected component
        label = morphology.label(label, background=0)
    
    label_id = np.unique(label)
    # print(np.unique(label_id))
    skeleton = np.zeros(label.shape, dtype=np.uint8)
    distance = np.zeros(label.shape, dtype=np.float32)
    Temp = np.zeros(label.shape, dtype=np.uint8)
    if len(label_id) == 1: # only one object within current volume
        if label_id[0] == 0:
            return distance, skeleton
        else:
            temp_id = label_id
    else:
        temp_id = label_id[1:]

    for idx in temp_id:
        temp1 = (label == idx)
        temp2 = morphology.remove_small_holes(temp1, 16, connectivity=1)

        #temp3 = erosion(temp2)
        temp3 = temp2.copy()
        Temp += temp3
        skeleton_mask = skeletonize(temp3).astype(np.uint8)
        skeleton += skeleton_mask

        skeleton_edt = ndimage.distance_transform_edt(
            1-skeleton_mask, resolution)
        dist_max = np.max(skeleton_edt*temp3)
        dist_max = np.clip(dist_max, a_min=2.0, a_max=None)
        skeleton_edt = skeleton_edt / (dist_max*alpha)
        skeleton_edt = skeleton_edt**(beta)

        reverse = 1.0-(skeleton_edt*temp3)
        distance += reverse*temp3

    # generate boundary
    distance[np.where(Temp == 0)] = -1.0

    return distance, skeleton 

def find_lines(lines_mask: np.ndarray) -> list:
    """
    Finds the longest central line for each connected component in the given binary mask.

    :param lines_mask: Binary mask of the detected line-areas
    :return: a list of Opencv-style polygonal lines (each contour encoded as [N,1,2] elements where each tuple is (x,y) )
    """
    # Make sure one-pixel wide 8-connected mask
    lines_mask = skeletonize(lines_mask)

    class MakeLineMCP(MCP_Connect):
        def __init__(self, *args, **kwargs):
            super().__init__(*args, **kwargs)
            self.connections = dict()
            self.scores = defaultdict(lambda: np.inf)

        def create_connection(self, id1, id2, pos1, pos2, cost1, cost2):
            k = (min(id1, id2), max(id1, id2))
            s = cost1 + cost2
            if self.scores[k] > s:
                self.connections[k] = (pos1, pos2, s)
                self.scores[k] = s

        def get_connections(self, subsample=5):
            results = dict()
            for k, (pos1, pos2, s) in self.connections.items():
                path = np.concatenate([self.traceback(pos1), self.traceback(pos2)[::-1]])
                results[k] = path[::subsample]
            return results

        def goal_reached(self, int_index, float_cumcost):
            if float_cumcost > 0:
                return 2
            else:
                return 0

    if np.sum(lines_mask) == 0:
        return []
    # Find extremities points
    end_points_candidates = np.stack(np.where((convolve2d(lines_mask, np.ones((3, 3)), mode='same') == 2) & lines_mask)).T
    connected_components = skimage_label(lines_mask, connectivity=2)
    # Group endpoint by connected components and keep only the two points furthest away
    d = defaultdict(list)
    for pt in end_points_candidates:
        d[connected_components[pt[0], pt[1]]].append(pt)
    end_points = []
    for pts in d.values():
        d = euclidean_distances(np.stack(pts), np.stack(pts))
        i, j = np.unravel_index(d.argmax(), d.shape)
        end_points.append(pts[i])
        end_points.append(pts[j])
    end_points = np.stack(end_points)

    mcp = MakeLineMCP(~lines_mask)
    mcp.find_costs(end_points)
    connections = mcp.get_connections()
    if not np.all(np.array(sorted([i for k in connections.keys() for i in k])) == np.arange(len(end_points))):
        print('Warning : find_lines seems weird')
    return [c[:, None, ::-1] for c in connections.values()]