diff --git a/brick/component/extension.py b/brick/component/extension.py index dc6c5562b8c257aca62e70ff9cd99ce734dd099b..ef55dd1b3e533a2cdcd8b11e63276ba64e93b1db 100644 --- a/brick/component/extension.py +++ b/brick/component/extension.py @@ -33,6 +33,7 @@ from __future__ import annotations import brick.processing.frangi3 as fg_ import brick.processing.map_labeling as ml_ +import brick.processing.input as in_ from brick.component.glial_cmp import glial_cmp_t from brick.general.type import array_t, site_h @@ -45,7 +46,8 @@ import skimage.morphology as mp_ from scipy import ndimage as im_ -min_area_c = 100 +ext_min_area_c = 100 +ext_selem_microns_c = 0.24050024 class extension_t(glial_cmp_t): @@ -174,10 +176,12 @@ class extension_t(glial_cmp_t): return enhanced_img, scale_map @staticmethod - def CoarseMap(image: array_t, low: float, high: float) -> array_t: + def CoarseMap(image: array_t, low: float, high: float, voxel_micron: array_t, ext_selem_microns=ext_selem_microns_c) -> array_t: # result = __HysterisisImage__(image, low, high) - result = __MorphologicalCleaning__(result) + + selem = mp_.disk(in_.ToPixels(ext_selem_microns, voxel_micron)) + result = __MorphologicalCleaning__(result, selem) return result @@ -188,7 +192,7 @@ class extension_t(glial_cmp_t): lmp = ms_.label(map_) for region in ms_.regionprops(lmp): - if region.area <= min_area_c: + if region.area <= ext_min_area_c: region_sites = (lmp == region.label).nonzero() result[region_sites] = 0 lmp[region_sites] = 0 @@ -232,11 +236,10 @@ def __HysterisisImage__(image: array_t, low: float, high: float) -> array_t: return result -def __MorphologicalCleaning__(image: array_t) -> array_t: +def __MorphologicalCleaning__(image: array_t, selem) -> array_t: # result = image.copy() - selem = mp_.disk(1) for dep in range(result.shape[0]): result[dep, :, :] = mp_.closing(result[dep, :, :], selem) result[dep, :, :] = mp_.opening(result[dep, :, :], selem) diff --git a/brick/general/type.py b/brick/general/type.py index a1820df690879682891079a5cc8c8fe54a38e8db..c9292db2066a7d4e22ee6cbebe58b8b4065089ad 100644 --- a/brick/general/type.py +++ b/brick/general/type.py @@ -29,12 +29,13 @@ # 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 typing import Tuple, Union +from typing import Tuple, Union, Type import numpy as np_ +# from Cython.Includes.numpy import ndarray +from numpy.core._multiarray_umath import ndarray - -array_t = np_.ndarray +array_t: Type[Union[ndarray, ndarray]] = np_.ndarray tintegers_h = Tuple[int, ...] diff --git a/brick/processing/input.py b/brick/processing/input.py index 56a081f7c94f9c49f86bd2a6ff35369f2e8d48db..cbe0996d43593cff370cc2c43fe67c0e241958d1 100644 --- a/brick/processing/input.py +++ b/brick/processing/input.py @@ -123,6 +123,19 @@ def VoxelDimensionInMicrons(data_path: str) -> array_t: return voxel_size_microns +def ToPixels(micron: float, voxel_size_microns: array_t) -> int: + # Conversion of pixels into microns. Used in morphomath (disk structuring element mp_.disk(n_pixel)). + # Assumes that the axes X and Y have the same ratio pixel:microns. + return round(micron/voxel_size_microns[0]) + + +def ToMicrons(pixel: int, voxel_size_microns: array_t) -> float: + # May not be used. if not used, delete. + # Conversion of microns into pixels. + # Assumes that the axes X and Y have the same ratio pixel:microns. + return float(pixel * voxel_size_microns[0]) + + def DijkstraCosts(image: array_t, som_map: array_t, ext_map: array_t) -> array_t: # # TODO: Ideally, the extension part should be dilated diff --git a/nutrimorph.py b/nutrimorph.py index 52d993cc9f79bfa0dc56ae930b038f6b0bc719e5..c51aaff421bdeb708f0a4817c84550ff0153f2b2 100644 --- a/nutrimorph.py +++ b/nutrimorph.py @@ -52,9 +52,8 @@ import time as tm_ import matplotlib.pyplot as pl_ import numpy as np_ import skimage.io as io_ -import skimage.measure as ms_ from skimage.segmentation import relabel_sequential -import exifread as xf_ +import skimage.morphology as mp_ print(sy_.argv, sy_.argv.__len__()) @@ -72,10 +71,10 @@ soma_low_c = None soma_high_c = None ext_low_c = None ext_high_c = None -soma_selem_c = None +soma_selem_microns_c = None max_straight_sq_dist_c = None max_weighted_length_c = None -min_area_c = None +soma_min_area_c = None in_parallel = None exec(open(sy_.argv[1]).read()) @@ -110,6 +109,7 @@ 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 = in_.IntensityNormalizedImage(image) image_for_ext = in_.IntensityNormalizedImage(image) print(f"NRM-IMG: t.{image_for_soma.dtype} m.{image_for_soma.min():.2f} M.{image_for_soma.max():.2f}") @@ -128,8 +128,10 @@ axes = {} # --- Somas print("\n--- Soma Detection") +soma_selem_c = mp_.disk(in_.ToPixels(soma_selem_microns_c, voxel_micron)) + som_nfo["map"] = soma_t.Map(image_for_soma, soma_low_c, soma_high_c, soma_selem_c) -som_nfo["map"], som_lmp = soma_t.FilteredMap(som_nfo["map"], min_area_c) +som_nfo["map"], som_lmp = soma_t.FilteredMap(som_nfo["map"], soma_min_area_c) # som_nfo["lmp"], n_somas= ms_.label(som_nfo["map"], return_num=True) som_nfo["lmp"] = relabel_sequential(som_lmp)[0] # Use relabel instead of label to optimize the algorithm. n_somas = som_nfo["lmp"].max() @@ -152,7 +154,7 @@ if with_plot: print("\n--- Extension Detection") enhanced_ext, ext_scales = extension_t.EnhancedForDetection(image_for_ext, in_parallel=in_parallel) -ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c) +ext_nfo["coarse_map"] = extension_t.CoarseMap(enhanced_ext, ext_low_c, ext_high_c, voxel_micron) ext_nfo["coarse_map"], ext_lmp = extension_t.FilteredCoarseMap(ext_nfo["coarse_map"]) ext_nfo["map"] = extension_t.FineMapFromCoarseMap(ext_nfo["coarse_map"]) ext_nfo["map"][som_nfo["map"] > 0] = 0 diff --git a/parameters.py b/parameters.py index d6af36fcf792fb031cbe7334589b1e816fad428c..ee09ffe1e2c7a6ee1bdb57403850d5e0796c2a4f 100644 --- a/parameters.py +++ b/parameters.py @@ -29,8 +29,6 @@ # 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 skimage.morphology as mp_ - data_path = "./data/DIO_6H_6_1.70bis_2.2_3.tif" channel = 'G' @@ -41,11 +39,15 @@ soma_high_c = 0.7126 ext_low_c = 0.2 # 0.02 # 0.2 # ext_low_c = 9.0e-4 ext_high_c = 0.6 # 0.04 # 0.6 # high_ext = 8.0e-3 -soma_selem_c = mp_.disk(2) +# soma_selem_c = mp_.disk(2) +soma_selem_microns_c = 0.48100048 +# ext_selem_c = mp_.disk(1) +ext_selem_microns_c = 0.24050024 # /!\ Not directly used in brick/component #TODO direct use of param ? max_straight_sq_dist_c = 30 ** 2 max_weighted_length_c = 20.0 -min_area_c = 1000 +soma_min_area_c = 1000 +ext_min_area_c = 100 # /!\ Not directly used in brick/component #TODO direct use of param ? in_parallel = False