From 8a5e0c7422e4f460a2e4700c06124e846f8e194c Mon Sep 17 00:00:00 2001 From: msorelli Date: Thu, 14 Nov 2024 17:20:23 +0100 Subject: [PATCH] Simplify code --- doc/source/usage.rst | 58 ++--- foa3d/__main__.py | 11 +- foa3d/frangi.py | 536 ++++++++++++++++++--------------------- foa3d/input.py | 551 +++++++++++++++++------------------------ foa3d/odf.py | 189 ++++++-------- foa3d/output.py | 91 ++++--- foa3d/pipeline.py | 535 ++++++++++++++++++++------------------- foa3d/preprocessing.py | 75 +++--- foa3d/printing.py | 169 ++++++------- foa3d/slicing.py | 521 +++++++++++++++++++++----------------- foa3d/spharm.py | 10 + foa3d/utils.py | 396 ++++++++++++++--------------- requirements.in | 2 +- requirements.txt | 123 +++++---- 14 files changed, 1578 insertions(+), 1689 deletions(-) diff --git a/doc/source/usage.rst b/doc/source/usage.rst index fb009a8..bd192ae 100644 --- a/doc/source/usage.rst +++ b/doc/source/usage.rst @@ -41,17 +41,15 @@ Usage Microscopy image formats ------------------------ -Foa3D supports 3D grayscale or RGB image stacks in TIF/TIFF or NumPy format. -Alternatively, a .yml stitch file generated by the ZetaStitcher tool for large volumetric stack alignment and stitching -[`ZetaStitcher GitHub `_] -may be also provided as input. This .yml file can be generated following the detailed documentation available at -[`ZetaStitcher Documentation `_] -from a collection of 3D stacks composing a tiled reconstruction of a brain tissue sample. -In detail, the Foa3D tool employs the 3D stack alignment information included in such file -to programmatically access and process basic image sub-volumes of suitable size, -thus enabling the analysis of high-resolution mesoscopic microscopy images -(e.g., 10¹ - 10³ GB) which exceed the typical memory available on low-resource machines. -The .yml and the image stack files must be located within the same directory. +Foa3D accepts 3D grayscale or RGB image stacks in TIFF format. +Alternatively, a YAML stitch file created by the ZetaStitcher tool for large volumetric stack alignment +and stitching [`ZetaStitcher GitHub `_] can be used as input. +To generate this stitch file from a collection of adjacent 3D stacks composing a tiled reconstruction of brain tissue, +refer to the documentation at [`ZetaStitcher GitHub `_]. +In particular, the Foa3D tool uses the 3D stack alignment information contained in such a file to programmatically +access and process basic image sub-volumes of appropriate size, allowing for the analysis of high-resolution mesoscopic +microscopy images that exceed the typical memory available on low-resource machines. +The YAML and image stack files must be in the same directory. .. code-block:: console @@ -61,8 +59,8 @@ The .yml and the image stack files must be located within the same directory. Microscopy image resolution --------------------------- -The lateral and longitudinal voxel size (in μm) must be specified via CLI, -along with the full width at half maximum of the point spread function of the employed microscopy apparatus: +The lateral and longitudinal voxel size (in μm) must be specified via the command line, +along with the 3D full width at half maximum of the point spread function of the employed microscopy apparatus: .. code-block:: console @@ -71,7 +69,7 @@ along with the full width at half maximum of the point spread function of the em This information is required at the preprocessing stage of the pipeline to properly isotropize the spatial resolution of the raw microscopy images. In detail, since two-photon scanning and light-sheet fluorescence microscopes are in general characterized by a poorer resolution along the direction of the optical axis, the XY-plane of the sliced -image sub-volumes tipically needs to be blurred. A tailored Gaussian smoothing kernel is used in this regard. +image sub-volumes typically needs to be blurred. A tailored Gaussian smoothing kernel is used in this regard. If not properly corrected, the residual anisotropy would otherwise introduce a systematic bias in the assessed 3D fiber orientations. @@ -82,7 +80,7 @@ Frangi filter configuration Fiber enhancement and masking is achieved via a multiscale 3D Frangi filter [`Frangi, et al., 1998 `_]. The spatial scales of the filter (in μm) can be provided via the ``-s/--scales`` option. As discussed in [`Sorelli, et al., 2023 `_], -the optimal scales which best preserve the original intensity +the optimal scales that best preserve the original intensity and cross-sectional size of the 3D tubular structures present in the analized images correspond to half of their expected radius. The response of the Frangi filter is also affected by three sensitivity parameters, α, β, and γ. @@ -103,16 +101,20 @@ cross-sectional diameter of 5 and 10 μm, with an automatic (local) contrast sen $ ... -a 0.00001 -b 0.1 -s 1.25 2.5 +Please keep in mind that the above automatic local adjustment of the γ sensitivity may produce discontinuities +between the fiber orientation vector fields resulting from adjacent image slices. + .. _parallelization: Parallelization --------------- -In order to speed up the fiber orientation analysis on large brain tissue sections, the Foa3D pipeline divides -the input image reconstruction into basic slices of suitable shape, and feeds them to separate concurrent workers. -By default, Foa3D will use all available logical cores, splitting the multiscale fiber enhancement among parallel -threads - e.g., 16 image slices will be simultaneously processed at 2 spatial scales of interest on a 32-core CPU. -The size of these slices is automatically set depending on the currently available RAM. -The ``--job`` and ``--ram`` options may be otherwise specified via CLI in order to limit the employed resources: +In order to speed up the fiber orientation analysis on large brain tissue sections, the Foa3D pipeline divides the input +image reconstruction into basic slices of appropriate shape and assigns them to separate concurrent workers. +By default, Foa3D will use all available logical cores—for instance, a batch of 32 image slices will be simultaneously +processed on a 32-core CPU. The multiscale Frangi filter, on the other hand, is not currently parallelized in order to +avoid unnecessary overhead and resource oversubscription within the nested parallel loop. The size of the basic image +slices is automatically determined by the available RAM. The ``--job`` and ``--ram`` options can be specified +differently via the command line to limit the employed resources. .. code-block:: console @@ -122,13 +124,11 @@ The ``--job`` and ``--ram`` options may be otherwise specified via CLI in order Soma rejection -------------- -A neuronal soma fluorescence channel may be optionally provided to Foa3D, -in order to improve the specificity of the resulting fiber orientation maps -achieved thanks to the inherent attenuation of non-tubular objects offered by the Frangi filter. -This is performed via a postprocessing step which further suppresses neuronal bodies -applying Yen's automatic thresholding algorithm to an optionally provided channel. -The enhanced neuronal body rejection may be activated via the ``-c/--cell-msk`` option modifying, -if required, the default channel related to the soma fluorescence: +A neuronal soma fluorescence channel may be optionally provided to Foa3D for improving the specificity of the resulting +fiber orientation maps, which is otherwise dependent on the inherent attenuation of non-tubular objects offered by the +Frangi filter. This is performed via a postprocessing step that further suppresses neuronal bodies by applying Yen's +automatic thresholding algorithm to an optionally provided channel. The enhanced neuronal body rejection may be +activated via the ``-c/--cell-msk`` option, modifying, if required, the default channel related to the soma fluorescence: .. code-block:: console @@ -191,7 +191,7 @@ by selecting the "export all" option (-e or --exp-all) via CLI. #. Orientation distribution functions (ODF) stage (one file for each super-voxel size requested via CLI): - * **ODF** (*path/to/save_dir/odf/odf_mrtrixview_\*cfg_sfx\, type: float32, format: NIfTI) + * **ODF** (*path/to/save_dir/odf/odf_mrtrixview_\*cfg_sfx\**, type: float32, format: NIfTI) * **ODF background** (*path/to/save_dir/odf/bg_mrtrixview_\*cfg_sfx\**, type: uint8, format: NIfTI) diff --git a/foa3d/__main__.py b/foa3d/__main__.py index ff736b9..0fd30a3 100644 --- a/foa3d/__main__.py +++ b/foa3d/__main__.py @@ -1,7 +1,7 @@ from foa3d.input import get_cli_parser, load_microscopy_image -from foa3d.pipeline import parallel_odf_over_scales, parallel_frangi_over_slices +from foa3d.pipeline import parallel_frangi_over_slices, parallel_odf_over_scales from foa3d.printing import print_pipeline_heading -from foa3d.utils import delete_tmp_folder +from foa3d.utils import delete_tmp_data def foa3d(cli_args): @@ -12,15 +12,14 @@ def foa3d(cli_args): # parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices out_img = parallel_frangi_over_slices(cli_args, save_dirs, in_img) - # generate 3D fiber ODF maps over the spatial scales of interest using concurrent workers - parallel_odf_over_scales(cli_args, save_dirs, out_img['vec'], out_img['iso'], out_img['px_sz'], in_img['name']) + # generate 3D fiber ODF maps at the spatial scales of interest using concurrent workers + parallel_odf_over_scales(cli_args, save_dirs, out_img, in_img['name']) # delete temporary folder - delete_tmp_folder(save_dirs['tmp'], (in_img, out_img)) + delete_tmp_data(save_dirs['tmp'], (in_img, out_img)) def main(): - # start Foa3D by terminal print_pipeline_heading() foa3d(cli_args=get_cli_parser()) diff --git a/foa3d/frangi.py b/foa3d/frangi.py index 39d9bcc..90aeba3 100644 --- a/foa3d/frangi.py +++ b/foa3d/frangi.py @@ -2,7 +2,6 @@ from itertools import combinations_with_replacement import numpy as np -from joblib import Parallel, delayed from scipy import ndimage as ndi from foa3d.utils import (create_background_mask, create_memory_map, @@ -11,9 +10,9 @@ def analyze_hessian_eigen(img, sigma, trunc=4): """ - Return the eigenvalues of the local Hessian matrices + Compute the eigenvalues of local Hessian matrices of the input image array, sorted by absolute value (in ascending order), - along with the related eigenvectors. + along with the related (dominant) eigenvectors. Parameters ---------- @@ -28,26 +27,22 @@ def analyze_hessian_eigen(img, sigma, trunc=4): Returns ------- - eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + eigval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) Hessian eigenvalues sorted by absolute value (ascending order) - dom_eigenvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + dom_eigvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) Hessian eigenvectors related to the dominant (minimum) eigenvalue """ - # compute scaled Hessian matrices hessian = compute_scaled_hessian(img, sigma=sigma, trunc=trunc) + eigval, dom_eigvec = compute_dominant_eigen(hessian) - # compute Hessian eigenvalues and orientation vectors - # related to the dominant one - eigenval, dom_eigenvec = compute_dominant_eigen(hessian) - - return eigenval, dom_eigenvec + return eigval, dom_eigvec def compute_dominant_eigen(hessian): """ Compute the eigenvalues (sorted by absolute value) - of the symmetrical Hessian matrix. + of symmetrical Hessian matrix, selecting the eigenvectors related to the dominant ones. Parameters ---------- @@ -56,21 +51,17 @@ def compute_dominant_eigen(hessian): Returns ------- - sorted_eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + srt_eigval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) Hessian eigenvalues sorted by absolute value (ascending order) - dom_eigenvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + dom_eigvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) Hessian eigenvectors related to the dominant (minimum) eigenvalue """ - # compute and sort the eigenvalues/eigenvectors - # of the image Hessian matrices eigenval, eigenvec = np.linalg.eigh(hessian) - sorted_eigenval, sorted_eigenvec = sort_eigen(eigenval, eigenvec) - - # select the eigenvectors related to dominant eigenvalues - dom_eigenvec = sorted_eigenvec[..., 0] + srt_eigval, srt_eigvec = sort_eigen(eigenval, eigenvec) + dom_eigvec = srt_eigvec[..., 0] - return sorted_eigenval, dom_eigenvec + return srt_eigval, dom_eigvec def compute_fractional_anisotropy(eigenval): @@ -244,7 +235,7 @@ def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma): return vesselness -def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, dark=False): +def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None): """ Compute fiber orientation vectors at the input spatial scale of interest @@ -265,104 +256,78 @@ def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, d gamma: float background score sensitivity - dark: bool - if True, enhance black 3D tubular structures - (i.e., negative contrast polarity) - Returns ------- frangi_img: numpy.ndarray (axis order=(Z,Y,X), dtype=float) Frangi's vesselness likelihood image - eigenvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + eigvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) 3D orientation map at the input spatial scale - eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + eigval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) Hessian eigenvalues sorted by absolute value (ascending order) """ - # Hessian matrix estimation and eigenvalue decomposition - eigenval, eigenvec = analyze_hessian_eigen(img, scale_px) + # compute local Hessian matrices and perform eigenvalue decomposition + eigval, eigvec = analyze_hessian_eigen(img, scale_px) - # compute Frangi's vesselness probability - eigen1, eigen2, eigen3 = eigenval + # compute Frangi's vesselness probability image + eigen1, eigen2, eigen3 = eigval vesselness = compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha=alpha, beta=beta, gamma=gamma) + frangi_img = reject_vesselness_background(vesselness, eigen2, eigen3) + eigval = np.stack(eigval, axis=-1) - # reject vesselness background - frangi_img = reject_vesselness_background(vesselness, eigen2, eigen3, dark) + return frangi_img, eigvec, eigval - # stack eigenvalues list - eigenval = np.stack(eigenval, axis=-1) - return frangi_img, eigenvec, eigenval - - -def compute_parall_scaled_orientation(scales_px, img, alpha=0.001, beta=1, gamma=None, dark=False): +def init_frangi_arrays(in_img, cfg, tmp_dir): """ - Compute fiber orientation vectors over the spatial scales of interest using concurrent workers. + Initialize the output datasets of the Frangi filter stage. Parameters ---------- - scales_px: numpy.ndarray (dtype=float) - Frangi filter scales [px] - - img: numpy.ndarray (axis order=(Z,Y,X)) - 3D microscopy image + in_img: dict + input image dictionary - alpha: float - plate-like score sensitivity + data: numpy.ndarray (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) + 3D microscopy image - beta: float - blob-like score sensitivity + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - gamma: float - background score sensitivity + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - dark: bool - if True, enhance black 3D tubular structures - (i.e., negative contrast polarity) + fb_ch: int + neuronal fibers channel - Returns - ------- - frangi_img: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - Frangi's vesselness likelihood image + bc_ch: int + brain cell soma channel - eigenvec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - 3D orientation map at the input spatial scale + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - Hessian eigenvalues sorted by absolute value (ascending order) - """ - n_scales = len(scales_px) - with Parallel(n_jobs=n_scales, prefer='threads', require='sharedmem') as parallel: - par_res = parallel( - delayed(compute_scaled_orientation)( - scales_px[i], img=img, - alpha=alpha, beta=beta, gamma=gamma, dark=dark) for i in range(n_scales)) + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - # unpack and stack results - enh_img_tpl, eigenvec_tpl, eigenval_tpl = zip(*par_res) - eigenval = np.stack(eigenval_tpl, axis=0) - eigenvec = np.stack(eigenvec_tpl, axis=0) - frangi = np.stack(enh_img_tpl, axis=0) + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - # get max scale-wise vesselness - best_idx = np.expand_dims(np.argmax(frangi, axis=0), axis=0) - frangi = np.take_along_axis(frangi, best_idx, axis=0).squeeze(axis=0) + name: str + name of the 3D microscopy image - # select fiber orientation vectors (and the associated eigenvalues) among different scales - best_idx = np.expand_dims(best_idx, axis=-1) - eigenval = np.take_along_axis(eigenval, best_idx, axis=0).squeeze(axis=0) - fbr_vec = np.take_along_axis(eigenvec, best_idx, axis=0).squeeze(axis=0) + is_vec: bool + vector field flag - return frangi, fbr_vec, eigenval + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape + shape_um: numpy.ndarray (shape=(3,), dtype=float) + total image shape [μm] -def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False): - """ - Initialize the output datasets of the Frangi filtering stage. + item_sz: int + image item size [B] - Parameters - ---------- cfg: dict Frangi filter configuration @@ -387,9 +352,6 @@ def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False): px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - z_rng: int - output z-range in [px] - bc_ch: int neuronal bodies channel @@ -406,102 +368,70 @@ def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False): exp_all: bool export all images - img_shp: numpy.ndarray (shape=(3,), dtype=int) - 3D image shape [px] - - slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [px] - - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio - tmp_dir: str - temporary file directory - - msk_bc: bool - if True, mask neuronal bodies - in the optionally provided image channel + path to temporary folder Returns ------- out_img: dict output image dictionary - vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - initialized fiber orientation 3D image + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) + initialized fiber orientation 3D image - clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) - initialized orientation colormap image + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + initialized orientation colormap image - fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized fractional anisotropy image + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized fractional anisotropy image - frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized Frangi-enhanced image + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized Frangi-enhanced image - iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized fiber image (isotropic resolution) + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized fiber image (isotropic resolution) - fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized fiber mask image + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized fiber mask image - bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized soma mask image - - z_sel: NumPy slice object - selected z-depth range + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized soma mask image """ - # shape copies - img_shp = img_shp.copy() - slc_shp = slc_shp.copy() - - # adapt output z-axis shape if required - z_min, z_max = cfg['z_rng'] - if z_min != 0 or z_max is not None: - if z_max is None: - z_max = slc_shp[0] - img_shp[0] = z_max - z_min - z_sel = slice(z_min, z_max, 1) - # output shape - img_dims = len(img_shp) + img_shp = in_img['shape'].copy() + img_shp[0] = cfg['z_out'].stop - cfg['z_out'].start + rsz_ratio = np.divide(in_img['px_sz'], cfg['px_sz']) tot_shp = tuple(np.ceil(rsz_ratio * img_shp).astype(int)) + vec_shp = tot_shp + (len(img_shp),) + cfg.update({'rsz': rsz_ratio}) # fiber channel arrays - iso_fbr_img = create_memory_map(tot_shp, dtype='uint8', name='iso_fiber', tmp_dir=tmp_dir) + iso_fbr = create_memory_map(tot_shp, dtype='uint8', name='iso', tmp=tmp_dir) if cfg['exp_all']: - frangi_img = create_memory_map(tot_shp, dtype='uint8', name='frangi', tmp_dir=tmp_dir) - fbr_msk = create_memory_map(tot_shp, dtype='uint8', name='fbr_msk', tmp_dir=tmp_dir) - fa_img = create_memory_map(tot_shp, dtype='float32', name='fa', tmp_dir=tmp_dir) + frangi = create_memory_map(tot_shp, dtype='uint8', name='frangi', tmp=tmp_dir) + fbr_msk = create_memory_map(tot_shp, dtype='uint8', name='fbr_msk', tmp=tmp_dir) + fa = create_memory_map(tot_shp, dtype='float32', name='fa', tmp=tmp_dir) # soma channel array - if msk_bc: - bc_msk = create_memory_map(tot_shp, dtype='uint8', name='bc_msk', tmp_dir=tmp_dir) + if in_img['msk_bc']: + bc_msk = create_memory_map(tot_shp, dtype='uint8', name='bc_msk', tmp=tmp_dir) else: bc_msk = None else: - frangi_img, fbr_msk, fa_img, bc_msk = (None, None, None, None) + frangi, fbr_msk, fa, bc_msk = 4 * (None,) # fiber orientation arrays - vec_shape = tot_shp + (img_dims,) - fbr_vec_img = create_memory_map(vec_shape, dtype='float32', name='fiber_vec', tmp_dir=tmp_dir) - fbr_vec_clr = create_memory_map(vec_shape, dtype='uint8', name='fiber_cmap', tmp_dir=tmp_dir) + fbr_vec = create_memory_map(vec_shp, dtype='float32', name='vec', tmp=tmp_dir) + fbr_clr = create_memory_map(vec_shp, dtype='uint8', name='clr', tmp=tmp_dir) # fill output image dictionary - out_img = dict() - out_img['vec'] = fbr_vec_img - out_img['clr'] = fbr_vec_clr - out_img['fa'] = fa_img - out_img['frangi'] = frangi_img - out_img['iso'] = iso_fbr_img - out_img['fbr_msk'] = fbr_msk - out_img['bc_msk'] = bc_msk + out_img = {'vec': fbr_vec, 'clr': fbr_clr, 'fa': fa, 'frangi': frangi, + 'iso': iso_fbr, 'fbr_msk': fbr_msk, 'bc_msk': bc_msk, 'px_sz': cfg['px_sz']} - return out_img, z_sel + return out_img -def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=False, hsv=False): +def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, hsv=False, _fa=False): """ Apply 3D Frangi filter to 3D microscopy image. @@ -522,74 +452,95 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=Fals gamma: float background score sensitivity (if None, gamma is automatically tailored) - dark: bool - if True, enhance black 3D tubular structures - (i.e., negative contrast polarity) - hsv: bool + generate an HSV colormap of 3D fiber orientations + + _fa: bool + compute fractional anisotropy Returns ------- - orient_slc: dict - slice orientation dictionary + out_slc: dict + slice output data - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - 3D fiber orientation field + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - Frangi's vesselness likelihood image + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + Frangi's vesselness likelihood image """ - # single-scale vesselness analysis - if len(scales_px) == 1: - frangi, fbr_vec, eigenval = \ - compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark) - - # parallel scaled vesselness analysis - else: - frangi, fbr_vec, eigenval = \ - compute_parall_scaled_orientation(scales_px, img, alpha=alpha, beta=beta, gamma=gamma, dark=dark) - - # generate fractional anisotropy image - fa = compute_fractional_anisotropy(eigenval) - - # generate fiber orientation color map + # single-scale or parallel multi-scale vesselness analysis + ns = len(scales_px) + frangi = np.zeros((ns,) + img.shape, dtype='float32') + eigvec = np.zeros((ns,) + img.shape + (3,), dtype='float32') + eigval = np.zeros((ns,) + img.shape + (3,), dtype='float32') + for s in range(ns): + frangi[s], eigvec[s], eigval[s] = \ + compute_scaled_orientation(scales_px[s], img, alpha=alpha, beta=beta, gamma=gamma) + + # get maximum response across the requested scales + max_idx = np.expand_dims(np.argmax(frangi, axis=0), axis=0) + frangi = np.take_along_axis(frangi, max_idx, axis=0).squeeze(axis=0) + + # select fiber orientation vectors (and the associated eigenvalues) among different scales + max_idx = np.expand_dims(max_idx, axis=-1) + eigval = np.take_along_axis(eigval, max_idx, axis=0).squeeze(axis=0) + fbr_vec = np.take_along_axis(eigvec, max_idx, axis=0).squeeze(axis=0) + + # compute fractional anisotropy image and fiber orientation color map + fa = compute_fractional_anisotropy(eigval) if _fa else None fbr_clr = hsv_orient_cmap(fbr_vec) if hsv else rgb_orient_cmap(fbr_vec) - # fill orientation dictionary - orient_slc = dict() - orient_slc['vec_slc'] = fbr_vec - orient_slc['clr_slc'] = fbr_clr - orient_slc['fa_slc'] = fa + # fill slice output dictionary + out_slc = {'vec': fbr_vec, 'clr': fbr_clr, 'fa': fa, 'frangi': frangi} - return orient_slc, frangi + return out_slc -def mask_background(img, orient, method='yen', invert=False, ts_msk=None): +def mask_background(out_slc, ref_img=None, method='yen', invert=False, ornt_keys=('vec', 'clr', 'fa')): """ - Mask fiber orientation arrays. + Mask fiber orientation data arrays. Parameters ---------- - img: numpy.ndarray (axis order=(Z,Y,X)) - fiber (or neuron) fluorescence 3D image + out_slc: dict + slice output dictionary + + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field + + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap + + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy + + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) + + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice - orient: dict - slice orientation dictionary + ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + tissue mask slice - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - fiber orientation vector slice + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask slice - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap slice + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + brain cell mask slice - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - fractional anisotropy slice + rng: NumPy slice object + output range + + ref_img: numpy.ndarray (axis order=(Z,Y,X) + reference image used for thresholding method: str thresholding method (refer to skimage.filters) @@ -597,48 +548,41 @@ def mask_background(img, orient, method='yen', invert=False, ts_msk=None): invert: bool mask inversion flag - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + ornt_keys: tuple + fiber orientation data keys Returns ------- - orient: dict - (masked) slice orientation dictionary - - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - fiber orientation vector slice (masked) - - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap slice (masked) - - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - fractional anisotropy slice (masked) + out_slc: dict bg: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) background mask """ # generate background mask - bg = create_background_mask(img, method=method) + if ref_img is None: + ref_img = out_slc['frangi'] + bg = create_background_mask(ref_img, method=method) # apply tissue reconstruction mask, when provided - if ts_msk is not None: - bg = np.logical_or(bg, np.logical_not(ts_msk)) + if out_slc['ts_msk'] is not None: + bg = np.logical_or(bg, np.logical_not(out_slc['ts_msk'])) # invert mask if invert: bg = np.logical_not(bg) # apply mask to input orientation data dictionary - for key in orient.keys(): - if orient[key].ndim == 3: - orient[key][bg] = 0 - else: - orient[key][bg, :] = 0 + for key in out_slc.keys(): + if key in ornt_keys and out_slc[key] is not None: + if out_slc[key].ndim == 3: + out_slc[key][bg] = 0 + else: + out_slc[key][bg, :] = 0 - return orient, bg + return out_slc, bg -def reject_vesselness_background(vesselness, eigen2, eigen3, dark): +def reject_vesselness_background(vesselness, eigen2, eigen3): """ Reject the fiber background, exploiting the sign of the "secondary" eigenvalues λ2 and λ3. @@ -654,33 +598,28 @@ def reject_vesselness_background(vesselness, eigen2, eigen3, dark): eigen3: numpy.ndarray (axis order=(Z,Y,X), dtype=float) highest Hessian eigenvalue - dark: bool - if True, enhance black 3D tubular structures - (i.e., negative contrast polarity) - Returns ------- vesselness: numpy.ndarray (axis order=(Z,Y,X), dtype=float) masked Frangi's vesselness likelihood image """ - bg_msk = np.logical_or(eigen2 < 0, eigen3 < 0) if dark else np.logical_or(eigen2 > 0, eigen3 > 0) - bg_msk = np.logical_or(bg_msk, np.isnan(vesselness)) + bg_msk = np.logical_or(np.logical_or(eigen2 > 0, eigen3 > 0), np.isnan(vesselness)) vesselness[bg_msk] = 0 return vesselness -def sort_eigen(eigenval, eigenvec, axis=-1): +def sort_eigen(eigval, eigvec, axis=-1): """ Sort eigenvalue and related eigenvector arrays by absolute value along the given axis. Parameters ---------- - eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + eigval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) original eigenvalue array - eigenvec: numpy.ndarray (axis order=(Z,Y,X,C,C), dtype=float) + eigvec: numpy.ndarray (axis order=(Z,Y,X,C,C), dtype=float) original eigenvector array axis: int @@ -688,102 +627,111 @@ def sort_eigen(eigenval, eigenvec, axis=-1): Returns ------- - sorted_eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + srt_eigval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) sorted eigenvalue array (ascending order) - sorted_eigenvec: numpy.ndarray (axis order=(Z,Y,X,C,C), dtype=float) + srt_eigvec: numpy.ndarray (axis order=(Z,Y,X,C,C), dtype=float) sorted eigenvector array """ # sort the eigenvalue array by absolute value (ascending order) - sorted_val_index = np.abs(eigenval).argsort(axis) - sorted_eigenval = np.take_along_axis(eigenval, sorted_val_index, axis) - sorted_eigenval = [np.squeeze(eigenval, axis=axis) - for eigenval in np.split(sorted_eigenval, sorted_eigenval.shape[axis], axis=axis)] + srt_val_idx = np.abs(eigval).argsort(axis) + srt_eigval = np.take_along_axis(eigval, srt_val_idx, axis) + srt_eigval = [np.squeeze(eigval, axis=axis) for eigval in np.split(srt_eigval, srt_eigval.shape[axis], axis=axis)] - # sort eigenvectors consistently - sorted_vec_index = sorted_val_index[:, :, :, np.newaxis, :] - sorted_eigenvec = np.take_along_axis(eigenvec, sorted_vec_index, axis) + # sort related eigenvectors consistently + srt_vec_idx = srt_val_idx[:, :, :, np.newaxis, :] + srt_eigvec = np.take_along_axis(eigvec, srt_vec_idx, axis) - return sorted_eigenval, sorted_eigenvec + return srt_eigval, srt_eigvec -def write_frangi_arrays(out_rng, z_sel, iso_slc, frangi_slc, fbr_msk_slc, bc_msk_slc, vec_slc, clr_slc, fa_slc, - vec, clr, fa, frangi, iso, fbr_msk, bc_msk): +def write_frangi_arrays(out_img, out_slc, rng, z_out=None): """ - Fill the memory-mapped output arrays of the Frangi filter stage. + Fill the output arrays of the Frangi filter stage. Parameters ---------- - out_rng: tuple - 3D slice output index range + out_img: dict + output image dictionary - z_sel: NumPy slice object - selected z-depth range + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - iso_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image slice + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - frangi_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - Frangi-enhanced image slice + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - fbr_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fiber mask image slice + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability image) - bc_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - soma mask image slice + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image - vec_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - fiber orientation vector image slice + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask image - clr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - orientation colormap image slice + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + soma mask image - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - fractional anisotropy image slice + px_sz: numpy.ndarray (shape=(3,), dtype=float) + output pixel size [μm] + + out_slc: dict + slice output dictionary + + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field + + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap + + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy - vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector field + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) - clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice - fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image + ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + tissue mask slice - frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - Frangi-enhanced image (fiber probability image) + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask slice - iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + brain cell mask slice - fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - fiber mask image + rng: NumPy slice object + output range - bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - soma mask image + z_out: NumPy slice object + output z-range Returns ------- None """ - # fill memory-mapped output arrays - vec_rng_out = tuple(np.append(out_rng, slice(0, 3, 1))) - vec[vec_rng_out] = vec_slc[z_sel, ...] - clr[vec_rng_out] = clr_slc[z_sel, ...] - iso[out_rng] = iso_slc[z_sel, ...].astype(np.uint8) - - # optional output images: Frangi filter response - if frangi is not None: - frangi[out_rng] = (255 * frangi_slc[z_sel, ...]).astype(np.uint8) + vec_rng = tuple(np.append(rng, slice(0, 3, 1))) + out_img['vec'][vec_rng] = out_slc['vec'][z_out, ...] + out_img['clr'][vec_rng] = out_slc['clr'][z_out, ...] + out_img['iso'][rng] = out_slc['iso'][z_out, ...].astype(np.uint8) # optional output images: fractional anisotropy - if fa is not None: - fa[out_rng] = fa_slc[z_sel, ...] + if out_img['fa'] is not None: + out_img['fa'][rng] = out_slc['fa'][z_out, ...] + + # optional output images: Frangi filter response + if out_img['frangi'] is not None: + out_img['frangi'][rng] = (255 * out_slc['frangi'][z_out, ...]).astype(np.uint8) # optional output images: fiber mask - if fbr_msk is not None: - fbr_msk[out_rng] = (255 * (1 - fbr_msk_slc[z_sel, ...])).astype(np.uint8) + if out_img['fbr_msk'] is not None: + out_img['fbr_msk'][rng] = (255 * (1 - out_slc['fbr_msk'][z_out, ...])).astype(np.uint8) # optional output images: neuronal soma mask - if bc_msk is not None: - bc_msk[out_rng] = (255 * bc_msk_slc[z_sel, ...]).astype(np.uint8) + if out_img['bc_msk'] is not None: + out_img['bc_msk'][rng] = (255 * out_slc['bc_msk'][z_out, ...]).astype(np.uint8) diff --git a/foa3d/input.py b/foa3d/input.py index a3f61e9..19f6d5d 100644 --- a/foa3d/input.py +++ b/foa3d/input.py @@ -54,7 +54,6 @@ def get_cli_parser(): '* supported formats:\n' ' - .tif .tiff (microscopy image or fiber orientation vectors)\n' ' - .yml (ZetaStitcher\'s stitch file of tiled microscopy reconstruction)\n' - ' - .npy (fiber orientation vectors)\n' '* image axes order:\n' ' - grayscale image: (Z, Y, X)\n' ' - RGB image: (Z, Y, X, C) or (Z, C, Y, X)\n' @@ -73,7 +72,7 @@ def get_cli_parser(): 'use one thread per logical core if None') cli_parser.add_argument('-r', '--ram', type=float, default=None, help='maximum RAM available to the Frangi filter stage [GB]: use all if None') - cli_parser.add_argument('--px-size-xy', type=float, default=0.878, help='lateral pixel size [μm]') + cli_parser.add_argument('--px-size-xy', type=float, default=1.0, help='lateral pixel size [μm]') cli_parser.add_argument('--px-size-z', type=float, default=1.0, help='longitudinal pixel size [μm]') cli_parser.add_argument('--psf-fwhm-x', type=float, default=1.0, help='PSF FWHM along horizontal x-axis [μm]') cli_parser.add_argument('--psf-fwhm-y', type=float, default=1.0, help='PSF FWHM along vertical y-axis [μm]') @@ -96,13 +95,10 @@ def get_cli_parser(): help='apply neuronal body mask (the optional channel of neuronal bodies must be available)') cli_parser.add_argument('-t', '--tissue-msk', action='store_true', default=False, help='apply tissue background mask') - cli_parser.add_argument('-v', '--vec', action='store_true', default=False, - help='fiber orientation vector image') cli_parser.add_argument('-e', '--exp-all', action='store_true', default=False, help='save the full range of images produced by the Frangi filter and ODF stages, ' 'e.g. for testing purposes (see documentation)') - # parse arguments cli_args = cli_parser.parse_args() return cli_args @@ -110,7 +106,7 @@ def get_cli_parser(): def get_image_size(in_img): """ - Get information on the input 3D microscopy image. + Get information on the size of the input 3D microscopy image. Parameters ---------- @@ -120,48 +116,12 @@ def get_image_size(in_img): data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) 3D microscopy image - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask - ch_ax: int RGB image channel axis (either 1, 3, or None for grayscale images) - fb_ch: int - neuronal fibers channel - - bc_ch: int - brain cell soma channel - - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel - - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - 3D FWHM of the PSF [μm] - - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] - - name: str - name of the 3D microscopy image - - is_vec: bool - vector field flag - - Returns - ------- - in_img: dict - input image dictionary (extended) - - data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) - 3D microscopy image - ts_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask - ch_ax: int - RGB image channel axis (either 1, 3, or None for grayscale images) - fb_ch: int neuronal fibers channel @@ -178,34 +138,35 @@ def get_image_size(in_img): px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] + path: str + path to the 3D microscopy image + name: str name of the 3D microscopy image - is_vec: bool - vector field flag + fmt: str + format of the 3D microscopy image - shape: numpy.ndarray (shape=(3,), dtype=int) - total image shape + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher - shape_um: numpy.ndarray (shape=(3,), dtype=float) - total image shape [μm] - - item_sz: int - image item size [B] + is_vec: bool + vector field flag + Returns + ------- + None """ # adapt channel axis - in_img['shape'] = np.asarray(in_img['data'].shape) + img_shp = np.asarray(in_img['data'].shape) if in_img['ch_ax'] is None: - in_img['fb_ch'] = None - in_img['msk_bc'] = False + in_img.update({'fb_ch': None, 'msk_bc': False}) else: - in_img['shape'] = np.delete(in_img['shape'], in_img['ch_ax']) - - in_img['shape_um'] = np.multiply(in_img['shape'], in_img['px_sz']) - in_img['item_sz'] = get_item_bytes(in_img['data']) + img_shp = np.delete(img_shp, in_img['ch_ax']) - return in_img + in_img.update({'shape': img_shp, + 'shape_um': np.multiply(img_shp, in_img['px_sz']), + 'item_sz': get_item_bytes(in_img['data'])}) def get_image_info(cli_args): @@ -219,40 +180,38 @@ def get_image_info(cli_args): Returns ------- - img_path: str - path to the 3D microscopy image - - img_name: str - name of the 3D microscopy image + in_img: dict + input image dictionary + fb_ch: int + neuronal fibers channel - img_fmt: str - format of the 3D microscopy image + bc_ch: int + brain cell soma channel - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - 3D FWHM of the PSF [μm] + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - is_tiled: bool - True for tiled reconstructions aligned using ZetaStitcher + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - is_vec: bool - True when pre-estimated fiber orientation vectors - are directly provided to the pipeline + path: str + path to the 3D microscopy image - mip_msk: bool - apply tissue reconstruction mask (binarized MIP) + name: str + name of the 3D microscopy image - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel + fmt: str + format of the 3D microscopy image - fb_ch: int - myelinated fibers channel + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher - bc_ch: int - brain cell soma channel + msk_mip: bool + apply tissue reconstruction mask (binarized MIP) """ # get microscopy image path and name img_path = cli_args.image_path @@ -262,31 +221,28 @@ def get_image_info(cli_args): # check image format if len(split_name) == 1: raise ValueError('Format must be specified for input volume images!') - else: - img_fmt = split_name[-1] - img_name = img_name.replace(f'.{img_fmt}', '') - is_tiled = True if img_fmt == 'yml' else False - is_vec = cli_args.vec - - # get image channels - fb_ch = cli_args.fb_ch - bc_ch = cli_args.bc_ch + img_fmt = split_name[-1] + img_name = img_name.replace(f'.{img_fmt}', '') + is_tiled = True if img_fmt == 'yml' else False - # apply tissue reconstruction mask (binarized MIP) + # apply tissue reconstruction mask (binarized MIP) and/or brain cell soma mask msk_mip = cli_args.tissue_msk - - # apply brain cell soma mask msk_bc = cli_args.cell_msk - # get Frangi filter configuration + # append Frangi filter configuration to input image name cfg_lbl = get_config_label(cli_args) img_name = f'{img_name}_{cfg_lbl}' - # get microscopy image resolution + # get microscopy image channels and resolution + fb_ch = cli_args.fb_ch + bc_ch = cli_args.bc_ch px_sz, psf_fwhm = get_resolution(cli_args) - return img_path, img_name, img_fmt, px_sz, psf_fwhm, \ - is_tiled, is_vec, msk_mip, msk_bc, fb_ch, bc_ch + # populate input image dictionary + in_img = {'fb_ch': fb_ch, 'bc_ch': bc_ch, 'msk_bc': msk_bc, 'psf_fwhm': psf_fwhm, 'px_sz': px_sz, + 'path': img_path, 'name': img_name, 'fmt': img_fmt, 'is_tiled': is_tiled} + + return in_img, msk_mip def get_frangi_config(cli_args, in_img): @@ -367,8 +323,8 @@ def get_frangi_config(cli_args, in_img): px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - z_rng: int - output z-range in [px] + fb_thr: float or skimage.filters thresholding method + Frangi filter probability response threshold bc_ch: int neuronal bodies channel @@ -386,34 +342,30 @@ def get_frangi_config(cli_args, in_img): exp_all: bool export all images - + z_out: NumPy slice object + output z-range """ - # initialize Frangi filter configuration dictionary - frangi_cfg = dict() - - # preprocessing configuration (in-plane smoothing) - frangi_cfg['smooth_sd'], px_sz_iso = config_anisotropy_correction(in_img['px_sz'], in_img['psf_fwhm']) - - # Frangi filter parameters - frangi_cfg['alpha'], frangi_cfg['beta'], frangi_cfg['gamma'] = (cli_args.alpha, cli_args.beta, cli_args.gamma) - frangi_cfg['scales_um'] = np.array(cli_args.scales) - frangi_cfg['scales_px'] = frangi_cfg['scales_um'] / px_sz_iso[0] + # preprocessing configuration (adaptive smoothing) + smooth_sd, out_px_sz = config_anisotropy_correction(in_img['px_sz'], in_img['psf_fwhm']) + scales_px = np.array(cli_args.scales) / np.max(out_px_sz) - # Frangi filter response threshold - frangi_cfg['fb_thr'] = cli_args.fb_thr + # adapted output z-axis range when required + z_min = max(0, int(np.floor(cli_args.z_min / np.max(in_img['px_sz'])))) + z_max = min(int(np.ceil(cli_args.z_max / np.max(in_img['px_sz']))), in_img['shape'][0]) \ + if cli_args.z_max is not None else in_img['shape'][0] - # fiber orientation colormap - frangi_cfg['hsv_cmap'] = cli_args.hsv + # get threshold applied to the Frangi filter response + fb_thr = cli_args.fb_thr + if fb_thr.replace('.', '', 1).isdigit(): + fb_thr = float(fb_thr) - # forced output z-range - z_min = int(np.floor(cli_args.z_min / np.max(in_img['px_sz']))) - z_max = int(np.ceil(cli_args.z_max / np.max(in_img['px_sz']))) if cli_args.z_max is not None else cli_args.z_max - frangi_cfg['z_rng'] = (z_min, z_max) + # compile Frangi filter configuration dictionary + frangi_cfg = {'alpha': cli_args.alpha, 'beta': cli_args.beta, 'gamma': cli_args.gamma, + 'scales_um': np.array(cli_args.scales), 'scales_px': scales_px, 'smooth_sd': smooth_sd, + 'px_sz': out_px_sz, 'fb_thr': fb_thr, 'hsv_cmap': cli_args.hsv, + 'exp_all': cli_args.exp_all, 'z_out': slice(z_min, z_max, 1)} - # export all flag - frangi_cfg['exp_all'] = cli_args.exp_all - - return frangi_cfg, px_sz_iso + return frangi_cfg def get_resolution(cli_args): @@ -427,20 +379,19 @@ def get_resolution(cli_args): Returns ------- - px_sz: numpy.ndarray (shape=(3,), dtype=float) + px_sz: tuple (shape=(3,), dtype=float) pixel size [μm] - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + psf_fwhm: tuple (shape=(3,), dtype=float) 3D PSF FWHM [μm] """ - # create pixel and psf size arrays - px_sz = np.array([cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy]) - psf_fwhm = np.array([cli_args.psf_fwhm_z, cli_args.psf_fwhm_y, cli_args.psf_fwhm_x]) + px_sz = (cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy) + psf_fwhm = (cli_args.psf_fwhm_z, cli_args.psf_fwhm_y, cli_args.psf_fwhm_x) return px_sz, psf_fwhm -def get_resource_config(cli_args): +def get_resource_config(cli_args, frangi_cfg): """ Retrieve resource usage configuration of the Foa3D tool. @@ -449,22 +400,59 @@ def get_resource_config(cli_args): cli_args: see ArgumentParser.parse_args populated namespace of command line arguments + frangi_cfg: dict + Frangi filter configuration + + alpha: float + plate-like score sensitivity + + beta: float + blob-like score sensitivity + + gamma: float + background score sensitivity + + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] + + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] + + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + bc_ch: int + neuronal bodies channel + + fb_ch: int + myelinated fibers channel + + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel + + hsv_cmap: bool + generate HSV colormap of 3D fiber orientations + + exp_all: bool + export all images + + rsz: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio + Returns ------- - max_ram: float - maximum RAM available to the Frangi filter stage [B] - - jobs: int - number of parallel jobs (threads) - used by the Frangi filter stage + None """ - # resource parameters (convert maximum RAM to bytes) jobs = cli_args.jobs - max_ram = cli_args.ram - if max_ram is not None: - max_ram *= 1024**3 + ram = cli_args.ram + if ram is not None: + ram *= 1024**3 - return max_ram, jobs + frangi_cfg.update({'jobs': jobs, 'ram': ram}) def load_microscopy_image(cli_args): @@ -484,60 +472,71 @@ def load_microscopy_image(cli_args): in_img: dict input image dictionary - img_data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) - 3D microscopy image + data: numpy.ndarray (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) + 3D microscopy image - ch_ax: int - RGB image channel axis (either 1, 3, or None for grayscale images) + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - img_name: str - name of the 3D microscopy image + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + + fb_ch: int + neuronal fibers channel - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + bc_ch: int + brain cell soma channel - is_vec: bool - vector field flag + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - save_dirs: dict - saving directories - ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) - """ - # retrieve input file information - img_path, img_name, img_fmt, px_sz, psf_fwhm, is_tiled, is_vec, \ - msk_mip, msk_bc, fb_ch, bc_ch = get_image_info(cli_args) + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - # create saving directory - save_dirs = create_save_dirs(img_path, img_name, cli_args, is_vec=is_vec) + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - # import fiber orientation vector data - tic = perf_counter() - if is_vec: - in_img = load_orient(img_path, img_name, img_fmt, tmp_dir=save_dirs['tmp']) + path: str + path to the 3D microscopy image - # import raw 3D microscopy image - else: - in_img = load_raw(img_path, img_name, img_fmt, psf_fwhm, - is_tiled=is_tiled, tmp_dir=save_dirs['tmp'], - msk_mip=msk_mip, msk_bc=msk_bc, fb_ch=fb_ch, bc_ch=bc_ch) + name: str + name of the 3D microscopy image - # add image pixel size to input data dictionary - in_img['px_sz'] = px_sz + fmt: str + format of the 3D microscopy image - # add image name to input data dictionary - in_img['name'] = img_name + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher - # add vector field flag - in_img['is_vec'] = is_vec + is_vec: bool + vector field flag - # get microscopy image size information - in_img = get_image_size(in_img) + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape - # print import time - print_import_time(tic) + shape_um: numpy.ndarray (shape=(3,), dtype=float) + total image shape [μm] + + item_sz: int + image item size [B] + + save_dirs: dict + saving directories + ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) + """ + # get input information and create saving directories + in_img, msk_mip = get_image_info(cli_args) + save_dirs = create_save_dirs(cli_args, in_img) + + # import fiber orientation vector data or raw 3D microscopy image + tic = perf_counter() + load_data(in_img, save_dirs['tmp'], msk_mip=msk_mip) - # print volume image shape and resolution - if not is_vec: + # print input data information + get_image_size(in_img) + print_import_time(tic) + if not in_img['is_vec']: print_image_info(in_img) else: print_flsh() @@ -545,167 +544,88 @@ def load_microscopy_image(cli_args): return in_img, save_dirs -def load_orient(img_path, img_name, img_fmt, tmp_dir=None): +def load_data(in_img, tmp_dir, msk_mip=False): """ - Load array of 3D fiber orientations. + Load 3D microscopy data. Parameters ---------- - img_path: str - path to the 3D microscopy image - - img_name: str - name of the 3D microscopy image - - img_fmt: str - format of the 3D microscopy image - - tmp_dir: str - temporary file directory - - Returns - ------- in_img: dict input image dictionary + fb_ch: int + neuronal fibers channel - data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) - 3D microscopy image - - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask - - ch_ax: int - RGB image channel axis (either 1, 3, or None for grayscale images) - """ - # print heading - print_flsh(color_text(0, 191, 255, "\nFiber Orientation Data Import\n")) - - # load fiber orientations - if img_fmt == 'npy': - vec_img = np.load(img_path, mmap_mode='r') - elif img_fmt == 'tif' or img_fmt == 'tiff': - vec_img = tiff.imread(img_path) - ch_ax = detect_ch_axis(vec_img) - if ch_ax != 3: - vec_img = np.moveaxis(vec_img, ch_ax, -1) - - # memory-map the input TIFF image - vec_img = create_memory_map(vec_img.shape, dtype=vec_img.dtype, name=img_name, - tmp_dir=tmp_dir, arr=vec_img, mmap_mode='r') - - # check array shape - if vec_img.ndim != 4: - raise ValueError('Invalid 3D fiber orientation dataset (ndim != 4)!') - else: - print_flsh(f"Loading {img_path} orientation vector field...\n") - - # populate input image dictionary - in_img = dict() - in_img['data'] = vec_img - in_img['ts_msk'] = None - in_img['ch_ax'] = None - - return in_img + bc_ch: int + brain cell soma channel + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel -def load_raw(img_path, img_name, img_fmt, psf_fwhm, - is_tiled=False, tmp_dir=None, msk_mip=False, msk_bc=False, fb_ch=1, bc_ch=0): - """ - Load 3D microscopy image. + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - Parameters - ---------- - img_path: str - path to the 3D microscopy image + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - img_name: str - name of the 3D microscopy image + path: str + path to the 3D microscopy image - img_fmt: str - format of the 3D microscopy image + name: str + name of the 3D microscopy image - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - 3D FWHM of the PSF [μm] + fmt: str + format of the 3D microscopy image - is_tiled: bool - True for tiled reconstructions aligned using ZetaStitcher + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher tmp_dir: str - temporary file directory + path to temporary folder msk_mip: bool apply tissue reconstruction mask (binarized MIP) - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel - - fb_ch: int - myelinated fibers channel - - bc_ch: int - brain cell soma channel - Returns ------- - in_img: dict - input image dictionary - - data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) - 3D microscopy image - - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask - - ch_ax: int - RGB image channel axis (either 1, 3, or None for grayscale images) - - fb_ch: int - neuronal fibers channel - - bc_ch: int - brain cell soma channel - - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel - - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - 3D FWHM of the PSF [μm] + None """ - # print heading print_flsh(color_text(0, 191, 255, "\nMicroscopy Image Import\n")) - # load microscopy tiled reconstruction (aligned using ZetaStitcher) - if is_tiled: - print_flsh(f"Loading {img_path} tiled reconstruction...\n") - img = VirtualFusedVolume(img_path) + # load tiled reconstruction (aligned using ZetaStitcher) + if in_img['is_tiled']: + print_flsh(f"Loading {in_img['path']} tiled reconstruction...\n") + img = VirtualFusedVolume(in_img['path']) + ch_ax = detect_ch_axis(img) + is_vec = False - # load microscopy z-stack + # load z-stack else: - print_flsh(f"Loading {img_path} z-stack...\n") - img_fmt = img_fmt.lower() - if img_fmt == 'tif' or img_fmt == 'tiff': - img = tiff.imread(img_path) - else: - raise ValueError('Unsupported image format!') + print_flsh(f"Loading {in_img['path']} z-stack...\n") - # create image memory map - img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r') + img_fmt = in_img['fmt'].lower() + if img_fmt in ('tif', 'tiff'): + img = tiff.imread(in_img['path']) + ch_ax = detect_ch_axis(img) - # detect channel axis (RGB images) - ch_ax = detect_ch_axis(img) + # detect vector field input + is_vec = img.ndim == 4 and img.dtype in (np.float32, float, 'float32') + if is_vec: + if ch_ax != 3: + img = np.moveaxis(img, ch_ax, -1) - # generate tissue reconstruction mask - if msk_mip: - dims = len(img.shape) + img = create_memory_map(img.shape, dtype=img.dtype, name=in_img['name'], + tmp=tmp_dir, arr=img, mmap_mode='r') + else: + raise ValueError('Unsupported image format!') - # grayscale image + # generate tissue background mask + if not is_vec and msk_mip: + dims = len(img.shape) if dims == 3: img_fbr = img - # RGB image elif dims == 4: - img_fbr = img[:, fb_ch, :, :] if ch_ax == 1 else img[..., fb_ch] + img_fbr = img[:, in_img['fb_ch'], :, :] if ch_ax == 1 else img[..., in_img['fb_ch']] else: raise ValueError('Invalid image (ndim != 3 and ndim != 4)!') @@ -719,14 +639,5 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm, else: ts_msk = None - # populate input image dictionary - in_img = dict() - in_img['data'] = img - in_img['ts_msk'] = ts_msk - in_img['ch_ax'] = ch_ax - in_img['fb_ch'] = fb_ch - in_img['bc_ch'] = bc_ch - in_img['msk_bc'] = msk_bc - in_img['psf_fwhm'] = psf_fwhm - - return in_img + # update input image dictionary + in_img.update({'data': img, 'ts_msk': ts_msk, 'ch_ax': ch_ax, 'is_vec': is_vec}) diff --git a/foa3d/odf.py b/foa3d/odf.py index 89a42fa..ca3bd9c 100644 --- a/foa3d/odf.py +++ b/foa3d/odf.py @@ -7,8 +7,7 @@ from foa3d.utils import create_memory_map, normalize_image, transform_axes -def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tensor_eigen, vxl_side, odf_norm, - odf_deg=6, vxl_thr=0.5, vec_thr=0.000001): +def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tnsr_eig, scale, norm, deg=6, vx_thr=0.5, vec_thr=1e-6): """ Compute the spherical harmonics coefficients iterating over super-voxels of fiber orientation vectors. @@ -22,7 +21,7 @@ def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tensor_eigen, vxl_si adjusted isotropic pixel size [μm] odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) - initialized array of ODF spherical harmonics coefficients + initialized array of ODF spherical harmonics coefficients odi: dict orientation dispersion dictionary @@ -39,21 +38,22 @@ def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tensor_eigen, vxl_si odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) orientation dispersion index anisotropy - fbr_dnst + fbr_dnst: NumPy memory-map object (axis order=(Z,Y,X), dtype=float) + initialized fiber density image - vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + vec_tnsr_eig: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) initialized array of orientation tensor eigenvalues - vxl_side: int + scale: int side of the ODF super-voxel [px] - odf_norm: numpy.ndarray (dtype: float) + norm: numpy.ndarray (dtype: float) 2D array of spherical harmonics normalization factors - odf_deg: int + deg: int degrees of the spherical harmonics series expansion - vxl_thr: float + vx_thr: float minimum relative threshold on the sliced voxel volume vec_thr: float @@ -63,57 +63,56 @@ def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tensor_eigen, vxl_si ------- odf: numpy.ndarray (axis order=(X,Y,Z,C), dtype=float32) volumetric map of real-valued spherical harmonics coefficients - """ - # iterate over ODF super-voxels - ref_vxl_size = min(vxl_side, fbr_vec.shape[0]) * vxl_side**2 - for z in range(0, fbr_vec.shape[0], vxl_side): - zmax = z + vxl_side - - for y in range(0, fbr_vec.shape[1], vxl_side): - ymax = y + vxl_side - - for x in range(0, fbr_vec.shape[2], vxl_side): - xmax = x + vxl_side - - # slice orientation voxel - vec_vxl = fbr_vec[z:zmax, y:ymax, x:xmax, :] - zerovec = np.count_nonzero(np.all(vec_vxl == 0, axis=-1)) - sli_vxl_size = np.prod(vec_vxl.shape[:-1]) - - # local fiber density estimate - zv, yv, xv = z // vxl_side, y // vxl_side, x // vxl_side - fbr_dnst[zv, yv, xv] = (sli_vxl_size - zerovec) / (sli_vxl_size * np.prod(px_sz)) - # skip boundary voxels and voxels without enough non-zero orientation vectors - if sli_vxl_size / ref_vxl_size > vxl_thr and 1 - zerovec / sli_vxl_size > vec_thr: - - # flatten orientation supervoxel - vec_arr = vec_vxl.ravel() - - # compute ODF and orientation tensor eigenvalues - odf[zv, yv, xv, :] = fiber_vectors_to_sph_harm(vec_arr, odf_deg, odf_norm) - vec_tensor_eigen[zv, yv, xv, :] = compute_vec_tensor_eigen(vec_arr) - - # compute slice-wise dispersion and anisotropy parameters - compute_orientation_dispersion(vec_tensor_eigen, **odi) + fbr_dnst: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + fiber density image + """ + # iterate over super-voxels + px_vol = np.prod(px_sz) + ref_vx_vol = min(scale, fbr_vec.shape[0]) * scale**2 + for z in range(0, fbr_vec.shape[0], scale): + z_max = z + scale + for y in range(0, fbr_vec.shape[1], scale): + y_max = y + scale + for x in range(0, fbr_vec.shape[2], scale): + x_max = x + scale + + # select super-voxel of fiber orientation data + vec_vx = fbr_vec[z:z_max, y:y_max, x:x_max, :] + zerovec = np.count_nonzero(np.all(vec_vx == 0, axis=-1)) + vx_vol = np.prod(vec_vx.shape[:-1]) + + # compute local fiber density + zv, yv, xv = z // scale, y // scale, x // scale + fbr_dnst[zv, yv, xv] = (vx_vol - zerovec) / (vx_vol * px_vol) + + # compute ODF and orientation tensor eigenvalues + # (skipping boundary voxels and voxels without enough data) + if vx_vol / ref_vx_vol > vx_thr and 1 - zerovec / vx_vol > vec_thr: + vec_arr = vec_vx.ravel() + odf[zv, yv, xv, :] = fiber_vectors_to_sph_harm(vec_arr, deg, norm) + vec_tnsr_eig[zv, yv, xv, :] = compute_vec_tensor_eigen(vec_arr) + + # compute dispersion and anisotropy parameters + compute_orientation_dispersion(vec_tnsr_eig, **odi) # set dispersion background to 0 - mask_orientation_dispersion(vec_tensor_eigen, odi) + mask_orientation_dispersion(vec_tnsr_eig, odi) - # transform axes + # manipulate ODF axes (for visualization in MRtrix3) odf = transform_axes(odf, swapped=(0, 2), flipped=(1, 2)) return odf, fbr_dnst @njit(cache=True) -def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis): +def compute_orientation_dispersion(vec_tnsr_eig, odi_pri, odi_sec, odi_tot, odi_anis): """ Compute orientation dispersion parameters. Parameters ---------- - vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + vec_tnsr_eig: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) orientation tensor eigenvalues computed from an ODF super-voxel @@ -133,13 +132,8 @@ def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, ------- None """ - # get difference between secondary tensor eigenvalues - diff = np.abs(vec_tensor_eigen[..., 1] - vec_tensor_eigen[..., 0]) - - # get absolute tensor eigenvalues - avte = np.abs(vec_tensor_eigen) - # primary dispersion (0.3183098861837907 = 1/π) + avte = np.abs(vec_tnsr_eig) if odi_pri is not None: odi_pri[:] = (1 - 0.3183098861837907 * np.arctan2(avte[..., 2], avte[..., 1])).astype(np.float32) @@ -148,6 +142,7 @@ def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_sec[:] = (1 - 0.3183098861837907 * np.arctan2(avte[..., 2], avte[..., 0])).astype(np.float32) # dispersion anisotropy + diff = np.abs(vec_tnsr_eig[..., 1] - vec_tnsr_eig[..., 0]) if odi_anis is not None: odi_anis[:] = (1 - 0.3183098861837907 * np.arctan2(avte[..., 2], diff)).astype(np.float32) @@ -183,42 +178,37 @@ def compute_vec_tensor_eigen(fbr_vec): return vec_tensor_eigen -def generate_odf_background(bg_mrtrix, fbr_vec, vxl_side, iso_fbr=None): +def generate_odf_background(bg_mrtrix, fbr_vec, scale, iso_fbr=None): """ - Generate the downsampled background image required + Generate the down-sampled background image required to visualize the 3D ODF map in MRtrix3. Parameters ---------- - bg_mrtrix: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8) + bg_mrtrix: numpy.ndarray (axis order=(X,Y,Z), dtype=uint8) initialized background array for ODF visualization in MRtrix3 - fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vectors - vxl_side: int + scale: int side of the ODF super-voxel [px] - iso_fbr: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + iso_fbr: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) isotropic fiber image Returns ------- None """ - # select background data - bg_img = fbr_vec if iso_fbr is None else iso_fbr - - # get shape of new downsampled array - new_shape = bg_mrtrix.shape[:-1] - # image normalization: get global minimum and maximum values + bg_img = fbr_vec if iso_fbr is None else iso_fbr if bg_img.ndim == 3: min_glob = np.min(bg_img) max_glob = np.max(bg_img) - # loop over z-slices, and resize them - for z in range(vxl_side // 2, bg_img.shape[0], vxl_side): + # loop over z-slices, rescale and resize them + for z in range(scale // 2, bg_img.shape[0], scale): if bg_img.ndim == 3: tmp_slice = normalize_image(bg_img[z], min_val=min_glob, max_val=max_glob) elif bg_img.ndim == 4: @@ -226,11 +216,11 @@ def generate_odf_background(bg_mrtrix, fbr_vec, vxl_side, iso_fbr=None): tmp_slice = np.where(tmp_slice <= 255.0, tmp_slice, 255.0) tmp_slice = transform_axes(tmp_slice, swapped=(0, 1), flipped=(0, 1)) - bg_mrtrix[..., z // vxl_side] = \ - resize(tmp_slice, output_shape=new_shape, anti_aliasing=True, preserve_range=True) + bg_mrtrix[..., z // scale] = \ + resize(tmp_slice, output_shape=bg_mrtrix.shape[:-1], anti_aliasing=True, preserve_range=True) -def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False): +def init_odf_arrays(vec_img_shp, tmp_dir, scale, deg=6, exp_all=False): """ Initialize the output datasets of the ODF analysis stage. @@ -240,12 +230,12 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False): vector volume shape [px] tmp_dir: str - temporary file directory + path to temporary folder - odf_scale: int + scale: int fiber ODF resolution (super-voxel side [px]) - odf_deg: int + deg: int degrees of the spherical harmonics series expansion exp_all: bool @@ -271,8 +261,8 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False): odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) initialized array of orientation dispersion anisotropy parameters - fbr_dnst: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - initialized fiber density map + fbr_dnst: NumPy memory-map object (axis order=(Z,Y,X), dtype=float) + initialized fiber density image bg_mrtrix: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8) initialized background for ODF visualization in MRtrix3 @@ -280,52 +270,34 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False): vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) initialized fiber orientation tensor eigenvalues """ - # ODI maps shape - odi_shp = tuple(np.ceil(np.divide(vec_img_shp, odf_scale)).astype(int)) - - # create downsampled background memory map - bg_shp = tuple(np.flip(odi_shp)) - bg_mrtrix = create_memory_map(bg_shp, dtype='uint8', name=f'bg_tmp{odf_scale}', tmp_dir=tmp_dir) - - # create ODF memory map - nc = get_sph_harm_ncoeff(odf_deg) - odf_shp = odi_shp + (nc,) - odf = create_memory_map(odf_shp, dtype='float32', name=f'odf_tmp{odf_scale}', tmp_dir=tmp_dir) - - # create orientation tensor memory map - vec_tensor_shape = odi_shp + (3,) - vec_tensor_eigen = create_memory_map(vec_tensor_shape, dtype='float32', - name=f'tensor_tmp{odf_scale}', tmp_dir=tmp_dir) - - # fiber density memory map - fbr_dnst = create_memory_map(odi_shp, dtype='float32', name=f'fbr_dnst_tmp{odf_scale}', tmp_dir=tmp_dir) - - # create ODI memory maps - odi_tot = create_memory_map(odi_shp, dtype='float32', name=f'odi_tot_tmp{odf_scale}', tmp_dir=tmp_dir) + # initialize ODF and orientation dispersion maps + nc = get_sph_harm_ncoeff(deg) + odi_shp = tuple(np.ceil(np.divide(vec_img_shp, scale)).astype(int)) + odf = create_memory_map(odi_shp + (nc,), dtype='float32', name=f'odf{scale}', tmp=tmp_dir) + odi_tot = create_memory_map(odi_shp, dtype='float32', name=f'odi_tot{scale}', tmp=tmp_dir) + bg_mrtrix = create_memory_map(tuple(np.flip(odi_shp)), dtype='uint8', name=f'bg{scale}', tmp=tmp_dir) if exp_all: - odi_pri = create_memory_map(odi_shp, dtype='float32', name=f'odi_pri_tmp{odf_scale}', tmp_dir=tmp_dir) - odi_sec = create_memory_map(odi_shp, dtype='float32', name=f'odi_sec_tmp{odf_scale}', tmp_dir=tmp_dir) - odi_anis = create_memory_map(odi_shp, dtype='float32', name=f'odi_anis_tmp{odf_scale}', tmp_dir=tmp_dir) + odi_pri = create_memory_map(odi_shp, dtype='float32', name=f'odi_pri{scale}', tmp=tmp_dir) + odi_sec = create_memory_map(odi_shp, dtype='float32', name=f'odi_sec{scale}', tmp=tmp_dir) + odi_anis = create_memory_map(odi_shp, dtype='float32', name=f'odi_anis{scale}', tmp=tmp_dir) else: odi_pri, odi_sec, odi_anis = (None, None, None) - # fill output image dictionary - odi = dict() - odi['odi_pri'] = odi_pri - odi['odi_sec'] = odi_sec - odi['odi_tot'] = odi_tot - odi['odi_anis'] = odi_anis + # initialize fiber orientation tensor, fiber density and orientation dispersion dictionary + vec_tnsr_eig = create_memory_map(odi_shp + (3,), dtype='float32', name=f'tensor{scale}', tmp=tmp_dir) + fbr_dnst = create_memory_map(odi_shp, dtype='float32', name=f'fbr_dnst{scale}', tmp=tmp_dir) + odi = {'odi_pri': odi_pri, 'odi_sec': odi_sec, 'odi_tot': odi_tot, 'odi_anis': odi_anis} - return odf, odi, fbr_dnst, bg_mrtrix, vec_tensor_eigen + return odf, odi, fbr_dnst, bg_mrtrix, vec_tnsr_eig -def mask_orientation_dispersion(vec_tensor_eigen, odi): +def mask_orientation_dispersion(vec_tnsr_eig, odi): """ Suppress orientation dispersion background. Parameters ---------- - vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + vec_tnsr_eig: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) orientation tensor eigenvalues computed from an ODF super-voxel @@ -348,8 +320,7 @@ def mask_orientation_dispersion(vec_tensor_eigen, odi): ------- None """ - # mask background - bg_msk = np.all(vec_tensor_eigen == 0, axis=-1) + bg_msk = np.all(vec_tnsr_eig == 0, axis=-1) for key in odi.keys(): if odi[key] is not None: odi[key][bg_msk] = 0 diff --git a/foa3d/output.py b/foa3d/output.py index c58467b..1387b2b 100644 --- a/foa3d/output.py +++ b/foa3d/output.py @@ -12,54 +12,75 @@ from foa3d.utils import get_item_size -def create_save_dirs(img_path, img_name, cli_args, is_vec=False): +def create_save_dirs(cli_args, in_img): """ Create saving directory. Parameters ---------- - img_path: str - path to input microscopy volume image - - img_name: str - name of the input volume image - cli_args: see ArgumentParser.parse_args updated namespace of command line arguments - is_vec: bool - True when fiber orientation vectors are provided as input - to the pipeline + in_img: dict + input image dictionary + fb_ch: int + neuronal fibers channel + + bc_ch: int + brain cell soma channel + + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel + + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + path: str + path to the 3D microscopy image + + name: str + name of the 3D microscopy image + + fmt: str + format of the 3D microscopy image + + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher + + is_vec: bool + vector field flag Returns ------- save_dirs: dict saving directories - ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) - """ - # get current time - time_stamp = datetime.now().strftime("%Y%m%d-%H%M%S") + frangi: Frangi filter + + odf: ODF analysis + + tmp: temporary data + """ # get output path out_path = cli_args.out if out_path is None: - out_path = path.dirname(img_path) + out_path = path.dirname(in_img['path']) # create saving directory - base_out_dir = path.join(out_path, f"Foa3D_{time_stamp}_{img_name}") + time_stamp = datetime.now().strftime("%Y%m%d-%H%M%S") + base_out_dir = path.join(out_path, f"Foa3D_{time_stamp}_{in_img['name']}") if not path.isdir(base_out_dir): makedirs(base_out_dir) - # initialize empty dictionary - save_dirs = dict() - # create Frangi filter output subdirectory - if not is_vec: - frangi_dir = path.join(base_out_dir, 'frangi') - mkdir(frangi_dir) - save_dirs['frangi'] = frangi_dir - else: - save_dirs['frangi'] = None + save_dirs = {} + frangi_dir = path.join(base_out_dir, 'frangi') + mkdir(frangi_dir) + save_dirs['frangi'] = frangi_dir # create ODF analysis output subdirectory if cli_args.odf_res is not None: @@ -70,8 +91,7 @@ def create_save_dirs(img_path, img_name, cli_args, is_vec=False): save_dirs['odf'] = None # create temporary directory - tmp_dir = tempfile.mkdtemp(dir=base_out_dir) - save_dirs['tmp'] = tmp_dir + save_dirs['tmp'] = tempfile.mkdtemp(dir=base_out_dir) return save_dirs @@ -113,7 +133,7 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): # check output format fmt = fmt.lower() - if fmt == 'tif' or fmt == 'tiff': + if fmt in ('tif', 'tiff'): # retrieve image pixel size px_sz_z, px_sz_y, px_sz_x = px_sz @@ -123,7 +143,7 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): nd_array = np.expand_dims(nd_array, 1) # adjust bigtiff optional argument - bigtiff = True if nd_array.itemsize * nd_array.size >= 4294967296 else False + bigtiff = nd_array.itemsize * nd_array.size >= 4294967296 metadata = {'axes': 'ZCYX', 'spacing': px_sz_z, 'unit': 'um'} out_name = f'{fname}.{fmt}' with TiffWriter(path.join(save_dir, out_name), bigtiff=bigtiff, append=True) as tif: @@ -144,7 +164,7 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): raise ValueError("Unsupported data format!!!") -def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None): +def save_frangi_arrays(save_dir, img_name, out_img, ram=None): """ Save the output arrays of the Frangi filter stage to TIF files. @@ -178,8 +198,8 @@ def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None): bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) neuron mask image - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size (Z,Y,X) [μm] + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size (Z,Y,X) [μm] ram: float maximum RAM available @@ -188,12 +208,11 @@ def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None): ------- None """ - # loop over output image dictionary fields and save to TIFF + # loop over output image dictionary fields and save to TIFF files for img_key in out_img.keys(): - if out_img[img_key] is not None: - save_array(f'{img_key}_{img_name}', save_dir, out_img[img_key], px_sz, ram=ram) + if isinstance(out_img[img_key], np.ndarray) and img_key not in (None, 'iso'): + save_array(f'{img_key}_{img_name}', save_dir, out_img[img_key], out_img['px_sz'], ram=ram) - # print output directory print_flsh(f"\nFrangi filter arrays saved to: {save_dir}\n") diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py index f79fed9..fe09bda 100644 --- a/foa3d/pipeline.py +++ b/foa3d/pipeline.py @@ -1,5 +1,4 @@ from time import perf_counter - import numpy as np from joblib import Parallel, delayed @@ -10,10 +9,11 @@ init_odf_arrays) from foa3d.output import save_frangi_arrays, save_odf_arrays from foa3d.preprocessing import correct_anisotropy -from foa3d.printing import (print_flsh, print_frangi_info, print_odf_info, - print_frangi_progress) -from foa3d.slicing import (check_background, config_frangi_batch, - generate_slice_lists, crop, crop_lst, slice_image) +from foa3d.printing import (print_flsh, print_frangi_config, + print_frangi_progress, print_odf_info) +from foa3d.slicing import (check_background, get_slicing_config, + generate_slice_ranges, crop, + crop_img_dict, slice_image) from foa3d.spharm import get_sph_harm_norm_factors from foa3d.utils import get_available_cores @@ -27,21 +27,27 @@ def parallel_frangi_over_slices(cli_args, save_dirs, in_img): cli_args: see ArgumentParser.parse_args populated namespace of command line arguments - save_dirs: str + save_dirs: dict saving directories + frangi: Frangi filter + + odf: ODF analysis + + tmp: temporary data + in_img: dict input image dictionary - data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) + data: NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) 3D microscopy image - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask - ch_ax: int RGB image channel axis (either 1, 3, or None for grayscale images) + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + fb_ch: int neuronal fibers channel @@ -58,9 +64,18 @@ def parallel_frangi_over_slices(cli_args, save_dirs, in_img): px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] + path: str + path to the 3D microscopy image + name: str name of the 3D microscopy image + fmt: str + format of the 3D microscopy image + + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher + is_vec: bool vector field flag @@ -78,100 +93,88 @@ def parallel_frangi_over_slices(cli_args, save_dirs, in_img): out_img: dict output image dictionary - fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector field + vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - fbr_clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image + clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image + fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - Frangi-enhanced image (fiber probability image) + frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability image) - iso_fbr_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image + iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image - fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - fiber mask image + fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fiber mask image - bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - soma mask image + bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + soma mask image - px_sz: numpy.ndarray (shape=(3,), dtype=float) - output pixel size [μm] + px_sz: numpy.ndarray (shape=(3,), dtype=float) + output pixel size [μm] """ - # skip Frangi filter stage if orientation vectors were directly provided as input + # skip the Frangi filter stage if an orientation vector field was provided as input if in_img['is_vec']: - out_img = dict() - out_img['fbr_vec'] = in_img['data'] - out_img['iso_fbr'] = None - out_img['px_sz'] = None + out_img = {'vec': in_img['data'], 'iso': None, 'px_sz': None} else: - - # get resources configuration - ram, jobs = get_resource_config(cli_args) - - # get Frangi filter configuration - frangi_cfg, px_sz_iso = get_frangi_config(cli_args, in_img) - - # configure the batch of basic image slices analyzed using parallel threads - batch_sz, slc_shp, slc_shp_um, rsz_ratio, ovlp, ovlp_rsz = \ - config_frangi_batch(in_img, frangi_cfg, px_sz_iso, ram=ram, jobs=jobs) - - # generate image range lists - slc_rng, out_slc_shp, tot_slc, batch_sz = \ - generate_slice_lists(slc_shp, in_img['shape'], batch_sz, rsz_ratio, ovlp, in_img['msk_bc']) - - # initialize output arrays - out_img, z_sel = init_frangi_arrays(frangi_cfg, in_img['shape'], out_slc_shp, rsz_ratio, save_dirs['tmp'], - msk_bc=in_img['msk_bc']) - - # print Frangi filter configuration - print_frangi_info(in_img, frangi_cfg, slc_shp_um, tot_slc) - - # parallel Frangi filter-based fiber orientation analysis of microscopy image slices - print_flsh(f"[Parallel(n_jobs={batch_sz})]: Using backend ThreadingBackend with {batch_sz} concurrent workers.") + # initialize Frangi filter stage output + frangi_cfg = get_frangi_config(cli_args, in_img) + out_img = init_frangi_arrays(in_img, frangi_cfg, save_dirs['tmp']) + + # get parallel processing configuration + get_resource_config(cli_args, frangi_cfg) + get_slicing_config(in_img, frangi_cfg) + + # conduct a Frangi-filter-based analysis of fiber orientations using concurrent workers + print_frangi_config(in_img, frangi_cfg) + slc_rng = generate_slice_ranges(in_img, frangi_cfg) + _fa = out_img['fa'] is not None t_start = perf_counter() - with Parallel(n_jobs=batch_sz, prefer='threads', require='sharedmem') as parallel: - parallel( - delayed(frangi_analysis)(s, in_img, out_img, frangi_cfg, ovlp_rsz, rsz_ratio, z_sel, **slc_rng, - print_info=(t_start, batch_sz, tot_slc)) - for s in range(tot_slc)) - - # save output arrays - save_frangi_arrays(save_dirs['frangi'], in_img['name'], out_img, px_sz_iso, ram=ram) + with Parallel(n_jobs=frangi_cfg['batch'], prefer='threads') as parallel: + parallel(delayed(frangi_analysis)(s, in_img, out_img, frangi_cfg, t_start, _fa=_fa) for s in slc_rng) - # add isotropic pixel size - out_img['px_sz'] = px_sz_iso + # save output arrays to file + save_frangi_arrays(save_dirs['frangi'], in_img['name'], out_img, ram=frangi_cfg['ram']) return out_img -def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_rng, out_rng, in_pad, - print_info=None, pad_mode='reflect', bc_thr='yen', ts_thr=0.0001, verbose=10): +def frangi_analysis(rng, in_img, out_img, cfg, t_start, _fa=False): """ Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image. Parameters ---------- - s: int - image slice index + rng: dict + in: np.ndarray + 3D input slice range [px] + + pad: np.ndarray + 3D slice padding range [px] + + out: np.ndarray + 3D output slice range [px] + + bc: np.ndarray + (optional) brain cell soma slice range in_img: dict - input image dictionary (extended) + input image dictionary - data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) + data: NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) 3D microscopy image - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask - ch_ax: int RGB image channel axis (either 1, 3, or None for grayscale images) + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + fb_ch: int neuronal fibers channel @@ -188,9 +191,18 @@ def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_ px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] + path: str + path to the 3D microscopy image + name: str name of the 3D microscopy image + fmt: str + format of the 3D microscopy image + + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher + is_vec: bool vector field flag @@ -254,14 +266,8 @@ def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_ px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - z_rng: int - output z-range in [px] - - bc_ch: int - neuronal bodies channel - - fb_ch: int - myelinated fibers channel + fb_thr: float or str + image thresholding applied to the Frangi filter response msk_bc: bool if True, mask neuronal bodies within @@ -273,78 +279,63 @@ def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_ exp_all: bool export all images - ovlp: int - overlapping range between slices along each axis [px] - - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio - - z_sel: NumPy slice object - selected z-depth range + rsz: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio - in_rng: NumPy slice object - input image range (fibers) + ram: float + maximum RAM available to the Frangi filter stage [B] - bc_rng: NumPy slice object - input image range (neurons) + jobs: int + number of parallel jobs (threads) + used by the Frangi filter stage - out_rng: NumPy slice object - output range + batch: int + slice batch size - in_pad: numpy.ndarray (axis order=(Z,Y,X)) - image padding range + slc_shp: numpy.ndarray (shape=(3,), dtype=int) + shape of the basic image slices + analyzed using parallel threads [px] - print_info: tuple - optional printed information + ovlp: int + overlapping range between image slices along each axis [px] - pad_mode: str - image padding mode + tot_slc: int + total number of image slices - bc_thr: str - thresholding method applied to the neuronal bodies channel + z_out: NumPy slice object + output z-range - ts_thr: float - relative threshold on non-zero tissue pixels + t_start: float + start time [s] - verbose: int - verbosity level (print progress every N=verbose image slices) + _fa: bool + compute fractional anisotropy Returns ------- None """ - # select fiber image slice - fbr_slc, ts_msk_slc = slice_image(in_img['data'], in_rng[s], in_img['ch_ax'], ch=in_img['fb_ch'], - ts_msk=in_img['ts_msk']) + # extract image slice + fbr_slc, ts_msk_slc = \ + slice_image(in_img['data'], rng['in'], in_img['ch_ax'], ch=in_img['fb_ch'], ts_msk=in_img['ts_msk']) # process non-background image slice - is_valid = check_background(fbr_slc, ts_msk=ts_msk_slc, ts_thr=ts_thr) - if is_valid: - orient_slc, frangi_slc, iso_slc, fbr_msk_slc = \ - analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng[s], in_pad[s], pad_mode=pad_mode, ts_msk=ts_msk_slc) - - # (optional) neuronal body masking - if in_img['msk_bc']: - - # get soma image slice - bc_slc, _ = slice_image(in_img['data'], bc_rng[s], in_img['ch_ax'], ch=in_img['bc_ch']) - - # suppress soma contribution - orient_slc, bc_msk_slc = reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng[s], bc_thr) + not_bg = check_background(fbr_slc, ts_msk=ts_msk_slc) + if not_bg: + out_slc = analyze_fibers(fbr_slc, cfg, rng['out'], rng['pad'], ts_msk=ts_msk_slc, _fa=_fa) - # soma mask not available - else: - bc_msk_slc = None + # (optional) brain cell soma masking + if rng['bc'] is not None: + bc_slc, _ = slice_image(in_img['data'], rng['bc'], in_img['ch_ax'], ch=in_img['bc_ch']) + out_slc = reject_brain_cells(bc_slc, out_slc, cfg['rsz'], rng['out']) - # fill memory-mapped output arrays - write_frangi_arrays(out_rng[s], z_sel, iso_slc, frangi_slc, fbr_msk_slc, bc_msk_slc, **orient_slc, **out_img) + # write memory-mapped output arrays + write_frangi_arrays(out_img, out_slc, rng['out'], z_out=cfg['z_out']) - # print progress - if print_info is not None: - print_frangi_progress(print_info, is_valid, verbose=verbose) + print_frangi_progress(t_start, cfg['batch'], cfg['tot_slc'], not_bg) -def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflect', ts_msk=None): +def analyze_fibers(fbr_slc, cfg, out_rng, pad_rng, ts_msk=None, _fa=False): """ Analyze 3D fiber orientations exploiting a Frangi-filter-based unsupervised enhancement of tubular structures. @@ -378,14 +369,8 @@ def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflec px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - z_rng: int - output z-range in [px] - - bc_ch: int - neuronal bodies channel - - fb_ch: int - myelinated fibers channel + fb_thr: float or str + image thresholding applied to the Frangi filter response msk_bc: bool if True, mask neuronal bodies within @@ -397,76 +382,99 @@ def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflec exp_all: bool export all images - ovlp: int - overlapping range between slices along each axis [px] + rsz: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio + ram: float + maximum RAM available to the Frangi filter stage [B] + + jobs: int + number of parallel jobs (threads) + used by the Frangi filter stage + + batch: int + slice batch size + + slc_shp: numpy.ndarray (shape=(3,), dtype=int) + shape of the basic image slices + analyzed using parallel threads [px] + + ovlp: int + overlapping range between image slices along each axis [px] + + tot_slc: int + total number of image slices + + z_out: NumPy slice object + output z-range out_rng: NumPy slice object output range - pad: numpy.ndarray (axis order=(Z,Y,X)) + pad_rng: numpy.ndarray (axis order=(Z,Y,X)) image padding range - pad_mode: str - image padding mode - ts_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask + _fa: bool + compute fractional anisotropy + Returns ------- - orient_slc: dict - slice orientation dictionary + out_slc: dict + slice output dictionary + + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - 3D fiber orientation field + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) - frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - Frangi-enhanced image slice (fiber probability image) + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice - iso_fbr_img: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image slice + ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + tissue mask slice - fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fiber mask image slice + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask slice """ - # preprocess fiber slice (lateral blurring and downsampling) - iso_fbr_slc, pad_rsz, ts_msk_slc_rsz = \ - correct_anisotropy(fbr_slc, rsz_ratio, sigma=cfg['smooth_sd'], pad=pad, ts_msk=ts_msk) + # preprocess fiber slice (adaptive blurring and downsampling) + iso_fbr_slc, pad_rsz, ts_msk_rsz = \ + correct_anisotropy(fbr_slc, cfg['rsz'], sigma=cfg['smooth_sd'], pad=pad_rng, ts_msk=ts_msk) # pad fiber slice if required if pad_rsz is not None: - iso_fbr_slc = np.pad(iso_fbr_slc, pad_rsz, mode=pad_mode) + iso_fbr_slc = np.pad(iso_fbr_slc, pad_rsz, mode='reflect') # pad tissue mask if available - if ts_msk_slc_rsz is not None: - ts_msk_slc_rsz = np.pad(ts_msk_slc_rsz, pad_rsz, mode='constant') + if ts_msk_rsz is not None: + ts_msk_rsz = np.pad(ts_msk_rsz, pad_rsz, mode='constant') # 3D Frangi filter - orient_slc, frangi_slc = \ - frangi_filter(iso_fbr_slc, scales_px=cfg['scales_px'], - alpha=cfg['alpha'], beta=cfg['beta'], gamma=cfg['gamma'], hsv=cfg['hsv_cmap']) + out_slc = frangi_filter(iso_fbr_slc, scales_px=cfg['scales_px'], + alpha=cfg['alpha'], beta=cfg['beta'], gamma=cfg['gamma'], hsv=cfg['hsv_cmap'], _fa=_fa) # crop resulting slices - orient_slc, frangi_slc, iso_fbr_slc, ts_msk_slc_rsz = \ - crop_lst([orient_slc, frangi_slc, iso_fbr_slc, ts_msk_slc_rsz], out_rng, ovlp) + out_slc.update({'iso': iso_fbr_slc, 'ts_msk': ts_msk_rsz}) + ovlp_rsz = np.multiply(cfg['ovlp'] * np.ones((3,)), cfg['rsz']).astype(int) + out_slc = crop_img_dict(out_slc, out_rng, ovlp_rsz) # remove Frangi filter background - orient_slc, fbr_msk_slc = \ - mask_background(frangi_slc, orient_slc, ts_msk=ts_msk_slc_rsz, method=cfg['fb_thr'], invert=False) + out_slc, msk = mask_background(out_slc, method=cfg['fb_thr'], invert=False) + out_slc['fbr_msk'] = msk - return orient_slc, frangi_slc, iso_fbr_slc, fbr_msk_slc + return out_slc -def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'): +def reject_brain_cells(bc_slc, out_slc, rsz_ratio, out_rng, bc_thr='yen'): """ Suppress soma contribution using the optional image channel of neuronal bodies. @@ -476,17 +484,29 @@ def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'): bc_slc: numpy.ndarray (axis order=(Z,Y,X)) brain cell soma image slice - orient_slc: dict - slice orientation dictionary + out_slc: dict + slice output dictionary - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) 3D fiber orientation field - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy + + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) + + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice + + ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + tissue mask slice + + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask slice rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) 3D image resize ratio @@ -495,38 +515,49 @@ def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'): output range bc_thr: str - thresholding method applied to the neuronal bodies channel + thresholding method applied to the channel of brain cell bodies Returns ------- - orient_slc: dict - (masked) slice orientation dictionary + out_slc: dict + (masked) slice output dictionary + + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field + + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - (masked) 3D fiber orientation field + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - (masked) orientation colormap image + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - (masked) fractional anisotropy image + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice - bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) - brain cell mask + ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + tissue mask slice + + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask slice + + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + brain cell mask """ - # resize soma slice (lateral downsampling) + # resize and crop brain cell soma slice iso_bc, _, _ = correct_anisotropy(bc_slc, rsz_ratio) - - # crop isotropic brain cell soma slice iso_bc = crop(iso_bc, out_rng) # mask neuronal bodies - orient_slc, bc_msk = mask_background(iso_bc, orient_slc, method=bc_thr, invert=True) + out_slc, msk = mask_background(out_slc, ref_img=iso_bc, method=bc_thr, invert=True) + out_slc['bc_msk'] = msk - return orient_slc, bc_msk + return out_slc -def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_name, backend='loky', verbose=100): +def parallel_odf_over_scales(cli_args, save_dirs, out_img, img_name): """ Iterate over the required spatial scales and apply the parallel ODF analysis implemented in parallel_odf_on_slices(). @@ -538,57 +569,65 @@ def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_n save_dirs: dict saving directories - ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) - fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector field + frangi: Frangi filter - iso_fbr: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image + odf: ODF analysis - px_sz: numpy.ndarray (axis order=(3,), dtype=float) - pixel size [μm] + tmp: temporary data - img_name: str - name of the input volume image + out_img: dict + output image dictionary - backend: str - backend module employed by joblib.Parallel + vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - verbose: int - joblib verbosity level + clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image + + fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image + + frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability image) + + iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image + + fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fiber mask image + + bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + soma mask image + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + output pixel size [μm] + + img_name: str + name of the input volume image Returns ------- None """ - # print ODF analysis heading - print_odf_info(cli_args.odf_res, cli_args.odf_deg) - - # export all images (optional) - exp_all = cli_args.exp_all + if cli_args.odf_res is not None: + print_odf_info(cli_args.odf_res, cli_args.odf_deg) - # compute spherical harmonics normalization factors (once for all scales) - odf_norm = get_sph_harm_norm_factors(cli_args.odf_deg) + # get pixel size from CLI arguments + # if a vector field was directly provided to Foa3D + px_sz = out_img['px_sz'] + if px_sz is None: + px_sz = (cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy) - # get number of logical cores - num_cpu = get_available_cores() + # parallel ODF analysis of fiber orientation vectors over the spatial scales of interest + batch_sz = min(len(cli_args.odf_res), get_available_cores()) + odf_norm = get_sph_harm_norm_factors(cli_args.odf_deg) + with Parallel(n_jobs=batch_sz, verbose=10, prefer='threads') as parallel: + parallel(delayed(odf_analysis)(out_img['vec'], out_img['iso'], px_sz, save_dirs, img_name, + odf_norm=odf_norm, odf_deg=cli_args.odf_deg, odf_scale_um=s, + exp_all=cli_args.exp_all) for s in cli_args.odf_res) - # get pixel size from CLI arguments if not provided - if px_sz is None: - px_sz = np.array([cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy]) - - # parallel ODF analysis of fiber orientation vectors - # over spatial scales of interest - n_scales = len(cli_args.odf_res) - batch_sz = min(num_cpu, n_scales) - with Parallel(n_jobs=batch_sz, backend=backend, verbose=verbose) as parallel: - parallel(delayed(odf_analysis)(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, - odf_norm=odf_norm, odf_deg=cli_args.odf_deg, odf_scale_um=s, exp_all=exp_all) - for s in cli_args.odf_res) - - # print output directory - print_flsh(f"\nODF and dispersion maps saved to: {save_dirs['odf']}\n") + print_flsh(f"\nODF and dispersion maps saved to: {save_dirs['odf']}\n") def odf_analysis(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_scale_um, odf_norm, odf_deg=6, exp_all=False): @@ -608,7 +647,7 @@ def odf_analysis(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_scale_um, odf save_dirs: dict saving directories - ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) + ('frangi': Frangi filter, 'odf': ODF analysis) img_name: str name of the 3D microscopy image @@ -629,18 +668,14 @@ def odf_analysis(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_scale_um, odf ------- None """ - # get the ODF kernel size in [px] + # initialize ODF stage output arrays odf_scale = int(np.ceil(odf_scale_um / px_sz[0])) - - # initialize ODF analysis output arrays odf, odi, dnst, bg_mrtrix, vec_tensor_eigen = \ - init_odf_arrays(fbr_vec.shape[:-1], save_dirs['tmp'], odf_scale=odf_scale, odf_deg=odf_deg, exp_all=exp_all) - - # generate downsampled background for MRtrix3 mrview - generate_odf_background(bg_mrtrix, fbr_vec, vxl_side=odf_scale, iso_fbr=iso_fbr) + init_odf_arrays(fbr_vec.shape[:-1], save_dirs['tmp'], scale=odf_scale, deg=odf_deg, exp_all=exp_all) - # compute ODF coefficients - odf, dnst = compute_odf_map(fbr_vec, px_sz, odf, odi, dnst, vec_tensor_eigen, odf_scale, odf_norm, odf_deg=odf_deg) + # generate ODF coefficients and down-sampled background for visualization in MRtrix3 + generate_odf_background(bg_mrtrix, fbr_vec, scale=odf_scale, iso_fbr=iso_fbr) + odf, dnst = compute_odf_map(fbr_vec, px_sz, odf, odi, dnst, vec_tensor_eigen, odf_scale, odf_norm, deg=odf_deg) - # save memory maps to TIFF or NIfTI files + # save output arrays to TIFF or NIfTI files save_odf_arrays(save_dirs['odf'], img_name, odf_scale_um, px_sz, odf, bg_mrtrix, dnst, **odi) diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py index 4f65851..9bfaf6e 100644 --- a/foa3d/preprocessing.py +++ b/foa3d/preprocessing.py @@ -2,8 +2,7 @@ from scipy.ndimage import gaussian_filter from skimage.transform import resize -from foa3d.printing import (print_blur, print_flsh, print_output_res, - print_prepro_heading) +from foa3d.printing import print_blur, print_flsh, print_prepro_heading from foa3d.utils import fwhm_to_sigma @@ -34,51 +33,42 @@ def config_anisotropy_correction(px_sz, psf_fwhm): px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) new isotropic pixel size [μm] """ - # set the isotropic pixel resolution equal to the z-sampling step - px_sz_iso = np.max(px_sz) * np.ones(shape=(3,)) + # set the isotropic pixel size to the maximum size along each axis + px_sz_iso = tuple(np.max(px_sz) * np.ones(shape=(3,))) - # detect preprocessing requirement - cndt_1 = not np.all(psf_fwhm == psf_fwhm[0]) - cndt_2 = not np.all(px_sz == px_sz[0]) + # detect preprocessing requirements + crt_1 = not np.all(psf_fwhm == psf_fwhm[0]) + crt_2 = not np.all(px_sz == px_sz[0]) smooth_sigma = None - if cndt_1 or cndt_2: - - # print preprocessing heading + if crt_1 or crt_2: print_prepro_heading() - # adjust PSF anisotropy via lateral Gaussian blurring - if cndt_1: + # anisotropic PSF + if crt_1: - # estimate the PSF variance from input FWHM values [μm²] + # compute the PSF variance from input FWHM values [μm²] + # and tailor the standard deviation of the smoothing Gaussian kernel [px] psf_var = np.square(fwhm_to_sigma(psf_fwhm)) - - # estimate the in-plane filter variance [μm²] - gauss_var = np.max(psf_var) - psf_var - - # ...and the corresponding standard deviation [px] - smooth_sigma_um = np.sqrt(gauss_var) + smooth_sigma_um = np.sqrt(np.max(psf_var) - psf_var) smooth_sigma = np.divide(smooth_sigma_um, px_sz) - # print preprocessing info print_blur(smooth_sigma_um, psf_fwhm) - # print pixel resize info - if cndt_2: - print_output_res(px_sz_iso) + # anisotropic pixel size + if crt_2: + print_flsh(f"Adjusted pixel size [μm]: ({px_sz_iso[0]:.3f}, {px_sz_iso[1]:.3f}, {px_sz_iso[2]:.3f})\n") else: print_flsh() - # skip line else: print_flsh() return smooth_sigma, px_sz_iso -def correct_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', - anti_alias=True, trunc=4, ts_msk=None): +def correct_anisotropy(img, rsz, sigma=None, pad=None, ts_msk=None): """ - Smooth and downsample the raw microscopy image so as to uniform the lateral sizes + Smooth and downsample the raw microscopy image to uniform the lateral sizes of the optical system's PSF and the original voxel size. Parameters @@ -86,7 +76,7 @@ def correct_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', img: numpy.ndarray (axis order=(Z,Y,X)) 3D microscopy image - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + rsz: numpy.ndarray (shape=(3,), dtype=float) 3D resize ratio sigma: numpy.ndarray (shape=(3,), dtype=int) @@ -96,16 +86,7 @@ def correct_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', pad: numpy.ndarray (shape=(3,2), dtype=int) image padding array to be resized - smooth_pad_mode: str - image padding mode adopted for the smoothing Gaussian filter - - anti_alias: bool - if True, apply an antialiasing filter when downsampling - - trunc: int - truncate the Gaussian kernel at this many standard deviations - - tissue_msk: numpy.ndarray (dtype=bool) + ts_msk: numpy.ndarray (dtype=bool) tissue binary mask Returns @@ -120,29 +101,29 @@ def correct_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', resized tissue binary mask """ # no resizing - if np.all(rsz_ratio == 1): + if np.all(rsz == 1): - # resize tissue mask, when available + # generate 3D tissue mask, when available if ts_msk is not None: ts_msk = ts_msk[np.newaxis, ...] ts_msk = np.repeat(ts_msk, img.shape[0], axis=0) return img, pad, ts_msk - # lateral blurring + # adaptive image blurring else: if sigma is not None: - img = gaussian_filter(img, sigma=sigma, mode=pad_mode, truncate=trunc, output=np.float32) + img = gaussian_filter(img, sigma=sigma, mode='reflect', output=np.float32) - # downsampling - iso_shp = np.ceil(np.multiply(np.asarray(img.shape), rsz_ratio)).astype(int) + # adaptive image downsampling + iso_shp = np.ceil(np.multiply(np.asarray(img.shape), rsz)).astype(int) iso_img = np.zeros(shape=iso_shp, dtype=img.dtype) for z in range(iso_shp[0]): iso_img[z, ...] = \ - resize(img[z, ...], output_shape=tuple(iso_shp[1:]), anti_aliasing=anti_alias, preserve_range=True) + resize(img[z, ...], output_shape=tuple(iso_shp[1:]), anti_aliasing=True, preserve_range=True) - # resize padding array accordingly - pad_rsz = np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad)).astype(int) \ + # resize padding range accordingly + pad_rsz = np.floor(np.multiply(np.array([rsz, rsz]).transpose(), pad)).astype(int) \ if pad is not None else None # resize tissue mask, when available diff --git a/foa3d/printing.py b/foa3d/printing.py index 5777e41..568f8e2 100644 --- a/foa3d/printing.py +++ b/foa3d/printing.py @@ -76,7 +76,7 @@ def print_blur(sigma_um, psf_fwhm): f"Adjusted PSF FWHM [μm]: ({psf_sz:.3f}, {psf_sz:.3f}, {psf_sz:.3f})") -def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num): +def print_frangi_config(in_img, cfg): """ Print Frangi filter stage heading. @@ -125,7 +125,7 @@ def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num): item_sz: int image item size [B] - frangi_cfg: dict + cfg: dict Frangi filter configuration alpha: float @@ -149,14 +149,8 @@ def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num): px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - z_rng: int - output z-range in [px] - - bc_ch: int - neuronal bodies channel - - fb_ch: int - myelinated fibers channel + fb_thr: float or str + image thresholding applied to the Frangi filter response msk_bc: bool if True, mask neuronal bodies within @@ -168,47 +162,76 @@ def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num): exp_all: bool export all images - in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=float) - shape of the analyzed image slices [μm] + rsz: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio + + ram: float + maximum RAM available to the Frangi filter stage [B] + + jobs: int + number of parallel jobs (threads) + used by the Frangi filter stage - tot_slc_num: int - total number of analyzed image slices + batch: int + slice batch size + + slc_shp: numpy.ndarray (shape=(3,), dtype=int) + shape of the basic image slices + analyzed using parallel threads [px] + + ovlp: int + overlapping range between image slices along each axis [px] + + tot_slc: int + total number of image slices + + z_out: NumPy slice object + output z-range Returns ------- None """ - # print Frangi filter sensitivity and scales - alpha = frangi_cfg['alpha'] - beta = frangi_cfg['beta'] - gamma = frangi_cfg['gamma'] - - scales_um = np.asarray(frangi_cfg['scales_um']) + # print Frangi filter sensitivity, scales and thresholding method + alpha = cfg['alpha'] + beta = cfg['beta'] + gamma = cfg['gamma'] + scales_um = np.asarray(cfg['scales_um']) if gamma is None: gamma = 'auto' - print_flsh(color_text(0, 191, 255, "\n\n\n3D Frangi Filter\n") + "\nSensitivity\n" + - f"• plate-like \u03B1: {alpha:.1e}\n• blob-like \u03B2: {beta:.1e}\n• background \u03B3: {gamma}\n\n" + - f"Enhanced scales [μm]: {scales_um}\nEnhanced diameters [μm]: {4 * scales_um}\n") + thr = cfg['fb_thr'] + thr_str = "Filter response threshold:" + thr_str = f"{thr_str} {thr:.2f} (global)\n" if isinstance(thr, float) else f"{thr_str} {thr.capitalize()} (local)\n" - # print iterative analysis information - print_slicing_info(in_img['shape_um'], in_slc_shp_um, tot_slc_num, in_img['px_sz'], in_img['item_sz']) + print_flsh(color_text(0, 191, 255, "\n3D Frangi Filter\n") + "\nSensitivity\n" + + f"• plate-like \u03B1: {alpha:.1e}\n• blob-like \u03B2: {beta:.1e}\n• background \u03B3: {gamma}\n\n" + + f"Enhanced scales [μm]: {scales_um}\nEnhanced diameters [μm]: {4 * scales_um}\n\n" + thr_str) - # print neuron masking info - print_soma_masking(in_img['msk_bc']) + # print parallel processing information + batch_sz = cfg['batch'] + slc_shp_um = np.multiply(cfg['px_sz'], cfg['slc_shp']) + print_slicing_info(in_img['shape_um'], slc_shp_um, in_img['px_sz'], in_img['item_sz'], in_img['msk_bc']) + print_flsh(f"[Parallel(n_jobs={batch_sz})]: Using backend ThreadingBackend with {batch_sz} concurrent workers.") -def print_frangi_progress(info, is_valid, verbose=1): +def print_frangi_progress(start_time, batch, tot, not_bg, verbose=5): """ Print Frangi filter progress. Parameters ---------- - info: tuple - info to be printed out + start_time: float + start time [s] + + batch: int + slice batch size + + tot: int + total number of image slices - is_valid: bool - image slice validity flag + not_bg: bool + foreground slice flag verbose: int verbosity level (print info only every "verbose" slices) @@ -221,12 +244,11 @@ def print_frangi_progress(info, is_valid, verbose=1): slc_cnt += 1 # print only every N=verbose image slices - start_time, batch_sz, tot_slc = info - if (slc_cnt % verbose == 0 and is_valid) or slc_cnt == tot_slc: - prog_prc = 100 * slc_cnt / tot_slc + if (slc_cnt % verbose == 0 and not_bg) or slc_cnt == tot: + prog = 100 * slc_cnt / tot _, hrs, mins, secs = elapsed_time(start_time) - print_flsh(f"[Parallel(n_jobs={batch_sz})]:\t{slc_cnt}/{tot_slc} done\t|\t" + - f"elapsed: {hrs} hrs {mins} min {int(secs)} s\t{prog_prc:.1f}%") + print_flsh( + f"[Parallel(n_jobs={batch})]:\t{slc_cnt}/{tot} done\t|\telapsed: {hrs} hr {mins} min {secs} s\t{prog:.1f}%") def print_image_info(in_img): @@ -237,19 +259,16 @@ def print_image_info(in_img): ---------- in_img: dict input image dictionary - ('img_data': image data, 'ts_msk': tissue sample mask, 'ch_ax': channel axis) Returns ------- None """ - # get pixel and PSF sizes + # get image info + ch_ax = in_img['ch_ax'] px_sz = in_img['px_sz'] psf_fwhm = in_img['psf_fwhm'] - # get channel axis (RGB image only) - ch_ax = in_img['ch_ax'] - # get image shape (ignore channel axis) img_shp = in_img['data'].shape if ch_ax is not None: @@ -294,25 +313,8 @@ def print_odf_info(odf_scales_um, odf_degrees): ------- None """ - print_flsh( - color_text(0, 191, 255, "\n3D ODF Analysis") + - f"\n\nResolution [μm]: {odf_scales_um}\nExpansion degrees: {odf_degrees}\n") - - -def print_output_res(px_sz_iso): - """ - Print the adjusted isotropic spatial resolution. - - Parameters - ---------- - px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) - new isotropic pixel size [μm] - - Returns - ------- - None - """ - print_flsh(f"Adjusted pixel size [μm]: ({px_sz_iso[0]:.3f}, {px_sz_iso[1]:.3f}, {px_sz_iso[2]:.3f})\n") + print_flsh(color_text(0, 191, 255, "\n3D ODF Analysis") + + f"\n\nResolution [μm]: {odf_scales_um}\nExpansion degrees: {odf_degrees}\n") def print_pipeline_heading(): @@ -338,7 +340,7 @@ def print_prepro_heading(): "\n Z Y X") -def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): +def print_slicing_info(img_shp_um, slc_shp_um, px_sz, item_sz, msk_bc): """ Print information on the slicing of the basic image sub-volumes processed by the Foa3D tool. @@ -350,15 +352,16 @@ def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): slc_shp_um: numpy.ndarray (shape=(3,), dtype=float) shape of the analyzed image slices [μm] - tot_slc_num: int - total number of analyzed image slices - px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - img_item_sz: int + item_sz: int image item size (in bytes) + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel + Returns ------- None @@ -367,37 +370,17 @@ def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): if np.any(img_shp_um < slc_shp_um): slc_shp_um = img_shp_um - # get image memory size - img_sz = img_item_sz * np.prod(np.divide(img_shp_um, px_sz)) + # get memory sizes + img_sz = item_sz * np.prod(np.divide(img_shp_um, px_sz)) + max_slc_sz = item_sz * np.prod(np.divide(slc_shp_um, px_sz)) - # get slice memory size - max_slc_sz = img_item_sz * np.prod(np.divide(slc_shp_um, px_sz)) - - # print info + # print total image and basic slices information print_flsh("\n Z Y X") print_flsh(f"Total image shape [μm]: ({img_shp_um[0]:.1f}, {img_shp_um[1]:.1f}, {img_shp_um[2]:.1f})\n" + f"Total image size [MB]: {np.ceil(img_sz / 1024**2).astype(int)}\n\n" + f"Image slice shape [μm]: ({slc_shp_um[0]:.1f}, {slc_shp_um[1]:.1f}, {slc_shp_um[2]:.1f})\n" + - f"Image slice size [MB]: {np.ceil(max_slc_sz / 1024**2).astype(int)}\n" + - f"Image slice number: {tot_slc_num}\n") - - -def print_soma_masking(msk_bc): - """ - Print information on the optional masking of neuronal bodies. - - Parameters - ---------- - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel - - Returns - ------- - None - """ - prt = 'Soma mask: ' + f"Image slice size [MB]: {np.ceil(max_slc_sz / 1024**2).astype(int)}\n") if msk_bc: - print_flsh(f'{prt}active\n') + print_flsh('Soma mask: active\n') else: - print(f'{prt}not active\n') + print_flsh('Soma mask: not active\n') diff --git a/foa3d/slicing.py b/foa3d/slicing.py index 32f5516..68f4830 100644 --- a/foa3d/slicing.py +++ b/foa3d/slicing.py @@ -4,7 +4,6 @@ import psutil from numba import njit -from foa3d.printing import print_flsh from foa3d.utils import get_available_cores @@ -45,10 +44,8 @@ def adjust_axis_range(ax_iter, img_shp, slc_per_ax, start, stop, ovlp=0): pad: numpy.ndarray (shape=(2,), dtype=int) lower and upper padding ranges [px] """ - # initialize slice padding array - pad = np.zeros((2,), dtype=np.int64) - # adjust start coordinate + pad = np.zeros((2,), dtype=np.int64) start -= ovlp if start < 0: pad[0] = -start @@ -68,20 +65,28 @@ def adjust_axis_range(ax_iter, img_shp, slc_per_ax, start, stop, ovlp=0): return start, stop, pad -def check_background(img, ts_msk=None, ts_thr=0.0001): +def check_background(img, ts_msk=None, ts_thr=1e-4): """ - Description + Check whether the input image predominantly includes background voxels. Parameters ---------- + img: numpy.ndarray (axis order=(Z,Y,X)) + + ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + tissue reconstruction binary mask + + ts_thr: float + relative threshold on non-zero voxels Returns ------- - is_valid: bool + not_bg: bool + foreground boolean flag """ - is_valid = np.count_nonzero(ts_msk) / np.prod(ts_msk.shape) > ts_thr if ts_msk is not None else np.max(img) != 0 + not_bg = np.count_nonzero(ts_msk) / np.prod(ts_msk.shape) > ts_thr if ts_msk is not None else np.max(img) != 0 - return is_valid + return not_bg @njit(cache=True) @@ -187,7 +192,7 @@ def compute_slice_range(ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip=Fal return rng, pad -def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2): +def compute_overlap(smooth_sigma, frangi_sigma_um, rsz_ratio=None, trunc=2): """ Compute lateral slice extension range for coping with smoothing-related boundary artifacts. @@ -200,7 +205,7 @@ def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2): frangi_sigma_um: numpy.ndarray (dtype=float) analyzed spatial scales [μm] - px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) 3D image resize ratio trunc: int @@ -211,22 +216,17 @@ def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2): ------- ovlp: int overlapping range between image slices along each axis [px] - - ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int) - resized overlapping range between image slices - (if px_rsz_ratio is not None) [px] """ - sigma = np.max(frangi_sigma_um) / px_rsz_ratio + sigma = np.max(frangi_sigma_um) / rsz_ratio if smooth_sigma is not None: sigma = np.concatenate((smooth_sigma, sigma)) ovlp = int(np.ceil(2 * trunc * np.max(sigma)) // 2) - ovlp_rsz = np.multiply(ovlp * np.ones((3,)), px_rsz_ratio).astype(np.int64) if px_rsz_ratio is not None else None - return ovlp, ovlp_rsz + return ovlp -def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0): +def compute_slice_shape(img_shp, item_sz, slc_sz=1e6, ovlp=0): """ Compute basic image chunk shape depending on its maximum size (in bytes). @@ -238,9 +238,6 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0): item_sz: int image item size [B] - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] - slc_sz: float maximum memory size of the basic image slices analyzed using parallel threads [B] @@ -253,11 +250,6 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0): slc_shp: numpy.ndarray (shape=(3,), dtype=int) shape of the basic image slices analyzed using parallel threads [px] - - slc_shp_um: numpy.ndarray (shape=(3,), dtype=float) - shape of the basic image slices - analyzed using parallel threads [μm] - (if px_size is provided) """ if len(img_shp) == 4: item_sz *= 3 @@ -267,129 +259,8 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0): side_z = min(img_shp[0] + tot_ovlp, side_xyz) side_xy = np.floor(np.sqrt(np.abs(slc_sz / side_z / item_sz))).astype(int) slc_shp = np.array([side_z, side_xy, side_xy]) - tot_ovlp - slc_shp_um = np.multiply(slc_shp, px_sz) if px_sz is not None else None - - return slc_shp, slc_shp_um - -def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1): - """ - Compute the size of the basic microscopy image slices fed to the Frangi filtering stage. - - Parameters - ---------- - max_ram: float - available RAM [B] - - mem_growth: float - empirical memory growth factor - - mem_fudge: float - memory fudge factor - - batch_sz: int - slice batch size - - ns: int - number of spatial scales - - Returns - ------- - slc_sz: float - memory size of the basic image slices - analyzed using parallel threads [B] - """ - slc_sz = max_ram / (batch_sz * mem_growth * mem_fudge * ns) - - return slc_sz - - -def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fudge=1.0, jobs=None, ram=None, shp_thr=7): - """ - Compute size and number of the batches of basic microscopy image slices - analyzed in parallel. - - Parameters - ---------- - - - px_sz_iso: int - isotropic pixel size [μm] - - mem_growth: float - empirical memory growth factor - of the Frangi filtering stage - - mem_fudge: float - memory fudge factor - - jobs: int - number of parallel jobs (threads) - used by the Frangi filtering stage - - ram: float - maximum RAM available to the Frangi filtering stage [B] - - shp_thr: int - minimum slice side [px] - - Returns - ------- - batch_sz: int - slice batch size - - in_slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [px] - - in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [μm] - - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio - - ovlp: int - overlapping range between image slices along each axis [px] - - ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int) - resized overlapping range between image slices [px] - """ - # maximum RAM not provided: use all - if ram is None: - ram = psutil.virtual_memory()[1] - - # number of logical cores - num_cpu = get_available_cores() - if jobs is None: - jobs = num_cpu - - # number of spatial scales - ns = len(frangi_cfg['scales_um']) - - # initialize slice batch size - batch_sz = np.min([jobs // ns, num_cpu]).astype(int) + 1 - - # get pixel resize ratio - rsz_ratio = np.divide(in_img['px_sz'], px_sz_iso) - - # compute slice overlap (boundary artifacts suppression) - ovlp, ovlp_rsz = compute_overlap(frangi_cfg['smooth_sd'], frangi_cfg['scales_um'], - px_rsz_ratio=rsz_ratio) - - # compute the shape of basic microscopy image slices - in_slc_shp = np.array([-1]) - while np.any(in_slc_shp < shp_thr): - batch_sz -= 1 - if batch_sz == 0: - raise ValueError( - "Basic image slices do not fit the available RAM: decrease spatial scales and/or maximum scale value!") - else: - slc_sz = compute_slice_size(ram, mem_growth, mem_fudge, batch_sz, ns) - in_slc_shp, in_slc_shp_um = compute_slice_shape(in_img['shape'], in_img['item_sz'], - px_sz=in_img['px_sz'], slc_sz=slc_sz, ovlp=ovlp) - - return batch_sz, in_slc_shp, in_slc_shp_um, rsz_ratio, ovlp, ovlp_rsz + return slc_shp def crop(img, rng, ovlp=None, flip=()): @@ -432,14 +303,14 @@ def crop(img, rng, ovlp=None, flip=()): return img_out -def crop_lst(img_lst, rng, ovlp=None, flip=()): +def crop_img_dict(img_dict, rng, ovlp=None, flip=()): """ Shrink list of image slices at total volume boundaries, for overall shape consistency. Parameters ---------- - img_lst: list - list of images to be cropped + img_dict: list + dictionary of image data to be cropped rng: tuple 3D slice index range @@ -452,114 +323,304 @@ def crop_lst(img_lst, rng, ovlp=None, flip=()): Returns ------- - img_lst: list - list of cropped image slices + img_dict: list + dictionary of cropped image slices """ - for s, img in enumerate(img_lst): - if img is not None: - if isinstance(img, dict): - for key in img.keys(): - img_lst[s][key] = crop(img[key], rng, ovlp=ovlp, flip=flip) - else: - img_lst[s] = crop(img, rng, ovlp=ovlp, flip=flip) + for key in img_dict.keys(): + if img_dict[key] is not None: + img_dict[key] = crop(img_dict[key], rng, ovlp=ovlp, flip=flip) - return img_lst + return img_dict -def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bc=False): +def generate_slice_ranges(in_img, cfg): """ Generate image slice ranges for the Frangi filtering stage. Parameters ---------- - in_slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [px] + in_img: dict + input image dictionary - img_shp: numpy.ndarray (shape=(3,), dtype=int) - total image shape [px] + data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) + 3D microscopy image - batch_sz: int - slice batch size + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - ovlp: int - overlapping range between image slices along each axis [px] + fb_ch: int + neuronal fibers channel - msk_bc: bool - if True, mask neuronal bodies - in the optionally provided image channel + bc_ch: int + brain cell soma channel - Returns - ------- - slc_rng: dict - in_rng: list - list of input slice index ranges + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel + + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + name: str + name of the 3D microscopy image + + is_vec: bool + vector field flag + + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape + + shape_um: numpy.ndarray (shape=(3,), dtype=float) + total image shape [μm] + + item_sz: int + image item size [B] + + cfg: dict + Frangi filter configuration + + alpha: float + plate-like score sensitivity + + beta: float + blob-like score sensitivity + + gamma: float + background score sensitivity + + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] + + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] + + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + bc_ch: int + neuronal bodies channel + + fb_ch: int + myelinated fibers channel + + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel + + hsv_cmap: bool + generate HSV colormap of 3D fiber orientations + + exp_all: bool + export all images + + rsz: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio - in_pad: list - list of slice padding arrays + ram: float + maximum RAM available to the Frangi filter stage [B] - out_rng: list - list of output slice index ranges + jobs: int + number of parallel jobs (threads) + used by the Frangi filter stage - bc_rng: list - (optional) list of neuronal body slice index ranges + batch: int + slice batch size - out_slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the processed image slices [px] + slc_shp: numpy.ndarray (shape=(3,), dtype=int) + shape of the basic image slices + analyzed using parallel threads [px] - tot_slc: int - total number of analyzed image slices + ovlp: int + overlapping range between image slices along each axis [px] - batch_sz: int - adjusted slice batch size + tot_slc: int + total number of image slices + + z_out: NumPy slice object + output z-range + + Returns + ------- + slc_rng: list of dictionaries + in: np.ndarray + 3D input slice range [px] + + pad: np.ndarray + 3D slice padding range [px] + + out: np.ndarray + 3D output slice range [px] + + bc: np.ndarray + (optional) brain cell soma slice range """ - # adjust output shapes - out_slc_shp = np.ceil(np.multiply(px_rsz_ratio, in_slc_shp)).astype(int) - out_img_shp = np.ceil(np.multiply(px_rsz_ratio, img_shp)).astype(int) - - # number of image slices along each axis - slc_per_dim = np.floor(np.divide(img_shp, in_slc_shp)).astype(int) - - # initialize empty range lists - # fill slice range dictionary - slc_rng = dict() - slc_rng['in_rng'] = list() - slc_rng['in_pad'] = list() - slc_rng['out_rng'] = list() - slc_rng['bc_rng'] = list() - tot_slc = np.prod(slc_per_dim) - for i, zyx in enumerate(product(range(slc_per_dim[0]), range(slc_per_dim[1]), range(slc_per_dim[2]))): - - # index ranges of analyzed image slices (with padding) - in_rng, pad = compute_slice_range(zyx, in_slc_shp, img_shp, slc_per_dim, ovlp=ovlp) - slc_rng['in_rng'].append(in_rng) - slc_rng['in_pad'].append(pad) - - # output index ranges + # compute i/o image slice index ranges + out_slc_shp = np.ceil(np.multiply(cfg['rsz'], cfg['slc_shp'])).astype(int) + out_img_shp = np.ceil(np.multiply(cfg['rsz'], in_img['shape'])).astype(int) + slc_per_dim = np.floor(np.divide(in_img['shape'], cfg['slc_shp'])).astype(int) + slc_rng = [] + for zyx in product(range(slc_per_dim[0]), range(slc_per_dim[1]), range(slc_per_dim[2])): + in_rng, pad = compute_slice_range(zyx, cfg['slc_shp'], in_img['shape'], slc_per_dim, ovlp=cfg['ovlp']) out_rng, _ = compute_slice_range(zyx, out_slc_shp, out_img_shp, slc_per_dim) - slc_rng['out_rng'].append(out_rng) # (optional) neuronal body channel - if msk_bc: - bc_rng, _ = compute_slice_range(zyx, in_slc_shp, img_shp, slc_per_dim) - slc_rng['bc_rng'].append(bc_rng) + if in_img['msk_bc']: + bc_rng, _ = compute_slice_range(zyx, cfg['slc_shp'], in_img['shape'], slc_per_dim) else: - slc_rng['bc_rng'].append(None) + bc_rng = None + + slc_rng.append({'in': in_rng, 'out': out_rng, 'pad': pad, 'bc': bc_rng}) + + return slc_rng + + +def get_slicing_config(in_img, frangi_cfg, mem_growth=149.7, shp_thr=7): + """ + Compute size and number of the batches of basic microscopy image slices + analyzed in parallel. + + Parameters + ---------- + in_img: dict + input image dictionary + + data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X)) + 3D microscopy image + + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) + + fb_ch: int + neuronal fibers channel + + bc_ch: int + brain cell soma channel + + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - # show progress - print_flsh(f"Generating image slice ranges:\t{100 * (i+1) / tot_slc:.1f}%", end='\r') + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - # total number of slices - tot_slc = len(slc_rng['in_rng']) + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + name: str + name of the 3D microscopy image + + is_vec: bool + vector field flag + + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape + + shape_um: numpy.ndarray (shape=(3,), dtype=float) + total image shape [μm] + + item_sz: int + image item size [B] + + frangi_cfg: dict + Frangi filter configuration + + alpha: float + plate-like score sensitivity + + beta: float + blob-like score sensitivity + + gamma: float + background score sensitivity + + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] + + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] + + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + bc_ch: int + neuronal bodies channel + + fb_ch: int + myelinated fibers channel + + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel + + hsv_cmap: bool + generate HSV colormap of 3D fiber orientations + + exp_all: bool + export all images + + rsz: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio + + ram: float + maximum RAM available to the Frangi filter stage [B] + + jobs: int + number of parallel jobs (threads) + used by the Frangi filter stage + + mem_growth: float + empirical memory growth factor + of the Frangi filtering stage + + shp_thr: int + minimum slice side [px] + + Returns + ------- + None + """ + # maximum RAM not provided: use all + ram = frangi_cfg['ram'] + if ram is None: + ram = psutil.virtual_memory()[1] + ram /= mem_growth + + # number of logical cores + jobs = frangi_cfg['jobs'] + num_cpu = get_available_cores() + if jobs is None: + jobs = num_cpu + + # compute the shape of basic microscopy image slices + slc_shp = np.array([-1]) + batch_sz = np.min([jobs, num_cpu]).astype(int) + 1 + ovlp = compute_overlap(frangi_cfg['smooth_sd'], frangi_cfg['scales_um'], rsz_ratio=frangi_cfg['rsz']) + while np.any(slc_shp < shp_thr): + batch_sz -= 1 + if batch_sz == 0: + raise ValueError( + "Basic image slices do not fit the available RAM: decrease spatial scales and/or maximum scale value!") - # adjust slice batch size - if batch_sz > tot_slc: - batch_sz = tot_slc + slc_shp = compute_slice_shape(in_img['shape'], in_img['item_sz'], slc_sz=ram/batch_sz, ovlp=ovlp) - return slc_rng, out_slc_shp, tot_slc, batch_sz + # update Frangi filter configuration dictionary + tot_slc = np.prod(np.floor(np.divide(in_img['shape'], slc_shp)).astype(int)) + frangi_cfg.update({'batch': min(batch_sz, tot_slc), 'slc_shp': slc_shp, 'ovlp': ovlp, 'tot_slc': tot_slc}) def slice_image(img, rng, ch_ax, ch, ts_msk=None): diff --git a/foa3d/spharm.py b/foa3d/spharm.py index 2e36dd4..b4e9a08 100644 --- a/foa3d/spharm.py +++ b/foa3d/spharm.py @@ -279,6 +279,8 @@ def sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm): return norm[1] * 3 * sin_theta * cos_theta * np.cos(phi) elif order == 2: return norm[2] * 3 * sin_theta**2 * np.cos(2 * phi) + else: + raise ValueError('Invalid spherical harmonics order!') @njit(cache=True) @@ -301,6 +303,8 @@ def sph_harm_degree_4(order, phi, sin_theta, cos_theta, norm): return norm[3] * 105 * sin_theta**3 * cos_theta * np.cos(3 * phi) elif order == 4: return norm[4] * 105 * sin_theta**4 * np.cos(4 * phi) + else: + raise ValueError('Invalid spherical harmonics order!') @njit(cache=True) @@ -333,6 +337,8 @@ def sph_harm_degree_6(order, phi, sin_theta, cos_theta, norm): return norm[5] * 10395 * sin_theta**5 * cos_theta * np.cos(5 * phi) elif order == 6: return norm[6] * 10395 * sin_theta**6 * np.cos(6 * phi) + else: + raise ValueError('Invalid spherical harmonics order!') @njit(cache=True) @@ -378,6 +384,8 @@ def sph_harm_degree_8(order, phi, sin_theta, cos_theta, norm): return norm[7] * 2027025 * sin_theta**7 * cos_theta * np.cos(7 * phi) elif order == 8: return norm[8] * 2027025 * sin_theta**8 * np.cos(8 * phi) + else: + raise ValueError('Invalid spherical harmonics order!') @njit(cache=True) @@ -442,3 +450,5 @@ def sph_harm_degree_10(order, phi, sin_theta, cos_theta, norm): return norm[9] * 654729075 * sin_theta**9 * cos_theta * np.cos(9 * phi) elif order == 10: return norm[10] * 654729075 * sin_theta**10 * np.cos(10 * phi) + else: + raise ValueError('Invalid spherical harmonics order!') diff --git a/foa3d/utils.py b/foa3d/utils.py index 0695187..2ad4821 100644 --- a/foa3d/utils.py +++ b/foa3d/utils.py @@ -1,5 +1,6 @@ import gc import tempfile + from multiprocessing import cpu_count from os import environ, path, unlink from shutil import rmtree @@ -14,6 +15,29 @@ threshold_yen) +def ceil_to_multiple(number, multiple): + """ + Round up number to the nearest multiple. + + Parameters + ---------- + number: + number to be rounded + + multiple: + the input number will be rounded + to the nearest multiple higher than this value + + Returns + ------- + rounded: + rounded up number + """ + rounded = multiple * np.ceil(number / multiple) + + return rounded + + def create_background_mask(img, method='yen', black_bg=False): """ Compute background mask. @@ -57,7 +81,7 @@ def create_background_mask(img, method='yen', black_bg=False): return bg_msk -def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mode='w+'): +def create_memory_map(shape, dtype, name='tmp', tmp=None, arr=None, mmap_mode='w+'): """ Create a memory-map to an array stored in a binary file on disk. @@ -72,7 +96,7 @@ def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mod name: str optional temporary filename - tmp_dir: str + tmp: str temporary file directory arr: numpy.ndarray @@ -86,12 +110,12 @@ def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mod mmap: NumPy memory map memory-mapped array """ - if tmp_dir is None: - tmp_dir = tempfile.mkdtemp() - mmap_path = path.join(tmp_dir, name + '.mmap') + if tmp is None: + tmp = tempfile.mkdtemp() + mmap_path = path.join(tmp, name + '.mmap') if path.exists(mmap_path): - unlink(mmap_path) + unlink(mmap_path) if arr is None: _ = open(mmap_path, mode='w+') @@ -106,83 +130,7 @@ def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mod return mmap -def detect_ch_axis(img): - """ - Detect image channel axis. - - Parameters - ---------- - img: numpy.ndarray - 3D microscopy image - - Returns - ------- - ch_ax: int - channel axis (either 1 or 3) - """ - if len(img.shape) < 4: - return None - else: - ch_ax = (np.array(img.shape) == 3).nonzero()[0] - ch_ax = ch_ax[np.logical_or(ch_ax == 1, ch_ax == 3)] - if len(ch_ax) != 1: - raise ValueError("Ambiguous image axes order: could not determine channel axis!") - else: - ch_ax = ch_ax[0] - - return ch_ax - - -def get_available_cores(): - """ - Return the number of available logical cores. - - Returns - ------- - num_cpu: int - number of available cores - """ - num_cpu = environ.pop('OMP_NUM_THREADS', default=None) - num_cpu = cpu_count() if num_cpu is None else int(num_cpu) - - return num_cpu - - -def get_item_size(dtype): - """ - Get the item size in bytes of a data type. - - Parameters - ---------- - dtype: str - data type identifier - - Returns - ------- - item_sz: int - item size in bytes - """ - # data type lists - lst_1 = ['uint8', 'int8'] - lst_2 = ['uint16', 'int16', 'float16', np.float16] - lst_3 = ['uint32', 'int32', 'float32', np.float32] - lst_4 = ['uint64', 'int64', 'float64', np.float64] - - if dtype in lst_1: - item_sz = 1 - elif dtype in lst_2: - item_sz = 2 - elif dtype in lst_3: - item_sz = 4 - elif dtype in lst_4: - item_sz = 8 - else: - raise ValueError("Unsupported data type!") - - return item_sz - - -def delete_tmp_folder(tmp_dir, tmp_data): +def delete_tmp_data(tmp_dir, tmp_data): """ Delete temporary folder. @@ -190,7 +138,7 @@ def delete_tmp_folder(tmp_dir, tmp_data): ---------- tmp_dir: str path to temporary folder to be removed - + tmp_data: tuple temporary data dictionaries @@ -209,6 +157,33 @@ def delete_tmp_folder(tmp_dir, tmp_data): pass +def detect_ch_axis(img): + """ + Detect image channel axis. + + Parameters + ---------- + img: numpy.ndarray + 3D microscopy image + + Returns + ------- + ch_ax: int + channel axis (either 1 or 3) + """ + if len(img.shape) < 4: + return None + + ch_ax = (np.array(img.shape) == 3).nonzero()[0] + ch_ax = ch_ax[np.logical_or(ch_ax == 1, ch_ax == 3)] + if len(ch_ax) != 1: + raise ValueError("Ambiguous image axes order: could not determine channel axis!") + + ch_ax = ch_ax[0] + + return ch_ax + + def divide_nonzero(nd_array1, nd_array2, new_val=1e-10): """ Divide two arrays handling zero denominator values. @@ -256,7 +231,7 @@ def elapsed_time(start_time): mins: int minutes - secs: float + secs: int seconds """ stop_time = perf_counter() @@ -266,7 +241,7 @@ def elapsed_time(start_time): hrs = int(secs // 3600) secs %= 3600 mins = int(secs // 60) - secs %= 60 + secs = int(secs % 60) return tot, hrs, mins, secs @@ -291,6 +266,43 @@ def fwhm_to_sigma(fwhm): return sigma +def get_available_cores(): + """ + Return the number of available logical cores. + + Returns + ------- + num_cpu: int + number of available cores + """ + num_cpu = environ.pop('OMP_NUM_THREADS', default=None) + num_cpu = cpu_count() if num_cpu is None else int(num_cpu) + + return num_cpu + + +def get_config_label(cli_args): + """ + Generate the output filename prefix including + pipeline configuration information. + + Parameters + ---------- + cli_args: see ArgumentParser.parse_args + updated namespace of command line arguments + + Returns + ------- + cfg_lbl: str + Frangi filter configuration label + """ + cfg_lbl = f'a{cli_args.alpha}_b{cli_args.beta}_g{cli_args.gamma}_t{cli_args.fb_thr}' + for s in cli_args.scales: + cfg_lbl += f'_s{s}' + + return cfg_lbl + + def get_item_bytes(data): """ Retrieve data item size in bytes. @@ -305,7 +317,6 @@ def get_item_bytes(data): bts: int item size in bytes """ - # type byte size try: bts = int(np.iinfo(data.dtype).bits / 8) except ValueError: @@ -314,27 +325,69 @@ def get_item_bytes(data): return bts -def get_config_label(cli_args): +def get_item_size(dtype): """ - Generate the output filename prefix including - pipeline configuration information. + Get the item size in bytes of a data type. Parameters ---------- - cli_args: see ArgumentParser.parse_args - updated namespace of command line arguments + dtype: str + data type identifier Returns ------- - cfg_lbl: str - Frangi filter configuration label + item_sz: int + item size in bytes """ - cfg_lbl = '_s' - for s in cli_args.scales: - cfg_lbl += f'{s}_' - cfg_lbl = f'a{cli_args.alpha}_b{cli_args.beta}_g{cli_args.gamma}_t{cli_args.fb_thr}{cfg_lbl}' + lst_1 = ['uint8', 'int8'] + lst_2 = ['uint16', 'int16', 'float16', np.float16] + lst_3 = ['uint32', 'int32', 'float32', np.float32] + lst_4 = ['uint64', 'int64', 'float64', np.float64] - return cfg_lbl + if dtype in lst_1: + item_sz = 1 + elif dtype in lst_2: + item_sz = 2 + elif dtype in lst_3: + item_sz = 4 + elif dtype in lst_4: + item_sz = 8 + else: + raise ValueError("Unsupported data type!") + + return item_sz + + +def hsv_orient_cmap(vec_img): + """ + Compute HSV colormap of vector orientations from 3D vector field. + + Parameters + ---------- + vec_img: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + orientation vectors + + Returns + ------- + rgb_map: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation color map + """ + # compute the in-plane versor length + vy, vx = (vec_img[..., 1], vec_img[..., 2]) + vxy_abs = np.sqrt(np.square(vx) + np.square(vy)) + vxy_abs = divide_nonzero(vxy_abs, np.max(vxy_abs)) + + # compute the in-plane angular orientation + vxy_ang = normalize_angle(np.arctan2(vy, vx), lower=0, upper=np.pi, dtype=np.float32) + vxy_ang = np.divide(vxy_ang, np.pi) + + # generate HSV orientation colormap + rgb_map = np.zeros(shape=tuple(list(vec_img.shape[:-1]) + [3]), dtype=np.uint8) + for z in range(vec_img.shape[0]): + hsv_map = np.stack((vxy_ang[z], vxy_abs[z], vxy_abs[z]), axis=-1) + rgb_map[z] = (255.0 * hsv_to_rgb(hsv_map)).astype(np.uint8) + + return rgb_map def normalize_angle(angle, lower=0.0, upper=360.0, dtype=None): @@ -343,8 +396,8 @@ def normalize_angle(angle, lower=0.0, upper=360.0, dtype=None): Parameters ---------- - angle: float or list (dtype=float) - angular value(s) to be normalized (in degrees) + angle: numpy.ndarray (dtype=float) + angular values to be normalized (in degrees) lower: float lower limit (default: 0.0) @@ -357,49 +410,33 @@ def normalize_angle(angle, lower=0.0, upper=360.0, dtype=None): Returns ------- - norm_angle: float or list (dtype=float) - angular value(s) (in degrees) normalized within [lower, upper) + angle: numpy.ndarray (dtype=float) + angular values (in degrees) normalized within [lower, upper) Raises ------ ValueError if lower >= upper """ - # convert to array if needed - is_list = False - if np.isscalar(angle): - angle = np.array(angle) - elif isinstance(angle, list): - angle = np.array(angle) - is_list = True - - # check limits + # check input variables if lower >= upper: - raise ValueError(f"Invalid lower and upper limits: ({lower}, {upper})") - - # view - norm_angle = angle + raise ValueError(f"Invalid lower and upper angular limits: ({lower}, {upper})") - # correction 1 - c1 = np.logical_or(angle > upper, angle == lower) - angle[c1] = lower + abs(angle[c1] + upper) % (abs(lower) + abs(upper)) + # apply corrections + dvsr = abs(lower) + abs(upper) + corr_1 = np.logical_or(angle > upper, angle == lower) + angle[corr_1] = lower + np.abs(angle[corr_1] + upper) % dvsr - # correction 2 - c2 = np.logical_or(angle < lower, angle == upper) - angle[c2] = upper - abs(angle[c2] - lower) % (abs(lower) + abs(upper)) + corr_2 = np.logical_or(angle < lower, angle == upper) + angle[corr_2] = upper - np.abs(angle[corr_2] - lower) % dvsr - # correction 3 angle[angle == upper] = lower # cast to desired data type if dtype is not None: - norm_angle = norm_angle.astype(dtype) - - # convert back to list - if is_list: - norm_angle = list(norm_angle) + angle = angle.astype(dtype) - return norm_angle + return angle def normalize_image(img, min_val=None, max_val=None, max_out_val=255.0, dtype=np.uint8): @@ -410,10 +447,10 @@ def normalize_image(img, min_val=None, max_val=None, max_out_val=255.0, dtype=np ---------- img: numpy.ndarray input image - + min_val: float minimum input value - + max_val: float maximum input value @@ -446,67 +483,41 @@ def normalize_image(img, min_val=None, max_val=None, max_out_val=255.0, dtype=np return norm_img -def hsv_orient_cmap(vec_img): +def rgb_orient_cmap(vec_img, minimum=0, stretch=1, q=8): """ - Compute HSV colormap of vector orientations from 3D vector field. + Compute RGB colormap of orientation vector components from 3D vector field. Parameters ---------- vec_img: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - orientation vectors + n-dimensional array of orientation vectors + + minimum: int + intensity that should be mapped to black (a scalar or array for R, G, B) + + stretch: int + linear stretch of the image + + q: int + asinh softening parameter Returns ------- rgb_map: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) orientation color map """ - # extract planar components - vy, vx = (vec_img[..., 1], vec_img[..., 2]) - - # compute the in-plane versor length - vxy_abs = np.sqrt(np.square(vx) + np.square(vy)) - vxy_abs = divide_nonzero(vxy_abs, np.max(vxy_abs)) - - # compute the in-plane angular orientation - vxy_ang = normalize_angle(np.arctan2(vy, vx), lower=0, upper=np.pi, dtype=np.float32) - vxy_ang = np.divide(vxy_ang, np.pi) - # initialize colormap - rgb_map = np.zeros(shape=tuple(list(vec_img.shape[:-1]) + [3]), dtype=np.uint8) + vec_img = np.abs(vec_img) + rgb_map = np.zeros(shape=vec_img.shape, dtype=np.uint8) for z in range(vec_img.shape[0]): - # generate HSV colormap slice by slice - hsv_map = np.stack((vxy_ang[z], vxy_abs[z], vxy_abs[z]), axis=-1) - - # convert to RGB map with 8-bit precision - rgb_map[z] = (255.0 * hsv_to_rgb(hsv_map)).astype(np.uint8) + # generate RGB orientation colormap + img_r, img_g, img_b = (vec_img[z, :, :, 2], vec_img[z, :, :, 1], vec_img[z, :, :, 0]) + rgb_map[z] = make_lupton_rgb(img_r, img_g, img_b, minimum=minimum, stretch=stretch, Q=q) return rgb_map -def ceil_to_multiple(number, multiple): - """ - Round up number to the nearest multiple. - - Parameters - ---------- - number: - number to be rounded - - multiple: - the input number will be rounded - to the nearest multiple higher than this value - - Returns - ------- - rounded: - rounded up number - """ - rounded = multiple * np.ceil(number / multiple) - - return rounded - - def transform_axes(nd_array, flipped=None, swapped=None, expand=None): """ Manipulate axes and dimensions of the input array. @@ -543,40 +554,3 @@ def transform_axes(nd_array, flipped=None, swapped=None, expand=None): nd_array = np.expand_dims(nd_array, axis=expand) return nd_array - - -def rgb_orient_cmap(vec_img, minimum=0, stretch=1, q=8): - """ - Compute RGB colormap of orientation vector components from 3D vector field. - - Parameters - ---------- - vec_img: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - n-dimensional array of orientation vectors - - minimum: int - intensity that should be mapped to black (a scalar or array for R, G, B) - - stretch: int - linear stretch of the image - - q: int - asinh softening parameter - - Returns - ------- - rgb_map: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation color map - """ - # take absolute value - vec_img = np.abs(vec_img) - - # initialize colormap - rgb_map = np.zeros(shape=vec_img.shape, dtype=np.uint8) - for z in range(vec_img.shape[0]): - - # generate colormap slice by slice - img_r, img_g, img_b = (vec_img[z, :, :, 2], vec_img[z, :, :, 1], vec_img[z, :, :, 0]) - rgb_map[z] = make_lupton_rgb(img_r, img_g, img_b, minimum=minimum, stretch=stretch, Q=q) - - return rgb_map diff --git a/requirements.in b/requirements.in index f89bf76..fe1f595 100644 --- a/requirements.in +++ b/requirements.in @@ -1,5 +1,5 @@ astropy -joblib +joblib==1.4.2 matplotlib nibabel numba diff --git a/requirements.txt b/requirements.txt index cd43c64..4f1022d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,68 +1,67 @@ # -# This file is autogenerated by pip-compile with Python 3.9 +# This file is autogenerated by pip-compile with Python 3.8 # by the following command: # # pip-compile requirements.in # -about-time==3.1.1 - # via alive-progress -astropy==5.1 +astropy==5.2.2 # via -r requirements.in -cachetools==5.3.1 +cachetools==5.5.0 # via zetastitcher -clarabel==0.6.0 +clarabel==0.9.0 # via cvxpy coloredlogs==15.0.1 # via zetastitcher -contourpy==1.0.5 +contourpy==1.1.1 # via matplotlib -cvxpy==1.4.1 +cvxpy==1.5.2 # via zetastitcher -cycler==0.11.0 +cycler==0.12.1 # via matplotlib -daqp==0.5.1 - # via qpsolvers -ecos==2.0.12 - # via - # cvxpy - # qpsolvers -fonttools==4.37.4 +ecos==2.0.14 + # via cvxpy +fonttools==4.55.0 # via matplotlib -grapheme==0.6.0 - # via alive-progress humanfriendly==10.0 # via coloredlogs -humanize==4.8.0 +humanize==4.10.0 # via zetastitcher -imageio==2.22.1 +imageio==2.35.1 # via # pims # scikit-image # zetastitcher -joblib==1.2.0 +importlib-metadata==8.5.0 + # via numba +importlib-resources==6.4.5 + # via + # matplotlib + # nibabel +joblib==1.4.2 # via -r requirements.in -kiwisolver==1.4.4 +kiwisolver==1.4.7 # via matplotlib -llvmlite==0.39.1 +lazy-loader==0.4 + # via scikit-image +llvmlite==0.41.1 # via numba -matplotlib==3.6.0 +matplotlib==3.7.5 # via -r requirements.in -networkx==2.8.7 +networkx==3.1 # via # scikit-image # zetastitcher -nibabel==4.0.2 +nibabel==5.2.1 # via -r requirements.in -numba==0.56.2 +numba==0.58.1 # via -r requirements.in -numpy==1.23.3 +numpy==1.24.4 # via # astropy # clarabel # contourpy # cvxpy # ecos - # h5py # imageio # matplotlib # nibabel @@ -73,91 +72,89 @@ numpy==1.23.3 # pims # pyerfa # pywavelets + # qdldl # qpsolvers # scikit-image # scipy # scs # tifffile # zetastitcher -opencv-python==4.8.1.78 +opencv-python==4.10.0.84 # via zetastitcher -osqp==0.6.3 - # via - # cvxpy - # qpsolvers -packaging==21.3 +osqp==0.6.7.post3 + # via cvxpy +packaging==24.2 # via # astropy + # lazy-loader # matplotlib # nibabel + # pims # scikit-image -pandas==2.1.1 +pandas==2.0.3 # via zetastitcher -pillow==9.2.0 +pillow==10.4.0 # via # imageio # matplotlib # scikit-image -pims==0.6.1 +pims==0.7 # via zetastitcher -psutil==5.9.4 +psutil==6.1.0 # via # -r requirements.in # zetastitcher -pybind11==2.11.1 - # via cvxpy -pyerfa==2.0.0.1 +pyerfa==2.0.0.3 # via astropy -pyparsing==3.0.9 - # via - # matplotlib - # packaging -pyreadline3==3.4.1 +pyparsing==3.1.4 + # via matplotlib +pyreadline3==3.5.4 # via humanfriendly -python-dateutil==2.8.2 +python-dateutil==2.9.0.post0 # via # matplotlib # pandas -pytz==2023.3.post1 +pytz==2024.2 # via pandas pywavelets==1.4.1 # via scikit-image -pyyaml==6.0 +pyyaml==6.0.2 # via # astropy # zetastitcher -qdldl==0.1.7.post0 +qdldl==0.1.7.post4 # via osqp -qpsolvers==4.0.0 +qpsolvers==4.4.0 # via zetastitcher -scikit-image==0.19.3 +scikit-image==0.21.0 # via -r requirements.in -scipy==1.9.1 +scipy==1.10.1 # via # clarabel # cvxpy # ecos # osqp + # qdldl # qpsolvers # scikit-image # scs # zetastitcher -scs==3.2.3 - # via - # cvxpy - # qpsolvers +scs==3.2.7 + # via cvxpy six==1.16.0 # via python-dateutil slicerator==1.1.0 # via pims -tifffile==2022.8.12 +tifffile==2023.7.10 # via + # pims # scikit-image # zetastitcher -tzdata==2023.3 +tzdata==2024.2 # via pandas zetastitcher==0.7.0.post1 # via -r requirements.in - -# The following packages are considered to be unsafe in a requirements file: -# setuptools +zipp==3.20.2 + # via + # importlib-metadata + # importlib-resources