Mentions légales du service

Skip to content
Snippets Groups Projects
cell.py 33.3 KiB
Newer Older
DEBREUVE Eric's avatar
DEBREUVE Eric committed
# 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.

import re as re_
from typing import Any, Dict, Tuple, Union

import numpy as nmpy
DEBREUVE Eric's avatar
DEBREUVE Eric committed
import pandas as pd_
import scipy.stats as st_
from logger_36 import LOGGER
DEBREUVE Eric's avatar
DEBREUVE Eric committed

import brick.task.ellipsoid_fit as bf_
import brick.task.unit
from brick.type.base import array_t
from brick.type.soma import soma_t

DEBREUVE Eric's avatar
DEBREUVE Eric committed

def FindGraphsRootWithEdges(
    soma: soma_t, ext_nfo: Dict[str, Union[array_t, Any]]
) -> dict:
DEBREUVE Eric's avatar
DEBREUVE Eric committed
    """
    Finds the soma roots of the graph extension.
    """
    # For a given soma, find the roots of the graphs
    root_nodes = {}

    # Finds the primary extensions
    primary_extension_uids = tuple(extension.uid for extension in soma.extensions)
    # print(primary_extension_uids, '\nn = ', len(primary_extension_uids))

    # List of the degree 1 nodes of the graph
    for node1_id, node2_id, edge_nfo in soma.skl_graph.edges.data("details"):
        if (soma.skl_graph.degree[node1_id] == 1) or (
            soma.skl_graph.degree[node2_id] == 1
        ):
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            # Find the pixels of the terminal extension
            sites = ext_nfo["lmp_soma"][edge_nfo.sites]
            ext_uid = nmpy.unique(sites)[-1]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            # sites > 0 because ext_nfo['lmp'] do not contain the connexions

            # Save the root node candidates (one-degree nodes)
            if ext_uid in primary_extension_uids:
                if soma.skl_graph.degree[node1_id] == 1:
                    root_node = node1_id
                else:
                    root_node = node2_id

                # Get the node coordinates and extend them to the 26 neighboring voxels
                root_node_coor = _GetNodesCoordinates((root_node,))[
                    0
                ]  # tuple('x-y-z') -> list[(x,y,z)]
DEBREUVE Eric's avatar
DEBREUVE Eric committed

                root_sites = set(
                    (
                        root_node_coor[0] + i,
                        root_node_coor[1] + j,
                        root_node_coor[2] + k,
                    )
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                    for i in (-1, 0, 1)
                    for j in (-1, 0, 1)
                    for k in (-1, 0, 1)
                    if i != 0 or j != 0 or k != 0
                )
DEBREUVE Eric's avatar
DEBREUVE Eric committed

                # Find the intersection between the extended root node candidate and the soma contour points
                intersections = set(soma.contour_points).intersection(root_sites)

                # if the graph root sites are included in the soma extensions sites (non-nul intersection):
                if len(intersections) > 0:
                    # Keep the info of the root node. Key = ext uid, Value = root node
                    root_nodes[ext_uid] = root_node
                    ## By construction, only one root node possible for an ext
DEBREUVE Eric's avatar
DEBREUVE Eric committed

    return (
        root_nodes  # TODO: find out why there are less root points than extensions !!
    )
DEBREUVE Eric's avatar
DEBREUVE Eric committed


def FindGraphsRootWithNodes(soma: soma_t) -> dict:
    """
    Find the roots of the {extension+connexion} graphs to be lined to the soma.
    Add a key "root" (bool) in the dict of nodes attributes.
    """

    # Find the nodes of degree == 1, and the str coordinates of the nodes
    node_degree_bool = tuple(degree == 1 for _, degree in soma.skl_graph.degree)
    node_coord = tuple(xyz for xyz, _ in soma.skl_graph.degree)

    root_nodes = {}

    # get the coordinates of the nodes (x,y,z)
    coordinates = _GetNodesCoordinates(node_coord)
