nutrimorph.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.

# Time profiling:
#     python -m cProfile -o runtime/profiling.log -s name main.py
# Memory profiling:
#     python -m memory_profiler main.py
#     or
#     mprof run main.py
#     mprof plot


import brick.component.connection as cn_
import brick.component.extension as xt_
import brick.component.soma as sm_
import brick.general.feedback as fb_
import brick.processing.input as in_
import brick.processing.plot as po_
import brick.processing.map_labeling as ml_
# import brick.processing.image_verification as iv_
from sklgraph.skl_fgraph import skl_graph_t
from sklgraph.skl_map import skl_map_t
import brick.processing.graph_feat_extraction as ge_
import features_analysis as fa_

import os as os_
import sys as sy_
import time as tm_
import pathlib
import psutil as mu_

import matplotlib.pyplot as pl_
import numpy as np_
import imageio as io_
import skimage.morphology as mp_
import skimage.measure as ms_
import skimage.exposure as ex_
import pandas as pd_
import networkx as nx_

print(sy_.argv, sy_.argv.__len__())

if sy_.argv.__len__() < 2:
    print("Missing parameter file argument")
    sy_.exit(0)
if not (os_.path.isfile(sy_.argv[1]) and os_.access(sy_.argv[1], os_.R_OK)):
    print("Wrong parameter file path or parameter file unreadable")
    sy_.exit(0)

data_path = None
channel = None
save_images = None
save_csv = None
dilatation_erosion = None
size_voxel_in_micron = None
statistical_analysis = None
soma_low_c = None
soma_high_c = None
soma_selem_micron_c = None
soma_min_area_c = None
soma_max_area_c = None
adapt_hist_equalization = None
ext_low_c = None
ext_high_c = None
ext_selem_micron_c = None
ext_min_area_c = None
ext_max_length_c = None
max_straight_sq_dist_c = None
max_weighted_length_c = None
scale_range = None
scale_step = None
alpha = None
beta = None
frangi_c = None
diff_mode = None
bright_on_dark = None
method = None
hist_min_length = None
hist_step_length = None
number_of_bins_length = None
hist_bins_borders_length = None
hist_min_curvature = None
hist_step_curvature = None
number_of_bins_curvature = None
hist_bins_borders_curvature = None
with_plot = None
in_parallel = None

exec(open(sy_.argv[1]).read())  # Only with search_parameters.py (stable version)


def Test_Tangency_ext_conn(soma):
    pb = []
    for extension in soma.extensions:
        for site in list(zip(*extension.sites)):
            close_end_pt = tuple(
                (site[0] + i, site[1] + j, site[2] + k)
                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)
            for ext in tuple(zip(soma.connection_path.values())):
                if ext[0] is not None:
                    for conn in list(zip(*np_.asarray(ext[0]))):
                        if (conn in close_end_pt) and (conn not in list(zip(*extension.sites))):
                            for extens in soma.extensions:
                                if conn not in zip(*extens.end_points):
                                    if conn not in pb:
                                        pb.append(conn)
                                        print("\nconn: ", conn, "\next_site: ", site, "\next_uid: ", extension.uid)


#########----------- Executable version with files AND repositories ------------###########

