soma.py

# Copyright CNRS/Inria/UNS
# Contributor(s): Eric Debreuve (since 2019), Morgane Nadal (2020)
#
# eric.debreuve@cnrs.fr
#
# This software is governed by the CeCILL  license under French law and
# abiding by the rules of distribution of free software.  You can  use,
# modify and/ or redistribute the software under the terms of the CeCILL
# license as circulated by CEA, CNRS and INRIA at the following URL
# "http://www.cecill.info".
#
# As a counterpart to the access to the source code and  rights to copy,
# modify and redistribute granted by the license, users are provided only
# with a limited warranty  and the software's author,  the holder of the
# economic rights,  and the successive licensors  have only  limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading,  using,  modifying and/or developing or reproducing the
# software by the user in light of its specific status of free software,
# that may mean  that it is complicated to manipulate,  and  that  also
# therefore means  that it is reserved for developers  and  experienced
# professionals having in-depth computer knowledge. Users are therefore
# encouraged to load and test the software's suitability as regards their
# requirements in conditions enabling the security of their systems and/or
# data to be ensured and,  more generally, to use and operate it in the
# same conditions as regards security.
#
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.

from __future__ import annotations

from brick.component.glial_cmp import glial_cmp_t
import brick.processing.map_labeling as ml_
from brick.general.type import array_t, py_array_picker_h, site_h

from collections import namedtuple as namedtuple_t
import sys as sy_
from typing import Optional, Sequence, Tuple

import numpy as np_
import scipy.ndimage as im_
import skimage.filters as fl_
import skimage.measure as ms_
import skimage.morphology as mp_
import matplotlib.pyplot as pl_


som_ext_path_t = namedtuple_t("som_ext_path_t", "extension length path idx")


class soma_t(glial_cmp_t):
    #
    __slots__ = ("contour_points", "centroid", "skl_graph", "ext_roots", "graph_roots", "volume_soma_micron", "axes_ellipsoid")

    def __init__(self):
        #
        super().__init__()
        for slot in self.__class__.__slots__:
            setattr(self, slot, None)
        self.ext_roots = []

    @classmethod
    def FromMap(cls, lmp: array_t, uid: int) -> soma_t:
        '''
        Initialize and create the soma object based on the labelled map.
        '''
        #
        instance = cls()

        # Create a boolean map keeping only the soma number 'uid'.
        bmp = lmp == uid
        # Initialize the object with its different fields
        instance.InitializeFromMap(bmp, uid)
        # Find the contour sites of the soma
        contour_map = cls.ContourMap(bmp)
        instance.contour_points = tuple(zip(*contour_map.nonzero()))
        # Find the centroid of the soma
        instance.centroid = np_.array([np_.median(instance.sites[i]) for i in range(3)])

        return instance

    def Extensions(self, max_level: int = sy_.maxsize) -> Tuple[glial_cmp_t]:
        #
        # max_level=1: directly connected extensions
        # ...
        soma_ext = self.extensions.copy()

        current_level_idx = 0
        for level in range(2, max_level + 1):
            next_level_idx = soma_ext.__len__()
            if next_level_idx == current_level_idx:
                break

            for extension in soma_ext[current_level_idx:]:
                soma_ext.extend(extension.extensions)

            current_level_idx = next_level_idx

        return tuple(soma_ext)

    def ContourPointsCloseTo(
        self, point: site_h, max_distance: float
    ) -> Tuple[Optional[Tuple[site_h, ...]], Optional[py_array_picker_h]]:
        '''
        Find the closest contour points of soma to an extension endpoint in the 3D map
        '''
        #
        # Find the points in the contour points of the soma which euclidian distance to the endpoint is inferior to the
        # maximum distance allowed.
        points = tuple(
            contour_point
            for contour_point in self.contour_points
            if (np_.subtract(point, contour_point) ** 2).sum() <= max_distance
        )
        # Return the reformatted closest points
        return __PointsCloseTo__(points)

    def ExtensionPointsCloseTo(
        self, point: site_h, max_distance: float
    ) -> Tuple[Optional[Tuple[site_h, ...]], Optional[py_array_picker_h]]:
        '''
        Find the closest extensions point to an extension endpoint in the 3D map.
        Leave connection paths apart because only detected extension pieces
        (as opposed to invented=connection paths) are considered valid connection targets.
        '''

        points = []

        # Find the points in the extensions sites which euclidian distance to the endpoint is inferior to the
        # maximum distance allowed.
        for extension in self.Extensions():
            ext_sites = tuple(zip(*extension.sites))
            points.extend(
                ext_point
                for ext_point in ext_sites
                if (np_.subtract(point, ext_point) ** 2).sum() <= max_distance
            )

        # Return the reformatted closest points
        return __PointsCloseTo__(tuple(points))

    def BackReferenceSoma(self, glial_cmp: glial_cmp_t) -> None:
        #
        raise NotImplementedError("Soma do not need to back-reference a soma")

    def __str__(self) -> str:
        #
        if self.extensions is None:
            n_extensions = 0
        else:
            n_extensions = self.extensions.__len__()

        return (
            f"Soma.{self.uid}, "
            f"area={self.sites[0].__len__()}, "
            f"contour points={self.contour_points.__len__()}, "
            f"extensions={n_extensions}"
        )

    @staticmethod
    def Map(image: array_t, low: float, high: float, selem: array_t) -> array_t:
        '''
        Perform hysteresis thresholding, dependant on the intensity of the input image, followed by closing/opening.
        '''

        result = __HysterisisImage__(image, low, high)
        # MaximumIntensityProjectionZ(result, output_image_file_name="D:\\MorganeNadal\\M2 report\\for the slides\\soma_hyst_mip.png")

        result = __MorphologicalCleaning__(image, selem)
        # MaximumIntensityProjectionZ(result, output_image_file_name="D:\\MorganeNadal\\M2 report\\for the slides\\soma_opclos_mip.png")

        return result

    @staticmethod
    def DeleteSmallSomas(map: array_t, min_area: int) -> array_t:

        result = map.copy()

        # Label the image region by measuring the connectivity
        lmp = ms_.label(map)

        # Measure properties of labeled image regions
        for region in ms_.regionprops(lmp):
            # Delete regions (set pixel to 0) which area is smaller that the allowed minimum area
            if region.area <= min_area:
                region_sites = (lmp == region.label).nonzero()
                result[region_sites] = 0
                lmp[region_sites] = 0

        return result, lmp

    @staticmethod
    def DeleteBigSomas(map: array_t, lmp: array_t, max_area: int) -> array_t:

        result = map.copy()

        # Measure properties of labeled image regions
        for region in ms_.regionprops(lmp):
            # Delete regions (set pixel to 0) which area is bigger that the allowed maximum area
            # These detected big somas are probably an aggregation of multiple somas
            if region.area >= max_area:
                region_sites = (lmp == region.label).nonzero()
                result[region_sites] = 0
                lmp[region_sites] = 0

        return result, lmp

    @staticmethod
    def ContourMap(map: array_t) -> array_t:
        '''
        Create a map containing only the contour points  of the soma.
        '''
        # Find the 26-connectivity of each voxel
        part_map = ml_.PartLMap(map)
        # The background is labeled with 27, and contour points have a connectivity of at most (<=) 25.
        result = part_map < 26

        return result.astype(np_.int8)

    @staticmethod
    def InfluenceMaps(map: array_t) -> Tuple[array_t, array_t]:
        '''
        Create a 3D influence map based on the somas' influence in the image.
        '''

        # Set background to 0
        background = (map == 0).astype(np_.int8)
        # Find the distance map and the influence map
        dist_map, idx_map = im_.morphology.distance_transform_edt(background, return_indices=True)

        # obj_map = np_.empty_like(map_)
        # for row in range(0, obj_map.shape[0]):
        #     for col in range(0, obj_map.shape[1]):
        #         for dep in range(0, obj_map.shape[2]):
        #             obj_map[row, col, dep] = map_[
        #                 idx_map[0, row, col, dep],
        #                 idx_map[1, row, col, dep],
        #                 idx_map[2, row, col, dep],
        #             ]

        return dist_map, np_.array(map[tuple(idx_map)])

    @staticmethod
    def SomasLMap(somas: Sequence[soma_t], with_extensions: bool = True) -> array_t:
        '''
        Create a map containing the labelled somas and their extensions and connection paths
        with the label of their respective soma.
        '''
        #
        # Initialization
        shape = somas[0].img_shape
        result = np_.zeros(shape, dtype=np_.int64)

        # for each soma, put into the map the soma label at the sites of soma, connexion paths and extensions.
        for soma in somas:
            result[soma.sites] = soma.uid

            # if the soma has extensions, add the soma label at the site of the extensions and of the connection path
            if with_extensions:
                # For each existing connexion path soma-ext, add it on the result map with soma label
                for connection_path in filter(
                    lambda path: path is not None, soma.connection_path.values()
                ):
                    result[connection_path] = soma.uid

                # For each existing connexion path ext-ext, add it on the result map with soma label
                for extension in soma.Extensions():
                    for connection_path in filter(
                        lambda path: path is not None,
                        extension.connection_path.values(),
                    ):
                        result[connection_path] = soma.uid

                    # Also add the extensions on the result map with soma label
                    result[extension.sites] = soma.uid

        return result


