# 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!")