def NutriMorph(data_path: str,
               channel=channel,
               dilatation_erosion=dilatation_erosion,
               size_voxel_in_micron=size_voxel_in_micron,
               soma_low_c=soma_low_c,
               soma_high_c=soma_high_c,
               soma_selem_micron_c=soma_selem_micron_c,
               soma_min_area_c=soma_min_area_c,
               soma_max_area_c=soma_max_area_c,
               adapt_hist_equalization=adapt_hist_equalization,
               ext_low_c=ext_low_c,
               ext_high_c=ext_high_c,
               ext_selem_micron_c=ext_selem_micron_c,
               ext_min_area_c=ext_min_area_c,
               max_straight_sq_dist_c=max_straight_sq_dist_c,
               max_weighted_length_c=max_weighted_length_c,
               scale_range=scale_range,
               scale_step=scale_step,
               alpha=alpha,
               beta=beta,
               frangi_c=frangi_c,
               diff_mode=diff_mode,
               bright_on_dark=bright_on_dark,
               method=method,
               hist_min_length=hist_min_length,
               hist_step_length=hist_step_length,
               number_of_bins_length=number_of_bins_length,
               hist_bins_borders_length=hist_bins_borders_length,
               hist_min_curvature=hist_min_curvature,
               hist_step_curvature=hist_step_curvature,
               number_of_bins_curvature=number_of_bins_curvature,
               hist_bins_borders_curvature=hist_bins_borders_curvature,
               with_plot=with_plot,
               in_parallel=in_parallel) -> pd_.DataFrame():

    soma_t = sm_.soma_t
    extension_t = xt_.extension_t

    print(f"STARTED: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}")
    start_time = tm_.time()

    # TODO add proper warning with warn module
    # --- Images

    # Find the name of the file, eq. to the biological tested condition number x
    name_file = os_.path.basename(data_path)
    print("FILE: ", name_file)
    name_file = name_file.replace(".tif", "")

    name_dir = os_.path.dirname(data_path)
    print("DIR: ", name_dir)

    # Find the dimension of the image voxels in micron
    # TODO do something more general, here too specific to one microscope metadata
    size_voxel_in_micron = in_.FindVoxelDimensionInMicron(data_path, size_voxel_in_micron=size_voxel_in_micron)

    # Read the image
    image = io_.volread(data_path)

    # Image size verification - simple version without user interface
    image = in_.ImageVerification(image, channel)

    # iv_.image_verification(image, channel)  # -> PySide2 user interface  # TODO: must return the modified image!
    # /!\ conflicts between some versions of PySide2 and Python3

    # image = image[:, 800:, 300:]
    img_shape = image.shape
    #
    print(f"IMAGE: s.{img_shape} t.{image.dtype} m.{image.min()} M.{image.max()}")

    # Intensity relative normalization (between 0 and 1).
    # Image for soma normalized for hysteresis
    image_for_soma = in_.IntensityNormalizedImage(image)
    # Image for extensions not normalized for frangi enhancement
    image_for_ext = image.copy()

    print(f"NRM-IMG: t.{image_for_soma.dtype} m.{image_for_soma.min():.2f} M.{image_for_soma.max():.2f}")

    # --- Initialization
    som_nfo = {}
    ext_nfo = {}
    axes = {}

    #
    # --- Somas
    print("\n--- Soma Detection")

    # Change the soma parameters from micron to pixel dimensions, using the previously determined voxel size
    soma_min_area_c = in_.ToPixel(soma_min_area_c, size_voxel_in_micron, dimension=(0, 1))
    soma_max_area_c = in_.ToPixel(soma_max_area_c, size_voxel_in_micron, dimension=(0, 1))
    soma_selem_c = mp_.disk(in_.ToPixel(soma_selem_micron_c, size_voxel_in_micron))

    # - Create the maps for enhancing and detecting the soma
    # Hysteresis thresholding and closing/opening
    print("Hysteresis thresholding: ", end="")
    som_nfo["map_small"] = soma_t.Map(image_for_soma, soma_low_c, soma_high_c, soma_selem_c)
    print("Done")
    # Deletion of too small elements
    print("Morphological cleaning: ", end="")
    som_nfo["map_small"], som_lmp_small = soma_t.DeleteSmallSomas(som_nfo["map_small"], soma_min_area_c)
    # Deletion of too big elements (probably aggregation of cells)
    # Separated from the deletion of small elements to avoid extensions to be detected where too big somas were just deleted
    som_nfo["map"], som_lmp = soma_t.DeleteBigSomas(som_nfo["map_small"], som_lmp_small, soma_max_area_c)
    print("Done")

    # Labelling of somas, find the number of somas
    print("Somas labelling: ", end="")
    som_nfo["lmp"], n_somas = ms_.label(som_nfo["map"], return_num=True)
    print("Done")

    # Use relabel instead of label to optimize the algorithm. /!\ But BUG.
    # som_nfo["lmp"] = relabel_sequential(som_lmp)[0]
    # n_somas = som_nfo["lmp"].max()

    # Create the distance and influence maps used next for connexions between somas-extensions and extensions-extensions
    print("Distance and Influence Map creation: ", end="")
    som_nfo["dist_to_closest"], som_nfo["influence_map"] = soma_t.InfluenceMaps(som_nfo["lmp"])
    print("Done")

    # Create the tuple somas, containing all the objects 'soma'
    print("Somas creation: ", end="")
    somas = tuple(soma_t().FromMap(som_nfo["lmp"], uid) for uid in range(1, n_somas + 1))
    print("Done")

    print(f"    n = {n_somas}")

    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")

    if with_plot:
        fb_.PlotSomas(somas, som_nfo, axes)

    #
    # -- Extentions
    print("\n--- Extension Detection")

    # Set to zeros the pixels belonging to the somas.
    # Take into account the aggregation of cells detected as one big soma to avoid extension creation there
    del_soma = (som_nfo["map_small"] > 0).nonzero()
    image_for_ext[del_soma] = 0

    # Change the extensions parameters from micron to pixel dimensions
    ext_min_area_c = in_.ToPixel(ext_min_area_c, size_voxel_in_micron, dimension=(0, 1))

    if ext_selem_micron_c == 0:
        ext_selem_pixel_c = None
    else:
        ext_selem_pixel_c = mp_.disk(in_.ToPixel(ext_selem_micron_c, size_voxel_in_micron))

    # - Perform frangi feature enhancement (via python or c - faster - implementation)
    enhanced_ext, ext_scales = extension_t.EnhancedForDetection(
        image_for_ext,
        scale_range,
        scale_step,
        alpha,
        beta,
        frangi_c,
        bright_on_dark,
        method,
        diff_mode=diff_mode,
        in_parallel=in_parallel)

    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"Elapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}\n")

    # - Creation of the enhanced maps
    # Enhanced the contrast in frangi output
    print('Enhancing Contrast: ', end='')
    if adapt_hist_equalization:
        # necessary but not sufficient, histogram equalisation (adaptative)
        enhanced_ext = ex_.equalize_adapthist(enhanced_ext)
    else:
        # necessary but not sufficient, histogram equalisation (global)
        enhanced_ext = ex_.equalize_hist(enhanced_ext)

    # Rescale the image by stretching the intensities between 0 and 255, necessary but not sufficient
    enhanced_ext = ex_.rescale_intensity(enhanced_ext, out_range=(0, 255))
    print('Done')
    # Hysteresis thersholding
    print('Hysteresis Thresholding: ', end='')
    ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c, ext_selem_pixel_c)
    print('Done')
    # Deletion of too small elements
    print('Morphological cleaning: ', end='')
    ext_nfo["coarse_map"], ext_lmp = extension_t.FilteredCoarseMap(ext_nfo["coarse_map"], ext_min_area_c)
    print('Done')
    # Creation of the extensions skeleton
    print('Skeletonization and Thinning: ', end='')
    ext_nfo["map"] = extension_t.FineMapFromCoarseMap(ext_nfo["coarse_map"])
    print('Done')
    # Deletion of extensions found within the somas
    print('Map Labelling: ', end='')
    ext_nfo["map"][som_nfo["map"] > 0] = 0
    # Labelling of extensions and number of extensions determined
    ext_nfo["lmp"], n_extensions = ms_.label(ext_nfo["map"], return_num=True)
    print('Done')

    # Use relabel instead of label to optimize the algorithm. BUT PROBLEM WITH THE NUMBER OF EXTENSIONS DETECTED !
    # ext_nfo["lmp"] = relabel_sequential(ext_lmp)[0]
    # n_extensions = ext_nfo["lmp"].max()

    # Create the tuple extensions, containing all the objects 'extension'
    print('Extensions creation: ', end='')
    extensions = tuple(
        extension_t().FromMap(ext_nfo["lmp"], ext_scales, uid)
        for uid in range(1, n_extensions + 1))
    print('Done')
    print(f"    n = {n_extensions}")

    # Create global end point map for extensions
    glob_ext_ep_map = xt_.EndPointMap(ext_nfo["lmp"] > 0)
    all_end_points = list(zip(*glob_ext_ep_map.nonzero()))  # updated in ValidateConnexion

    #
    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")

    if with_plot:
        fb_.PlotExtensions(extensions, ext_nfo, img_shape)

    # -- Preparation for Connections
    # Calculate the COSTS of each pixel on the image and DILATE the existing extensions to
    # avoid tangency of connexions to extensions
    dijkstra_costs = in_.DijkstraCosts(image, som_nfo["map"], ext_nfo["map"])

    if dilatation_erosion:
        # Dilate all extensions
        for extension in extensions:
            cn_.Dilate(extension.sites, dijkstra_costs)

    # -- Soma-Extention
    print("\n--- Soma <-> Extension")

    # Change the extensions parameters from micron to pixel dimensions
    max_straight_sq_dist_c = in_.ToPixel(max_straight_sq_dist_c, size_voxel_in_micron)
    max_weighted_length_c = in_.ToPixel(max_weighted_length_c, size_voxel_in_micron)

    # Find the candidate extensions, with their endpoints, for connexion to the somas
    candidate_conn_nfo = cn_.CandidateConnections(
        somas,
        som_nfo["influence_map"],
        som_nfo["dist_to_closest"],
        extensions,
        max_straight_sq_dist_c,
    )

    # For each candidate, verify that it is valid based on the shortest path and the maximum allowed straight distance
    for ep_idx, soma, extension, end_point in candidate_conn_nfo:
        # Only try to connect if the extension is not already connected
        if extension.is_unconnected:
            print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
            if dilatation_erosion:
                # Erode the end_point of the extension in the costs map
                cn_.Erode((end_point,), dijkstra_costs, extension, extensions, image, all_end_points=all_end_points)

            # Use of Dijkstra shortest weighted path
            path, length = cn_.ShortestPathFromToN(
                image=image,
                point=end_point,
                extension=extension,
                costs=dijkstra_costs,
                candidate_points_fct=soma.ContourPointsCloseTo,
                max_straight_sq_dist=max_straight_sq_dist_c,
                erode_path=False,
            )

            # Keep the connexion only if inferior to the allowed max weighted distance
            if path.__len__() > 0:
                # length in pixel and max_weighted_length_c too
                if length <= max_weighted_length_c:
                    # Validate and update all the fields + dilate again the whole extension
                    cn_.ValidateConnection(soma, extension, end_point, path, dijkstra_costs, all_ep=all_end_points)
                    if dilatation_erosion:
                        # Dilate extensions
                        cn_.Dilate(extension.sites, dijkstra_costs)
                    print(": Made")
                else:
                    if dilatation_erosion:
                        cn_.Dilate(extension.sites, dijkstra_costs)
                    print("")
            else:
                if dilatation_erosion:
                    cn_.Dilate(extension.sites, dijkstra_costs)
                print("")

    fb_.PrintConnectedExtensions(extensions)

    #
    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")

    if with_plot:
        fb_.PlotSomasWithExtensions(somas, som_nfo, "all")

    # -- Extention-Extention
    print("\n--- Extension <-> Extension")

    should_look_for_connections = True

    # Update the soma maps with next Ext-Ext connexions
    while should_look_for_connections:

        if dilatation_erosion:
            # Update the costs by dilating everything again plus the contour points of the soma
            for soma in somas:
                cn_.Dilate(soma.contour_points, dijkstra_costs)
                # Below probably not necessary but security
                for extension in soma.extensions:
                    cn_.Dilate(extension.sites, dijkstra_costs)

        # Update the final somas + ext map by adding the new extensions and their connexion paths
        som_nfo["soma_w_ext_lmp"] = soma_t.SomasLMap(somas)

        # Update the influence map
        som_nfo["dist_to_closest"], som_nfo["influence_map"] = soma_t.InfluenceMaps(som_nfo["soma_w_ext_lmp"])

        # Find the candidate extensions for connexion to the primary extensions
        candidate_conn_nfo = cn_.CandidateConnections(
            somas,
            som_nfo["influence_map"],
            som_nfo["dist_to_closest"],
            extensions,
            max_straight_sq_dist_c,
        )

        should_look_for_connections = False

        # For each candidate, verify if valid based on the shortest path and the maximum allowed straight distance
        for ep_idx, soma, extension, end_point in candidate_conn_nfo:
            # Only try to connect if the extension is not already connected
            if extension.is_unconnected:
                print(f"    Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
                if dilatation_erosion:
                    # Erode the end_point of the extension in the costs map
                    cn_.Erode((end_point,), dijkstra_costs, extension, extensions, image, all_end_points=all_end_points)
                # Use of Dijkstra shortest weighted path
                path, length = cn_.ShortestPathFromToN(
                    image=image,
                    point=end_point,
                    extension=extension,
                    costs=dijkstra_costs,
                    candidate_points_fct=soma.ExtensionPointsCloseTo,
                    extensions=extensions,
                    all_end_points=all_end_points,
                    max_straight_sq_dist=max_straight_sq_dist_c,
                    erode_path=dilatation_erosion,
                )

                # Keep the connexion only if inferior to the allowed max weighted distance
                if path.__len__() > 0:
                    if length <= max_weighted_length_c:
                        # Verify that the last element of the connexion path is contained in the extension to be connected
                        # If so, returns the extensions containing the path last site.
                        tgt_extenstion = extension_t.ExtensionContainingSite(extensions, path[-1])
                        # Validate connexion: update the glial component objects and store more info in their fields
                        cn_.ValidateConnection(tgt_extenstion, extension, end_point, path, dijkstra_costs, all_ep=all_end_points)
                        if dilatation_erosion:
                            # Dilate extensions
                            cn_.Dilate(extension.sites, dijkstra_costs)

                        should_look_for_connections = True

                        print(": Made")
                    else:
                        if dilatation_erosion:
                            cn_.Dilate(extension.sites, dijkstra_costs)
                        print("")
                else:
                    if dilatation_erosion:
                        cn_.Dilate(extension.sites, dijkstra_costs)
                    print("")

    # -- Create an extension map containing all the ext + connexions, without the somas.
    # Used for the graphs extraction
    ext_nfo["lmp_soma"] = som_nfo['soma_w_ext_lmp'].copy()
    ext_nfo["lmp_soma"][som_nfo["map"] > 0] = 0
    #

    fb_.PrintConnectedExtensions(extensions)

    #
    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")

    if with_plot:
        fb_.PlotSomasWithExtensions(somas, som_nfo, "with_ext_of_ext")

    # -- Summary
    print("\n")
    for soma in somas:
        print(soma)

    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")

    if with_plot:
        pl_.show()

    if save_images is not None:
        po_.MaximumIntensityProjectionZ(som_nfo['soma_w_ext_lmp'],
                                        block=False,
                                        output_image_file_name=f"{save_images}\\MIP_{name_file}.png")

    # --- Extract all the extensions of all somas as a graph
    print('\n--- Graph extraction')

    # Create the graphs with SKLGraph module (available on Eric Debreuve Gitlab)
    for soma in somas:
        print(f" Soma {soma.uid}", end="")
        # Create SKLGraph skeletonized map
        ext_map = skl_map_t.FromShapeMap(ext_nfo['lmp_soma'] == soma.uid,
                                         store_widths=True,
                                         skeletonize=True,
                                         do_post_thinning=True)

        # Create the graph from the SKLGaph skeletonized map
        soma.skl_graph = skl_graph_t.FromSkeleton(ext_map, end_point, size_voxel=size_voxel_in_micron)

        # --- Find the root of the {ext+conn} graphs.
        # Roots are the nodes of degree 1 that are to be linked to the soma
        soma.graph_roots = ge_.FindGraphsRootWithNodes(soma)

        # Add a node "soma" and link it to the root nodes
        soma_node = f"S-{int(soma.centroid[0])}-{int(soma.centroid[1])}-{int(soma.centroid[2])}"
        soma.skl_graph.add_node(soma_node, soma=True, soma_nfo=soma)

        for node in soma.graph_roots.values():
            soma.skl_graph.add_edge(node, soma_node, root=True)

        if save_images is not None:
            nx_.draw_networkx(soma.skl_graph)
            pl_.savefig(f"{save_images}\\graph_{name_file}_soma{soma.uid}.png")
            # pl_.show(block=True)
            pl_.close()

        print(": Done")

    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")

    # --- Extract features
    print('\n--- Features Extraction\n')

    # Parameters
    if hist_bins_borders_length is None:
        number_of_bins_length = int(number_of_bins_length)
        bins_length = np_.linspace(hist_min_length, hist_min_length + hist_step_length * number_of_bins_length,
                                   num=number_of_bins_length)
        bins_length[-1] = np_.inf
    else:
        bins_length = np_.array(hist_bins_borders_length)
        bins_length[-1] = np_.inf

    if hist_bins_borders_curvature is None:
        number_of_bins_curvature = int(number_of_bins_curvature)
        bins_curvature = np_.linspace(hist_min_curvature,
                                      hist_min_curvature + hist_step_curvature * number_of_bins_curvature,
                                      num=number_of_bins_curvature)
        bins_curvature[-1] = np_.inf
    else:
        bins_curvature = np_.array(hist_bins_borders_curvature)
        bins_curvature[-1] = np_.inf

    # Pandas dataframe creation with all the measured features
    features_df = ge_.ExtractFeaturesInDF(name_file,
                                          somas,
                                          size_voxel_in_micron,
                                          bins_length,
                                          bins_curvature,
                                          ext_scales,
                                          )

    # Save the pandas df into .csv
    if save_csv is not None:
        features_df.to_csv(f"{save_csv}\\{name_file}.csv")
    else:
        features_df.to_csv(f"{name_dir}\\{name_file}.csv")


    elapsed_time = tm_.gmtime(tm_.time() - start_time)
    print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
    print(f"DONE: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}\n")

    return features_df


if __name__ == '__main__':

    # --- Extract cell graphs and features from microscope images using NutriMorph function.
    #
    # Differentiate between path to a tiff file or to a repository
    if pathlib.Path(data_path).is_file():
        # Perform NutriMorph algorithm on the file entered in parameters
        print("WARNING: Will not perform features analysis on a single image.\n For features analysis, "
              "give a directory path.\n")
        features_df = NutriMorph(data_path)

    elif pathlib.Path(data_path).is_dir():
        # Keep the directory to the repository
        name_dir = os_.path.dirname(data_path)
        # Initialize the future concatenated features
        concatenated_features_df = pd_.DataFrame()
        # Find all the tiff files in the parent and child repositories
        for path in pathlib.Path(data_path).glob("**/*.tif"):
            if path.is_file():
                name_file = os_.path.basename(path)
                # Perform NutriMorph algorithm
                features_df = NutriMorph(path)
                # Concatenate all the dataframes
                concatenated_features_df = concatenated_features_df.append(features_df)
                # If some rows (ie. somas0) have NaN, drop them
                # -- due to best fitting ellipsoid algo (JTJ is singular due to convex hull being flat)
                concatenated_features_df = concatenated_features_df.dropna(axis=0, how="any")

        # Save to .csv in the parent repository
        if save_csv is None:
            concatenated_features_df.to_csv(f"{data_path}\\features.csv")
        else:
            concatenated_features_df.to_csv(f"{save_csv}\\features.csv")

        # Clustering with this df and module features_analysis.py
        if statistical_analysis:
            os_.system(f"feature_analysis.py {save_csv}\\features.csv")

    else:
        raise ImportError("Not a valid data path!")