From a61a4f456da75116371c1755355b1a6e1736b7b3 Mon Sep 17 00:00:00 2001 From: msorelli Date: Thu, 14 Nov 2024 17:20:23 +0100 Subject: [PATCH] Use generators and simplify code --- foa3d/__main__.py | 8 +- foa3d/frangi.py | 521 +++++++++++++++++--------------------- foa3d/input.py | 537 ++++++++++++++++++++------------------- foa3d/odf.py | 141 +++++------ foa3d/output.py | 63 +++-- foa3d/pipeline.py | 555 ++++++++++++++++++++++------------------- foa3d/preprocessing.py | 51 ++-- foa3d/printing.py | 88 ++++--- foa3d/slicing.py | 403 ++++++++++++++++++------------ foa3d/spharm.py | 10 + foa3d/utils.py | 156 +++++------- requirements.in | 2 +- requirements.txt | 123 +++++---- 13 files changed, 1354 insertions(+), 1304 deletions(-) diff --git a/foa3d/__main__.py b/foa3d/__main__.py index ff736b9..96a25d5 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): @@ -13,10 +13,10 @@ def foa3d(cli_args): 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']) + 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(): diff --git a/foa3d/frangi.py b/foa3d/frangi.py index 39d9bcc..fbc0196 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) + srt_eigval, srt_eigvec = sort_eigen(eigenval, eigenvec) + dom_eigvec = srt_eigvec[..., 0] - # select the eigenvectors related to dominant eigenvalues - dom_eigenvec = sorted_eigenvec[..., 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) - - # stack eigenvalues list - eigenval = np.stack(eigenval, axis=-1) + return frangi_img, eigvec, eigval - 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] + in_img: dict + input image dictionary - img: numpy.ndarray (axis order=(Z,Y,X)) - 3D microscopy image - - 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 - - clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) - initialized orientation colormap image + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) + initialized fiber orientation 3D image - fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized fractional anisotropy image + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + initialized orientation colormap image - frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized Frangi-enhanced image + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized fractional anisotropy image - iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized fiber image (isotropic resolution) + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized Frangi-enhanced image - fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized fiber mask image + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized fiber image (isotropic resolution) - bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized soma mask image + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + initialized fiber 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)) + 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_img = create_memory_map(tot_shp, dtype='uint8', name='iso', dir=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_img = create_memory_map(tot_shp, dtype='uint8', name='frangi', dir=tmp_dir) + fbr_msk = create_memory_map(tot_shp, dtype='uint8', name='fbr_msk', dir=tmp_dir) + fa_img = create_memory_map(tot_shp, dtype='float32', name='fa', dir=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', dir=tmp_dir) else: bc_msk = None else: frangi_img, fbr_msk, fa_img, bc_msk = (None, None, None, 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) + vec_shp = tot_shp + (len(img_shp),) + fbr_vec_img = create_memory_map(vec_shp, dtype='float32', name='vec', dir=tmp_dir) + fbr_vec_clr = create_memory_map(vec_shp, dtype='uint8', name='clr', dir=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_img, 'clr': fbr_vec_clr, 'fa': fa_img, 'frangi': frangi_img, + 'iso': iso_fbr_img, '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): """ Apply 3D Frangi filter to 3D microscopy image. @@ -522,74 +452,92 @@ 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 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) 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 - orient: dict - slice orientation dictionary + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy - vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - fiber orientation vector slice + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) - clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap slice + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - fractional anisotropy 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 + + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + brain cell mask 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 +545,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: + 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,17 +595,12 @@ 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 @@ -707,83 +643,96 @@ def sort_eigen(eigenval, eigenvec, axis=-1): return sorted_eigenval, sorted_eigenvec -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, 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 + + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - z_sel: NumPy slice object - selected z-depth range + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - iso_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image slice + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - frangi_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - Frangi-enhanced image slice + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability image) - fbr_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fiber mask image slice + iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image + + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask image + + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + soma mask image + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + output pixel size [μm] - bc_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - soma mask image slice + out_slc: dict + slice output dictionary - vec_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - fiber orientation vector image slice + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - clr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - orientation colormap image slice + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - fractional anisotropy image slice + 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) - - # optional output images: fractional anisotropy - if fa is not None: - fa[out_rng] = fa_slc[z_sel, ...] - - # 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) - - # 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) + # write valid slice output + if out_slc is not None: + out_rng = out_slc['rng'] + vec_out_rng = tuple(np.append(out_rng, slice(0, 3, 1))) + out_img['vec'][vec_out_rng] = out_slc['vec'][z_out, ...] + out_img['clr'][vec_out_rng] = out_slc['clr'][z_out, ...] + out_img['iso'][out_rng] = out_slc['iso'][z_out, ...].astype(np.uint8) + + # optional output images: fractional anisotropy + if out_img['fa'] is not None: + out_img['fa'][out_rng] = out_slc['fa'][z_out, ...] + + # optional output images: Frangi filter response + if out_img['frangi'] is not None: + out_img['frangi'][out_rng] = (255 * out_slc['frangi'][z_out, ...]).astype(np.uint8) + + # optional output images: fiber mask + if out_img['fbr_msk'] is not None: + out_img['fbr_msk'][out_rng] = (255 * (1 - out_slc['fbr_msk'][z_out, ...])).astype(np.uint8) + + # optional output images: neuronal soma mask + if out_img['bc_msk'] is not None: + out_img['bc_msk'][out_rng] = (255 * out_slc['bc_msk'][z_out, ...]).astype(np.uint8) diff --git a/foa3d/input.py b/foa3d/input.py index a3f61e9..9ecbf85 100644 --- a/foa3d/input.py +++ b/foa3d/input.py @@ -73,7 +73,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]') @@ -110,7 +110,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 +120,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 +142,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 - - shape: numpy.ndarray (shape=(3,), dtype=int) - total image shape + fmt: str + format of the 3D microscopy image - shape_um: numpy.ndarray (shape=(3,), dtype=float) - total image shape [μm] + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher - 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 +184,41 @@ def get_image_info(cli_args): Returns ------- - img_path: str - path to the 3D microscopy image + in_img: dict + input image dictionary + fb_ch: int + neuronal fibers channel - img_name: str - name of the 3D microscopy image + bc_ch: int + brain cell soma channel - img_fmt: str - format of the 3D microscopy image + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] + 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] + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - is_tiled: bool - True for tiled reconstructions aligned using ZetaStitcher + path: str + path to the 3D microscopy image - is_vec: bool - True when pre-estimated fiber orientation vectors - are directly provided to the pipeline + name: str + name of the 3D microscopy image - mip_msk: bool - apply tissue reconstruction mask (binarized MIP) + fmt: str + format of the 3D microscopy image - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher - fb_ch: int - myelinated fibers channel + is_vec: bool + vector field flag - 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 @@ -268,25 +234,26 @@ def get_image_info(cli_args): 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 - - # 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, 'is_vec': is_vec} + + return in_img, msk_mip def get_frangi_config(cli_args, in_img): @@ -367,8 +334,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 +353,27 @@ 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']) - # Frangi filter response threshold - frangi_cfg['fb_thr'] = cli_args.fb_thr + # spatial scales + scales_px = np.array(cli_args.scales) / np.max(out_px_sz) - # fiber orientation colormap - frangi_cfg['hsv_cmap'] = cli_args.hsv + # 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] - # 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': cli_args.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): @@ -434,13 +394,13 @@ def get_resolution(cli_args): 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 +409,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 +481,77 @@ 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) + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - img_name: str - name of the 3D microscopy image + 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 + + 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 + + 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] 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) - - # create saving directory - save_dirs = create_save_dirs(img_path, img_name, cli_args, is_vec=is_vec) + # 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 + # import fiber orientation vector data or raw 3D microscopy image tic = perf_counter() - if is_vec: - in_img = load_orient(img_path, img_name, img_fmt, tmp_dir=save_dirs['tmp']) - - # import raw 3D microscopy image + if in_img['is_vec']: + load_orient(in_img, save_dirs['tmp']) 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) - - # add image pixel size to input data dictionary - in_img['px_sz'] = px_sz - - # add image name to input data dictionary - in_img['name'] = img_name - - # add vector field flag - in_img['is_vec'] = is_vec + load_raw(in_img, save_dirs['tmp'], msk_mip=msk_mip) # get microscopy image size information - in_img = get_image_size(in_img) + get_image_size(in_img) - # print import time + # print input data information print_import_time(tic) - - # print volume image shape and resolution - if not is_vec: + if not in_img['is_vec']: print_image_info(in_img) else: print_flsh() @@ -545,121 +559,85 @@ 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_orient(in_img, tmp_dir): """ Load array of 3D fiber orientations. Parameters ---------- - img_path: str - path to the 3D microscopy image + in_img: dict + input image dictionary + fb_ch: int + neuronal fibers channel - img_name: str - name of the 3D microscopy image + bc_ch: int + brain cell soma channel - img_fmt: str - format of the 3D microscopy image + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - tmp_dir: str - temporary file directory + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - Returns - ------- - in_img: dict - input image dictionary + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - 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 + path: str + path to the 3D microscopy image - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + name: str + name of the 3D microscopy image - ch_ax: int - RGB image channel axis (either 1, 3, or None for grayscale images) + fmt: str + format of the 3D microscopy image + + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher + + is_vec: bool + vector field flag + + tmp_dir: str + path to temporary folder + + Returns + ------- + None """ - # print heading print_flsh(color_text(0, 191, 255, "\nFiber Orientation Data Import\n")) # load fiber orientations + img_fmt = in_img['fmt'].lower() 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) + vec_img = np.load(in_img['path'], mmap_mode='r') + elif img_fmt in ('tif', 'tiff'): + vec_img = tiff.imread(in_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') + vec_img = create_memory_map(vec_img.shape, dtype=vec_img.dtype, name=in_img['name'], 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") + print_flsh(f"Loading {in_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 + # update input image dictionary + in_img.update({'data': vec_img, 'ts_msk': None, 'ch_ax': None}) - return in_img - -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): +def load_raw(in_img, tmp_dir, msk_mip=False): """ Load 3D microscopy image. 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 - - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - 3D FWHM of the PSF [μm] - - is_tiled: bool - True for tiled reconstructions aligned using ZetaStitcher - - tmp_dir: str - temporary file directory - - 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 @@ -672,40 +650,62 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm, 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 + + tmp_dir: str + path to temporary folder + + msk_mip: bool + apply tissue reconstruction mask (binarized MIP) + + Returns + ------- + 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) + if in_img['is_tiled']: + print_flsh(f"Loading {in_img['path']} tiled reconstruction...\n") + img = VirtualFusedVolume(in_img['path']) # load microscopy 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) + print_flsh(f"Loading {in_img['path']} z-stack...\n") + img_fmt = in_img['fmt'].lower() + if img_fmt in ('tif', 'tiff'): + img = tiff.imread(in_img['path']) else: raise ValueError('Unsupported image format!') - # 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') - - # detect channel axis (RGB images) - ch_ax = detect_ch_axis(img) + img = create_memory_map(img.shape, dtype=img.dtype, name=in_img['name'], + dir=tmp_dir, arr=img, mmap_mode='r') # generate tissue reconstruction mask + ch_ax = detect_ch_axis(img) if msk_mip: dims = len(img.shape) - - # grayscale image 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 +719,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}) diff --git a/foa3d/odf.py b/foa3d/odf.py index 89a42fa..70cc5c5 100644 --- a/foa3d/odf.py +++ b/foa3d/odf.py @@ -7,8 +7,8 @@ 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, odf_scale, odf_norm, + odf_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 +22,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,7 +39,8 @@ 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) initialized array of orientation tensor eigenvalues @@ -63,44 +64,43 @@ 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 + 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(odf_scale, fbr_vec.shape[0]) * odf_scale**2 + for z in range(0, fbr_vec.shape[0], odf_scale): + z_max = z + odf_scale + for y in range(0, fbr_vec.shape[1], odf_scale): + y_max = y + odf_scale + for x in range(0, fbr_vec.shape[2], odf_scale): + x_max = x + odf_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 // odf_scale, y // odf_scale, x // odf_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, odf_deg, odf_norm) - vec_tensor_eigen[zv, yv, xv, :] = compute_vec_tensor_eigen(vec_arr) + vec_tnsr_eig[zv, yv, xv, :] = compute_vec_tensor_eigen(vec_arr) - # compute slice-wise dispersion and anisotropy parameters - compute_orientation_dispersion(vec_tensor_eigen, **odi) + # 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 @@ -190,34 +190,29 @@ def generate_odf_background(bg_mrtrix, fbr_vec, vxl_side, iso_fbr=None): 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 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 + # loop over z-slices, rescale and resize them for z in range(vxl_side // 2, bg_img.shape[0], vxl_side): if bg_img.ndim == 3: tmp_slice = normalize_image(bg_img[z], min_val=min_glob, max_val=max_glob) @@ -227,7 +222,7 @@ def generate_odf_background(bg_mrtrix, fbr_vec, vxl_side, iso_fbr=None): 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) + 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): @@ -240,7 +235,7 @@ 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 fiber ODF resolution (super-voxel side [px]) @@ -271,8 +266,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,43 +275,27 @@ 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 + # initialize ODF and orientation dispersion maps 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) + odi_shp = tuple(np.ceil(np.divide(vec_img_shp, odf_scale)).astype(int)) + odf = create_memory_map(odi_shp + (nc,), dtype='float32', name=f'odf{odf_scale}', dir=tmp_dir) + odi_tot = create_memory_map(odi_shp, dtype='float32', name=f'odi_tot{odf_scale}', dir=tmp_dir) + bg_mrtrix = create_memory_map(tuple(np.flip(odi_shp)), dtype='uint8', name=f'bg{odf_scale}', dir=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{odf_scale}', dir=tmp_dir) + odi_sec = create_memory_map(odi_shp, dtype='float32', name=f'odi_sec{odf_scale}', dir=tmp_dir) + odi_anis = create_memory_map(odi_shp, dtype='float32', name=f'odi_anis{odf_scale}', dir=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 and density arrays + vec_tnsr_eig = create_memory_map(odi_shp + (3,), dtype='float32', name=f'tensor{odf_scale}', dir=tmp_dir) + fbr_dnst = create_memory_map(odi_shp, dtype='float32', name=f'fbr_dnst{odf_scale}', dir=tmp_dir) + + # initialize output dispersion dictionary + 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): diff --git a/foa3d/output.py b/foa3d/output.py index c58467b..08b3391 100644 --- a/foa3d/output.py +++ b/foa3d/output.py @@ -12,30 +12,53 @@ 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) + ('frangi': Frangi filter, 'odf': ODF analysis) """ # get current time time_stamp = datetime.now().strftime("%Y%m%d-%H%M%S") @@ -43,10 +66,10 @@ def create_save_dirs(img_path, img_name, cli_args, is_vec=False): # 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}") + 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) @@ -54,7 +77,7 @@ def create_save_dirs(img_path, img_name, cli_args, is_vec=False): save_dirs = dict() # create Frangi filter output subdirectory - if not is_vec: + if not in_img['is_vec']: frangi_dir = path.join(base_out_dir, 'frangi') mkdir(frangi_dir) save_dirs['frangi'] = frangi_dir @@ -113,7 +136,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 @@ -144,7 +167,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 +201,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 @@ -190,8 +213,8 @@ def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None): """ # loop over output image dictionary fields and save to TIFF 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 out_img[img_key] is not None: + 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..67ab7b8 100644 --- a/foa3d/pipeline.py +++ b/foa3d/pipeline.py @@ -1,7 +1,6 @@ -from time import perf_counter - import numpy as np from joblib import Parallel, delayed +from time import perf_counter from foa3d.frangi import (frangi_filter, init_frangi_arrays, mask_background, write_frangi_arrays) @@ -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,22 @@ 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) 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 +59,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 +88,92 @@ 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: + # 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 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) + # get parallel processing configuration + get_resource_config(cli_args, frangi_cfg) + get_slicing_config(in_img, frangi_cfg) - # 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.") + # 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) 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)) + with Parallel(n_jobs=frangi_cfg['batch'], prefer='threads', require='sharedmem', + return_as="generator_unordered") as parallel: + out_gen = parallel(delayed(frangi_analysis)(s, in_img, frangi_cfg, t_start) for s in slc_rng) - # save output arrays - save_frangi_arrays(save_dirs['frangi'], in_img['name'], out_img, px_sz_iso, ram=ram) + # write results to output arrays once available + for out_slc in out_gen: + write_frangi_arrays(out_img, out_slc, z_out=frangi_cfg['z_out']) - # 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, cfg, t_start): """ 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 +190,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 @@ -203,33 +214,6 @@ def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_ item_sz: int image item size [B] - out_img: dict - output image dictionary - - vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector field - - 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] - cfg: dict Frangi filter configuration @@ -254,9 +238,6 @@ 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 @@ -273,78 +254,91 @@ 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: 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] - z_sel: NumPy slice object - selected z-depth range + jobs: int + number of parallel jobs (threads) + used by the Frangi filter stage - in_rng: NumPy slice object - input image range (fibers) + batch: int + slice batch size - bc_rng: NumPy slice object - input image range (neurons) + slc_shp: numpy.ndarray (shape=(3,), dtype=int) + shape of the basic image slices + analyzed using parallel threads [px] - out_rng: NumPy slice object - output range + ovlp: int + overlapping range between image slices along each axis [px] - in_pad: numpy.ndarray (axis order=(Z,Y,X)) - image padding range + tot_slc: int + total number of image slices - print_info: tuple - optional printed information + z_out: NumPy slice object + output z-range - pad_mode: str - image padding mode + t_start: float + start time [s] - bc_thr: str - thresholding method applied to the neuronal bodies channel + Returns + ------- + out_slc: dict + slice output dictionary - ts_thr: float - relative threshold on non-zero tissue pixels + vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - verbose: int - verbosity level (print progress every N=verbose image slices) + clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap - 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']) + fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy - # 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) + 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 + + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + brain cell mask slice - # (optional) neuronal body masking - if in_img['msk_bc']: + rng: NumPy slice object + output range + """ + # 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']) - # get soma image slice - bc_slc, _ = slice_image(in_img['data'], bc_rng[s], in_img['ch_ax'], ch=in_img['bc_ch']) + # process non-background image slice + 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) - # suppress soma contribution - orient_slc, bc_msk_slc = reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng[s], bc_thr) + # (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']) - # soma mask not available - else: - bc_msk_slc = None + out_slc['rng'] = rng['out'] + else: + out_slc = None - # 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) + print_frangi_progress(t_start, cfg['batch'], cfg['tot_slc'], not_bg) - # print progress - if print_info is not None: - print_frangi_progress(print_info, is_valid, verbose=verbose) + return out_slc -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): """ Analyze 3D fiber orientations exploiting a Frangi-filter-based unsupervised enhancement of tubular structures. @@ -378,9 +372,6 @@ 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 @@ -397,76 +388,96 @@ 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 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']) # 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 +487,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 +518,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 +572,60 @@ 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) + ('frangi': Frangi filter, 'odf': ODF analysis) - fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector field + out_img: dict + output image dictionary - iso_fbr: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image + vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - px_sz: numpy.ndarray (axis order=(3,), dtype=float) - pixel size [μm] + clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - img_name: str - name of the input volume image + fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - backend: str - backend module employed by joblib.Parallel + frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability image) - verbose: int - joblib verbosity level + 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 - - # compute spherical harmonics normalization factors (once for all scales) - odf_norm = get_sph_harm_norm_factors(cli_args.odf_deg) + if cli_args.odf_res is not None: + print_odf_info(cli_args.odf_res, cli_args.odf_deg) - # get number of logical cores - num_cpu = get_available_cores() + # 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 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 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) - # 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 +645,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 +666,16 @@ 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] - odf_scale = int(np.ceil(odf_scale_um / px_sz[0])) - # initialize ODF analysis output arrays + odf_scale = int(np.ceil(odf_scale_um / px_sz[0])) 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 down-sampled background for visualization in MRtrix3 generate_odf_background(bg_mrtrix, fbr_vec, vxl_side=odf_scale, iso_fbr=iso_fbr) # 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) - # save memory maps to TIFF or NIfTI files + # save 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..a95f49e 100644 --- a/foa3d/preprocessing.py +++ b/foa3d/preprocessing.py @@ -34,51 +34,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: + # anisotropic pixel size + if crt_2: print_output_res(px_sz_iso) 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_ratio, sigma=None, pad=None, pad_mode='reflect', anti_alias=True, trunc=4, 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 @@ -96,7 +87,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 + pad_mode: str image padding mode adopted for the smoothing Gaussian filter anti_alias: bool @@ -105,7 +96,7 @@ def correct_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', 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 @@ -122,26 +113,26 @@ def correct_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', # no resizing if np.all(rsz_ratio == 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) - # downsampling + # adaptive image downsampling iso_shp = np.ceil(np.multiply(np.asarray(img.shape), rsz_ratio)).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 padding array accordingly + # resize padding range accordingly pad_rsz = np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad)).astype(int) \ if pad is not None else None diff --git a/foa3d/printing.py b/foa3d/printing.py index 5777e41..d9a9005 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,9 +149,6 @@ 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 @@ -168,47 +165,73 @@ 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 + + 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_num: int - total number of analyzed image slices + 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']) + 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" + + 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") - # 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 neuron masking info + # 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']) print_soma_masking(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] - is_valid: bool - image slice validity flag + batch: int + slice batch size + + tot: int + total number of image slices + + 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): @@ -338,7 +360,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, img_item_sz): """ Print information on the slicing of the basic image sub-volumes processed by the Foa3D tool. @@ -350,9 +372,6 @@ 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] @@ -378,8 +397,7 @@ def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): 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") + f"Image slice size [MB]: {np.ceil(max_slc_sz / 1024**2).astype(int)}\n") def print_soma_masking(msk_bc): diff --git a/foa3d/slicing.py b/foa3d/slicing.py index 32f5516..260115e 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 @@ -187,7 +186,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 +199,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 +210,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 +232,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 +244,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,12 +253,11 @@ 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 + return slc_shp -def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1): +def compute_slice_size(max_ram, mem_growth, batch_sz): """ Compute the size of the basic microscopy image slices fed to the Frangi filtering stage. @@ -284,112 +269,160 @@ def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1): 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) + slc_sz = max_ram / (batch_sz * mem_growth) 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): +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 - px_sz_iso: int - isotropic pixel size [μm] + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - mem_growth: float - empirical memory growth factor - of the Frangi filtering stage + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - mem_fudge: float - memory fudge factor + fb_ch: int + neuronal fibers channel - jobs: int - number of parallel jobs (threads) - used by the Frangi filtering stage + bc_ch: int + brain cell soma channel - ram: float - maximum RAM available to the Frangi filtering stage [B] + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - shp_thr: int - minimum slice side [px] + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - Returns - ------- - batch_sz: int - slice batch size + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - in_slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [px] + name: str + name of the 3D microscopy image - in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [μm] + is_vec: bool + vector field flag - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape - ovlp: int - overlapping range between image slices along each axis [px] + 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] - ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int) - resized overlapping range between image slices [px] + Returns + ------- + None """ # maximum RAM not provided: use all + ram = frangi_cfg['ram'] if ram is None: ram = psutil.virtual_memory()[1] # number of logical cores + jobs = frangi_cfg['jobs'] 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) + ovlp = compute_overlap(frangi_cfg['smooth_sd'], frangi_cfg['scales_um'], rsz_ratio=frangi_cfg['rsz']) # compute the shape of basic microscopy image slices - in_slc_shp = np.array([-1]) - while np.any(in_slc_shp < shp_thr): + batch_sz = np.min([jobs, num_cpu]).astype(int) + 1 + slc_shp = np.array([-1]) + 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!") 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) + slc_sz = compute_slice_size(ram, mem_growth, batch_sz) + slc_shp = compute_slice_shape(in_img['shape'], in_img['item_sz'], slc_sz=slc_sz, ovlp=ovlp) - return batch_sz, in_slc_shp, in_slc_shp_um, rsz_ratio, ovlp, ovlp_rsz + # 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 crop(img, rng, ovlp=None, flip=()): @@ -432,14 +465,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 +485,164 @@ 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 - in_pad: list - list of slice padding arrays + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - out_rng: list - list of output slice index ranges + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - bc_rng: list - (optional) list of neuronal body slice index ranges + name: str + name of the 3D microscopy image - out_slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the processed image slices [px] + is_vec: bool + vector field flag - tot_slc: int - total number of analyzed image slices + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape - batch_sz: int - adjusted slice batch size + 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 + + 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 + + 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 = list() + 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) - - # show progress - print_flsh(f"Generating image slice ranges:\t{100 * (i+1) / tot_slc:.1f}%", end='\r') - - # total number of slices - tot_slc = len(slc_rng['in_rng']) + bc_rng = None - # adjust slice batch size - if batch_sz > tot_slc: - batch_sz = tot_slc + slc_rng.append({'in': in_rng, 'out': out_rng, 'pad': pad, 'bc': bc_rng}) - return slc_rng, out_slc_shp, tot_slc, batch_sz + return slc_rng 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..9f62e7b 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 @@ -57,7 +58,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', dir=None, arr=None, mmap_mode='w+'): """ Create a memory-map to an array stored in a binary file on disk. @@ -86,12 +87,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 dir is None: + dir = tempfile.mkdtemp() + mmap_path = path.join(dir, name + '.mmap') if path.exists(mmap_path): - unlink(mmap_path) + unlink(mmap_path) if arr is None: _ = open(mmap_path, mode='w+') @@ -106,6 +107,33 @@ def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mod return mmap +def delete_tmp_data(tmp_dir, tmp_data): + """ + Delete temporary folder. + + Parameters + ---------- + tmp_dir: str + path to temporary folder to be removed + + tmp_data: tuple + temporary data dictionaries + + Returns + ------- + None + """ + try: + for data in tmp_data: + for k in list(data.keys()): + del data[k] + + rmtree(tmp_dir) + + except OSError: + pass + + def detect_ch_axis(img): """ Detect image channel axis. @@ -122,15 +150,15 @@ def detect_ch_axis(img): """ 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 + 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 get_available_cores(): @@ -162,7 +190,6 @@ def get_item_size(dtype): 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] @@ -182,33 +209,6 @@ def get_item_size(dtype): return item_sz -def delete_tmp_folder(tmp_dir, tmp_data): - """ - Delete temporary folder. - - Parameters - ---------- - tmp_dir: str - path to temporary folder to be removed - - tmp_data: tuple - temporary data dictionaries - - Returns - ------- - None - """ - try: - for data in tmp_data: - for k in list(data.keys()): - del data[k] - - rmtree(tmp_dir) - - except OSError: - pass - - def divide_nonzero(nd_array1, nd_array2, new_val=1e-10): """ Divide two arrays handling zero denominator values. @@ -256,7 +256,7 @@ def elapsed_time(start_time): mins: int minutes - secs: float + secs: int seconds """ stop_time = perf_counter() @@ -266,7 +266,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 @@ -305,7 +305,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: @@ -329,10 +328,9 @@ def get_config_label(cli_args): cfg_lbl: str Frangi filter configuration label """ - cfg_lbl = '_s' + 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}_' - cfg_lbl = f'a{cli_args.alpha}_b{cli_args.beta}_g{cli_args.gamma}_t{cli_args.fb_thr}{cfg_lbl}' + cfg_lbl += f'_s{s}' return cfg_lbl @@ -343,8 +341,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 +355,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) + angle = angle.astype(dtype) - # convert back to list - if is_list: - norm_angle = list(norm_angle) - - 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 +392,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 @@ -460,10 +442,8 @@ def hsv_orient_cmap(vec_img): 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 + 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)) @@ -471,14 +451,10 @@ def hsv_orient_cmap(vec_img): 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 + # 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]): - - # 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) return rgb_map @@ -568,14 +544,12 @@ def rgb_orient_cmap(vec_img, minimum=0, stretch=1, q=8): 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 + 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 colormap slice by slice + # 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) 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