Source code for ClearMap.Analysis.Measurements.maxima_detection

# -*- coding: utf-8 -*-
"""
MaximaDetection
===============

Collection of routines to detect maxima in scalar images.

Used for finding cells or intensity peaks.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__ = 'GPLv3 - GNU General Pulic License v3 (see LICENSE)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__ = 'https://idisco.info'
__download__ = 'https://www.github.com/ChristophKirst/ClearMap2'

import numpy as np

import scipy.ndimage.measurements as ndm
import scipy.ndimage.filters as ndf

import ClearMap.ImageProcessing.GreyReconstruction as gr

import ClearMap.Utils.Timer as tmr
import ClearMap.Utils.HierarchicalDict as hdict


##############################################################################
# Basic Transforms
##############################################################################
   


def h_max_transform(source, h_max):
    """
    H-maximum transform of an array.

    Arguments
    ---------
    source : array
        Input image.
    h_max : float or None
        H parameter of h-max transform, if None return source.

    Returns
    -------
    transformed : array
        H-max transformed image if h is not None.
    """
    if h_max is None:
        return source
    else:
        return gr.reconstruct(source - h_max, source)





def local_max(source, shape=5):
    """
    Calculates local maxima of an image

    Arguments
    ---------
    source : array
        Input image.
    shape : int, tuple or None
        Shape of the volume to search for maxima.

    Returns
    -------
    local_max : array
        Mask that is True at local maxima.
    """

    if shape is None:
        return source

    if not isinstance(shape, tuple):
        shape = (shape,) * source.ndim

    return ndf.maximum_filter(source, size=shape) == source

    



def extended_max(source, h_max=0):
    """
    Calculates extended h-maxima of an image

    Arguments
    ---------
    source : array
        Input image.
    h_max : float or None
        H parameter of h-max transform. If None, calculate local maxima only.

    Returns
    -------
    maxima : array
        Extended maxima of the image.

    Note
    ----
    Extended maxima are the local maxima of the h-max transform.
    """
    # h max transform
    if h_max is not None and h_max > 0:
        source = h_max_transform(source, h_max)

    # max
    return local_max(source)



##############################################################################
# Maxima Detection 
##############################################################################




def find_maxima(source, h_max=None, shape=5, threshold=None, verbose=None):
    """
    Find local and extended maxima in an image.

    Arguments
    ---------
    source : array
        The source data.
    h_max : float or None
        H parameter for the initial h-Max transform.
        If None, do not perform a h-max transform.
    shape : int or tuple
        Shape for the structure element for the local maxima filter.
    threshold : float or None
        If float, include only maxima larger than this threshold.
    verbose : bool
        Print progress info.

    Returns
    -------
    maxima : array
        Binary image with True pixel at extended maxima.

    Notes
    -----
    This routine performs a h-max transform, followed by a local maxima search
    and thresholding of the maxima.

    See also
    --------
    :func:`h_max_transform`, :func:`local_max`
    """
    if verbose:
        timer = tmr.Timer()
        hdict.pprint(head='Find Maxima:', h_max=h_max, shape=shape, threshold=threshold)

    # extended maxima
    maxima = h_max_transform(source, h_max=h_max)

    # local maxima
    maxima = local_max(maxima, shape=shape)

    # thresholding
    if threshold is not None:
        maxima = np.logical_and(maxima, source >= threshold)

    if verbose:
        timer.print_elapsed_time(head='Find Maxima')

    return maxima





def find_center_of_maxima(source, maxima=None, label=None, verbose=False):
    """
    Find center of detected maxima weighted by intensity

    Arguments
    ---------
    source : array
        Intensity image data.
    maxima : array or None
        Binary array indicating the maxima. If label is not None this can be None.
    label : array or None
        Labeled image of the shapes of the maxima.
        If None, determined from maxima.
    verbose : bool
        Print progress info.

    Returns
    -------
    coordinates : array
        Coordinates of the n centers of maxima as (n,d)-array.
    """
    if verbose:
        timer = tmr.Timer()
        print('Center of Maxima initialized!')

    # center of maxima
    if label is None:
        label, n_label = ndm.label(maxima)
    else:
        n_label = label.max()

    if n_label > 0:
        centers = np.array(ndm.center_of_mass(source, label, index=np.arange(1, n_label)))
    else:
        centers = np.zeros((0, source.ndim))

    if verbose:
        timer.print_elapsed_time(f'Center of Maxima: {centers.shape[0]} maxima detected')

    return centers