Source code for pywsi.segmentation.max_clustering

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
#import pyximport; pyximport.install()
import numpy as np
import skimage.measure

from ._max_clustering_cython import _max_clustering_cython
"""
This and _max_clustering_cython.pyx moudles are courtesy of
https://github.com/DigitalSlideArchive/HistomicsTK/tree/master/histomicstk
(Apache Licensed, see: https://github.com/DigitalSlideArchive/HistomicsTK/blob/dcd0d53f0cdc1424e1bedc0f0c0871a7a0e9fd70/LICENSE)
"""


def max_clustering(im_response, im_fgnd_mask, r=10):
    """Local max clustering pixel aggregation for nuclear segmentation.
    Takes as input a constrained log or other filtered nuclear image, a binary
    nuclear mask, and a clustering radius. For each pixel in the nuclear mask,
    the local max is identified. A hierarchy of local maxima is defined, and
    the root nodes used to define the label image. References [#]_ and [#]_

    Parameters
    ----------
    im_response : array_like
        A filtered-smoothed image where the maxima correspond to nuclear
        center. Typically obtained by constrained-LoG filtering on a
        hematoxylin intensity image obtained from ColorDeconvolution.
    im_fgnd_mask : array_like
        A binary mask of type boolean where nuclei pixels have value
        'True', and non-nuclear pixels have value 'False'.
    r : float
        A scalar defining the clustering radius. Default value = 10.

    Returns
    -------
    im_label : array_like
        im_label image where positive values correspond to foreground pixels that
        share mutual sinks.
    seeds : array_like
        An N x 2 array defining the (x,y) coordinates of nuclei seeds.
    max_response : array_like
        An N x 1 array containing the maximum response value corresponding to
        seeds.


    References
    ----------
    .. [#] XW. Wu et al "The local maximum clustering method and its
       application in microarray gene expression data analysis,"
       EURASIP J. Appl. Signal Processing, volume 2004, no.1, pp.53-63,
       2004.
    .. [#] Y. Al-Kofahi et al "Improved Automatic Detection and Segmentation
       of Cell Nuclei in Histopathology Images" in IEEE Transactions on
       Biomedical Engineering,vol.57,no.4,pp.847-52, 2010.

    """

    # find local maxima of all foreground pixels
    mval, mind = _max_clustering_cython(im_response,
                                        im_fgnd_mask.astype(np.int32), r)
    # identify connected regions of local maxima and define their seeds
    im_label = skimage.measure.label(im_fgnd_mask & (im_response == mval))

    # compute normalized response
    min_resp = np.nanmin(im_response)
    max_resp = np.nanmin(im_response)
    resp_range = max_resp - min_resp

    im_response_nmzd = (im_response - min_resp) / resp_range

    # compute object properties
    obj_props = skimage.measure.regionprops(im_label, im_response_nmzd)

    obj_props = [
        prop for prop in obj_props
        if np.isfinite(prop.weighted_centroid).all()
    ]

    num_labels = len(obj_props)

    # extract object seeds
    seeds = np.array(
        [obj_props[i].weighted_centroid for i in range(num_labels)])
    seeds = np.round(seeds).astype(np.int)

    # fix seeds outside the object region - happens for non-convex objects
    for i in range(num_labels):

        sy = seeds[i, 0]
        sx = seeds[i, 1]

        if im_label[sy, sx] == obj_props[i].label:
            continue

        # find object point with closest manhattan distance to center of mass
        pts = obj_props[i].coords

        ydist = np.abs(pts[:, 0] - sy)
        xdist = np.abs(pts[:, 1] - sx)

        seeds[i, :] = pts[np.argmin(xdist + ydist), :]

        assert im_label[seeds[i, 0], seeds[i, 1]] == obj_props[i].label

    # get seed responses

    assert im_response.shape[1] >= 1, 'im_response.shape: {}'.format(
        im_response.shape)

    try:
        max_response = im_response[seeds[:, 0], seeds[:, 1]]
    except IndexError:
        # No seeds found!
        max_response = np.array([])
    # set label of each foreground pixel to the label of its nearest peak
    im_label_flat = im_label.ravel()

    pind = np.flatnonzero(im_fgnd_mask)

    mind_flat = mind.ravel()

    im_label_flat[pind] = im_label_flat[mind_flat[pind]]

    # return
    return im_label, seeds, max_response