DEBREUVE Eric's avatar
DEBREUVE Eric committed

    # get a list with elements = (extension_uid, root coordinates), which length is the number of primary extensions
    roots = soma.ext_roots

    # for each node in the graph, search among the degree 1 nodes the nodes that are roots (linked to soma)
    for node in range(len(coordinates)):
        if node_degree_bool[node]:
            # compare the coor with end points
            for ext_root in roots:
                if ext_root[1] == coordinates[node]:
                    root_nodes[ext_root[0]] = node_coord[node]

    if root_nodes.__len__() != roots.__len__():
        # raise ValueError("Number of extensions roots not equal to number of graph roots.")
            f"\nNumber of extensions roots: {root_nodes.__len__()} not equal to number of graph roots: {roots.__len__()}."
        )
DEBREUVE Eric's avatar
DEBREUVE Eric committed

    return root_nodes


def _GetNodesCoordinates(node_coord: Tuple[str, ...]) -> list:
DEBREUVE Eric's avatar
DEBREUVE Eric committed
    """
    Input: nodes attributes -> Tuple('x1-y1-z1', 'x2-y2-z2', ...) .
    Output: coordinates -> List[Tuple(x1,y1,z1), Tuple(x2,y2,z2), ...]
    """
    coord = []
    for c in node_coord:
        coord.append(c)

    for node in range(len(node_coord)):
        coord_node = coord[node]
        pattern = r"\d+"
        coord_node = re_.findall(pattern, coord_node)
        coor = []
        for i in range(3):
            coor.append(int(coord_node[i]))
        coor = tuple(coor)
        coord[node] = coor

    return coord


