cytometer/utils.py
def segment_dmap_contour(dmap, contour=None,
                         sigma=10, min_seed_object_size=50, border_dilation=0, boundary_threshold=0.1,
                         median_size=11, closing_size=11, contour_seed_threshold=0,
                         version=2):
    """
    Segment cells from a distance transformation image, and optionally, a contour estimate image.

    This function computes the normal curvature of the dmap seen as a Monge patch. The "valleys" in
    the dmap (the cell contours) correspond to higher normal curvature values.

    If provided, the normal curvature is element-wise multiplied by a prior estimate of where the
    cell contours are. This prior estimate can be computed e.g. with a convolutional neural
    network classifier.

    Then a threshold <= 0 is applied to the curvatures (or weighted curvatures). Those pixels are
    considered to belong to the inside of cells. This usually produces separate "islands", one per
    cell, that don't extend all the way to the cells edges. "Islands" with less that
    min_seed_object_size pixels are removed as noise.

    Next, a watershed algorithm is initialised with the "islands", and run on the (weighted) mean
    curvature image to extend the segmentations all the way to the cell borders. This produces
    one different label per cell.

    The borders of each cell are then extracted, and optionally, dilated with a 3x3 kernel of ones,
    "border_dilation" times.

    :param dmap: numpy.ndarray matrix with distance transformation, distance range image,
    topographic map, etc.
    :param contour: (def None) numpy.ndarray matrix, same size as dmap, with an estimate of where the
    contours are.
    :param sigma: (def 10) Standard deviation in pixels of Gaussian blurring applied to the dmap
    before computing the normal curvature. Because normal curvature is a function of second
    derivatives, slight noise in the dmap gets amplified, and produces poor curvature estimates.
    Gaussian blurring removes part of that noise.
    :param min_seed_object_size: (def 50). Objects with fewer pixels than this value will be discarded.
    :param border_dilation: (def 0) Optional dilation of the watershed borders.
    :param boundary_threshold: (def 0.1) Threshold for the 90-percentile of contour values along kissing points
    between adjacent labels. If the 90-percentile < boundary_threshold, we assume that there's no membrane between
    both labels, and thus they belong to the same cell and should be merged.
    :param median_size: (def 11) Kernel size for the median filter to reduce noise in the output of the contour CNN.
    Larger values may blur contours too much, so there won't be seeds inside of the cells. Smaller values may leave
    fragmented seeds, that split cells into two or more labels.
    :param closing_size: (def 11) Kernel size for the dilation and erosion operation that reduces noise in the
    segmentation seeds (ideally, you get a connected component seed per cell). Larger values may merge neighbour cells
    together, whereas a low value may fragment a cell into two or more labels.
    :param contour_seed_threshold: (def 0) The output of the contour CNN are contour values. Segmentations seeds inside
    of cells are computed as contour <= contour_seed_threshold. Usually, 0 works well. This value can be
    increased a bit if some cells fail to get a seed. However, contour values on membranes are often as low as 0.1.
    Thus, even small values of contour_seed_threshold can create seeds that overflow the cell.
    :param version: (def 2) Implementation version of the segmentation algorithm. version=1 corresponds to a deprecated
    algorithm, kept only for historical comparisons.
    :return: labels, labels_borders
    """

    # auxiliary functions for merge_hierarchical() taken from
    # https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_rag_merge.html#sphx-glr-auto-examples-segmentation-plot-rag-merge-py
    #
    # For the purpose of this function, the details of what they do are not very important
    def _weight_mean_color(graph, src, dst, n):
        diff = graph.node[dst]['mean color'] - graph.node[n]['mean color']
        diff = np.linalg.norm(diff)
        return {'weight': diff}

    def merge_mean_color(graph, src, dst):
        graph.node[dst]['total color'] += graph.node[src]['total color']
        graph.node[dst]['pixel count'] += graph.node[src]['pixel count']
        graph.node[dst]['mean color'] = (graph.node[dst]['total color'] / graph.node[dst]['pixel count'])

    # check size of inputs
    if dmap.ndim != 2:
        raise ValueError('dmap array must have 2 dimensions')
    if contour is not None and contour.shape != dmap.shape:
        raise ValueError('if provided, contour must have the same shape as dmap')

    if version == 1:

        # compute mean curvature from dmap
        _, mean_curvature, _, _ = principal_curvatures_range_image(dmap, sigma=sigma)

        # multiply mean curvature by estimated contours, and clip negative values
        if contour is not None:
            contour *= mean_curvature
        else:
            contour = mean_curvature.copy()
        contour[contour < 0] = 0

        # rough segmentation of inner areas
        labels = (mean_curvature <= 0).astype('uint8')

        # label areas with a different label per connected area
        labels = measure.label(labels)

        # remove very small labels (noise)
        labels_prop = measure.regionprops(labels)
        for j in range(1, np.max(labels)):
            # label of region under consideration is not the same as index j
            lab = labels_prop[j]['label']
            if labels_prop[j]['area'] < min_seed_object_size:
                labels[labels == lab] = 0

        # extend labels using watershed
        labels = watershed(mean_curvature, labels)

    elif version == 2:

        if DEBUG:
            plt.clf()
            plt.subplot(221)
            plt.imshow(contour)

        # remove noise from contour image
        # Note: scipy.signal.medfilt is 9x slower for (21, 21) kernel compared to scipy.ndimage.median_filter
        if median_size > 1:
            contour = median_filter(contour, size=(median_size, median_size))

        if DEBUG:
            plt.subplot(222)
            plt.imshow(contour)

        # threshold contours to find seeds (inside areas of cells)
        seg = (contour <= contour_seed_threshold).astype(np.uint8)

        if DEBUG:
            plt.subplot(223)
            plt.imshow(seg)

        # dilate and erode to reduce the segmentation noise
        seg = binary_closing(seg.astype(np.bool), selem=np.ones(shape=(closing_size, closing_size))).astype(np.uint8)

        if DEBUG:
            plt.subplot(224)
            plt.imshow(seg)

        # remove small holes from the segmentation of insides of cells
        seg = remove_small_holes(seg.astype(np.bool), area_threshold=10e3).astype(np.uint8)

        if DEBUG:
            plt.subplot(223)
            plt.imshow(seg)

        # assign different label to each connected components
        # Note: cv2.connectedComponentsWithStats is 10x faster than skimage.measure.label
        nlabels, labels, stats, centroids = cv2.connectedComponentsWithStats(seg)

        if DEBUG:
            plt.subplot(224)
            plt.cla()
            plt.imshow(labels)

        # use watershed to expand the seeds
        labels = watershed(-dmap, labels, watershed_line=False)

        if DEBUG:
            plt.subplot(224)
            plt.imshow(labels)

            plt.subplot(222)
            plt.imshow(contour)
            plt.contour(labels, colors='r', levels=range(np.max(labels)))

        # For every pair of adjacent labels, we want to find out whether both belong to different cells (i.e.
        # they have a membrane between them), or they belong to the same cell and should be merged.
        #
        # To do this efficiently, we are going to use a Region Adjacency Graph (RAG) that will create a graph
        # where each node corresponds to one label, and edges indicate that two labels are adjacent.
        #
        # The edge weight = norm([la, la, la] - [lb, lb, lb]), where la, lb are the two labels. This is how it's
        # implemented in rag_mean_color(). We don't care about the value of the weight as long as it's > 0.
        #
        # Then, we are going to check pairs of adjacent labels. If they should be merge, we make the edge weight = 0.
        # The criterion to decide whether labels should be merged is to check whether on the kissing points between
        # both labels we have a membrane or not.
        #
        # Once all edges have been examined, the graph can be simplified in one go with merge_hierarchical().

        if len(np.unique(labels)) > 1:

            # compute the Region Adjacency Graph using mean colours (the weights from the mean colours will be ignored)
            rag = rag_mean_color(labels, labels, mode='distance')

            # iterate pairs of adjacent labels (note that if rag.edges has (1,2), it doesn't also have (2,1))
            for lab_a, lab_b in rag.edges:

                # dilate each label so that they overlap wherever they are adjacent
                labels_a = (labels == lab_a).astype(np.uint8)
                labels_b = (labels == lab_b).astype(np.uint8)
                labels_a = cv2.dilate(labels_a, kernel=np.ones(shape=(3, 3), dtype=np.uint8))
                labels_b = cv2.dilate(labels_b, kernel=np.ones(shape=(3, 3), dtype=np.uint8))

                # boundary pixels that are common to both labels (kissing points)
                common = np.logical_and(labels_a, labels_b)

                # if the 90-percentile of the contour values along the kissing points is small, that means that the
                # separation between the two labels is spurious, because there's probably no membrane between them
                if np.percentile(contour[common], 90) <= boundary_threshold:
                    if DEBUG:
                        print('Merging ' + str((lab_a, lab_b)))
                    nx.set_edge_attributes(rag, {(lab_a, lab_b): {'weight': 0}})

            # update the labels by greedy merging
            labels = merge_hierarchical(labels, rag, thresh=0.1, rag_copy=False, in_place_merge=True,
                                        merge_func=merge_mean_color, weight_func=_weight_mean_color)

        if DEBUG:
            plt.subplot(224)
            plt.contour(labels, levels=np.unique(labels), colors='r')

    else:

        raise ValueError('You are asking for a version of this function that does not exist')

    # extract borders of watershed regions for plots
    labels_borders = borders(labels)

    # dilate borders for easier visualization
    if border_dilation > 0:
        kernel = np.ones((3, 3), np.uint8)
        labels_borders = cv2.dilate(labels_borders.astype(np.uint8), kernel=kernel, iterations=border_dilation) > 0

    return labels, labels_borders

    Segment cells from a distance transformation image, and optionally, a contour estimate image.

    This function computes the normal curvature of the dmap seen as a Monge patch. The "valleys" in
    the dmap (the cell contours) correspond to higher normal curvature values.

    If provided, the normal curvature is element-wise multiplied by a prior estimate of where the
    cell contours are. This prior estimate can be computed e.g. with a convolutional neural
    network classifier.

    Then a threshold <= 0 is applied to the curvatures (or weighted curvatures). Those pixels are
    considered to belong to the inside of cells. This usually produces separate "islands", one per
    cell, that don't extend all the way to the cells edges. "Islands" with less that
    min_seed_object_size pixels are removed as noise.

    Next, a watershed algorithm is initialised with the "islands", and run on the (weighted) mean
    curvature image to extend the segmentations all the way to the cell borders. This produces
    one different label per cell.

    The borders of each cell are then extracted, and optionally, dilated with a 3x3 kernel of ones,
    "border_dilation" times.

    :param dmap: numpy.ndarray matrix with distance transformation, distance range image,
    topographic map, etc.
    :param contour: (def None) numpy.ndarray matrix, same size as dmap, with an estimate of where the
    contours are.
    :param sigma: (def 10) Standard deviation in pixels of Gaussian blurring applied to the dmap
    before computing the normal curvature. Because normal curvature is a function of second
    derivatives, slight noise in the dmap gets amplified, and produces poor curvature estimates.
    Gaussian blurring removes part of that noise.
    :param min_seed_object_size: (def 50). Objects with fewer pixels than this value will be discarded.
    :param border_dilation: (def 0) Optional dilation of the watershed borders.
    :param boundary_threshold: (def 0.1) Threshold for the 90-percentile of contour values along kissing points
    between adjacent labels. If the 90-percentile < boundary_threshold, we assume that there's no membrane between
    both labels, and thus they belong to the same cell and should be merged.
    :param median_size: (def 11) Kernel size for the median filter to reduce noise in the output of the contour CNN.
    Larger values may blur contours too much, so there won't be seeds inside of the cells. Smaller values may leave
    fragmented seeds, that split cells into two or more labels.
    :param closing_size: (def 11) Kernel size for the dilation and erosion operation that reduces noise in the
    segmentation seeds (ideally, you get a connected component seed per cell). Larger values may merge neighbour cells
    together, whereas a low value may fragment a cell into two or more labels.
    :param contour_seed_threshold: (def 0) The output of the contour CNN are contour values. Segmentations seeds inside
    of cells are computed as contour <= contour_seed_threshold. Usually, 0 works well. This value can be
    increased a bit if some cells fail to get a seed. However, contour values on membranes are often as low as 0.1.
    Thus, even small values of contour_seed_threshold can create seeds that overflow the cell.
    :param version: (def 2) Implementation version of the segmentation algorithm. version=1 corresponds to a deprecated
    algorithm, kept only for historical comparisons.
    :return: labels, labels_borders
    