def __PointsCloseTo__(
    points: Tuple[site_h, ...]
) -> Tuple[Optional[Tuple[site_h, ...]], Optional[py_array_picker_h]]:
    '''
    Reformat the 3D points coordinates.
    '''
    #
    if points.__len__() > 0:
        points_as_picker = tuple(zip(*points))
    else:
        points = None
        points_as_picker = None

    return points, points_as_picker


def __HysterisisImage__(image: array_t, low: float, high: float) -> array_t:

    # Find the image max value
    max_image = image.max()
    # Find the image minimum non zero value
    nonzero_sites = image.nonzero()
    nonzero_values = image[nonzero_sites]
    min_image = nonzero_values.min()

    # Adapt the hysteresis thresholds to the max and min of the image
    low = low * (max_image - min_image) + min_image
    high = high * (max_image - min_image) + min_image

    # Perform hysteresis
    result = fl_.apply_hysteresis_threshold(image, low, high)
    result = result.astype(np_.int8, copy=False)

    return result


def __MorphologicalCleaning__(image: array_t, selem) -> array_t:

    result = image.copy()

    for dep in range(image.shape[0]):
        result[dep, :, :] = mp_.closing(result[dep, :, :], selem)
        result[dep, :, :] = mp_.opening(result[dep, :, :], selem)

    return result


def MaximumIntensityProjectionZ(img: array_t, cmap: str ='tab20', axis: int = 0, output_image_file_name: str = None) -> None:
    """ Maximum Image Projection on the Z axis. """
    #
    xy = np_.amax(img, axis=axis)
    pl_.imshow(xy, cmap=cmap)
    pl_.show(block=True)
    if output_image_file_name is not None:
        pl_.imsave(output_image_file_name, xy, cmap=cmap)
        print('Image saved in', output_image_file_name)