def ExtractFeaturesInDF(
    name_file,
    somas,
    voxel_size_in_micron: array_t,
    bins_length: array_t,
    bins_curvature: array_t,
    scale_map: array_t,
    decimals: int = 4,
    condition=None,
    duration=None,
):
DEBREUVE Eric's avatar
DEBREUVE Eric committed
    """
    Extract the features from somas and graphs.
    Returns a pandas dataframe, without NaN.
    """
    # Store the condition and duration of the image
    # if (condition is None) or (duration is None):
    # Condition = re_.findall(r"[A-Z]{3}", name_file)[0]
    # Duration = re_.findall(r"\dH", name_file)
    # Find our whether the duration is in hours or in weeks
    # if Duration.__len__() == 0:
    #     Duration = re_.findall(r"\dW", name_file)[0]
    # else:
    #     Duration = Duration[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed

    somas_features_dict = {}  # Dict{soma 1: [features], soma 2: [features], ...}
    columns = [
        "condition",
        "duration",
        "soma uid",
        "coef_V_soma__V_convex_hull",
        "coef_axes_ellips_b__a",
        "coef_axes_ellips_c__a",
        "spherical_angles_eva",
        "spherical_angles_evb",
        #
        "N_nodes",
        "N_ext",
        "N_primary_ext",
        "N_sec_ext",
        "min_degree",
        "mean_degree",
        "median_degree",
        "max_degree",
        "std_degree",
        #
        "total_ext_length",
        "min_length",
        "mean_length",
        "median_length",
        "max_length",
        "std_lengths",
        "entropy_lengths",
        "hist_lengths",
        "min_thickness",
        "mean_thickness",
        "median_thickness",
        "max_thickness",
        "std_thickness",
        "entropy_thickness",
        "min_volume",
        "mean_volume",
        "median_volume",
        "max_volume",
        "std_volume",
        "entropy_volume",
        "min_curvature",
        "max_curvature",
        "mean_curvature",
        "median_curvature",
        "std_curvature",
        "entropy_curvature",
        "hist_curvature",
        "min_torsion",
        "max_torsion",
        "mean_torsion",
        "median_torsion",
        "std_torsion",
        "entropy_torsion",
        #
        "total_ext_length_P",
        "min_length_P",
        "mean_length_P",
        "median_length_P",
        "max_length_P",
        "std_lengths_P",
        "entropy_lengths_P",
        "hist_lengths_P",
        "min_thickness_P",
        "mean_thickness_P",
        "median_thickness_P",
        "max_thickness_P",
        "std_thickness_P",
        "entropy_thickness_P",
        "min_volume_P",
        "mean_volume_P",
        "median_volume_P",
        "max_volume_P",
        "std_volume_P",
        "entropy_volume_P",
        "min_curvature_P",
        "max_curvature_P",
        "mean_curvature_P",
        "median_curvature_P",
        "std_curvature_P",
        "entropy_curvature_P",
        "hist_curvature_P",
        "min_torsion_P",
        "max_torsion_P",
        "mean_torsion_P",
        "median_torsion_P",
        "std_torsion_P",
        "entropy_torsion_P",
        #
        "total_ext_length_S",
        "min_length_S",
        "mean_length_S",
        "median_length_S",
        "max_length_S",
        "std_lengths_S",
        "entropy_lengths_S",
        "hist_lengths_S",
        "min_thickness_S",
        "mean_thickness_S",
        "median_thickness_S",
        "max_thickness_S",
        "std_thickness_S",
        "entropy_thickness_S",
        "min_volume_S",
        "mean_volume_S",
        "median_volume_S",
        "max_volume_S",
        "std_volume_S",
        "entropy_volume_S",
        "min_curvature_S",
        "max_curvature_S",
        "mean_curvature_S",
        "median_curvature_S",
        "std_curvature_S",
        "entropy_curvature_S",
        "hist_curvature_S",
        "min_torsion_S",
        "max_torsion_S",
        "mean_torsion_S",
        "median_torsion_S",
        "std_torsion_S",
        "entropy_torsion_S",
    ]

    Condition = (
        Duration
    ) = (
        min_degree
    ) = (
        mean_degree
    ) = (
        median_degree
    ) = (
        max_degree
    ) = (
        std_degree
    ) = (
        total_ext_length_S
    ) = (
        min_length_S
    ) = (
        mean_length_S
    ) = (
        median_length_S
    ) = (
        max_length_S
    ) = (
        std_lengths_S
    ) = (
        entropy_lengths_S
    ) = (
        hist_lengths_S
    ) = (
        min_thickness_S
    ) = (
        mean_thickness_S
    ) = (
        median_thickness_S
    ) = (
        max_thickness_S
    ) = (
        std_thickness_S
    ) = (
        entropy_thickness_S
    ) = (
        min_volume_S
    ) = (
        mean_volume_S
    ) = (
        median_volume_S
    ) = (
        max_volume_S
    ) = (
        std_volume_S
    ) = (
        entropy_volume_S
    ) = (
        min_curvature_S
    ) = (
        max_curvature_S
    ) = (
        mean_curvature_S
    ) = (
        median_curvature_S
    ) = (
        std_curvature_S
    ) = (
        entropy_curvature_S
    ) = (
        hist_curvature_S
    ) = (
        min_torsion_S
    ) = (
        max_torsion_S
    ) = mean_torsion_S = median_torsion_S = std_torsion_S = entropy_torsion_S = None

DEBREUVE Eric's avatar
DEBREUVE Eric committed
    for soma in somas:
        # -- Soma features

        # Axes of the best fitting ellipsoid
        # a > b > c
        (
            _,
            _,
            soma.axes_ellipsoid,
            _,
            spherical_coor,
            _,
            volume_convex_hull,
        ) = bf_.FindBestFittingEllipsoid3D(soma)
DEBREUVE Eric's avatar
DEBREUVE Eric committed

        # These ratios give info about the shape of the soma. ex: rather flat, rather patatoide, rather spherical...