Segment cells from a distance transformation image, and optionally, a contour estimate image.
This function computes the normal curvature of the dmap seen as a Monge patch. The "valleys" in
the dmap (the cell contours) correspond to higher normal curvature values.

If provided, the normal curvature is element-wise multiplied by a prior estimate of where the
cell contours are. This prior estimate can be computed e.g. with a convolutional neural
network classifier.

Then a threshold <= 0 is applied to the curvatures (or weighted curvatures). Those pixels are
considered to belong to the inside of cells. This usually produces separate "islands", one per
cell, that don't extend all the way to the cells edges. "Islands" with less that
min_seed_object_size pixels are removed as noise.

Next, a watershed algorithm is initialised with the "islands", and run on the (weighted) mean
curvature image to extend the segmentations all the way to the cell borders. This produces
one different label per cell.

The borders of each cell are then extracted, and optionally, dilated with a 3x3 kernel of ones,
"border_dilation" times.
Segment cells from a distance transformation image and optionally a contour estimate image . This function computes the normal curvature of the dmap seen as a Monge patch . The " valleys " in the dmap ( the cell contours ) correspond to higher normal curvature values . If provided the normal curvature is element - wise multiplied by a prior estimate of where the cell contours are . This prior estimate can be computed e . g . with a convolutional neural network classifier . Then a threshold < = 0 is applied to the curvatures ( or weighted curvatures ) . Those pixels are considered to belong to the inside of cells . This usually produces separate " islands " one per cell that don ' t extend all the way to the cells edges . " Islands " with less that min_seed_object_size pixels are removed as noise . Next a watershed algorithm is initialised with the " islands " and run on the ( weighted ) mean curvature image to extend the segmentations all the way to the cell borders . This produces one different label per cell . The borders of each cell are then extracted and optionally dilated with a 3x3 kernel of ones " border_dilation " times .
def segment_dmap_contour(dmap, contour=None,
                         sigma=10, min_seed_object_size=50, border_dilation=0, boundary_threshold=0.1,
                         median_size=11, closing_size=11, contour_seed_threshold=0,
                         version=2):
     For the purpose of this function, the details of what they do are not very important
    def _weight_mean_color(graph, src, dst, n):
        diff = graph.node[dst]['mean color'] - graph.node[n]['mean color']
        diff = np.linalg.norm(diff)
        return {'weight': diff}
    def merge_mean_color(graph, src, dst):
        graph.node[dst]['total color'] += graph.node[src]['total color']
        graph.node[dst]['pixel count'] += graph.node[src]['pixel count']
        graph.node[dst]['mean color'] = (graph.node[dst]['total color'] / graph.node[dst]['pixel count'])
     check size of inputs
    if dmap.ndim != 2:
        raise ValueError('dmap array must have 2 dimensions')
    if contour is not None and contour.shape != dmap.shape:
        raise ValueError('if provided, contour must have the same shape as dmap')
    if version == 1:
         compute mean curvature from dmap
        _, mean_curvature, _, _ = principal_curvatures_range_image(dmap, sigma=sigma)
         multiply mean curvature by estimated contours, and clip negative values
        if contour is not None:
            contour *= mean_curvature
        else:
            contour = mean_curvature.copy()
        contour[contour < 0] = 0
         rough segmentation of inner areas
        labels = (mean_curvature <= 0).astype('uint8')
         label areas with a different label per connected area
        labels = measure.label(labels)
         remove very small labels (noise)
        labels_prop = measure.regionprops(labels)
        for j in range(1, np.max(labels)):
             label of region under consideration is not the same as index j
            lab = labels_prop[j]['label']
            if labels_prop[j]['area'] < min_seed_object_size:
                labels[labels == lab] = 0
         extend labels using watershed
        labels = watershed(mean_curvature, labels)
    elif version == 2:
        if DEBUG:
            plt.clf()
            plt.subplot(221)
            plt.imshow(contour)
         remove noise from contour image
         Note: scipy.signal.medfilt is 9x slower for (21, 21) kernel compared to scipy.ndimage.median_filter
        if median_size > 1:
            contour = median_filter(contour, size=(median_size, median_size))
        if DEBUG:
            plt.subplot(222)
            plt.imshow(contour)
         threshold contours to find seeds (inside areas of cells)
        seg = (contour <= contour_seed_threshold).astype(np.uint8)
        if DEBUG:
            plt.subplot(223)
            plt.imshow(seg)
         dilate and erode to reduce the segmentation noise
        seg = binary_closing(seg.astype(np.bool), selem=np.ones(shape=(closing_size, closing_size))).astype(np.uint8)
        if DEBUG:
            plt.subplot(224)
            plt.imshow(seg)
         remove small holes from the segmentation of insides of cells
        seg = remove_small_holes(seg.astype(np.bool), area_threshold=10e3).astype(np.uint8)
        if DEBUG:
            plt.subplot(223)
            plt.imshow(seg)
         assign different label to each connected components
         Note: cv2.connectedComponentsWithStats is 10x faster than skimage.measure.label
        nlabels, labels, stats, centroids = cv2.connectedComponentsWithStats(seg)
        if DEBUG:
            plt.subplot(224)
            plt.cla()
            plt.imshow(labels)
         use watershed to expand the seeds
        labels = watershed(-dmap, labels, watershed_line=False)
        if DEBUG:
            plt.subplot(224)
            plt.imshow(labels)
            plt.subplot(222)
            plt.imshow(contour)
            plt.contour(labels, colors='r', levels=range(np.max(labels)))
         For every pair of adjacent labels, we want to find out whether both belong to different cells (i.e.
         they have a membrane between them), or they belong to the same cell and should be merged.
         To do this efficiently, we are going to use a Region Adjacency Graph (RAG) that will create a graph
         where each node corresponds to one label, and edges indicate that two labels are adjacent.
         The edge weight = norm([la, la, la] - [lb, lb, lb]), where la, lb are the two labels. This is how it's
         implemented in rag_mean_color(). We don't care about the value of the weight as long as it's > 0.
         Then, we are going to check pairs of adjacent labels. If they should be merge, we make the edge weight = 0.
         The criterion to decide whether labels should be merged is to check whether on the kissing points between
         both labels we have a membrane or not.
         Once all edges have been examined, the graph can be simplified in one go with merge_hierarchical().
        if len(np.unique(labels)) > 1:
             compute the Region Adjacency Graph using mean colours (the weights from the mean colours will be ignored)
            rag = rag_mean_color(labels, labels, mode='distance')
             iterate pairs of adjacent labels (note that if rag.edges has (1,2), it doesn't also have (2,1))
            for lab_a, lab_b in rag.edges:
                 dilate each label so that they overlap wherever they are adjacent
                labels_a = (labels == lab_a).astype(np.uint8)
                labels_b = (labels == lab_b).astype(np.uint8)
                labels_a = cv2.dilate(labels_a, kernel=np.ones(shape=(3, 3), dtype=np.uint8))
                labels_b = cv2.dilate(labels_b, kernel=np.ones(shape=(3, 3), dtype=np.uint8))
                 boundary pixels that are common to both labels (kissing points)
                common = np.logical_and(labels_a, labels_b)
                 if the 90-percentile of the contour values along the kissing points is small, that means that the
                 separation between the two labels is spurious, because there's probably no membrane between them
                if np.percentile(contour[common], 90) <= boundary_threshold:
                    if DEBUG:
                        print('Merging ' + str((lab_a, lab_b)))
                    nx.set_edge_attributes(rag, {(lab_a, lab_b): {'weight': 0}})
             update the labels by greedy merging
            labels = merge_hierarchical(labels, rag, thresh=0.1, rag_copy=False, in_place_merge=True,
                                        merge_func=merge_mean_color, weight_func=_weight_mean_color)
        if DEBUG:
            plt.subplot(224)
            plt.contour(labels, levels=np.unique(labels), colors='r')
    else:
        raise ValueError('You are asking for a version of this function that does not exist')
     extract borders of watershed regions for plots
    labels_borders = borders(labels)
     dilate borders for easier visualization
    if border_dilation > 0:
        kernel = np.ones((3, 3), np.uint8)
        labels_borders = cv2.dilate(labels_borders.astype(np.uint8), kernel=kernel, iterations=border_dilation) > 0
    return labels, labels_borders
def segment_dmap_contour ( dmap , contour = None , sigma = 10 , min_seed_object_size = 50 , border_dilation = 0 , boundary_threshold = 0.1 , median_size = 11 , closing_size = 11 , contour_seed_threshold = 0 , version = 2 ) : def _weight_mean_color ( graph , src , dst , n ) : diff = graph . node [ dst ] [ 'mean color' ] - graph . node [ n ] [ 'mean color' ] diff = np . linalg . norm ( diff ) return { 'weight' : diff } def merge_mean_color ( graph , src , dst ) : graph . node [ dst ] [ 'total color' ] += graph . node [ src ] [ 'total color' ] graph . node [ dst ] [ 'pixel count' ] += graph . node [ src ] [ 'pixel count' ] graph . node [ dst ] [ 'mean color' ] = ( graph . node [ dst ] [ 'total color' ] / graph . node [ dst ] [ 'pixel count' ] ) if dmap . ndim != 2 : raise ValueError ( 'dmap array must have 2 dimensions' ) if contour is not None and contour . shape != dmap . shape : raise ValueError ( 'if provided, contour must have the same shape as dmap' ) if version == 1 : _ , mean_curvature , _ , _ = principal_curvatures_range_image ( dmap , sigma = sigma ) if contour is not None : contour *= mean_curvature else : contour = mean_curvature . copy ( ) contour [ contour < 0 ] = 0 labels = ( mean_curvature <= 0 ) . astype ( 'uint8' ) labels = measure . label ( labels ) labels_prop = measure . regionprops ( labels ) for j in range ( 1 , np . max ( labels ) ) : lab = labels_prop [ j ] [ 'label' ] if labels_prop [ j ] [ 'area' ] < min_seed_object_size : labels [ labels == lab ] = 0 labels = watershed ( mean_curvature , labels ) elif version == 2 : if DEBUG : plt . clf ( ) plt . subplot ( 221 ) plt . imshow ( contour ) if median_size > 1 : contour = median_filter ( contour , size = ( median_size , median_size ) ) if DEBUG : plt . subplot ( 222 ) plt . imshow ( contour ) seg = ( contour <= contour_seed_threshold ) . astype ( np . uint8 ) if DEBUG : plt . subplot ( 223 ) plt . imshow ( seg ) seg = binary_closing ( seg . astype ( np . bool ) , selem = np . ones ( shape = ( closing_size , closing_size ) ) ) . astype ( np . uint8 ) if DEBUG : plt . subplot ( 224 ) plt . imshow ( seg ) seg = remove_small_holes ( seg . astype ( np . bool ) , area_threshold = 10e3 ) . astype ( np . uint8 ) if DEBUG : plt . subplot ( 223 ) plt . imshow ( seg ) nlabels , labels , stats , centroids = cv2 . connectedComponentsWithStats ( seg ) if DEBUG : plt . subplot ( 224 ) plt . cla ( ) plt . imshow ( labels ) labels = watershed ( - dmap , labels , watershed_line = False ) if DEBUG : plt . subplot ( 224 ) plt . imshow ( labels ) plt . subplot ( 222 ) plt . imshow ( contour ) plt . contour ( labels , colors = 'r' , levels = range ( np . max ( labels ) ) ) if len ( np . unique ( labels ) ) > 1 : rag = rag_mean_color ( labels , labels , mode = 'distance' ) for lab_a , lab_b in rag . edges : labels_a = ( labels == lab_a ) . astype ( np . uint8 ) labels_b = ( labels == lab_b ) . astype ( np . uint8 ) labels_a = cv2 . dilate ( labels_a , kernel = np . ones ( shape = ( 3 , 3 ) , dtype = np . uint8 ) ) labels_b = cv2 . dilate ( labels_b , kernel = np . ones ( shape = ( 3 , 3 ) , dtype = np . uint8 ) ) common = np . logical_and ( labels_a , labels_b ) if np . percentile ( contour [ common ] , 90 ) <= boundary_threshold : if DEBUG : print ( 'Merging ' + str ( ( lab_a , lab_b ) ) ) nx . set_edge_attributes ( rag , { ( lab_a , lab_b ) : { 'weight' : 0 } } ) labels = merge_hierarchical ( labels , rag , thresh = 0.1 , rag_copy = False , in_place_merge = True , merge_func = merge_mean_color , weight_func = _weight_mean_color ) if DEBUG : plt . subplot ( 224 ) plt . contour ( labels , levels = np . unique ( labels ) , colors = 'r' ) else : raise ValueError ( 'You are asking for a version of this function that does not exist' ) labels_borders = borders ( labels ) if border_dilation > 0 : kernel = np . ones ( ( 3 , 3 ) , np . uint8 ) labels_borders = cv2 . dilate ( labels_borders . astype ( np . uint8 ) , kernel = kernel , iterations = border_dilation ) > 0 return labels , labels_borders
Segment cells from a distance transformation image, and optionally, a contour estimate image.
Segment cells from a distance transformation image and optionally a contour estimate image .
"""
    Segment cells from a distance transformation image, and optionally, a contour estimate image.

    This function computes the normal curvature of the dmap seen as a Monge patch. The "valleys" in
    the dmap (the cell contours) correspond to higher normal curvature values.

    If provided, the normal curvature is element-wise multiplied by a prior estimate of where the
    cell contours are. This prior estimate can be computed e.g. with a convolutional neural
    network classifier.

    Then a threshold <= 0 is applied to the curvatures (or weighted curvatures). Those pixels are
    considered to belong to the inside of cells. This usually produces separate "islands", one per
    cell, that don't extend all the way to the cells edges. "Islands" with less that
    min_seed_object_size pixels are removed as noise.

    Next, a watershed algorithm is initialised with the "islands", and run on the (weighted) mean
    curvature image to extend the segmentations all the way to the cell borders. This produces
    one different label per cell.

    The borders of each cell are then extracted, and optionally, dilated with a 3x3 kernel of ones,
    "border_dilation" times.

    :param dmap: numpy.ndarray matrix with distance transformation, distance range image,
    topographic map, etc.
    :param contour: (def None) numpy.ndarray matrix, same size as dmap, with an estimate of where the
    contours are.
    :param sigma: (def 10) Standard deviation in pixels of Gaussian blurring applied to the dmap
    before computing the normal curvature. Because normal curvature is a function of second
    derivatives, slight noise in the dmap gets amplified, and produces poor curvature estimates.
    Gaussian blurring removes part of that noise.
    :param min_seed_object_size: (def 50). Objects with fewer pixels than this value will be discarded.
    :param border_dilation: (def 0) Optional dilation of the watershed borders.
    :param boundary_threshold: (def 0.1) Threshold for the 90-percentile of contour values along kissing points
    between adjacent labels. If the 90-percentile < boundary_threshold, we assume that there's no membrane between
    both labels, and thus they belong to the same cell and should be merged.
    :param median_size: (def 11) Kernel size for the median filter to reduce noise in the output of the contour CNN.
    Larger values may blur contours too much, so there won't be seeds inside of the cells. Smaller values may leave
    fragmented seeds, that split cells into two or more labels.
    :param closing_size: (def 11) Kernel size for the dilation and erosion operation that reduces noise in the
    segmentation seeds (ideally, you get a connected component seed per cell). Larger values may merge neighbour cells
    together, whereas a low value may fragment a cell into two or more labels.
    :param contour_seed_threshold: (def 0) The output of the contour CNN are contour values. Segmentations seeds inside
    of cells are computed as contour <= contour_seed_threshold. Usually, 0 works well. This value can be
    increased a bit if some cells fail to get a seed. However, contour values on membranes are often as low as 0.1.
    Thus, even small values of contour_seed_threshold can create seeds that overflow the cell.
    :param version: (def 2) Implementation version of the segmentation algorithm. version=1 corresponds to a deprecated
    algorithm, kept only for historical comparisons.
    :return: labels, labels_borders
    """ # auxiliary functions for merge_hierarchical() taken from # https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_rag_merge.html#sphx-glr-auto-examples-segmentation-plot-rag-merge-py # # For the purpose of this function, the details of what they do are not very important # check size of inputs # compute mean curvature from dmap # multiply mean curvature by estimated contours, and clip negative values # rough segmentation of inner areas # label areas with a different label per connected area # remove very small labels (noise) # label of region under consideration is not the same as index j # extend labels using watershed # remove noise from contour image # Note: scipy.signal.medfilt is 9x slower for (21, 21) kernel compared to scipy.ndimage.median_filter # threshold contours to find seeds (inside areas of cells) # dilate and erode to reduce the segmentation noise # remove small holes from the segmentation of insides of cells # assign different label to each connected components # Note: cv2.connectedComponentsWithStats is 10x faster than skimage.measure.label # use watershed to expand the seeds # For every pair of adjacent labels, we want to find out whether both belong to different cells (i.e. # they have a membrane between them), or they belong to the same cell and should be merged. # # To do this efficiently, we are going to use a Region Adjacency Graph (RAG) that will create a graph # where each node corresponds to one label, and edges indicate that two labels are adjacent. # # The edge weight = norm([la, la, la] - [lb, lb, lb]), where la, lb are the two labels. This is how it's # implemented in rag_mean_color(). We don't care about the value of the weight as long as it's > 0. # # Then, we are going to check pairs of adjacent labels. If they should be merge, we make the edge weight = 0. # The criterion to decide whether labels should be merged is to check whether on the kissing points between # both labels we have a membrane or not. # # Once all edges have been examined, the graph can be simplified in one go with merge_hierarchical(). # compute the Region Adjacency Graph using mean colours (the weights from the mean colours will be ignored) # iterate pairs of adjacent labels (note that if rag.edges has (1,2), it doesn't also have (2,1)) # dilate each label so that they overlap wherever they are adjacent # boundary pixels that are common to both labels (kissing points) # if the 90-percentile of the contour values along the kissing points is small, that means that the # separation between the two labels is spurious, because there's probably no membrane between them # update the labels by greedy merging # extract borders of watershed regions for plots # dilate borders for easier visualization
None
dmap
numpy.ndarray matrix with distance transformation, distance range image,
topographic map, etc.
numpy . ndarray matrix with distance transformation distance range image topographic map etc .
contour
(def None) numpy.ndarray matrix, same size as dmap, with an estimate of where the
contours are.
( def None ) numpy . ndarray matrix same size as dmap with an estimate of where the contours are .
sigma
(def 10) Standard deviation in pixels of Gaussian blurring applied to the dmap
before computing the normal curvature. Because normal curvature is a function of second
derivatives, slight noise in the dmap gets amplified, and produces poor curvature estimates.
Gaussian blurring removes part of that noise.
( def 10 ) Standard deviation in pixels of Gaussian blurring applied to the dmap before computing the normal curvature . Because normal curvature is a function of second derivatives slight noise in the dmap gets amplified and produces poor curvature estimates . Gaussian blurring removes part of that noise .
min_seed_object_size
None
border_dilation
(def 0) Optional dilation of the watershed borders.
( def 0 ) Optional dilation of the watershed borders .
boundary_threshold
(def 0.1) Threshold for the 90-percentile of contour values along kissing points
between adjacent labels. If the 90-percentile < boundary_threshold, we assume that there's no membrane between
both labels, and thus they belong to the same cell and should be merged.
( def 0 . 1 ) Threshold for the 90 - percentile of contour values along kissing points between adjacent labels . If the 90 - percentile < boundary_threshold we assume that there ' s no membrane between both labels and thus they belong to the same cell and should be merged .
median_size
(def 11) Kernel size for the median filter to reduce noise in the output of the contour CNN.
Larger values may blur contours too much, so there won't be seeds inside of the cells. Smaller values may leave
fragmented seeds, that split cells into two or more labels.
( def 11 ) Kernel size for the median filter to reduce noise in the output of the contour CNN . Larger values may blur contours too much so there won ' t be seeds inside of the cells . Smaller values may leave fragmented seeds that split cells into two or more labels .
closing_size
(def 11) Kernel size for the dilation and erosion operation that reduces noise in the
segmentation seeds (ideally, you get a connected component seed per cell). Larger values may merge neighbour cells
together, whereas a low value may fragment a cell into two or more labels.
( def 11 ) Kernel size for the dilation and erosion operation that reduces noise in the segmentation seeds ( ideally you get a connected component seed per cell ) . Larger values may merge neighbour cells together whereas a low value may fragment a cell into two or more labels .
contour_seed_threshold
(def 0) The output of the contour CNN are contour values. Segmentations seeds inside
of cells are computed as contour <= contour_seed_threshold. Usually, 0 works well. This value can be
increased a bit if some cells fail to get a seed. However, contour values on membranes are often as low as 0.1.
Thus, even small values of contour_seed_threshold can create seeds that overflow the cell.
( def 0 ) The output of the contour CNN are contour values . Segmentations seeds inside of cells are computed as contour < = contour_seed_threshold . Usually 0 works well . This value can be increased a bit if some cells fail to get a seed . However contour values on membranes are often as low as 0 . 1 . Thus even small values of contour_seed_threshold can create seeds that overflow the cell .
version
(def 2) Implementation version of the segmentation algorithm. version=1 corresponds to a deprecated
algorithm, kept only for historical comparisons.
( def 2 ) Implementation version of the segmentation algorithm . version = 1 corresponds to a deprecated algorithm kept only for historical comparisons .