DEBREUVE Eric's avatar
DEBREUVE Eric committed
        if type(soma.axes_ellipsoid[0]) is str:
            Coef_axes_ellips_b__a = None
            Coef_axes_ellips_c__a = None
            spherical_angles_eva = None
            spherical_angles_evb = None
            soma.volume_soma_micron = None
            Coef_V_soma__V_convex_hull = None
        else:
            Coef_axes_ellips_b__a = soma.axes_ellipsoid[0] / soma.axes_ellipsoid[2]
            Coef_axes_ellips_c__a = soma.axes_ellipsoid[1] / soma.axes_ellipsoid[2]

            # Spherical angles give the orientation of the somas in the 3D volume
            spherical_angles_eva = (spherical_coor[0][1], spherical_coor[0][2])
            spherical_angles_evb = (spherical_coor[1][1], spherical_coor[1][2])

            # Volume of the  in micron**3
            soma.volume_soma_micron = brick.task.unit.ToMicron(
                voxel_size_in_micron,
                dimension=(0, 1, 2),
                decimals=2,
            )
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            # Calculates volume of soma's convex hull in voxel volume
            # Take into account anisotropy of the 3D space ( volume = x * y * z with z > x=y)
            volume_convex_hull = (
                volume_convex_hull * voxel_size_in_micron[2] / voxel_size_in_micron[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed

            # Volume of the soma / Volume of its convex hull gives the info about the convexity of the soma
            # If close to 0, the soma has a lot of invaginations, if close to 1, it is smooth and convex
            Coef_V_soma__V_convex_hull = len(soma.sites[0]) / round(
                volume_convex_hull + 0.5
            )
DEBREUVE Eric's avatar
DEBREUVE Eric committed

        # -- Extension features
        # Graph features

        # number of nodes except the constructed ones from node soma to the roots
        N_nodes = soma.skl_graph.n_nodes - len(soma.graph_roots)
        # number of edges except the constructed ones from node soma to the roots
        N_ext = soma.skl_graph.n_edges - len(soma.graph_roots)
        # number of primary edges = linked to the soma except the constructed ones from node soma to the roots
        N_primary_ext = len(soma.graph_roots)
        # number of secondary edges = not linked to the soma.
        N_sec_ext = N_ext - N_primary_ext

        print(
            f" Soma {soma.uid}\n"
            f"N nodes = {N_nodes}\n"
            f"N edges = {N_ext}\n"
            f"N primary extensions = {N_primary_ext}\n"
            f"N secondary extensions = {N_sec_ext}\n"
        )

        if N_primary_ext > 0:
            ext_lengths = [_elm[-1] for _elm in soma.skl_graph.lengths]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            for idx, length in enumerate(ext_lengths):
                ext_lengths[idx] = brick.task.unit.ToMicron(
                    length, voxel_size_in_micron, decimals=decimals
            total_ext_length = brick.task.unit.ToMicron(
                soma.skl_graph.length, voxel_size_in_micron, decimals=decimals
DEBREUVE Eric's avatar
DEBREUVE Eric committed

            # Lengths histogram
            hist_lengths = nmpy.histogram(ext_lengths, bins_length)[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            #
            # min, mean, median, max and standard deviation of the ALL extensions
            min_length = brick.task.unit.ToMicron(
                min(ext_lengths), voxel_size_in_micron, decimals=decimals
            mean_length = brick.task.unit.ToMicron(
                nmpy.mean(ext_lengths).item(), voxel_size_in_micron, decimals=decimals
            median_length = brick.task.unit.ToMicron(
                nmpy.median(ext_lengths).item(), voxel_size_in_micron, decimals=decimals
            max_length = brick.task.unit.ToMicron(
                max(ext_lengths), voxel_size_in_micron, decimals=decimals
            std_lengths = nmpy.std(ext_lengths)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(ext_lengths) > 0:
                entropy_lengths = st_.entropy(ext_lengths)
            else:
                entropy_lengths = 0

            # Find the thickness of the extensions
            for ___, ___, edge in soma.skl_graph.edges.data("details"):
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                if edge is not None:
                    edge.widths = scale_map[edge.sites] * voxel_size_in_micron[1]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            mean_widths = soma.skl_graph.edge_reduced_widths()
            ext_thickness = nmpy.array(mean_widths) ** 2
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            min_thickness = min(ext_thickness)
            mean_thickness = nmpy.mean(ext_thickness)
            median_thickness = nmpy.median(ext_thickness)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            max_thickness = max(ext_thickness)
            std_thickness = nmpy.std(ext_thickness)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(ext_lengths) > 0:
                entropy_thickness = st_.entropy(ext_thickness)
            else:
                entropy_thickness = 0
            #
            # Find the volume of the extensions
            ext_volume = nmpy.array(ext_lengths) * ext_thickness
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            min_volume = min(ext_volume)
            mean_volume = nmpy.mean(ext_volume)
            median_volume = nmpy.median(ext_volume)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            max_volume = max(ext_volume)
            std_volume = nmpy.std(ext_volume)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(ext_volume) > 0:
                entropy_volume = st_.entropy(ext_volume)
            else:
                entropy_volume = 0
            #
            # Curvature and Torsion
            curvatures, torsions = soma.skl_graph.curvature_and_torsion(
                size_voxel=voxel_size_in_micron
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            #
            min_curvature = min(curvatures)
            max_curvature = max(curvatures)
            mean_curvature = nmpy.mean(curvatures)
            median_curvature = nmpy.median(curvatures)
            std_curvature = nmpy.std(curvatures)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(curvatures) > 0:
                entropy_curvature = st_.entropy(curvatures)
            else:
                entropy_curvature = 0
            hist_curvature = nmpy.histogram(curvatures, bins_curvature)[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            #
            min_torsion = min(torsions)
            max_torsion = max(torsions)
            mean_torsion = nmpy.mean(torsions)
            median_torsion = nmpy.median(torsions)
            std_torsion = nmpy.std(torsions)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(torsions) > 0:
                entropy_torsion = st_.entropy(torsions)
            else:
                entropy_torsion = 0

            # PRIMARY extensions
            ext_lengths_P = list(soma.skl_graph.primary_edge_lengths(soma))
            for idx, length in enumerate(ext_lengths_P):
                ext_lengths_P[idx] = brick.task.unit.ToMicron(
                    length, voxel_size_in_micron, decimals=decimals
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            total_ext_length_P = sum(ext_lengths_P)
            #
            # Lengths histogram
            hist_lengths_P = nmpy.histogram(ext_lengths_P, bins_length)[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            #
            # min, mean, median, max and standard deviation of the PRIMARY extensions
            min_length_P = min(ext_lengths_P)
            mean_length_P = nmpy.mean(ext_lengths_P)
            median_length_P = nmpy.median(ext_lengths_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            max_length_P = max(ext_lengths_P)
            std_lengths_P = nmpy.std(ext_lengths_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(ext_lengths_P) > 0:
                entropy_lengths_P = st_.entropy(ext_lengths_P)
            else:
                entropy_lengths_P = 0

            #
            mean_widths_P = soma.skl_graph.P_edge_reduced_widths(soma)
            ext_thickness_P = nmpy.array(mean_widths_P) ** 2
DEBREUVE Eric's avatar
DEBREUVE Eric committed

            min_thickness_P = min(ext_thickness_P)
            mean_thickness_P = nmpy.mean(ext_thickness_P)
            median_thickness_P = nmpy.median(ext_thickness_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            max_thickness_P = max(ext_thickness_P)
            std_thickness_P = nmpy.std(ext_thickness_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(ext_thickness_P) > 0:
                entropy_thickness_P = st_.entropy(ext_thickness_P)
            else:
                entropy_thickness_P = 0

            #
            #
            ext_volume_P = nmpy.array(ext_lengths_P) * ext_thickness_P
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            min_volume_P = min(ext_volume_P)
            mean_volume_P = nmpy.mean(ext_volume_P)
            median_volume_P = nmpy.median(ext_volume_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            max_volume_P = max(ext_volume_P)
            std_volume_P = nmpy.std(ext_volume_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(ext_volume_P) > 0:
                entropy_volume_P = st_.entropy(ext_volume_P)
            else:
                entropy_volume_P = 0

            #
            # Curvature and Torsion
            curvatures_P, torsions_P = soma.skl_graph.P_curvature_and_torsion(
                size_voxel=voxel_size_in_micron, soma=soma
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            #
            min_curvature_P = min(curvatures_P)
            max_curvature_P = max(curvatures_P)
            mean_curvature_P = nmpy.mean(curvatures_P)
            median_curvature_P = nmpy.median(curvatures_P)
            std_curvature_P = nmpy.std(curvatures_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(curvatures_P) > 0:
                entropy_curvature_P = st_.entropy(curvatures_P)
            else:
                entropy_curvature_P = 0

            hist_curvature_P = nmpy.histogram(curvatures_P, bins_curvature)[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            #
            min_torsion_P = min(torsions_P)
            max_torsion_P = max(torsions_P)
            mean_torsion_P = nmpy.mean(torsions_P)
            median_torsion_P = nmpy.median(torsions_P)
            std_torsion_P = nmpy.std(torsions_P)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            if any(torsions_P) > 0:
                entropy_torsion_P = st_.entropy(torsions_P)
            else:
                entropy_torsion_P = 0
            #
            # Secondary extensions
            if N_sec_ext > 0:
                # min, mean, median, max and standard deviation of the degrees of non-leaves nodes
                min_degree = soma.skl_graph.min_degree_except_leaves_and_roots
                mean_degree = soma.skl_graph.mean_degree_except_leaves_and_roots
                median_degree = soma.skl_graph.median_degree_except_leaves_and_roots
                max_degree = soma.skl_graph.max_degree_except_leaves_an_roots
                std_degree = soma.skl_graph.std_degree_except_leaves_and_roots

                # SECONDARY extensions length
                ext_lengths_S = list(soma.skl_graph.secondary_edge_lengths(soma))
                for idx, length in enumerate(ext_lengths_S):
                    ext_lengths_S[idx] = brick.task.unit.ToMicron(
                        length, voxel_size_in_micron, decimals=decimals
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                total_ext_length_S = sum(ext_lengths_S)
                #
                # Lengths histogram
                hist_lengths_S = nmpy.histogram(ext_lengths_S, bins_length)[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                #
                # min, mean, median, max and standard deviation of the PRIMARY extensions
                min_length_S = min(ext_lengths_S)
                mean_length_S = nmpy.mean(ext_lengths_S)
                median_length_S = nmpy.median(ext_lengths_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                max_length_S = max(ext_lengths_S)
                std_lengths_S = nmpy.std(ext_lengths_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                if any(ext_lengths_S) > 0:
                    entropy_lengths_S = st_.entropy(ext_lengths_S)
                else:
                    entropy_lengths_S = 0

                #
                mean_widths_S = soma.skl_graph.S_edge_reduced_widths(soma)
                ext_thickness_S = nmpy.array(mean_widths_S) ** 2
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                min_thickness_S = min(ext_thickness_S)
                mean_thickness_S = nmpy.mean(ext_thickness_S)
                median_thickness_S = nmpy.median(ext_thickness_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                max_thickness_S = max(ext_thickness_S)
                std_thickness_S = nmpy.std(ext_thickness_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                if any(ext_thickness_S) > 0:
                    entropy_thickness_S = st_.entropy(ext_thickness_S)
                else:
                    entropy_thickness_S = 0

                #
                ext_volume_S = nmpy.array(ext_lengths_S) * ext_thickness_S
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                min_volume_S = min(ext_volume_S)
                mean_volume_S = nmpy.mean(ext_volume_S)
                median_volume_S = nmpy.median(ext_volume_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                max_volume_S = max(ext_volume_S)
                std_volume_S = nmpy.std(ext_volume_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                if any(ext_volume_S) > 0:
                    entropy_volume_S = st_.entropy(ext_volume_S)
                else:
                    entropy_volume_S = 0

                #
                # Curvature and Torsion
                curvatures_S, torsions_S = soma.skl_graph.S_curvature_and_torsion(
                    size_voxel=voxel_size_in_micron, soma=soma
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                #
                min_curvature_S = min(curvatures_S)
                max_curvature_S = max(curvatures_S)
                mean_curvature_S = nmpy.mean(curvatures_S)
                median_curvature_S = nmpy.median(curvatures_S)
                std_curvature_S = nmpy.std(curvatures_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                if any(curvatures_S) > 0:
                    entropy_curvature_S = st_.entropy(curvatures_S)
                else:
                    entropy_curvature_S = 0

                hist_curvature_S = nmpy.histogram(ext_lengths_S, bins_curvature)[0]
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                #
                min_torsion_S = min(torsions_S)
                max_torsion_S = max(torsions_S)
                mean_torsion_S = nmpy.mean(torsions_S)
                median_torsion_S = nmpy.median(torsions_S)
                std_torsion_S = nmpy.std(torsions_S)
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                if any(torsions_S) > 0:
                    entropy_torsion_S = st_.entropy(torsions_S)
                else:
                    entropy_torsion_S = 0
                #
            # If no secondary extensions, give certain value to parameters
            if N_sec_ext == 0:
                # min, mean, median, max and standard deviation of the degrees of non-leaves nodes
DEBREUVE Eric's avatar
DEBREUVE Eric committed
                min_degree = 1
                mean_degree = 1
                median_degree = 1
                max_degree = 1
                std_degree = 0

                total_ext_length_S = 0
                min_length_S = 0
                mean_length_S = 0
                median_length_S = 0
                max_length_S = 0
                std_lengths_S = 0
                entropy_lengths_S = 0
                hist_lengths_S = 0
                #
                min_thickness_S = 0
                mean_thickness_S = 0
                median_thickness_S = 0
                max_thickness_S = 0
                std_thickness_S = 0
                entropy_thickness_S = 0
                #
                min_volume_S = 0
                mean_volume_S = 0
                median_volume_S = 0
                max_volume_S = 0
                std_volume_S = 0
                entropy_volume_S = 0
                #
                min_curvature_S = -1
                max_curvature_S = -1
                mean_curvature_S = -1
                median_curvature_S = -1
                std_curvature_S = 0
                entropy_curvature_S = 0
                hist_curvature_S = 0
                #
                min_torsion_S = -1
                max_torsion_S = -1
                mean_torsion_S = -1
                median_torsion_S = -1
                std_torsion_S = 0
                entropy_torsion_S = 0

        else:
            min_degree = 0
            mean_degree = 0
            median_degree = 0
            max_degree = 0
            std_degree = 0
            #
            total_ext_length = 0
            min_length = 0
            mean_length = 0
            median_length = 0
            max_length = 0
            std_lengths = 0
            entropy_lengths = 0
            hist_lengths = 0
            min_thickness = 0
            mean_thickness = 0
            median_thickness = 0
            max_thickness = 0
            std_thickness = 0
            entropy_thickness = 0
            min_volume = 0
            mean_volume = 0
            median_volume = 0
            max_volume = 0
            std_volume = 0
            entropy_volume = 0
            min_curvature = -1
            max_curvature = -1
            mean_curvature = -1
            median_curvature = -1
            std_curvature = 0
            entropy_curvature = 0
            hist_curvature = 0
            min_torsion = -1
DEBREUVE Eric's avatar
DEBREUVE Eric committed
            mean_torsion = -1
            median_torsion = -1
            std_torsion = 0
            entropy_torsion = 0
            #
            total_ext_length_P = 0
            min_length_P = 0
            mean_length_P = 0
            median_length_P = 0
            max_length_P = 0
            std_lengths_P = 0
            entropy_lengths_P = 0
            hist_lengths_P = 0
            min_thickness_P = 0
            mean_thickness_P = 0
            median_thickness_P = 0
            max_thickness_P = 0
            std_thickness_P = 0
            entropy_thickness_P = 0
            min_volume_P = 0
            mean_volume_P = 0
            median_volume_P = 0
            max_volume_P = 0
            std_volume_P = 0
            entropy_volume_P = 0
            min_curvature_P = -1
            max_curvature_P = -1
            mean_curvature_P = -1
            median_curvature_P = -1
            std_curvature_P = 0
            entropy_curvature_P = 0
            hist_curvature_P = 0
            min_torsion_P = -1
            max_torsion_P = -1
            mean_torsion_P = -1
            median_torsion_P = -1
            std_torsion_P = 0
            entropy_torsion_P = 0
            #
            total_ext_length_S = 0
            min_length_S = 0
            mean_length_S = 0
            median_length_S = 0
            max_length_S = 0
            std_lengths_S = 0
            entropy_lengths_S = 0
            hist_lengths_S = 0
            min_thickness_S = 0
            mean_thickness_S = 0
            median_thickness_S = 0
            max_thickness_S = 0
            std_thickness_S = 0
            entropy_thickness_S = 0
            min_volume_S = 0
            mean_volume_S = 0
            median_volume_S = 0
            max_volume_S = 0
            std_volume_S = 0
            entropy_volume_S = 0
            min_curvature_S = -1
            max_curvature_S = -1
            mean_curvature_S = -1
            median_curvature_S = -1
            std_curvature_S = 0
            entropy_curvature_S = 0
            hist_curvature_S = 0
            min_torsion_S = -1
            max_torsion_S = -1
            mean_torsion_S = -1
            median_torsion_S = -1
            std_torsion_S = 0
            entropy_torsion_S = 0

        # Create a dictionary with all the features for every somas
        somas_features_dict[name_file + f" soma {soma.uid}"] = [
            Condition,
            Duration,
            soma.uid,
            Coef_V_soma__V_convex_hull,
            Coef_axes_ellips_b__a,
            Coef_axes_ellips_c__a,
            spherical_angles_eva,
            spherical_angles_evb,
            N_nodes,
            N_ext,
            N_primary_ext,
            N_sec_ext,
            min_degree,
            mean_degree,
            median_degree,
            max_degree,
            std_degree,
            #
            total_ext_length,
            min_length,
            mean_length,
            median_length,
            max_length,
            std_lengths,
            entropy_lengths,
            hist_lengths,
            min_thickness,
            mean_thickness,
            median_thickness,
            max_thickness,
            std_thickness,
            entropy_thickness,
            min_volume,
            mean_volume,
            median_volume,
            max_volume,
            std_volume,
            entropy_volume,
            min_curvature,
            max_curvature,
            mean_curvature,
            median_curvature,
            std_curvature,
            entropy_curvature,
            hist_curvature,
            min_torsion,
            max_torsion,
            mean_torsion,
            median_torsion,
            std_torsion,
            entropy_torsion,
            #
            total_ext_length_P,
            min_length_P,
            mean_length_P,
            median_length_P,
            max_length_P,
            std_lengths_P,
            entropy_lengths_P,
            hist_lengths_P,
            min_thickness_P,
            mean_thickness_P,
            median_thickness_P,
            max_thickness_P,
            std_thickness_P,
            entropy_thickness_P,
            min_volume_P,
            mean_volume_P,
            median_volume_P,
            max_volume_P,
            std_volume_P,
            entropy_volume_P,
            min_curvature_P,
            max_curvature_P,
            mean_curvature_P,
            median_curvature_P,
            std_curvature_P,
            entropy_curvature_P,
            hist_curvature_P,
            min_torsion_P,
            max_torsion_P,
            mean_torsion_P,
            median_torsion_P,
            std_torsion_P,
            entropy_torsion_P,
            #
            total_ext_length_S,
            min_length_S,
            mean_length_S,
            median_length_S,
            max_length_S,
            std_lengths_S,
            entropy_lengths_S,
            hist_lengths_S,
            min_thickness_S,
            mean_thickness_S,
            median_thickness_S,
            max_thickness_S,
            std_thickness_S,
            entropy_thickness_S,
            min_volume_S,
            mean_volume_S,
            median_volume_S,
            max_volume_S,
            std_volume_S,
            entropy_volume_S,
            min_curvature_S,
            max_curvature_S,
            mean_curvature_S,
            median_curvature_S,
            std_curvature_S,
            entropy_curvature_S,
            hist_curvature_S,
            min_torsion_S,
            max_torsion_S,
            mean_torsion_S,
            median_torsion_S,
            std_torsion_S,
            entropy_torsion_S,
        ]

    # Convert the dictionary into pandas dataframe
    features_df = pd_.DataFrame.from_dict(
        somas_features_dict, orient="index", columns=columns
    )
DEBREUVE Eric's avatar
DEBREUVE Eric committed

    return features_df