From 37c516ad6d51ae02452dacb76ef39257d48fa5c2 Mon Sep 17 00:00:00 2001 From: msorelli Date: Fri, 18 Oct 2024 16:43:56 +0200 Subject: [PATCH] Improve function parameters and reorganize modules --- doc/source/modules.rst | 1 + doc/source/spharm.rst | 4 + foa3d/__main__.py | 23 +- foa3d/frangi.py | 502 ++++++++++++++++--- foa3d/input.py | 512 +++++++++++++------- foa3d/odf.py | 648 +++++-------------------- foa3d/output.py | 171 ++++++- foa3d/pipeline.py | 1044 +++++++++++++++------------------------- foa3d/preprocessing.py | 13 +- foa3d/printing.py | 319 +++++++----- foa3d/slicing.py | 105 ++-- foa3d/spharm.py | 444 +++++++++++++++++ foa3d/utils.py | 16 +- 13 files changed, 2136 insertions(+), 1666 deletions(-) create mode 100644 doc/source/spharm.rst create mode 100644 foa3d/spharm.py diff --git a/doc/source/modules.rst b/doc/source/modules.rst index b7f689b..40585b5 100644 --- a/doc/source/modules.rst +++ b/doc/source/modules.rst @@ -12,4 +12,5 @@ Python module index preprocessing printing slicing + spharm utils \ No newline at end of file diff --git a/doc/source/spharm.rst b/doc/source/spharm.rst new file mode 100644 index 0000000..ca54ad0 --- /dev/null +++ b/doc/source/spharm.rst @@ -0,0 +1,4 @@ +foa3d.spharm +------------ +.. automodule:: foa3d.spharm + :members: \ No newline at end of file diff --git a/foa3d/__main__.py b/foa3d/__main__.py index bf6b4de..a318280 100644 --- a/foa3d/__main__.py +++ b/foa3d/__main__.py @@ -1,29 +1,22 @@ from foa3d.input import get_cli_parser, load_microscopy_image -from foa3d.pipeline import (parallel_odf_at_scales, parallel_frangi_on_slices) +from foa3d.pipeline import parallel_odf_over_scales, parallel_frangi_over_slices from foa3d.printing import print_pipeline_heading from foa3d.utils import delete_tmp_folder def foa3d(cli_args): - # load 3D microscopy image or 4D array of fiber orientation vectors - img, ts_msk, ch_ax, is_fovec, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args) + # load 3D grayscale or RGB microscopy image or 4D array of fiber orientation vectors + in_img, save_dirs = load_microscopy_image(cli_args) - # conduct parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices - if not is_fovec: - fbr_vec_img, iso_fbr_img, px_sz \ - = parallel_frangi_on_slices(img, ch_ax, cli_args, save_dir[0], tmp_dir, img_name, ts_msk=ts_msk) + # parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices + out_img = parallel_frangi_over_slices(cli_args, save_dirs, in_img) - # skip Frangi filter stage if orientation vectors were directly provided as input - else: - fbr_vec_img, iso_fbr_img, px_sz = (img, None, None) - - # estimate 3D fiber ODF maps at the spatial scales of interest using concurrent workers - if cli_args.odf_res: - parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir[1], tmp_dir, img_name) + # 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']) # delete temporary folder - delete_tmp_folder(tmp_dir) + delete_tmp_folder(save_dirs['tmp']) def main(): diff --git a/foa3d/frangi.py b/foa3d/frangi.py index 417b4fe..39d9bcc 100644 --- a/foa3d/frangi.py +++ b/foa3d/frangi.py @@ -5,7 +5,8 @@ from joblib import Parallel, delayed from scipy import ndimage as ndi -from foa3d.utils import divide_nonzero +from foa3d.utils import (create_background_mask, create_memory_map, + divide_nonzero, hsv_orient_cmap, rgb_orient_cmap) def analyze_hessian_eigen(img, sigma, trunc=4): @@ -33,7 +34,6 @@ def analyze_hessian_eigen(img, sigma, trunc=4): dom_eigenvec: 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) @@ -62,7 +62,6 @@ def compute_dominant_eigen(hessian): dom_eigenvec: 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) @@ -74,6 +73,30 @@ def compute_dominant_eigen(hessian): return sorted_eigenval, dom_eigenvec +def compute_fractional_anisotropy(eigenval): + """ + Compute structure tensor fractional anisotropy + as in Schilling et al. (2018). + + Parameters + ---------- + eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + structure tensor eigenvalues (at the best local spatial scale) + + Returns + ------- + fa: numpy.ndarray (shape=(3,), dtype=float) + fractional anisotropy + """ + fa = np.sqrt(0.5 * divide_nonzero( + np.square((eigenval[..., 0] - eigenval[..., 1])) + + np.square((eigenval[..., 0] - eigenval[..., 2])) + + np.square((eigenval[..., 1] - eigenval[..., 2])), + np.sum(eigenval ** 2, axis=-1))) + + return fa + + def compute_frangi_features(eigen1, eigen2, eigen3, gamma): """ Compute the basic image features employed by the Frangi filter. @@ -139,7 +162,6 @@ def compute_scaled_hessian(img, sigma=1, trunc=4): hessian: numpy.ndarray (axis order=(Z,Y,X,C,C), dtype=float) Hessian matrix of image second derivatives """ - # scale selection scaled_img = ndi.gaussian_filter(img, sigma=sigma, output=np.float32, truncate=trunc) @@ -168,9 +190,63 @@ def compute_scaled_hessian(img, sigma=1, trunc=4): return hessian +def compute_plate_like_score(ra, alpha): + return 1 - np.exp(np.divide(np.negative(np.square(ra)), 2 * np.square(alpha))) + + +def compute_blob_like_score(rb, beta): + return np.exp(np.divide(np.negative(np.square(rb)), 2 * np.square(beta))) + + +def compute_background_score(s, gamma): + return 1 - np.exp(np.divide(np.negative(np.square(s)), 2 * np.square(gamma))) + + +def compute_structureness(eigen1, eigen2, eigen3): + return np.sqrt(np.square(eigen1) + np.square(eigen2) + np.square(eigen3)) + + +def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma): + """ + Estimate Frangi's vesselness probability. + + Parameters + ---------- + eigen1: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + lowest Hessian eigenvalue (i.e., the dominant eigenvalue) + + eigen2: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + middle Hessian eigenvalue + + eigen3: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + highest Hessian eigenvalue + + alpha: float + plate-like score sensitivity + + beta: float + blob-like score sensitivity + + gamma: float + background score sensitivity + + Returns + ------- + vesselness: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + Frangi's vesselness likelihood image + """ + ra, rb, s, gamma = compute_frangi_features(eigen1, eigen2, eigen3, gamma) + plate = compute_plate_like_score(ra, alpha) + blob = compute_blob_like_score(rb, beta) + background = compute_background_score(s, gamma) + vesselness = plate * blob * background + + return vesselness + + def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, dark=False): """ - Estimate fiber orientation vectors at the input spatial scale. + Compute fiber orientation vectors at the input spatial scale of interest Parameters ---------- @@ -204,7 +280,6 @@ def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, d eigenval: 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) @@ -213,7 +288,7 @@ def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, d vesselness = compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha=alpha, beta=beta, gamma=gamma) # reject vesselness background - frangi_img = reject_fiber_background(vesselness, eigen2, eigen3, dark) + frangi_img = reject_vesselness_background(vesselness, eigen2, eigen3, dark) # stack eigenvalues list eigenval = np.stack(eigenval, axis=-1) @@ -221,20 +296,17 @@ def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, d return frangi_img, eigenvec, eigenval -def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma): +def compute_parall_scaled_orientation(scales_px, img, alpha=0.001, beta=1, gamma=None, dark=False): """ - Estimate Frangi's vesselness probability. + Compute fiber orientation vectors over the spatial scales of interest using concurrent workers. Parameters ---------- - eigen1: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - lowest Hessian eigenvalue (i.e., the dominant eigenvalue) - - eigen2: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - middle Hessian eigenvalue + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] - eigen3: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - highest Hessian eigenvalue + img: numpy.ndarray (axis order=(Z,Y,X)) + 3D microscopy image alpha: float plate-like score sensitivity @@ -245,22 +317,191 @@ def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma): gamma: float background score sensitivity + 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) + 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) + 3D orientation map at the input spatial scale + + 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)) - ra, rb, s, gamma = compute_frangi_features(eigen1, eigen2, eigen3, gamma) - plate = compute_plate_like_score(ra, alpha) - blob = compute_blob_like_score(rb, beta) - background = compute_background_score(s, gamma) - vesselness = plate * blob * background + # 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) - return vesselness + # 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) + # 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) -def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True): + return frangi, fbr_vec, eigenval + + +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. + + Parameters + ---------- + 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] + + z_rng: int + output z-range in [px] + + 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 + + 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 + + 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 + + fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + initialized fractional anisotropy image + + frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + initialized Frangi-enhanced image + + iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + initialized fiber image (isotropic resolution) + + fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + initialized fiber mask image + + bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + initialized soma mask image + + z_sel: NumPy slice object + selected z-depth range + """ + # 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) + tot_shp = tuple(np.ceil(rsz_ratio * img_shp).astype(int)) + + # fiber channel arrays + iso_fbr_img = create_memory_map(tot_shp, dtype='uint8', name='iso_fiber', tmp_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) + + # soma channel array + if msk_bc: + bc_msk = create_memory_map(tot_shp, dtype='uint8', name='bc_msk', tmp_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) + + # 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 + + return out_img, z_sel + + +def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=False, hsv=False): """ Apply 3D Frangi filter to 3D microscopy image. @@ -285,52 +526,119 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True if True, enhance black 3D tubular structures (i.e., negative contrast polarity) + hsv: bool + Returns ------- - frangi_img: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - Frangi's vesselness likelihood image + orient_slc: dict + slice orientation dictionary - fbr_vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - 3D fiber orientation map + vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - Hessian eigenvalues sorted by absolute value (ascending order) - """ + clr_slc: 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 + + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + Frangi's vesselness likelihood image + """ # single-scale vesselness analysis - n_scales = len(scales_px) - if n_scales == 1: - frangi_img, fbr_vec, eigenval \ - = compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark) + 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: - 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)) + 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 + fbr_clr = hsv_orient_cmap(fbr_vec) if hsv else rgb_orient_cmap(fbr_vec) - # unpack and stack results - enhanced_img_tpl, eigenvectors_tpl, eigenvalues_tpl = zip(*par_res) - eigenval = np.stack(eigenvalues_tpl, axis=0) - eigenvec = np.stack(eigenvectors_tpl, axis=0) - frangi_img = np.stack(enhanced_img_tpl, axis=0) + # fill orientation dictionary + orient_slc = dict() + orient_slc['vec_slc'] = fbr_vec + orient_slc['clr_slc'] = fbr_clr + orient_slc['fa_slc'] = fa - # get max scale-wise vesselness - best_idx = np.argmax(frangi_img, axis=0) - best_idx = np.expand_dims(best_idx, axis=0) - frangi_img = np.take_along_axis(frangi_img, best_idx, axis=0).squeeze(axis=0) + return orient_slc, frangi - # 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) - return frangi_img, fbr_vec, eigenval +def mask_background(img, orient, method='yen', invert=False, ts_msk=None): + """ + Mask fiber orientation arrays. + + Parameters + ---------- + img: numpy.ndarray (axis order=(Z,Y,X)) + fiber (or neuron) fluorescence 3D image + + orient: dict + slice orientation dictionary + + vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + fiber orientation vector slice + + clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap slice + + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + fractional anisotropy slice + + method: str + thresholding method (refer to skimage.filters) + + invert: bool + mask inversion flag + + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + + 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) -def reject_fiber_background(vesselness, eigen2, eigen3, dark): + bg: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + background mask + """ + # generate background mask + bg = create_background_mask(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)) + + # 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 + + return orient, bg + + +def reject_vesselness_background(vesselness, eigen2, eigen3, dark): """ Reject the fiber background, exploiting the sign of the "secondary" eigenvalues λ2 and λ3. @@ -355,7 +663,6 @@ def reject_fiber_background(vesselness, eigen2, eigen3, dark): 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)) vesselness[bg_msk] = 0 @@ -387,7 +694,6 @@ def sort_eigen(eigenval, eigenvec, axis=-1): sorted_eigenvec: numpy.ndarray (axis order=(Z,Y,X,C,C), dtype=float) sorted eigenvector array """ - # sort the eigenvalue array by absolute value (ascending order) sorted_val_index = np.abs(eigenval).argsort(axis) sorted_eigenval = np.take_along_axis(eigenval, sorted_val_index, axis) @@ -401,17 +707,83 @@ def sort_eigen(eigenval, eigenvec, axis=-1): return sorted_eigenval, sorted_eigenvec -def compute_plate_like_score(ra, alpha): - return 1 - np.exp(np.divide(np.negative(np.square(ra)), 2 * np.square(alpha))) +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): + """ + Fill the memory-mapped output arrays of the Frangi filter stage. + Parameters + ---------- + out_rng: tuple + 3D slice output index range -def compute_blob_like_score(rb, beta): - return np.exp(np.divide(np.negative(np.square(rb)), 2 * np.square(beta))) + z_sel: NumPy slice object + selected z-depth range + iso_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice -def compute_background_score(s, gamma): - return 1 - np.exp(np.divide(np.negative(np.square(s)), 2 * np.square(gamma))) + frangi_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + Frangi-enhanced image slice + fbr_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask image slice -def compute_structureness(eigen1, eigen2, eigen3): - return np.sqrt(np.square(eigen1) + np.square(eigen2) + np.square(eigen3)) + bc_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + soma mask image slice + + vec_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + fiber orientation vector image slice + + clr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + orientation colormap image slice + + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + fractional anisotropy image slice + + 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 + + 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) diff --git a/foa3d/input.py b/foa3d/input.py index afeceea..a3f61e9 100644 --- a/foa3d/input.py +++ b/foa3d/input.py @@ -13,10 +13,10 @@ from foa3d.output import create_save_dirs from foa3d.preprocessing import config_anisotropy_correction -from foa3d.printing import (color_text, print_flushed, print_image_shape, - print_import_time, print_native_res) -from foa3d.utils import (create_background_mask, create_memory_map, detect_ch_axis, - get_item_bytes, get_config_label) +from foa3d.printing import (color_text, print_flsh, print_image_info, + print_import_time) +from foa3d.utils import (create_background_mask, create_memory_map, + detect_ch_axis, get_item_bytes, get_config_label) class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawTextHelpFormatter): @@ -73,21 +73,19 @@ 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('-m', '--mmap', action='store_true', default=False, - help='create a memory-mapped array of input microscopy or fiber orientation images') cli_parser.add_argument('--px-size-xy', type=float, default=0.878, 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=0.692, help='PSF FWHM along horizontal x-axis [μm]') - cli_parser.add_argument('--psf-fwhm-y', type=float, default=0.692, help='PSF FWHM along vertical x-axis [μm]') - cli_parser.add_argument('--psf-fwhm-z', type=float, default=2.612, help='PSF FWHM along depth z-axis [μ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]') + cli_parser.add_argument('--psf-fwhm-z', type=float, default=1.0, help='PSF FWHM along depth z-axis [μm]') cli_parser.add_argument('--fb-ch', type=int, default=1, help='neuronal fibers channel') - cli_parser.add_argument('--bc-ch', type=int, default=0, help='neuronal bodies channel') + cli_parser.add_argument('--bc-ch', type=int, default=0, help='brain cell soma channel') cli_parser.add_argument('--fb-thr', default='li', type=str, help='Frangi filter probability response threshold (t ∈ [0, 1] or skimage.filters method)') cli_parser.add_argument('--z-min', type=float, default=0, help='forced minimum output z-depth [μm]') cli_parser.add_argument('--z-max', type=float, default=None, help='forced maximum output z-depth [μm]') cli_parser.add_argument('--hsv', action='store_true', default=False, - help='generate HSV colormap for 3D fiber orientations') + help='generate HSV colormap of 3D fiber orientations') cli_parser.add_argument('--odf-res', nargs='+', type=float, help='side of the fiber ODF super-voxels: ' 'do not generate ODFs if None [μm]') cli_parser.add_argument('--odf-deg', type=int, default=6, @@ -110,62 +108,107 @@ def get_cli_parser(): return cli_args -def get_image_info(img, px_sz, msk_bc, fb_ch, ch_ax): +def get_image_size(in_img): """ Get information on the input 3D microscopy image. Parameters ---------- - img: numpy.ndarray - 3D microscopy image + 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 - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - fb_ch: int - myelinated fibers channel + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) + + fb_ch: int + neuronal fibers channel - ch_ax: int - RGB image channel axis (either 1 or 3, or None for grayscale images) + 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 ------- - img_shp: numpy.ndarray (shape=(3,), dtype=int) - volume image shape [px] + in_img: dict + input image dictionary (extended) - img_shp_um: numpy.ndarray (shape=(3,), dtype=float) - volume image shape [μ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 - img_item_sz: int - image item size [B] + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - fb_ch: int - myelinated fibers channel + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel - """ + 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 + + 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] + + """ # adapt channel axis - img_shp = np.asarray(img.shape) - if ch_ax is None: - fb_ch = None - msk_bc = False + in_img['shape'] = np.asarray(in_img['data'].shape) + if in_img['ch_ax'] is None: + in_img['fb_ch'] = None + in_img['msk_bc'] = False else: - img_shp = np.delete(img_shp, ch_ax) + in_img['shape'] = np.delete(in_img['shape'], in_img['ch_ax']) - img_shp_um = np.multiply(img_shp, px_sz) - img_item_sz = get_item_bytes(img) + in_img['shape_um'] = np.multiply(in_img['shape'], in_img['px_sz']) + in_img['item_sz'] = get_item_bytes(in_img['data']) - return img_shp, img_shp_um, img_item_sz, fb_ch, msk_bc + return in_img -def get_file_info(cli_args): +def get_image_info(cli_args): """ Get microscopy image file path and format. @@ -185,27 +228,33 @@ def get_file_info(cli_args): img_fmt: str format of the 3D microscopy image + 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] + is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher - is_fovec: bool + is_vec: bool True when pre-estimated fiber orientation vectors are directly provided to the pipeline - is_mmap: bool - create a memory-mapped array of the 3D microscopy image, - increasing the parallel processing performance - (the image will be preliminarily loaded to RAM) - mip_msk: 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 + """ # get microscopy image path and name - is_mmap = cli_args.mmap img_path = cli_args.image_path img_name = path.basename(img_path) split_name = img_name.split('.') @@ -217,20 +266,30 @@ def get_file_info(cli_args): img_fmt = split_name[-1] img_name = img_name.replace(f'.{img_fmt}', '') is_tiled = True if img_fmt == 'yml' else False - is_fovec = cli_args.vec + 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) msk_mip = cli_args.tissue_msk - fb_ch = cli_args.fb_ch + + # apply brain cell soma mask + msk_bc = cli_args.cell_msk # get Frangi filter configuration cfg_lbl = get_config_label(cli_args) img_name = f'{img_name}_{cfg_lbl}' - return img_path, img_name, img_fmt, is_tiled, is_fovec, is_mmap, msk_mip, fb_ch + # get microscopy image resolution + 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 -def get_frangi_config(cli_args): + +def get_frangi_config(cli_args, in_img): """ Get Frangi filter configuration. @@ -239,83 +298,122 @@ def get_frangi_config(cli_args): cli_args: see ArgumentParser.parse_args populated namespace of command line arguments + 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 + + 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 + + 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] + Returns ------- - alpha: float - plate-like score sensitivity + frangi_cfg: dict + Frangi filter configuration - beta: float - blob-like score sensitivity + alpha: float + plate-like score sensitivity - gamma: float - background score sensitivity + beta: float + blob-like score sensitivity - scales_px: numpy.ndarray (dtype=float) - Frangi filter scales [px] + gamma: float + background score sensitivity - scales_um: numpy.ndarray (dtype=float) - Frangi filter scales [μm] + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] - smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the smoothing Gaussian filter [px] + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] - px_sz_iso: int - isotropic pixel size [μm] + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - z_rng: int - output z-range in [px] + z_rng: int + output z-range in [px] - bc_ch: int - neuronal bodies channel + bc_ch: int + neuronal bodies channel - fb_ch: int - myelinated fibers channel + fb_ch: int + myelinated fibers channel - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - hsv_vec_cmap: bool + hsv_cmap: bool + generate HSV colormap of 3D fiber orientations - exp_all: bool - export all images - """ + exp_all: bool + export all images - # microscopy image pixel and PSF size - px_sz, psf_fwhm = get_resolution(cli_args) + + """ + # initialize Frangi filter configuration dictionary + frangi_cfg = dict() # preprocessing configuration (in-plane smoothing) - smooth_sigma, px_sz_iso = config_anisotropy_correction(px_sz, psf_fwhm) + frangi_cfg['smooth_sd'], px_sz_iso = config_anisotropy_correction(in_img['px_sz'], in_img['psf_fwhm']) # Frangi filter parameters - alpha, beta, gamma = (cli_args.alpha, cli_args.beta, cli_args.gamma) - scales_um = np.array(cli_args.scales) - scales_px = scales_um / px_sz_iso[0] - - # image channels - fb_ch = cli_args.fb_ch - bc_ch = cli_args.bc_ch - msk_bc = cli_args.cell_msk + 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] # Frangi filter response threshold - fb_thr = cli_args.fb_thr + frangi_cfg['fb_thr'] = cli_args.fb_thr # fiber orientation colormap - hsv_vec_cmap = cli_args.hsv + frangi_cfg['hsv_cmap'] = cli_args.hsv # forced output z-range - z_min = int(np.floor(cli_args.z_min / px_sz[0])) - z_max = int(np.ceil(cli_args.z_max / px_sz[0])) if cli_args.z_max is not None else cli_args.z_max - z_rng = (z_min, z_max) + 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) # export all flag - exp_all = cli_args.exp_all + frangi_cfg['exp_all'] = cli_args.exp_all - return alpha, beta, gamma, scales_px, scales_um, smooth_sigma, \ - px_sz, px_sz_iso, z_rng, bc_ch, fb_ch, fb_thr, msk_bc, hsv_vec_cmap, exp_all + return frangi_cfg, px_sz_iso def get_resolution(cli_args): @@ -335,14 +433,10 @@ def get_resolution(cli_args): psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) 3D PSF FWHM [μm] """ - # create pixel and psf size arrays px_sz = np.array([cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy]) psf_fwhm = np.array([cli_args.psf_fwhm_z, cli_args.psf_fwhm_y, cli_args.psf_fwhm_x]) - # print resolution info - print_native_res(px_sz, psf_fwhm) - return px_sz, psf_fwhm @@ -364,7 +458,6 @@ def get_resource_config(cli_args): number of parallel jobs (threads) used by the Frangi filter stage """ - # resource parameters (convert maximum RAM to bytes) jobs = cli_args.jobs max_ram = cli_args.ram @@ -376,8 +469,8 @@ def get_resource_config(cli_args): def load_microscopy_image(cli_args): """ - Load 3D microscopy image from TIFF, NumPy or ZetaStitcher .yml file. - Alternatively, the processing pipeline accepts as input NumPy or HDF5 + Load 3D microscopy image from TIFF, or ZetaStitcher .yml file. + Alternatively, the processing pipeline accepts as input TIFF or NumPy files of fiber orientation vector data: in this case, the Frangi filter stage will be skipped. @@ -388,59 +481,71 @@ def load_microscopy_image(cli_args): Returns ------- - img: numpy.ndarray or NumPy memory-map object - 3D microscopy image or array of fiber orientation vectors + in_img: dict + input image dictionary - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + 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 - ch_ax: int - RGB image channel axis (either 1 or 3, or None for grayscale images) + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - is_fovec: bool - True when pre-estimated fiber orientation vectors - are directly provided to the pipeline + img_name: str + name of the 3D microscopy image - save_dir: list (dtype=str) - saving subdirectory string paths + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - tmp_dir: str - temporary file directory + is_vec: bool + vector field flag - img_name: str - microscopy image filename + save_dirs: dict + saving directories + ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) """ - # retrieve input file information - img_path, img_name, img_fmt, is_tiled, is_fovec, is_mmap, msk_mip, fb_ch = get_file_info(cli_args) + 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_dir, tmp_dir = create_save_dirs(img_path, img_name, cli_args, is_fovec=is_fovec) + save_dirs = create_save_dirs(img_path, img_name, cli_args, is_vec=is_vec) # import fiber orientation vector data tic = perf_counter() - if is_fovec: - img = load_orient(img_path, img_name, img_fmt, is_mmap=is_mmap, tmp_dir=tmp_dir) - ts_msk = None + if is_vec: + in_img = load_orient(img_path, img_name, img_fmt, tmp_dir=save_dirs['tmp']) # import raw 3D microscopy image else: - img, ts_msk, ch_ax = load_raw(img_path, img_name, img_fmt, - is_tiled=is_tiled, is_mmap=is_mmap, tmp_dir=tmp_dir, msk_mip=msk_mip, fb_ch=fb_ch) + 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 + + # get microscopy image size information + in_img = get_image_size(in_img) # print import time print_import_time(tic) - # print volume image shape - if not is_fovec: - print_image_shape(cli_args, img, ch_ax) + # print volume image shape and resolution + if not is_vec: + print_image_info(in_img) else: - print_flushed() + print_flsh() - return img, ts_msk, ch_ax, is_fovec, save_dir, tmp_dir, img_name + return in_img, save_dirs -def load_orient(img_path, img_name, img_fmt, is_mmap=False, tmp_dir=None): +def load_orient(img_path, img_name, img_fmt, tmp_dir=None): """ Load array of 3D fiber orientations. @@ -455,46 +560,56 @@ def load_orient(img_path, img_name, img_fmt, is_mmap=False, tmp_dir=None): img_fmt: str format of the 3D microscopy image - is_mmap: bool - create a memory-mapped array of the 3D microscopy image, - increasing the parallel processing performance - (the image will be preliminarily loaded to RAM) - tmp_dir: str temporary file directory Returns ------- - img: numpy.ndarray - 3D fiber orientation vectors - """ + 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) + """ # print heading - print_flushed(color_text(0, 191, 255, "\nFiber Orientation Data Import\n")) + print_flsh(color_text(0, 191, 255, "\nFiber Orientation Data Import\n")) # load fiber orientations if img_fmt == 'npy': - img = np.load(img_path, mmap_mode='r') - if is_mmap: - img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r') + vec_img = np.load(img_path, mmap_mode='r') elif img_fmt == 'tif' or img_fmt == 'tiff': - img = tiff.imread(img_path) - ch_ax = detect_ch_axis(img) - if ch_ax == 1: - img = np.moveaxis(img, 1, -1) - if is_mmap: - img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r') + vec_img = tiff.imread(img_path) + ch_ax = detect_ch_axis(vec_img) + if ch_ax != 3: + vec_img = np.moveaxis(vec_img, ch_ax, -1) + + # memory-map the input TIFF image + vec_img = create_memory_map(vec_img.shape, dtype=vec_img.dtype, name=img_name, + tmp_dir=tmp_dir, arr=vec_img, mmap_mode='r') # check array shape - if img.ndim != 4: + if vec_img.ndim != 4: raise ValueError('Invalid 3D fiber orientation dataset (ndim != 4)!') else: - print_flushed(f"Loading {img_path} orientation vector field...\n") + print_flsh(f"Loading {img_path} orientation vector field...\n") + + # populate input image dictionary + in_img = dict() + in_img['data'] = vec_img + in_img['ts_msk'] = None + in_img['ch_ax'] = None - return img + return in_img -def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None, msk_mip=False, fb_ch=1): +def load_raw(img_path, img_name, img_fmt, psf_fwhm, + is_tiled=False, tmp_dir=None, msk_mip=False, msk_bc=False, fb_ch=1, bc_ch=0): """ Load 3D microscopy image. @@ -509,57 +624,74 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir 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 - is_mmap: bool - create a memory-mapped array of the 3D microscopy image, - increasing the parallel processing performance - (the image will be preliminarily loaded to RAM) - tmp_dir: str temporary file directory - mip_msk: bool + 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 ------- - img: numpy.ndarray or NumPy memory-map object - 3D microscopy image + in_img: dict + input image dictionary - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + 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 - ch_ax: int - RGB image channel axis (either 1 or 3) - """ + 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] + """ # print heading - print_flushed(color_text(0, 191, 255, "\nMicroscopy Image Import\n")) + print_flsh(color_text(0, 191, 255, "\nMicroscopy Image Import\n")) # load microscopy tiled reconstruction (aligned using ZetaStitcher) if is_tiled: - print_flushed(f"Loading {img_path} tiled reconstruction...\n") + print_flsh(f"Loading {img_path} tiled reconstruction...\n") img = VirtualFusedVolume(img_path) # load microscopy z-stack else: - print_flushed(f"Loading {img_path} z-stack...\n") + print_flsh(f"Loading {img_path} z-stack...\n") img_fmt = img_fmt.lower() - if img_fmt == 'npy': - img = np.load(img_path) - elif img_fmt == 'tif' or img_fmt == 'tiff': + if img_fmt == 'tif' or img_fmt == 'tiff': img = tiff.imread(img_path) else: raise ValueError('Unsupported image format!') # create image memory map - if is_mmap: - img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r') + 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) @@ -574,6 +706,8 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir # RGB image elif dims == 4: img_fbr = img[:, fb_ch, :, :] if ch_ax == 1 else img[..., fb_ch] + else: + raise ValueError('Invalid image (ndim != 3 and ndim != 4)!') # compute MIP (naive for loop to minimize the required RAM) ts_mip = np.zeros(img_fbr.shape[1:], dtype=img_fbr.dtype) @@ -585,4 +719,14 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir else: ts_msk = None - return img, ts_msk, ch_ax + # 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 diff --git a/foa3d/odf.py b/foa3d/odf.py index 679a4dc..fd0d27a 100644 --- a/foa3d/odf.py +++ b/foa3d/odf.py @@ -3,65 +3,11 @@ from numba import njit from skimage.transform import resize -from foa3d.utils import normalize_image, transform_axes +from foa3d.spharm import fiber_vectors_to_sph_harm, get_sph_harm_ncoeff +from foa3d.utils import create_memory_map, normalize_image, transform_axes -def compute_disarray(fbr_vec): - """ - Compute local angular disarray, i.e. the misalignment degree of nearby orientation vectors - with respect to the mean direction (Giardini et al., Front. Physiol. 2021). - - Parameters - ---------- - fbr_vec: numpy.ndarray (shape=(N,3), dtype=float) - array of fiber orientation vectors - (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) - - Returns - ------- - disarray: float32 - local angular disarray - """ - fbr_vec = np.delete(fbr_vec, np.all(fbr_vec == 0, axis=-1), axis=0) - disarray = (1 - np.linalg.norm(np.mean(fbr_vec, axis=0))).astype(np.float32) - - return disarray - - -@njit(cache=True) -def compute_fiber_angles(fbr_vec, norm): - """ - Estimate the spherical coordinates (φ azimuth and θ polar angles) - of the fiber orientation vectors returned by the Frangi filtering stage - (all-zero background vectors are excluded). - - Parameters - ---------- - fbr_vec: numpy.ndarray (shape=(N,3), dtype=float) - array of fiber orientation vectors - (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) - - norm: numpy.ndarray (shape=(N,), dtype=float) - 2-norm of fiber orientation vectors - - Returns - ------- - phi: numpy.ndarray (shape=(N,), dtype=float) - fiber azimuth angle [rad] - - theta: numpy.ndarray (shape=(N,), dtype=float) - fiber polar angle [rad] - """ - - fbr_vec = fbr_vec[norm > 0, :] - phi = np.arctan2(fbr_vec[:, 1], fbr_vec[:, 2]) - theta = np.arccos(fbr_vec[:, 0] / norm[norm > 0]) - - return phi, theta - - -def compute_odf_map(fbr_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, vec_tensor_eigen, vxl_side, odf_norm, - odf_deg=6, vxl_thr=0.5, vec_thr=0.000001): +def compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, vxl_side, odf_norm, odf_deg=6, vxl_thr=0.5, vec_thr=0.000001): """ Compute the spherical harmonics coefficients iterating over super-voxels of fiber orientation vectors. @@ -72,22 +18,22 @@ def compute_odf_map(fbr_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, fiber orientation vectors 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_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - primary orientation dispersion index + odi: dict + orientation dispersion dictionary - odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - secondary orientation dispersion index + odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + primary orientation dispersion index - odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - total orientation dispersion index + odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + secondary orientation dispersion index - odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - orientation dispersion index anisotropy + odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + total orientation dispersion index - disarray: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - local angular disarray + odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + orientation dispersion index anisotropy vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) initialized array of orientation tensor eigenvalues @@ -112,7 +58,6 @@ def compute_odf_map(fbr_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, 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): @@ -140,15 +85,11 @@ def compute_odf_map(fbr_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, 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) - # (optional) compute angular disarray - if disarray is not None: - disarray[zv, yv, xv] = compute_disarray(vec_arr) - # compute slice-wise dispersion and anisotropy parameters - compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis) + compute_orientation_dispersion(vec_tensor_eigen, **odi) - # keep a black background - mask_orientation_dispersion([odi_pri, odi_sec, odi_tot, odi_anis], vec_tensor_eigen) + # set dispersion background to 0 + mask_orientation_dispersion(vec_tensor_eigen, odi) # transform axes odf = transform_axes(odf, swapped=(0, 2), flipped=(1, 2)) @@ -167,20 +108,6 @@ def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, orientation tensor eigenvalues computed from an ODF super-voxel - odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - primary orientation dispersion index - - odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - secondary orientation dispersion index - - odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - total orientation dispersion index - - odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - orientation dispersion index anisotropy - - Returns - ------- odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) primary orientation dispersion index @@ -192,77 +119,32 @@ def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) orientation dispersion index anisotropy - """ - k1 = np.abs(vec_tensor_eigen[..., 2]) - k2 = np.abs(vec_tensor_eigen[..., 1]) - k3 = np.abs(vec_tensor_eigen[..., 0]) + Returns + ------- + None + """ + # get difference between secondary tensor eigenvalues diff = np.abs(vec_tensor_eigen[..., 1] - vec_tensor_eigen[..., 0]) + # get absolute tensor eigenvalues + avte = np.abs(vec_tensor_eigen) + # primary dispersion (0.3183098861837907 = 1/π) if odi_pri is not None: - odi_pri[:] = (1 - 0.3183098861837907 * np.arctan2(k1, k2)).astype(np.float32) + odi_pri[:] = (1 - 0.3183098861837907 * np.arctan2(avte[..., 2], avte[..., 1])).astype(np.float32) # secondary dispersion if odi_sec is not None: - odi_sec[:] = (1 - 0.3183098861837907 * np.arctan2(k1, k3)).astype(np.float32) + odi_sec[:] = (1 - 0.3183098861837907 * np.arctan2(avte[..., 2], avte[..., 0])).astype(np.float32) # dispersion anisotropy if odi_anis is not None: - odi_anis[:] = (1 - 0.3183098861837907 * np.arctan2(k1, diff)).astype(np.float32) + odi_anis[:] = (1 - 0.3183098861837907 * np.arctan2(avte[..., 2], diff)).astype(np.float32) # total dispersion - odi_tot[:] = (1 - 0.3183098861837907 * np.arctan2(k1, np.sqrt(np.abs(np.multiply(k2, k3))))).astype(np.float32) - - -@njit(cache=True) -def compute_real_sph_harm(degree, order, phi, sin_theta, cos_theta, norm_factors): - """ - Estimate the coefficients of the real spherical harmonics series expansion - as described by Alimi et al. (Medical Image Analysis, 2020). - - Parameters - ---------- - degree: int - degree index of the spherical harmonics expansion - - order: int - order index of the spherical harmonics expansion - - phi: float - azimuth angle [rad] - - sin_theta: float - polar angle sine - - cos_theta: float - polar angle cosine - - norm_factors: numpy.ndarray (dtype: float) - normalization factors - - Returns - ------- - real_sph_harm: float - real-valued spherical harmonic coefficient - """ - - if degree == 0: - real_sph_harm = norm_factors[0, 0] - elif degree == 2: - real_sph_harm = sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm_factors[1, :]) - elif degree == 4: - real_sph_harm = sph_harm_degree_4(order, phi, sin_theta, cos_theta, norm_factors[2, :]) - elif degree == 6: - real_sph_harm = sph_harm_degree_6(order, phi, sin_theta, cos_theta, norm_factors[3, :]) - elif degree == 8: - real_sph_harm = sph_harm_degree_8(order, phi, sin_theta, cos_theta, norm_factors[4, :]) - elif degree == 10: - real_sph_harm = sph_harm_degree_10(order, phi, sin_theta, cos_theta, norm_factors[5, :]) - else: - raise ValueError("\n Invalid degree of the spherical harmonics series expansion!!!") - - return real_sph_harm + odi_tot[:] = (1 - 0.3183098861837907 * + np.arctan2(avte[..., 2], np.sqrt(np.abs(np.multiply(avte[..., 1], avte[..., 0]))))).astype(np.float32) def compute_vec_tensor_eigen(fbr_vec): @@ -281,7 +163,6 @@ def compute_vec_tensor_eigen(fbr_vec): vec_tensor_eigen: numpy.ndarray (shape=(3,), dtype=float32) orientation tensor eigenvalues in ascending order """ - fbr_vec = np.delete(fbr_vec, np.all(fbr_vec == 0, axis=-1), axis=0) fbr_vec.shape = (-1, 3) t = np.zeros((3, 3)) @@ -293,144 +174,34 @@ def compute_vec_tensor_eigen(fbr_vec): return vec_tensor_eigen -factorial_lut = np.array([ - 1, 1, 2, 6, 24, 120, 720, 5040, 40320, - 362880, 3628800, 39916800, 479001600, - 6227020800, 87178291200, 1307674368000, - 20922789888000, 355687428096000, 6402373705728000, - 121645100408832000, 2432902008176640000], dtype=np.double) - - -@njit(cache=True) -def factorial(n): - """ - Retrieve factorial using pre-computed LUT. - - Parameters - ---------- - n: int - integer number (max: 20) - - Returns - ------- - f: int - factorial - """ - if n > 20: - raise ValueError - - return factorial_lut[n] - - -@njit(cache=True) -def fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff): - """ - Generate the real-valued symmetric spherical harmonics series expansion - from fiber φ azimuth and θ polar angles, - i.e. the spherical coordinates of the fiber orientation vectors. - - Parameters - ---------- - phi: numpy.ndarray (shape=(N,), dtype=float) - fiber azimuth angles [rad] - (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) - - theta: numpy.ndarray (shape=(N,), dtype=float) - fiber polar angle [rad] - (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) - - degrees: int - degrees of the spherical harmonics expansion - - norm_factors: numpy.ndarray (dtype: float) - normalization factors - - ncoeff: int - number of spherical harmonics coefficients - - Returns - ------- - real_sph_harm: numpy.ndarray (shape=(ncoeff,), dtype=float) - array of real-valued spherical harmonics coefficients - building the spherical harmonics series expansion - """ - - real_sph_harm = np.zeros(ncoeff) - i = 0 - for n in range(0, degrees + 1, 2): - for m in range(-n, n + 1, 1): - for p, t in zip(phi, theta): - real_sph_harm[i] += compute_real_sph_harm(n, m, p, np.sin(t), np.cos(t), norm_factors) - i += 1 - - real_sph_harm /= phi.size - - return real_sph_harm - - -def fiber_vectors_to_sph_harm(fbr_vec, degrees, norm_factors): - """ - Generate the real-valued symmetric spherical harmonics series expansion - from the fiber orientation vectors returned by the Frangi filter stage. - - Parameters - ---------- - fbr_vec: numpy.ndarray (shape=(N,3), dtype=float) - array of fiber orientation vectors - (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) - - degrees: int - degrees of the spherical harmonics expansion - - norm_factors: numpy.ndarray (dtype: float) - normalization factors - - Returns - ------- - real_sph_harm: numpy.ndarray (shape=(ncoeff,), dtype=float) - real-valued spherical harmonics coefficients - """ - - fbr_vec.shape = (-1, 3) - ncoeff = get_sph_harm_ncoeff(degrees) - - norm = np.linalg.norm(fbr_vec, axis=-1) - if np.sum(norm) < np.sqrt(fbr_vec.shape[0]): - return np.zeros(ncoeff) - - phi, theta = compute_fiber_angles(fbr_vec, norm) - - real_sph_harm = fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff) - - return real_sph_harm - - -def generate_odf_background(bg_img, bg_mrtrix_mmap, vxl_side): +def generate_odf_background(bg_mrtrix, fbr_vec, vxl_side, iso_fbr=None): """ Generate the downsampled background image required to visualize the 3D ODF map in MRtrix3. Parameters ---------- - bg_img: NumPy memory-map object or HDF5 dataset - (axis order=(Z,Y,X), dtype=uint8; or axis order=(Z,Y,X,C), dtype=float32) - fiber volume image or orientation vector volume - to be used as the MRtrix3 background image - - bg_mrtrix_mmap: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z), dtype=uint8) + bg_mrtrix: NumPy memory-map object (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) + 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) + isotropic fiber image + Returns ------- - bg_mrtrix_mmap: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z), dtype=uint8) - downsampled background for ODF visualization in MRtrix3 + 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_mmap.shape[:-1] + new_shape = bg_mrtrix.shape[:-1] # image normalization: get global minimum and maximum values if bg_img.ndim == 3: @@ -438,295 +209,132 @@ def generate_odf_background(bg_img, bg_mrtrix_mmap, vxl_side): max_glob = np.max(bg_img) # loop over z-slices, and resize them - for z in range(0, bg_img.shape[0], vxl_side): + 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:z + vxl_side, ...], min_val=min_glob, max_val=max_glob) - tmp_slice = np.mean(tmp_slice, axis=0) + tmp_slice = normalize_image(bg_img[z], min_val=min_glob, max_val=max_glob) elif bg_img.ndim == 4: tmp_slice = 255.0 * np.sum(np.abs(bg_img[z, ...]), axis=-1) tmp_slice = np.where(tmp_slice <= 255.0, tmp_slice, 255.0) tmp_slice = transform_axes(tmp_slice, swapped=(0, 1), flipped=(0, 1)) - bg_mrtrix_mmap[..., z // vxl_side] = \ + bg_mrtrix[..., z // vxl_side] = \ resize(tmp_slice, output_shape=new_shape, anti_aliasing=True, preserve_range=True) -@njit(cache=True) -def get_sph_harm_ncoeff(degrees): +def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False): """ - Get the number of coefficients of the real spherical harmonics series - expansion. + Initialize the output datasets of the ODF analysis stage. Parameters ---------- - degrees: int - degrees of the spherical harmonics series expansion - - Returns - ------- - ncoeff: int - number of spherical harmonics coefficients - """ - ncoeff = (2 * (degrees // 2) + 1) * ((degrees // 2) + 1) - - return ncoeff + vec_img_shp: tuple + vector volume shape [px] + tmp_dir: str + temporary file directory -@njit(cache=True) -def get_sph_harm_norm_factors(degrees): - """ - Estimate the normalization factors of the real spherical harmonics series - expansion. + odf_scale: int + fiber ODF resolution (super-voxel side [px]) - Parameters - ---------- - degrees: int + odf_deg: int degrees of the spherical harmonics series expansion + exp_all: bool + export all images + Returns ------- - norm_factors: numpy.ndarray (dtype: float) - 2D array of spherical harmonics normalization factors - """ + odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) + initialized array of ODF spherical harmonics coefficients - norm_factors = np.zeros(shape=(degrees + 1, 2 * degrees + 1)) - for n in range(0, degrees + 1, 2): - for m in range(0, n + 1, 1): - norm_factors[n, m] = norm_factor(n, m) + odi: dict + dispersion image dictionary - norm_factors = norm_factors[::2] + odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + initialized array of primary orientation dispersion parameters - return norm_factors + odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + initialized array of secondary orientation dispersion parameters + odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + initialized array of total orientation dispersion parameters -def mask_orientation_dispersion(odi_lst, vec_tensor_eigen): - """ - Keep a black background where no fiber orientation vector - is available. + odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + initialized array of orientation dispersion anisotropy parameters - Parameters - ---------- - odi_lst: list - list including the ODI maps estimated - in compute_orientation_dispersion() + bg_mrtrix: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8) + initialized background for ODF visualization in MRtrix3 vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - orientation tensor eigenvalues - computed from an ODF super-voxel - - Returns - ------- - None + initialized fiber orientation tensor eigenvalues """ + # ODI maps shape + odi_shp = tuple(np.ceil(np.divide(vec_img_shp, odf_scale)).astype(int)) + + # create downsampled background memory map + bg_shp = tuple(np.flip(odi_shp)) + bg_mrtrix = create_memory_map(bg_shp, dtype='uint8', name=f'bg_tmp{odf_scale}', tmp_dir=tmp_dir) + + # create ODF memory map + nc = get_sph_harm_ncoeff(odf_deg) + odf_shp = odi_shp + (nc,) + odf = create_memory_map(odf_shp, dtype='float32', name=f'odf_tmp{odf_scale}', tmp_dir=tmp_dir) + + # create orientation tensor memory map + vec_tensor_shape = odi_shp + (3,) + vec_tensor_eigen = create_memory_map(vec_tensor_shape, dtype='float32', + name=f'tensor_tmp{odf_scale}', tmp_dir=tmp_dir) + + # create ODI memory maps + odi_tot = create_memory_map(odi_shp, dtype='float32', name=f'odi_tot_tmp{odf_scale}', tmp_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) + else: + odi_pri, odi_sec, odi_anis = (None, None, None) - bg_msk = np.all(vec_tensor_eigen == 0, axis=-1) + # 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 - for odi in odi_lst: - if odi is not None: - odi[bg_msk] = 0 + return odf, odi, bg_mrtrix, vec_tensor_eigen -@njit(cache=True) -def norm_factor(n, m): +def mask_orientation_dispersion(vec_tensor_eigen, odi): """ - Compute the normalization factor of the term of degree n and order m - of the real-valued spherical harmonics series expansion. + Suppress orientation dispersion background. Parameters ---------- - n: int - degree index - - m: int - order index - - Returns - ------- - nf: float - normalization factor - """ - - if m == 0: - nf = np.sqrt((2 * n + 1) / (4 * np.pi)) - else: - nf = (-1)**m * np.sqrt(2) * np.sqrt(((2 * n + 1) / (4 * np.pi) * - (factorial(n - np.abs(m)) / factorial(n + np.abs(m))))) - - return nf - - -@njit(cache=True) -def sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm): - if order == -2: - return norm[2] * 3 * sin_theta**2 * np.sin(2 * phi) - elif order == -1: - return norm[1] * 3 * sin_theta * cos_theta * np.sin(phi) - elif order == 0: - return norm[0] * 0.5 * (3 * cos_theta ** 2 - 1) - elif order == 1: - 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) - - -@njit(cache=True) -def sph_harm_degree_4(order, phi, sin_theta, cos_theta, norm): - if order == -4: - return norm[4] * 105 * sin_theta**4 * np.sin(4 * phi) - elif order == -3: - return norm[3] * 105 * sin_theta**3 * cos_theta * np.sin(3 * phi) - elif order == -2: - return norm[2] * 7.5 * sin_theta**2 * (7 * cos_theta ** 2 - 1) * np.sin(2 * phi) - elif order == -1: - return norm[1] * 2.5 * sin_theta * (7 * cos_theta ** 3 - 3 * cos_theta) * np.sin(phi) - elif order == 0: - return norm[0] * 0.125 * (35 * cos_theta ** 4 - 30 * cos_theta ** 2 + 3) - elif order == 1: - return norm[1] * 2.5 * sin_theta * (7 * cos_theta ** 3 - 3 * cos_theta) * np.cos(phi) - elif order == 2: - return norm[2] * 7.5 * sin_theta**2 * (7 * cos_theta ** 2 - 1) * np.cos(2 * phi) - elif order == 3: - 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) + vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + orientation tensor eigenvalues + computed from an ODF super-voxel + odi: dict + orientation dispersion dictionary -@njit(cache=True) -def sph_harm_degree_6(order, phi, sin_theta, cos_theta, norm): - if order == -6: - return norm[6] * 10395 * sin_theta**6 * np.sin(6 * phi) - elif order == -5: - return norm[5] * 10395 * sin_theta**5 * cos_theta * np.sin(5 * phi) - elif order == -4: - return norm[4] * 472.5 * sin_theta**4 * (11 * cos_theta ** 2 - 1) * np.sin(4 * phi) - elif order == -3: - return norm[3] * 157.5 * sin_theta**3 * (11 * cos_theta ** 3 - 3 * cos_theta) * np.sin(3 * phi) - elif order == -2: - return norm[2] * 13.125 * sin_theta**2 * (33 * cos_theta ** 4 - 18 * cos_theta ** 2 + 1) * np.sin(2 * phi) - elif order == -1: - return norm[1] * 2.625 * sin_theta \ - * (33 * cos_theta**5 - 30 * cos_theta**3 + 5 * cos_theta) * np.sin(phi) - elif order == 0: - return norm[0] * 0.0625 * (231 * cos_theta ** 6 - 315 * cos_theta ** 4 + 105 * cos_theta ** 2 - 5) - elif order == 1: - return norm[1] * 2.625 * sin_theta \ - * (33 * cos_theta**5 - 30 * cos_theta**3 + 5 * cos_theta) * np.cos(phi) - elif order == 2: - return norm[2] * 13.125 * sin_theta**2 * (33 * cos_theta ** 4 - 18 * cos_theta ** 2 + 1) * np.cos(2 * phi) - elif order == 3: - return norm[3] * 157.5 * sin_theta**3 * (11 * cos_theta ** 3 - 3 * cos_theta) * np.cos(3 * phi) - elif order == 4: - return norm[4] * 472.5 * sin_theta**4 * (11 * cos_theta ** 2 - 1) * np.cos(4 * phi) - elif order == 5: - 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) + odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + array of primary orientation dispersion parameters + odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + array of secondary orientation dispersion parameters -@njit(cache=True) -def sph_harm_degree_8(order, phi, sin_theta, cos_theta, norm): - if order == -8: - return norm[8] * 2027025 * sin_theta**8 * np.sin(8 * phi) - elif order == -7: - return norm[7] * 2027025 * sin_theta**7 * cos_theta * np.sin(7 * phi) - elif order == -6: - return norm[6] * 67567.5 * sin_theta**6 * (15 * cos_theta ** 2 - 1) * np.sin(6 * phi) - elif order == -5: - return norm[5] * 67567.5 * sin_theta**5 * (5 * cos_theta ** 3 - cos_theta) * np.sin(5 * phi) - elif order == -4: - return norm[4] * 1299.375 * sin_theta**4 * (65 * cos_theta ** 4 - 26 * cos_theta ** 2 + 1) * np.sin(4 * phi) - elif order == -3: - return norm[3] * 433.125 * sin_theta**3 \ - * (39 * cos_theta**5 - 26 * cos_theta**3 + 3 * cos_theta) * np.sin(3 * phi) - elif order == -2: - return norm[2] * 19.6875 * sin_theta**2 \ - * (143 * cos_theta**6 - 143 * cos_theta**4 + 33 * cos_theta**2 - 1) * np.sin(2 * phi) - elif order == -1: - return norm[1] * 0.5625 * sin_theta \ - * (715 * cos_theta**7 - 1001 * cos_theta**5 + 385 * cos_theta**3 - 35 * cos_theta) * np.sin(phi) - elif order == 0: - return norm[0] * 0.0078125 \ - * (6435 * cos_theta**8 - 12012 * cos_theta**6 + 6930 * cos_theta**4 - 1260 * cos_theta**2 + 35) - elif order == 1: - return norm[1] * 0.5625 * sin_theta \ - * (715 * cos_theta**7 - 1001 * cos_theta**5 + 385 * cos_theta**3 - 35 * cos_theta) * np.cos(phi) - elif order == 2: - return norm[2] * 19.6875 * sin_theta**2 \ - * (143 * cos_theta**6 - 143 * cos_theta**4 + 33 * cos_theta**2 - 1) * np.cos(2 * phi) - elif order == 3: - return norm[3] * 433.125 * sin_theta**3 \ - * (39 * cos_theta**5 - 26 * cos_theta**3 + 3 * cos_theta) * np.cos(3 * phi) - elif order == 4: - return norm[4] * 1299.375 * sin_theta**4 * (65 * cos_theta ** 4 - 26 * cos_theta ** 2 + 1) * np.cos(4 * phi) - elif order == 5: - return norm[5] * 67567.5 * sin_theta**5 * (5 * cos_theta ** 3 - cos_theta) * np.cos(5 * phi) - elif order == 6: - return norm[6] * 67567.5 * sin_theta**6 * (15 * cos_theta ** 2 - 1) * np.cos(6 * phi) - elif order == 7: - 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) + odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + array of total orientation dispersion parameters + odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + array of orientation dispersion anisotropy parameters -@njit(cache=True) -def sph_harm_degree_10(order, phi, sin_theta, cos_theta, norm): - if order == -10: - return norm[10] * 654729075 * sin_theta**10 * np.sin(10 * phi) - elif order == -9: - return norm[9] * 654729075 * sin_theta**9 * cos_theta * np.sin(9 * phi) - elif order == -8: - return norm[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta ** 2 - 1) * np.sin(8 * phi) - elif order == -7: - return norm[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta ** 3 - 3 * cos_theta) * np.sin(7 * phi) - elif order == -6: - return norm[6] * 84459.375 * sin_theta**6 \ - * (323 * cos_theta**4 - 102 * cos_theta**2 + 3) * np.sin(6 * phi) - elif order == -5: - return norm[5] * 16891.875 * sin_theta**5 \ - * (323 * cos_theta**5 - 170 * cos_theta**3 + 15 * cos_theta) * np.sin(5 * phi) - elif order == -4: - return norm[4] * 2815.3125 * sin_theta**4 \ - * (323 * cos_theta**6 - 255 * cos_theta**4 + 45 * cos_theta**2 - 1) * np.sin(4 * phi) - elif order == -3: - return norm[3] * 402.1875 * sin_theta**3 \ - * (323 * cos_theta**7 - 357 * cos_theta**5 + 105 * cos_theta**3 - 7 * cos_theta) * np.sin(3 * phi) - elif order == -2: - return norm[2] * 3.8671875 * sin_theta**2 \ - * (4199 * cos_theta**8 - 6188 * cos_theta**6 + 2730 * cos_theta**4 - 364 * cos_theta**2 + 7) \ - * np.sin(2 * phi) - elif order == -1: - return norm[1] * 0.4296875 * sin_theta \ - * (4199 * cos_theta**9 - 7956 * cos_theta**7 + 4914 * cos_theta**5 - 1092 * cos_theta**3 + 63 * cos_theta) \ - * np.sin(phi) - elif order == 0: - return norm[0] * 0.00390625 \ - * (46189 * cos_theta**10 - 109395 * cos_theta**8 + 90090 * cos_theta**6 - - 30030 * cos_theta**4 + 3465 * cos_theta**2 - 63) - elif order == 1: - return norm[1] * 0.4296875 * sin_theta \ - * (4199 * cos_theta**9 - 7956 * cos_theta**7 + 4914 * cos_theta**5 - - 1092 * cos_theta**3 + 63 * cos_theta) * np.cos(phi) - elif order == 2: - return norm[2] * 3.8671875 * sin_theta**2 \ - * (4199 * cos_theta**8 - 6188 * cos_theta**6 + 2730 * cos_theta**4 - 364 * cos_theta**2 + 7) \ - * np.cos(2 * phi) - elif order == 3: - return norm[3] * 402.1875 * sin_theta**3 \ - * (323 * cos_theta**7 - 357 * cos_theta**5 + 105 * cos_theta**3 - 7 * cos_theta) * np.cos(3 * phi) - elif order == 4: - return norm[4] * 2815.3125 * sin_theta**4 \ - * (323 * cos_theta**6 - 255 * cos_theta**4 + 45 * cos_theta**2 - 1) * np.cos(4 * phi) - elif order == 5: - return norm[5] * 16891.875 * sin_theta**5 \ - * (323 * cos_theta**5 - 170 * cos_theta**3 + 15 * cos_theta) * np.cos(5 * phi) - elif order == 6: - return norm[6] * 84459.375 * sin_theta**6 \ - * (323 * cos_theta**4 - 102 * cos_theta**2 + 3) * np.cos(6 * phi) - elif order == 7: - return norm[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta ** 3 - 3 * cos_theta) * np.cos(7 * phi) - elif order == 8: - return norm[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta ** 2 - 1) * np.cos(8 * phi) - elif order == 9: - 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) + Returns + ------- + None + """ + # mask background + bg_msk = np.all(vec_tensor_eigen == 0, axis=-1) + for key in odi.keys(): + if odi[key] is not None: + odi[key][bg_msk] = 0 diff --git a/foa3d/output.py b/foa3d/output.py index 06b9721..c45fb5e 100644 --- a/foa3d/output.py +++ b/foa3d/output.py @@ -2,16 +2,17 @@ import tempfile from datetime import datetime -from os import makedirs, path +from os import makedirs, mkdir, path import nibabel as nib import numpy as np from tifffile import TiffWriter +from foa3d.printing import print_flsh from foa3d.utils import get_item_size -def create_save_dirs(img_path, img_name, cli_args, is_fovec=False): +def create_save_dirs(img_path, img_name, cli_args, is_vec=False): """ Create saving directory. @@ -26,19 +27,16 @@ def create_save_dirs(img_path, img_name, cli_args, is_fovec=False): cli_args: see ArgumentParser.parse_args updated namespace of command line arguments - is_fovec: bool + is_vec: bool True when fiber orientation vectors are provided as input to the pipeline Returns ------- - save_subdirs: list (dtype=str) - saving subdirectory string paths - - tmp_dir: str - temporary directory (for memory-map objects) + save_dirs: dict + saving directories + ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) """ - # get current time time_stamp = datetime.now().strftime("%Y%m%d-%H%M%S") @@ -49,30 +47,33 @@ def create_save_dirs(img_path, img_name, cli_args, is_fovec=False): # create saving directory base_out_dir = path.join(out_path, f"Foa3D_{time_stamp}_{img_name}") - save_dir_lst = list() if not path.isdir(base_out_dir): makedirs(base_out_dir) + # initialize empty dictionary + save_dirs = dict() + # create Frangi filter output subdirectory - if not is_fovec: + if not is_vec: frangi_dir = path.join(base_out_dir, 'frangi') - makedirs(frangi_dir) - save_dir_lst.append(frangi_dir) + mkdir(frangi_dir) + save_dirs['frangi'] = frangi_dir else: - save_dir_lst.append(None) + save_dirs['frangi'] = None # create ODF analysis output subdirectory if cli_args.odf_res is not None: odf_dir = path.join(base_out_dir, 'odf') - makedirs(odf_dir) - save_dir_lst.append(odf_dir) + mkdir(odf_dir) + save_dirs['odf'] = odf_dir else: - save_dir_lst.append(None) + save_dirs['odf'] = None # create temporary directory - tmp_dir = tempfile.mkdtemp(dir=out_path) + tmp_dir = tempfile.mkdtemp(dir=base_out_dir) + save_dirs['tmp'] = tmp_dir - return save_dir_lst, tmp_dir + return save_dirs def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): @@ -103,7 +104,6 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): ------- None """ - # get maximum RAM and initialized array memory size if ram is None: ram = psutil.virtual_memory()[1] @@ -118,18 +118,20 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): # retrieve image pixel size px_sz_z, px_sz_y, px_sz_x = px_sz + # adjust metadata + metadata = {'spacing': px_sz_z, 'unit': 'um'} if nd_array.ndim == 4 \ + else {'axes': 'ZYX', 'spacing': px_sz_z, 'unit': 'um'} + # adjust bigtiff optional argument bigtiff = True if nd_array.itemsize * nd_array.size >= 4294967296 else False out_name = f'{fname}.{fmt}' with TiffWriter(path.join(save_dir, out_name), bigtiff=bigtiff, append=True) as tif: for z in range(nz): zs = z * dz - tif.write(nd_array[zs:zs + dz, ...], contiguous=True, resolution=(1 / px_sz_x, 1 / px_sz_y), - metadata={'spacing': px_sz_z, 'unit': 'um'}) - - # save array to NumPy file - elif fmt == 'npy': - np.save(path.join(save_dir, fname + '.npy'), nd_array) + tif.write(nd_array[zs:zs + dz, ...], + contiguous=True, + resolution=(1 / px_sz_x, 1 / px_sz_y), + metadata=metadata) # save array to NIfTI file elif fmt == 'nii': @@ -139,3 +141,120 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None): # raise error else: raise ValueError("Unsupported data format!!!") + + +def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None): + """ + Save the output arrays of the Frangi filter stage to TIF files. + + Parameters + ---------- + save_dir: str + saving directory string path + + img_name: str + name of the input microscopy image + + out_img: dict + fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field + + fbr_vec_clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image + + fa_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image + + frangi_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability) + + iso_fbr: 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) + neuron mask image + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size (Z,Y,X) [μm] + + ram: float + maximum RAM available + + Returns + ------- + 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) + + # print output directory + print_flsh(f"\nFrangi filter arrays saved to: {save_dir}\n") + + +def save_odf_arrays(save_dir, img_name, odf_scale_um, px_sz, odf, bg, odi_pri, odi_sec, odi_tot, odi_anis): + """ + Save the output arrays of the ODF analysis stage to TIF and Nifti files. + Arrays tagged with 'mrtrixview' are preliminarily transformed + so that ODF maps viewed in MRtrix3 are spatially consistent + with the analyzed microscopy volume, and the output TIF files. + + Parameters + ---------- + save_dir: str + saving directory string path + + img_name: str + name of the 3D microscopy image + + odf_scale_um: float + fiber ODF resolution (super-voxel side [μm]) + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size (Z,Y,X) [μm] + + odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) + ODF spherical harmonics coefficients + + bg: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8) + background for ODF visualization in MRtrix3 + + odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + primary orientation dispersion parameter + + odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + secondary orientation dispersion parameter + + odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + total orientation dispersion parameter + + odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) + orientation dispersion anisotropy parameter + + Returns + ------- + None + """ + # save ODF image with background to NIfTI files (adjusted view for MRtrix3) + sbfx = f'{odf_scale_um}_{img_name}' + save_array(f'bg_mrtrixview_sv{sbfx}', save_dir, bg, fmt='nii') + save_array(f'odf_mrtrixview_sv{sbfx}', save_dir, odf, fmt='nii') + + # save total orientation dispersion + save_array(f'odi_tot_sv{sbfx}', save_dir, odi_tot, px_sz) + + # save primary orientation dispersion + if odi_pri is not None: + save_array(f'odi_pri_sv{sbfx}', save_dir, odi_pri, px_sz) + + # save secondary orientation dispersion + if odi_sec is not None: + save_array(f'odi_sec_sv{sbfx}', save_dir, odi_sec, px_sz) + + # save orientation dispersion anisotropy + if odi_anis is not None: + save_array(f'odi_anis_sv{sbfx}', save_dir, odi_anis, px_sz) diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py index ac69484..710943e 100644 --- a/foa3d/pipeline.py +++ b/foa3d/pipeline.py @@ -3,315 +3,296 @@ import numpy as np from joblib import Parallel, delayed -from foa3d.frangi import frangi_filter -from foa3d.input import get_image_info, get_frangi_config, get_resource_config +from foa3d.frangi import (frangi_filter, init_frangi_arrays, mask_background, + write_frangi_arrays) +from foa3d.input import get_frangi_config, get_resource_config from foa3d.odf import (compute_odf_map, generate_odf_background, - get_sph_harm_ncoeff, get_sph_harm_norm_factors) -from foa3d.output import save_array -from foa3d.preprocessing import correct_image_anisotropy -from foa3d.printing import print_flushed, print_frangi_info, print_odf_info -from foa3d.slicing import (config_frangi_batch, generate_slice_lists, - crop, crop_lst, slice_image) -from foa3d.utils import (create_background_mask, create_memory_map, - divide_nonzero, elapsed_time, - get_available_cores, hsv_orient_cmap, - rgb_orient_cmap) - - -slc_cnt = 0 - - -def compute_fractional_anisotropy(eigenval): + 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.spharm import get_sph_harm_norm_factors +from foa3d.utils import get_available_cores + + +def parallel_frangi_over_slices(cli_args, save_dirs, in_img): """ - Compute structure tensor fractional anisotropy - as in Schilling et al. (2018). + Apply 3D Frangi filtering to basic TPFM image slices using parallel threads. Parameters ---------- - eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - structure tensor eigenvalues (best local spatial scale) + cli_args: see ArgumentParser.parse_args + populated namespace of command line arguments - Returns - ------- - fa: numpy.ndarray (shape=(3,), dtype=float) - fractional anisotropy - """ + save_dirs: str + saving directories - fa = np.sqrt(0.5 * divide_nonzero( - np.square((eigenval[..., 0] - eigenval[..., 1])) + - np.square((eigenval[..., 0] - eigenval[..., 2])) + - np.square((eigenval[..., 1] - eigenval[..., 2])), - np.sum(eigenval ** 2, axis=-1))) + in_img: dict + input image dictionary - return fa + 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 -def init_frangi_volumes(img_shp, slc_shp, rsz_ratio, tmp_dir, z_rng=(0, None), msk_bc=False, exp_all=False): - """ - Initialize the output datasets of the Frangi filtering stage. + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - Parameters - ---------- - img_shp: numpy.ndarray (shape=(3,), dtype=int) - volume image shape [px] + fb_ch: int + neuronal fibers channel - slc_shp: numpy.ndarray (shape=(3,), dtype=int) - shape of the basic image slices - analyzed using parallel threads [px] + bc_ch: int + brain cell soma channel - rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio + 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] - z_rng: int - output z-range in [px] + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - msk_bc: bool - if True, mask neuronal bodies - in the optionally provided image channel + name: str + name of the 3D microscopy image - exp_all: bool - export all images + 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] Returns ------- - fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) - initialized fiber orientation volume image + out_img: dict + output image dictionary - fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) - initialized orientation colormap image + fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - fa_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - initialized fractional anisotropy image + fbr_clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - frangi_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - initialized Frangi-enhanced image + fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - initialized fiber mask 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 or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - initialized fiber image (isotropic resolution) + iso_fbr_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image - bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - initialized soma mask image + fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fiber mask image - z_sel: NumPy slice object - selected z-depth range + 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] """ + # skip Frangi filter stage if orientation vectors were directly 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 - # shape copies - img_shp = img_shp.copy() - slc_shp = slc_shp.copy() - - # adapt output z-axis shape if required - z_min, z_max = 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) - tot_shp = tuple(np.ceil(rsz_ratio * img_shp).astype(int)) - - # fiber channel arrays - iso_fbr_img = create_memory_map(tot_shp, dtype='uint8', name='iso_fiber', tmp_dir=tmp_dir) - if 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) - - # soma channel array - if msk_bc: - bc_msk = create_memory_map(tot_shp, dtype='uint8', name='bc_msk', tmp_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) + # get resources configuration + ram, jobs = get_resource_config(cli_args) + + # get Frangi filter configuration + frangi_cfg, px_sz_iso = get_frangi_config(cli_args, in_img) + + # configure the batch of basic image slices analyzed using parallel threads + batch_sz, slc_shp, slc_shp_um, rsz_ratio, ovlp, ovlp_rsz = \ + config_frangi_batch(in_img, frangi_cfg, px_sz_iso, ram=ram, jobs=jobs) + + # generate image range lists + slc_rng, out_slc_shp, tot_slc, batch_sz = \ + generate_slice_lists(slc_shp, in_img['shape'], batch_sz, rsz_ratio, ovlp, in_img['msk_bc']) + + # initialize output arrays + out_img, z_sel = init_frangi_arrays(frangi_cfg, in_img['shape'], out_slc_shp, rsz_ratio, save_dirs['tmp'], + msk_bc=in_img['msk_bc']) + + # print Frangi filter configuration + print_frangi_info(in_img, frangi_cfg, slc_shp_um, tot_slc) + + # parallel Frangi filter-based fiber orientation analysis of microscopy image slices + print_flsh(f"[Parallel(n_jobs={batch_sz})]: Using backend ThreadingBackend with {batch_sz} concurrent workers.") + 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)) - return fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, fbr_msk, iso_fbr_img, bc_msk, z_sel + # save output arrays + save_frangi_arrays(save_dirs['frangi'], in_img['name'], out_img, px_sz_iso, ram=ram) + # add isotropic pixel size + out_img['px_sz'] = px_sz_iso -def init_odf_volumes(vec_img_shp, tmp_dir, odf_scale, odf_degrees=6, exp_all=False): + 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): """ - Initialize the output datasets of the ODF analysis stage. + Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image. Parameters ---------- - vec_img_shp: tuple - vector volume shape [px] + s: int + image slice index - tmp_dir: str - temporary file directory + in_img: dict + input image dictionary (extended) - odf_scale: int - fiber ODF resolution (super-voxel side [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 - odf_degrees: int - degrees of the spherical harmonics series expansion + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - exp_all: bool - export all images + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - Returns - ------- - odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) - initialized array of ODF spherical harmonics coefficients + fb_ch: int + neuronal fibers channel - bg_mrtrix: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8) - initialized background for ODF visualization in MRtrix3 + bc_ch: int + brain cell soma channel - odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of primary orientation dispersion parameters + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of secondary orientation dispersion parameters + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] - odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of total orientation dispersion parameters + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of orientation dispersion anisotropy parameters + name: str + name of the 3D microscopy image - disarray: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - local angular disarray + is_vec: bool + vector field flag - vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) - initialized array of fiber orientation tensor eigenvalues - """ + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape - # ODI maps shape - odi_shp = tuple(np.ceil(np.divide(vec_img_shp, odf_scale)).astype(int)) - - # create downsampled background memory map - bg_shp = tuple(np.flip(odi_shp)) - bg_mrtrix = create_memory_map(bg_shp, dtype='uint8', name=f'bg_tmp{odf_scale}', tmp_dir=tmp_dir) - - # create ODF memory map - nc = get_sph_harm_ncoeff(odf_degrees) - 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}'.format(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) - 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) - disarray = create_memory_map(odi_shp, dtype='float32', name=f'disarray{odf_scale}', tmp_dir=tmp_dir) - else: - odi_pri, odi_sec, odi_anis, disarray = (None, None, None, None) + shape_um: numpy.ndarray (shape=(3,), dtype=float) + total image shape [μm] - return odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray, vec_tensor_eigen + item_sz: int + image item size [B] + out_img: dict + output image dictionary -def fiber_analysis(img, in_rng, bc_rng, out_rng, pad, ovlp, smooth_sigma, scales_px, px_rsz_ratio, z_sel, - fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, ts_msk=None, - bc_ch=0, fb_ch=1, ch_ax=None, alpha=0.05, beta=1, gamma=100, dark=False, msk_bc=False, - print_info=None, hsv_vec_cmap=False, pad_mode='reflect', fb_thr='li', bc_thr='yen', ts_thr=0.0001, - verbose=10): - """ - Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image. + vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - Parameters - ---------- - img: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X)) - fiber fluorescence volume image + clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - in_rng: NumPy slice object - input image range (fibers) + fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - bc_rng: NumPy slice object - input image range (neurons) + frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image (fiber probability image) - out_rng: NumPy slice object - output range + iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image - pad: numpy.ndarray (axis order=(Z,Y,X)) - image padding array + fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + fiber mask image - ovlp: int - overlapping range between slices along each axis [px] + bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + soma mask image - smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the smoothing Gaussian filter [px] + px_sz: numpy.ndarray (shape=(3,), dtype=float) + output pixel size [μm] - scales_px: numpy.ndarray (dtype=int) - spatial scales [px] + cfg: dict + Frangi filter configuration - px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) - 3D image resize ratio + alpha: float + plate-like score sensitivity - z_sel: NumPy slice object - selected z-depth range + beta: float + blob-like score sensitivity - fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector image + gamma: float + background score sensitivity - fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] - fa_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] - frangi_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - Frangi-enhanced volume image (fiber probability volume) + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] - iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - fiber mask image + z_rng: int + output z-range in [px] - bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - soma mask image + bc_ch: int + neuronal bodies channel - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + fb_ch: int + myelinated fibers channel - bc_ch: int - neuronal bodies channel + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - fb_ch: int - myelinated fibers channel + hsv_cmap: bool + generate HSV colormap of 3D fiber orientations - ch_ax: int - RGB image channel axis (either 1 or 3, or None for grayscale images) + exp_all: bool + export all images - alpha: float - plate-like score sensitivity + ovlp: int + overlapping range between slices along each axis [px] - beta: float - blob-like score sensitivity + rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio - gamma: float - background score sensitivity + z_sel: NumPy slice object + selected z-depth range - dark: bool - if True, enhance black 3D tubular structures - (i.e., negative contrast polarity) + in_rng: NumPy slice object + input image range (fibers) - hsv_vec_cmap: bool - if True, generate an HSV orientation color map based on XY-plane orientation angles - (instead of an RGB map using the cartesian components of the estimated vectors) + bc_rng: NumPy slice object + input image range (neurons) + + out_rng: NumPy slice object + output range - msk_bc: bool - if True, mask neuronal bodies exploiting the autofluorescence - signal of lipofuscin pigments + in_pad: numpy.ndarray (axis order=(Z,Y,X)) + image padding range print_info: tuple optional printed information @@ -319,12 +300,6 @@ def fiber_analysis(img, in_rng, bc_rng, out_rng, pad, ovlp, smooth_sigma, scales pad_mode: str image padding mode - is_tiled: bool - must be True for tiled reconstructions aligned using ZetaStitcher - - fb_thr: str - thresholding method applied to the myelinated fibers channel - bc_thr: str thresholding method applied to the neuronal bodies channel @@ -338,402 +313,241 @@ def fiber_analysis(img, in_rng, bc_rng, out_rng, pad, ovlp, smooth_sigma, scales ------- 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']) - # slice fiber image slice - fbr_slc, ts_msk_slc = slice_image(img, in_rng, fb_ch, ch_ax, ts_msk=ts_msk) - - # skip background slice - crt = np.count_nonzero(ts_msk_slc) / np.prod(ts_msk_slc.shape) > ts_thr \ - if ts_msk_slc is not None else np.max(fbr_slc) != 0 - if crt: - - # preprocess fiber slice - iso_fbr_slc, pad_rsz, ts_msk_slc_rsz = \ - correct_image_anisotropy(fbr_slc, px_rsz_ratio, sigma=smooth_sigma, pad=pad, ts_msk=ts_msk_slc) - - # pad fiber slice if required - if pad_rsz is not None: - iso_fbr_slc = np.pad(iso_fbr_slc, pad_rsz, mode=pad_mode) - - # 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') - - # 3D Frangi filter - frangi_slc, fbr_vec_slc, eigenval_slc = \ - frangi_filter(iso_fbr_slc, scales_px=scales_px, alpha=alpha, beta=beta, gamma=gamma, dark=dark) - - # crop resulting slices - iso_fbr_slc, frangi_slc, fbr_vec_slc, eigenval_slc, ts_msk_slc_rsz = \ - crop_lst([iso_fbr_slc, frangi_slc, fbr_vec_slc, eigenval_slc, ts_msk_slc_rsz], out_rng, ovlp) - - # generate fractional anisotropy image - fa_slc = compute_fractional_anisotropy(eigenval_slc) - - # generate fiber orientation color map - fbr_clr_slc = hsv_orient_cmap(fbr_vec_slc) if hsv_vec_cmap else rgb_orient_cmap(fbr_vec_slc) - - # remove background - fbr_vec_slc, fbr_clr_slc, fa_slc, fbr_msk_slc = \ - mask_background(frangi_slc, fbr_vec_slc, fbr_clr_slc, fa_slc, - ts_msk=ts_msk_slc_rsz, method=fb_thr, invert=False) + # process non-background image slice + is_valid = check_background(fbr_slc, ts_msk=ts_msk_slc, ts_thr=ts_thr) + if is_valid: + orient_slc, frangi_slc, iso_slc, fbr_msk_slc = \ + analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng[s], in_pad[s], pad_mode=pad_mode, ts_msk=ts_msk_slc) # (optional) neuronal body masking - if msk_bc: + if in_img['msk_bc']: # get soma image slice - bc_slc, _ = slice_image(img, bc_rng, bc_ch, ch_ax) - - # resize soma slice (lateral blurring and downsampling) - iso_bc_slc, _, _ = correct_image_anisotropy(bc_slc, px_rsz_ratio) + bc_slc, _ = slice_image(in_img['data'], bc_rng[s], in_img['ch_ax'], ch=in_img['bc_ch']) - # crop isotropized soma slice - iso_bc_slc = crop(iso_bc_slc, out_rng) - - # mask neuronal bodies - fbr_vec_slc, fbr_clr_slc, fa_slc, bc_msk_slc = \ - mask_background(iso_bc_slc, fbr_vec_slc, fbr_clr_slc, fa_slc, method=bc_thr, invert=True) + # suppress soma contribution + orient_slc, bc_msk_slc = reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng[s], bc_thr) + # soma mask not available else: bc_msk_slc = None # fill memory-mapped output arrays - fill_frangi_volumes(fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, - fbr_vec_slc, fbr_clr_slc, fa_slc, frangi_slc, iso_fbr_slc, - fbr_msk_slc, bc_msk_slc, out_rng, z_sel) + write_frangi_arrays(out_rng[s], z_sel, iso_slc, frangi_slc, fbr_msk_slc, bc_msk_slc, **orient_slc, **out_img) # print progress if print_info is not None: - global slc_cnt - slc_cnt += 1 - - # print only every N=verbose image slices - start_time, batch_sz, tot_slc_num = print_info - if slc_cnt % verbose == 0 or slc_cnt == tot_slc_num: - prog_prc = 100 * slc_cnt / tot_slc_num - _, hrs, mins, secs = elapsed_time(start_time) - print_flushed(f"[Parallel(n_jobs={batch_sz})]:\t{slc_cnt}/{tot_slc_num} done\t|\telapsed: {hrs} hrs {mins} min {secs:.1f} s\t{prog_prc:.1f}%") + print_frangi_progress(print_info, is_valid, verbose=verbose) -def fill_frangi_volumes(fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, fbr_vec_slc, - fbr_clr_slc, fa_slc, frangi_slc, iso_fbr_slc, fbr_msk_slc, bc_msk_slc, out_rng, z_sel): +def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflect', ts_msk=None): """ - Fill the memory-mapped output arrays of the Frangi filter stage. + Analyze 3D fiber orientations exploiting a Frangi-filter-based + unsupervised enhancement of tubular structures. Parameters ---------- - fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector image - - fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image - - fa_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image - - frangi_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - Frangi-enhanced image (fiber probability image) + fbr_slc: numpy.ndarray (axis order=(Z,Y,X)) + fiber image slice - iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image + cfg: dict + Frangi filter configuration - fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - fiber mask image + alpha: float + plate-like score sensitivity - bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - soma mask image + beta: float + blob-like score sensitivity - fbr_vec_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - fiber orientation vector image slice + gamma: float + background score sensitivity - fbr_clr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - orientation colormap image slice + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - fractional anisotropy image slice + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] - frangi_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - Frangi-enhanced image slice + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] - iso_fbr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image slice + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] - fbr_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - fiber mask image slice + z_rng: int + output z-range in [px] - bc_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) - soma mask image slice + bc_ch: int + neuronal bodies channel - out_rng: tuple - 3D slice output index range + fb_ch: int + myelinated fibers channel - z_sel: NumPy slice object - selected z-depth range + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - Returns - ------- - None - """ - - # fill memory-mapped output arrays - vec_rng_out = tuple(np.append(out_rng, slice(0, 3, 1))) - fbr_vec_img[vec_rng_out] = fbr_vec_slc[z_sel, ...] - fbr_vec_clr[vec_rng_out] = fbr_clr_slc[z_sel, ...] - iso_fbr_img[out_rng] = iso_fbr_slc[z_sel, ...].astype(np.uint8) - - # optional output images: Frangi filter response - if frangi_img is not None: - frangi_img[out_rng] = (255 * frangi_slc[z_sel, ...]).astype(np.uint8) - - # optional output images: fractional anisotropy - if fa_img is not None: - fa_img[out_rng] = fa_slc[z_sel, ...] + hsv_cmap: bool + generate HSV colormap of 3D fiber orientations - # 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) + exp_all: bool + export all images - # 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) - - -def mask_background(img, fbr_vec_slc, fbr_clr_slc, fa_slc=None, method='yen', invert=False, ts_msk=None): - """ - Mask fiber orientation arrays. - - Parameters - ---------- - img: numpy.ndarray (axis order=(Z,Y,X)) - fiber (or neuron) fluorescence volume image - - fbr_vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - fiber orientation vector slice + ovlp: int + overlapping range between slices along each axis [px] - fbr_clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - fiber orientation colormap slice + rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - fractional anisotropy slice + out_rng: NumPy slice object + output range - method: str - thresholding method (refer to skimage.filters) + pad: numpy.ndarray (axis order=(Z,Y,X)) + image padding range - invert: bool - mask inversion flag + pad_mode: str + image padding mode ts_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask Returns ------- - fbr_vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - orientation vector patch (masked) - - orientcol_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap patch (masked) + orient_slc: dict + slice orientation dictionary - fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) - fractional anisotropy patch (masked) + vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - bg: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) - background mask - """ - - # generate background mask - bg = create_background_mask(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)) - - # invert mask - if invert: - bg = np.logical_not(bg) + clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - # apply mask to input arrays - fbr_vec_slc[bg, :] = 0 - fbr_clr_slc[bg, :] = 0 - fa_slc[bg] = 0 + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - return fbr_vec_slc, fbr_clr_slc, fa_slc, bg + frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + Frangi-enhanced image slice (fiber probability image) + iso_fbr_img: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image slice -def odf_analysis(fbr_vec_img, iso_fbr_img, px_sz_iso, save_dir, tmp_dir, img_name, odf_scale_um, odf_norm, - odf_deg=6, exp_all=False): + fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fiber mask image slice """ - Estimate 3D fiber ODFs from basic orientation data chunks using parallel threads. + # 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) - Parameters - ---------- - fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vectors - - iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber volume - - px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) - adjusted isotropic pixel size [μm] + # pad fiber slice if required + if pad_rsz is not None: + iso_fbr_slc = np.pad(iso_fbr_slc, pad_rsz, mode=pad_mode) - save_dir: str - saving directory string path + # 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') - tmp_dir: str - temporary file directory - - img_name: str - name of the input volume image + # 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']) - odf_scale_um: float - fiber ODF resolution (super-voxel side [μm]) + # 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) - odf_norm: numpy.ndarray (dtype: float) - 2D array of spherical harmonics normalization factors + # 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) - odf_deg: int - degrees of the spherical harmonics series expansion + return orient_slc, frangi_slc, iso_fbr_slc, fbr_msk_slc - exp_all: bool - export all images - Returns - ------- - None +def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'): """ - - # derive the ODF kernel size in [px] - odf_scale = int(np.ceil(odf_scale_um / px_sz_iso[0])) - - # initialize ODF analysis output volumes - odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray, tensor \ - = init_odf_volumes(fbr_vec_img.shape[:-1], tmp_dir, odf_scale=odf_scale, odf_degrees=odf_deg, exp_all=exp_all) - - # generate downsampled background for MRtrix3 mrview - bg_img = fbr_vec_img if iso_fbr_img is None else iso_fbr_img - generate_odf_background(bg_img, bg_mrtrix, vxl_side=odf_scale) - - # compute ODF coefficients - odf = compute_odf_map(fbr_vec_img, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, tensor, - odf_scale, odf_norm, odf_deg=odf_deg) - - # save memory maps to file - save_odf_arrays(odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray, - px_sz_iso, save_dir, img_name, odf_scale_um) - - -def parallel_frangi_on_slices(img, ch_ax, cli_args, save_dir, tmp_dir, img_name, ts_msk=None): - """ - Apply 3D Frangi filtering to basic TPFM image slices using parallel threads. + Suppress soma contribution using the optional image channel + of neuronal bodies. Parameters ---------- - img: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X)) - microscopy volume image + bc_slc: numpy.ndarray (axis order=(Z,Y,X)) + brain cell soma image slice - ch_ax: int - RGB image channel axis (either 1 or 3, or None for grayscale images) + orient_slc: dict + slice orientation dictionary - cli_args: see ArgumentParser.parse_args - populated namespace of command line arguments + vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 3D fiber orientation field - save_dir: str - saving directory string path + clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + orientation colormap image - tmp_dir: str - temporary file directory + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fractional anisotropy image - img_name: str - name of the microscopy volume image + rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio - ts_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + out_rng: NumPy slice object + output range + + bc_thr: str + thresholding method applied to the neuronal bodies channel Returns ------- - fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector image - - iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image - - px_sz_iso: int - isotropic pixel size [μm] - - img_name: str - name of the input microscopy image - """ - - # get resources configuration - ram, jobs = get_resource_config(cli_args) - - # get Frangi filter configuration - alpha, beta, gamma, frangi_sigma, frangi_sigma_um, smooth_sigma, px_sz, px_sz_iso, \ - z_rng, bc_ch, fb_ch, fb_thr, msk_bc, hsv_vec_cmap, exp_all = get_frangi_config(cli_args) - - # get info about the input microscopy image - img_shp, img_shp_um, img_item_sz, fb_ch, msk_bc = get_image_info(img, px_sz, msk_bc, fb_ch, ch_ax) + orient_slc: dict + (masked) slice orientation dictionary - # configure the batch of basic image slices analyzed using parallel threads - batch_sz, in_slc_shp, in_slc_shp_um, px_rsz_ratio, ovlp, ovlp_rsz = \ - config_frangi_batch(px_sz, px_sz_iso, img_shp, img_item_sz, smooth_sigma, frangi_sigma_um, ram=ram, jobs=jobs) + vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + (masked) 3D fiber orientation field - # get info about the processed image slices - in_rng_lst, in_pad_lst, out_rng_lst, bc_rng_lst, out_slc_shp, tot_slc_num, batch_sz = \ - generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp, msk_bc=msk_bc, jobs=jobs) + clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + (masked) orientation colormap image - # initialize output arrays - fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, fbr_msk, iso_fbr_img, bc_msk, z_sel = \ - init_frangi_volumes(img_shp, out_slc_shp, px_rsz_ratio, tmp_dir, z_rng=z_rng, msk_bc=msk_bc, exp_all=exp_all) + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + (masked) fractional anisotropy image - # print Frangi filter configuration - print_frangi_info(alpha, beta, gamma, frangi_sigma_um, img_shp_um, in_slc_shp_um, tot_slc_num, - px_sz, img_item_sz, msk_bc) + bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + brain cell mask + """ + # resize soma slice (lateral downsampling) + iso_bc, _, _ = correct_anisotropy(bc_slc, rsz_ratio) - # parallel Frangi filter-based fiber orientation analysis of microscopy image slices - print_flushed(f"[Parallel(n_jobs={batch_sz})]: Using backend ThreadingBackend with {batch_sz} concurrent workers.") - start_time = perf_counter() - with Parallel(n_jobs=batch_sz, prefer='threads', require='sharedmem') as parallel: - parallel( - delayed(fiber_analysis)(img, in_rng_lst[i], bc_rng_lst[i], out_rng_lst[i], in_pad_lst[i], ovlp_rsz, - smooth_sigma, frangi_sigma, px_rsz_ratio, z_sel, fbr_vec_img, fbr_vec_clr, - fa_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, ts_msk=ts_msk, bc_ch=bc_ch, - fb_ch=fb_ch, fb_thr=fb_thr, ch_ax=ch_ax, alpha=alpha, beta=beta, gamma=gamma, - msk_bc=msk_bc, hsv_vec_cmap=hsv_vec_cmap, - print_info=(start_time, batch_sz, tot_slc_num)) - for i in range(tot_slc_num)) + # crop isotropic brain cell soma slice + iso_bc = crop(iso_bc, out_rng) - # save output arrays - save_frangi_arrays(fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, fbr_msk, bc_msk, px_sz_iso, - save_dir, img_name, ram=ram) + # mask neuronal bodies + orient_slc, bc_msk = mask_background(iso_bc, orient_slc, method=bc_thr, invert=True) - return fbr_vec_img, iso_fbr_img, px_sz_iso + return orient_slc, bc_msk -def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir, tmp_dir, img_name, - backend='loky', verbose=100): +def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_name, backend='loky', verbose=100): """ Iterate over the required spatial scales and apply the parallel ODF analysis implemented in parallel_odf_on_slices(). Parameters ---------- - fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) - fiber orientation vector image - - iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - isotropic fiber image - cli_args: see ArgumentParser.parse_args populated namespace of command line arguments - px_sz: numpy.ndarray (axis order=(3,), dtype=float) - pixel size [μm] + save_dirs: dict + saving directories + ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) + + fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fiber orientation vector field - save_dir: str - saving directory string path + iso_fbr: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image - tmp_dir: str - temporary file directory + px_sz: numpy.ndarray (axis order=(3,), dtype=float) + pixel size [μm] img_name: str name of the input volume image @@ -741,9 +555,6 @@ def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir, backend: str backend module employed by joblib.Parallel - ram: float - maximum RAM available to the ODF analysis stage [B] - verbose: int joblib verbosity level @@ -751,10 +562,9 @@ def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir, ------- 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 @@ -764,157 +574,73 @@ def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir, # get number of logical cores num_cpu = get_available_cores() - # generate pixel size if not provided + # 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 required spatial scales + # 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_img, iso_fbr_img, px_sz, save_dir, tmp_dir, img_name, + 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_flushed(f"\nODF and dispersion maps saved to: {save_dir}\n") + print_flsh(f"\nODF and dispersion maps saved to: {save_dirs['odf']}\n") -def save_frangi_arrays(fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, fbr_msk, bc_msk, px_sz, save_dir, img_name, - ram=None): +def odf_analysis(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_scale_um, odf_norm, odf_deg=6, exp_all=False): """ - Save the output arrays of the Frangi filter stage to TIF files. + Estimate 3D fiber ODFs from basic orientation data chunks using parallel threads. Parameters ---------- - fbr_vec_img: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vector field - fbr_vec_clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap image - - fa_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - fractional anisotropy image - - frangi_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - Frangi-enhanced image (fiber probability) - - 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) - neuron mask image + iso_fbr: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) + isotropic fiber image px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size (Z,Y,X) [μm] + adjusted isotropic pixel size [μm] - save_dir: str - saving directory string path + save_dirs: dict + saving directories + ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files) img_name: str - name of the input microscopy image - - ram: float - maximum RAM available - - Returns - ------- - None - """ - - # save orientation vectors to TIFF - save_array(f'fiber_vec_{img_name}', save_dir, np.moveaxis(fbr_vec_img, -1, 1), px_sz, ram=ram) - - # save orientation color map to TIFF - save_array(f'fiber_cmap_{img_name}', save_dir, fbr_vec_clr, px_sz, ram=ram) - - # save fractional anisotropy map to TIFF - if fa_img is not None: - save_array(f'frac_anis_{img_name}', save_dir, fa_img, px_sz, ram=ram) - - # save Frangi-enhanced fiber volume to TIFF - if frangi_img is not None: - save_array(f'frangi_{img_name}', save_dir, frangi_img, px_sz, ram=ram) - - # save masked fiber volume to TIFF - if fbr_msk is not None: - save_array(f'fiber_msk_{img_name}', save_dir, fbr_msk, px_sz, ram=ram) - - # save masked soma volume to TIFF - if bc_msk is not None: - save_array(f'soma_msk_{img_name}', save_dir, bc_msk, px_sz, ram=ram) + name of the 3D microscopy image - # print output directory - print_flushed(f"\nFrangi filter arrays saved to: {save_dir}\n") - - -def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, disarray, px_sz, save_dir, img_name, odf_scale_um): - """ - Save the output arrays of the ODF analysis stage to TIF and Nifti files. - Arrays tagged with 'mrtrixview' are preliminarily transformed - so that ODF maps viewed in MRtrix3 are spatially consistent - with the analyzed microscopy volume, and the output TIF files. - - Parameters - ---------- - odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) - ODF spherical harmonics coefficients - - bg: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8) - background for ODF visualization in MRtrix3 - - odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - primary orientation dispersion parameter - - odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - secondary orientation dispersion parameter - - odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - total orientation dispersion parameter - - odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - orientation dispersion anisotropy parameter - - disarray: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32) - - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size (Z,Y,X) [μm] + odf_scale_um: float + fiber ODF resolution (super-voxel side [μm]) - save_dir: str - saving directory string path + odf_norm: numpy.ndarray (dtype: float) + 2D array of spherical harmonic normalization factors - img_name: str - name of the input volume image + odf_deg: int + degrees of the spherical harmonic series expansion - odf_scale_um: float - fiber ODF resolution (super-voxel side [μm]) + exp_all: bool + export all images Returns ------- None """ + # get the ODF kernel size in [px] + odf_scale = int(np.ceil(odf_scale_um / px_sz[0])) - # save ODF image with background to NIfTI files (adjusted view for MRtrix3) - sbfx = f'{odf_scale_um}_{img_name}' - save_array(f'bg_mrtrixview_sv{sbfx}', save_dir, bg, fmt='nii') - save_array(f'odf_mrtrixview_sv{sbfx}', save_dir, odf, fmt='nii') - - # save total orientation dispersion - save_array(f'odi_tot_sv{sbfx}', save_dir, odi_tot, px_sz) + # initialize ODF analysis output arrays + odf, odi, 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) - # save primary orientation dispersion - if odi_pri is not None: - save_array(f'odi_pri_sv{sbfx}', save_dir, odi_pri, px_sz) - - # save secondary orientation dispersion - if odi_sec is not None: - save_array(f'odi_sec_sv{sbfx}', save_dir, odi_sec, px_sz) + # generate downsampled background for MRtrix3 mrview + generate_odf_background(bg_mrtrix, fbr_vec, vxl_side=odf_scale, iso_fbr=iso_fbr) - # save orientation dispersion anisotropy - if odi_anis is not None: - save_array(f'odi_anis_sv{sbfx}', save_dir, odi_anis, px_sz) + # compute ODF coefficients + odf = compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, odf_scale, odf_norm, odf_deg=odf_deg) - # save angular disarray - if disarray is not None: - save_array(f'disarray_sv{sbfx}', save_dir, disarray, px_sz) + # save memory maps to TIFF or NIfTI files + save_odf_arrays(save_dirs['odf'], img_name, odf_scale_um, px_sz, odf, bg_mrtrix, **odi) diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py index 53e4105..4f65851 100644 --- a/foa3d/preprocessing.py +++ b/foa3d/preprocessing.py @@ -2,8 +2,8 @@ from scipy.ndimage import gaussian_filter from skimage.transform import resize -from foa3d.printing import (print_blur, print_flushed, - print_new_res, print_prepro_heading) +from foa3d.printing import (print_blur, print_flsh, print_output_res, + print_prepro_heading) from foa3d.utils import fwhm_to_sigma @@ -34,7 +34,6 @@ 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,)) @@ -65,18 +64,18 @@ def config_anisotropy_correction(px_sz, psf_fwhm): # print pixel resize info if cndt_2: - print_new_res(px_sz_iso) + print_output_res(px_sz_iso) else: - print_flushed() + print_flsh() # skip line else: - print_flushed() + print_flsh() return smooth_sigma, px_sz_iso -def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', +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 diff --git a/foa3d/printing.py b/foa3d/printing.py index 131f1d0..84a0c19 100644 --- a/foa3d/printing.py +++ b/foa3d/printing.py @@ -1,14 +1,10 @@ -import os -from platform import system - import numpy as np from foa3d.utils import elapsed_time -# adjust ANSI escape sequence -# decoding to Windows OS -if system == 'Windows': - os.system("color") + +# slice progress counter +slc_cnt = 0 def color_text(r, g, b, text): @@ -39,7 +35,7 @@ def color_text(r, g, b, text): return clr_text -def print_flushed(string_to_print=""): +def print_flsh(string_to_print=""): """ Print string and flush output data buffer. @@ -55,27 +51,119 @@ def print_flushed(string_to_print=""): print(string_to_print, flush=True) -def print_frangi_info(alpha, beta, gamma, scales_um, img_shp_um, in_slc_shp_um, - tot_slc_num, px_sz, img_item_sz, msk_bc): +def print_blur(sigma_um, psf_fwhm): + """ + Print the standard deviation of the smoothing Gaussian filter. + + Parameters + ---------- + sigma_um: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [μm] + (resolution anisotropy correction) + + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] + + Returns + ------- + None + """ + psf_sz = np.max(psf_fwhm) + print_flsh(f"Gaussian blur σ [μm]: ({sigma_um[0]:.3f}, {sigma_um[1]:.3f}, {sigma_um[2]:.3f})\n" + + 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): """ Print Frangi filter stage heading. Parameters ---------- - alpha: float - plate-like score sensitivity + in_img: dict + input image dictionary (extended) - beta: float - blob-like score sensitivity + 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 - gamma: float - background score sensitivity + ts_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask - scales_um: list (dtype=float) - analyzed spatial scales [μm] + ch_ax: int + RGB image channel axis (either 1, 3, or None for grayscale images) - img_shp_um: numpy.ndarray (shape=(3,), dtype=float) - volume image shape [μm] + 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 + + shape: numpy.ndarray (shape=(3,), dtype=int) + total image shape + + shape_um: numpy.ndarray (shape=(3,), dtype=float) + total image shape [μm] + + item_sz: int + image item size [B] + + frangi_cfg: dict + Frangi filter configuration + + alpha: float + plate-like score sensitivity + + beta: float + blob-like score sensitivity + + gamma: float + background score sensitivity + + scales_px: numpy.ndarray (dtype=float) + Frangi filter scales [px] + + scales_um: numpy.ndarray (dtype=float) + Frangi filter scales [μm] + + smooth_sd: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the smoothing Gaussian filter [px] + + px_sz: numpy.ndarray (shape=(3,), dtype=float) + pixel size [μm] + + z_rng: int + output z-range in [px] + + 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 in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=float) shape of the analyzed image slices [μm] @@ -83,56 +171,91 @@ def print_frangi_info(alpha, beta, gamma, scales_um, img_shp_um, in_slc_shp_um, tot_slc_num: int total number of analyzed image slices - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] - - img_item_sz: int - image item size (in bytes) - - msk_bc: bool - if True, mask neuronal bodies within - the optionally provided channel - Returns ------- None """ + # print Frangi filter sensitivity and scales + alpha = frangi_cfg['alpha'] + beta = frangi_cfg['beta'] + gamma = frangi_cfg['gamma'] - scales_um = np.asarray(scales_um) + scales_um = np.asarray(frangi_cfg['scales_um']) if gamma is None: gamma = 'auto' - print_flushed(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_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(img_shp_um, in_slc_shp_um, tot_slc_num, px_sz, img_item_sz) + 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_soma_masking(msk_bc) + print_soma_masking(in_img['msk_bc']) -def print_blur(sigma_um, psf_fwhm): +def print_frangi_progress(info, is_valid, verbose=1): """ - Print the standard deviation of the smoothing Gaussian filter. + Print Frangi filter progress. Parameters ---------- - sigma_um: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the smoothing Gaussian filter [μm] - (resolution anisotropy correction) + info: tuple + info to be printed out - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - 3D FWHM of the PSF [μm] + is_valid: bool + image slice validity flag + + verbose: int + verbosity level (print info only every "verbose" slices) Returns ------- None """ - psf_sz = np.max(psf_fwhm) - print_flushed(f"Gaussian blur σ [μm]: ({sigma_um[0]:.3f}, {sigma_um[1]:.3f}, {sigma_um[2]:.3f})\n" + \ - f"Adjusted PSF FWHM [μm]: ({psf_sz:.3f}, {psf_sz:.3f}, {psf_sz:.3f})") + global slc_cnt + 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 + _, 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}%") + + +def print_image_info(in_img): + """ + Print information on the input microscopy image (shape, voxel size, PSF size). + + Parameters + ---------- + in_img: dict + input image dictionary + ('img_data': image data, 'ts_msk': tissue sample mask, 'ch_ax': channel axis) + + Returns + ------- + None + """ + # get pixel and PSF sizes + px_sz = in_img['px_sz'] + psf_fwhm = in_img['psf_fwhm'] + + # get channel axis (RGB image only) + ch_ax = in_img['ch_ax'] + + # get image shape (ignore channel axis) + img_shp = in_img['data'].shape + if ch_ax is not None: + img_shp = np.delete(img_shp, ch_ax) + + print_flsh("\n Z Y X\nImage shape [μm]: " + + f"({img_shp[0] * px_sz[0]:.1f}, {img_shp[1] * px_sz[1]:.1f}, {img_shp[2] * px_sz[2]:.1f})\n" + + f"Pixel size [μm]: ({px_sz[0]:.3f}, {px_sz[1]:.3f}, {px_sz[2]:.3f})\n" + + f"PSF FWHM [μm]: ({psf_fwhm[0]:.3f}, {psf_fwhm[1]:.3f}, {psf_fwhm[2]:.3f})") def print_import_time(start_time): @@ -149,7 +272,7 @@ def print_import_time(start_time): None """ _, _, mins, secs = elapsed_time(start_time) - print_flushed(f"Image loaded in: {mins} min {secs:3.1f} s") + print_flsh(f"Image loaded in: {mins} min {secs:3.1f} s") def print_odf_info(odf_scales_um, odf_degrees): @@ -168,68 +291,48 @@ def print_odf_info(odf_scales_um, odf_degrees): ------- None """ - print_flushed(color_text(0, 191, 255, "\n3D ODF Analysis\n") + \ - f"\nResolution [μm]: {odf_scales_um}\n" + \ - f"Expansion degrees: {odf_degrees}\n") + print_flsh( + color_text(0, 191, 255, + f"\n3D ODF Analysis\n\nResolution [μm]: {odf_scales_um}\nExpansion degrees: {odf_degrees}\n")) -def print_pipeline_heading(): - """ - Print Foa3D tool heading. - - Returns - ------- - None +def print_output_res(px_sz_iso): """ - print_flushed(color_text(0, 250, 154, "\n3D Fiber Orientation Analysis")) - + Print the adjusted isotropic spatial resolution. -def print_prepro_heading(): - """ - Print preprocessing heading. + Parameters + ---------- + px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) + new isotropic pixel size [μm] Returns ------- None """ - print_flushed(color_text(0, 191, 255, "\n\nMicroscopy Image Preprocessing\n") + \ - "\n Z Y X") + print_flsh(f"Adjusted pixel size [μm]: ({px_sz_iso[0]:.3f}, {px_sz_iso[1]:.3f}, {px_sz_iso[2]:.3f})\n") -def print_native_res(px_sz, psf_fwhm): +def print_pipeline_heading(): """ - Print the native pixel size and optical resolution of the microscopy system. - - Parameters - ---------- - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] - - psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) - PSF 3D FWHM [μm] + Print Foa3D tool heading. Returns ------- None """ - print_flushed(f"Pixel size [μm]: ({px_sz[0]:.3f}, {px_sz[1]:.3f}, {px_sz[2]:.3f})\n" + \ - f"PSF FWHM [μm]: ({psf_fwhm[0]:.3f}, {psf_fwhm[1]:.3f}, {psf_fwhm[2]:.3f})") + print_flsh(color_text(0, 250, 154, "\n3D Fiber Orientation Analysis")) -def print_new_res(px_sz_iso): +def print_prepro_heading(): """ - Print the adjusted isotropic spatial resolution. - - Parameters - ---------- - px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) - new isotropic pixel size [μm] + Print preprocessing heading. Returns ------- None """ - print_flushed(f"Adjusted pixel size [μm]: ({px_sz_iso[0]:.3f}, {px_sz_iso[1]:.3f}, {px_sz_iso[2]:.3f})\n") + print_flsh(color_text(0, 191, 255, "\n\nMicroscopy Image Preprocessing\n") + + "\n Z Y X") def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): @@ -257,7 +360,6 @@ def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): ------- None """ - # adjust slice shape if np.any(img_shp_um < slc_shp_um): slc_shp_um = img_shp_um @@ -269,12 +371,12 @@ def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): max_slc_sz = img_item_sz * np.prod(np.divide(slc_shp_um, px_sz)) # print info - print_flushed("\n Z Y X") - print_flushed(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") + print_flsh("\n Z Y X") + print_flsh(f"Total image shape [μm]: ({img_shp_um[0]:.1f}, {img_shp_um[1]:.1f}, {img_shp_um[2]:.1f})\n" + + f"Total image size [MB]: {np.ceil(img_sz / 1024**2).astype(int)}\n\n" + + f"Image slice shape [μm]: ({slc_shp_um[0]:.1f}, {slc_shp_um[1]:.1f}, {slc_shp_um[2]:.1f})\n" + + f"Image slice size [MB]: {np.ceil(max_slc_sz / 1024**2).astype(int)}\n" + + f"Image slice number: {tot_slc_num}\n") def print_soma_masking(msk_bc): @@ -293,39 +395,6 @@ def print_soma_masking(msk_bc): """ prt = 'Soma mask: ' if msk_bc: - print_flushed(f'{prt}active\n') + print_flsh(f'{prt}active\n') else: print(f'{prt}not active\n') - - -def print_image_shape(cli_args, img, ch_ax): - """ - Print 3D microscopy image shape. - - Parameters - ---------- - cli_args: see ArgumentParser.parse_args - populated namespace of command line arguments - - img: numpy.ndarray (shape=(Z,Y,X)) - microscopy volume image - - ch_ax: int - RGB image channel axis (either 1 or 3) - - Returns - ------- - None - """ - - # get pixel size - px_sz_z = cli_args.px_size_z - px_sz_xy = cli_args.px_size_xy - - # get image shape (ignore channel axis) - img_shp = img.shape - if ch_ax is not None: - img_shp = np.delete(img_shp, ch_ax) - - print_flushed("\n Z Y X\nImage shape [μm]: " + \ - f"({img_shp[0] * px_sz_z:.1f}, {img_shp[1] * px_sz_xy:.1f}, {img_shp[2] * px_sz_xy:.1f})") diff --git a/foa3d/slicing.py b/foa3d/slicing.py index ea97bd2..148caa6 100644 --- a/foa3d/slicing.py +++ b/foa3d/slicing.py @@ -44,7 +44,6 @@ def adjust_axis_range(ax_iter, img_shp, slc_per_ax, start, stop, ovlp=0): pad: numpy.ndarray (shape=(2,), dtype=int) lower and upper padding ranges [px] """ - # initialize slice padding array pad = np.zeros((2,), dtype=np.int64) @@ -68,6 +67,22 @@ def adjust_axis_range(ax_iter, img_shp, slc_per_ax, start, stop, ovlp=0): return start, stop, pad +def check_background(img, ts_msk=None, ts_thr=0.0001): + """ + Description + + Parameters + ---------- + + Returns + ------- + is_valid: bool + """ + is_valid = np.count_nonzero(ts_msk) / np.prod(ts_msk.shape) > ts_thr if ts_msk is not None else np.max(img) != 0 + + return is_valid + + @njit(cache=True) def compute_axis_range(ax, ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip=False): """ @@ -108,7 +123,6 @@ def compute_axis_range(ax, ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip= pad: numpy.ndarray (shape=(2,), dtype=int) lower and upper padding ranges [px] """ - # compute start and stop coordinates start = ax_iter[ax] * slc_shp[ax] stop = start + slc_shp[ax] @@ -158,7 +172,6 @@ def compute_slice_range(ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip=Fal pad: np.ndarray (shape=(3,2), dtype=int) slice padding range [px] """ - # generate axis range and padding array dim = len(ax_iter) slc = tuple() @@ -290,32 +303,18 @@ def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1): return slc_sz -def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi_sigma_um, - mem_growth=149.7, mem_fudge=1.0, jobs=None, ram=None, shp_thr=7): +def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fudge=1.0, jobs=None, ram=None, shp_thr=7): """ Compute size and number of the batches of basic microscopy image slices analyzed in parallel. Parameters ---------- - px_sz: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] + px_sz_iso: int isotropic pixel size [μm] - img_shp: numpy.ndarray (shape=(3,), dtype=int) - total image shape [px] - - item_sz: int - image item size [B] - - smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the smoothing Gaussian filter [px] - - frangi_sigma_um: numpy.ndarray (dtype=float) - Frangi filter scales [μm] - mem_growth: float empirical memory growth factor of the Frangi filtering stage @@ -346,7 +345,7 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi shape of the basic image slices analyzed using parallel threads [μm] - px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) 3D image resize ratio ovlp: int @@ -365,16 +364,17 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi jobs = num_cpu # number of spatial scales - ns = len(frangi_sigma_um) + 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 - px_rsz_ratio = np.divide(px_sz, px_sz_iso) + rsz_ratio = np.divide(in_img['px_sz'], px_sz_iso) # compute slice overlap (boundary artifacts suppression) - ovlp, ovlp_rsz = compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=px_rsz_ratio) + ovlp, ovlp_rsz = compute_overlap(frangi_cfg['smooth_sd'], frangi_cfg['scales_um'], + px_rsz_ratio=rsz_ratio) # compute the shape of basic microscopy image slices in_slc_shp = np.array([-1]) @@ -385,9 +385,10 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi "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(img_shp, item_sz, px_sz=px_sz, slc_sz=slc_sz, ovlp=ovlp) + in_slc_shp, in_slc_shp_um = compute_slice_shape(in_img['shape'], in_img['item_sz'], + px_sz=in_img['px_sz'], slc_sz=slc_sz, ovlp=ovlp) - return batch_sz, in_slc_shp, in_slc_shp_um, px_rsz_ratio, ovlp, ovlp_rsz + return batch_sz, in_slc_shp, in_slc_shp_um, rsz_ratio, ovlp, ovlp_rsz def crop(img, rng, ovlp=None, flip=()): @@ -413,7 +414,6 @@ def crop(img, rng, ovlp=None, flip=()): img_out: numpy.ndarray cropped microscopy image """ - # delete overlapping boundaries if ovlp is not None: img = img[ovlp[0]:img.shape[0] - ovlp[0], ovlp[1]:img.shape[1] - ovlp[1], ovlp[2]:img.shape[2] - ovlp[2]] @@ -454,15 +454,18 @@ def crop_lst(img_lst, rng, ovlp=None, flip=()): img_lst: list list of cropped image slices """ - for s, img in enumerate(img_lst): if img is not None: - img_lst[s] = crop(img, rng, ovlp=ovlp, flip=flip) + 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) return img_lst -def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bc=False, jobs=None): +def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bc=False): """ Generate image slice ranges for the Frangi filtering stage. @@ -488,33 +491,30 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms if True, mask neuronal bodies in the optionally provided image channel - jobs: int - number of parallel jobs (threads) - Returns ------- - in_rng_lst: list - list of input slice index ranges + slc_rng: dict + in_rng: list + list of input slice index ranges - in_pad_lst: list - list of slice padding arrays + in_pad: list + list of slice padding arrays - out_rng_lst: list - list of output slice index ranges + out_rng: list + list of output slice index ranges - bc_rng_lst: list - (optional) list of neuronal body slice index ranges + bc_rng: list + (optional) list of neuronal body slice index ranges out_slc_shp: numpy.ndarray (shape=(3,), dtype=int) shape of the processed image slices [px] - tot_slc_num: int + tot_slc: int total number of analyzed image slices batch_sz: int adjusted slice batch size """ - # 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) @@ -522,10 +522,6 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms # number of image slices along each axis slc_per_dim = np.floor(np.divide(img_shp, in_slc_shp)).astype(int) - # number of logical cores - if jobs is None: - jobs = get_available_cores() - # initialize empty range lists in_pad_lst = list() in_rng_lst = list() @@ -550,16 +546,23 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms bc_rng_lst.append(None) # total number of slices - tot_slc_num = len(in_rng_lst) + tot_slc = len(in_rng_lst) # adjust slice batch size - if batch_sz > tot_slc_num: - batch_sz = tot_slc_num + if batch_sz > tot_slc: + batch_sz = tot_slc + + # fill slice range dictionary + slc_rng = dict() + slc_rng['in_rng'] = in_rng_lst + slc_rng['in_pad'] = in_pad_lst + slc_rng['out_rng'] = out_rng_lst + slc_rng['bc_rng'] = bc_rng_lst - return in_rng_lst, in_pad_lst, out_rng_lst, bc_rng_lst, out_slc_shp, tot_slc_num, batch_sz + return slc_rng, out_slc_shp, tot_slc, batch_sz -def slice_image(img, rng, ch, ch_ax, ts_msk=None): +def slice_image(img, rng, ch_ax, ch, ts_msk=None): """ Slice desired channel from input image volume. diff --git a/foa3d/spharm.py b/foa3d/spharm.py new file mode 100644 index 0000000..2e36dd4 --- /dev/null +++ b/foa3d/spharm.py @@ -0,0 +1,444 @@ +import numpy as np +from numba import njit + + +@njit(cache=True) +def compute_fiber_angles(fbr_vec, norm): + """ + Estimate the spherical coordinates (φ azimuth and θ polar angles) + of the fiber orientation vectors returned by the Frangi filtering stage + (all-zero background vectors are excluded). + + Parameters + ---------- + fbr_vec: numpy.ndarray (shape=(N,3), dtype=float) + array of fiber orientation vectors + (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) + + norm: numpy.ndarray (shape=(N,), dtype=float) + 2-norm of fiber orientation vectors + + Returns + ------- + phi: numpy.ndarray (shape=(N,), dtype=float) + fiber azimuth angle [rad] + + theta: numpy.ndarray (shape=(N,), dtype=float) + fiber polar angle [rad] + """ + fbr_vec = fbr_vec[norm > 0, :] + phi = np.arctan2(fbr_vec[:, 1], fbr_vec[:, 2]) + theta = np.arccos(fbr_vec[:, 0] / norm[norm > 0]) + + return phi, theta + + +@njit(cache=True) +def compute_real_sph_harm(degree, order, phi, sin_theta, cos_theta, norm_factors): + """ + Estimate the coefficients of the real spherical harmonics series expansion + as described by Alimi et al. (Medical Image Analysis, 2020). + + Parameters + ---------- + degree: int + degree index of the spherical harmonics expansion + + order: int + order index of the spherical harmonics expansion + + phi: float + azimuth angle [rad] + + sin_theta: float + polar angle sine + + cos_theta: float + polar angle cosine + + norm_factors: numpy.ndarray (dtype: float) + normalization factors + + Returns + ------- + real_sph_harm: float + real-valued spherical harmonic coefficient + """ + if degree == 0: + real_sph_harm = norm_factors[0, 0] + elif degree == 2: + real_sph_harm = sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm_factors[1, :]) + elif degree == 4: + real_sph_harm = sph_harm_degree_4(order, phi, sin_theta, cos_theta, norm_factors[2, :]) + elif degree == 6: + real_sph_harm = sph_harm_degree_6(order, phi, sin_theta, cos_theta, norm_factors[3, :]) + elif degree == 8: + real_sph_harm = sph_harm_degree_8(order, phi, sin_theta, cos_theta, norm_factors[4, :]) + elif degree == 10: + real_sph_harm = sph_harm_degree_10(order, phi, sin_theta, cos_theta, norm_factors[5, :]) + else: + raise ValueError("\n Invalid degree of the spherical harmonics series expansion!!!") + + return real_sph_harm + + +factorial_lut = np.array([ + 1, 1, 2, 6, 24, 120, 720, 5040, 40320, + 362880, 3628800, 39916800, 479001600, + 6227020800, 87178291200, 1307674368000, + 20922789888000, 355687428096000, 6402373705728000, + 121645100408832000, 2432902008176640000], dtype=np.double) + + +@njit(cache=True) +def factorial(n): + """ + Retrieve factorial using pre-computed LUT. + + Parameters + ---------- + n: int + integer number (max: 20) + + Returns + ------- + f: int + factorial + """ + if n > 20: + raise ValueError + + return factorial_lut[n] + + +@njit(cache=True) +def fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff): + """ + Generate the real-valued symmetric spherical harmonics series expansion + from fiber φ azimuth and θ polar angles, + i.e. the spherical coordinates of the fiber orientation vectors. + + Parameters + ---------- + phi: numpy.ndarray (shape=(N,), dtype=float) + fiber azimuth angles [rad] + (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) + + theta: numpy.ndarray (shape=(N,), dtype=float) + fiber polar angle [rad] + (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) + + degrees: int + degrees of the spherical harmonics expansion + + norm_factors: numpy.ndarray (dtype: float) + normalization factors + + ncoeff: int + number of spherical harmonics coefficients + + Returns + ------- + real_sph_harm: numpy.ndarray (shape=(ncoeff,), dtype=float) + array of real-valued spherical harmonics coefficients + building the spherical harmonics series expansion + """ + real_sph_harm = np.zeros(ncoeff) + i = 0 + for n in range(0, degrees + 1, 2): + for m in range(-n, n + 1, 1): + for p, t in zip(phi, theta): + real_sph_harm[i] += compute_real_sph_harm(n, m, p, np.sin(t), np.cos(t), norm_factors) + i += 1 + + real_sph_harm /= phi.size + + return real_sph_harm + + +def fiber_vectors_to_sph_harm(fbr_vec, degrees, norm_factors): + """ + Generate the real-valued symmetric spherical harmonics series expansion + from the fiber orientation vectors returned by the Frangi filter stage. + + Parameters + ---------- + fbr_vec: numpy.ndarray (shape=(N,3), dtype=float) + array of fiber orientation vectors + (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) + + degrees: int + degrees of the spherical harmonics expansion + + norm_factors: numpy.ndarray (dtype: float) + normalization factors + + Returns + ------- + real_sph_harm: numpy.ndarray (shape=(ncoeff,), dtype=float) + real-valued spherical harmonics coefficients + """ + fbr_vec.shape = (-1, 3) + ncoeff = get_sph_harm_ncoeff(degrees) + + norm = np.linalg.norm(fbr_vec, axis=-1) + if np.sum(norm) < np.sqrt(fbr_vec.shape[0]): + return np.zeros(ncoeff) + + phi, theta = compute_fiber_angles(fbr_vec, norm) + + real_sph_harm = fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff) + + return real_sph_harm + + +@njit(cache=True) +def get_sph_harm_ncoeff(degrees): + """ + Get the number of coefficients of the real spherical harmonics series + expansion. + + Parameters + ---------- + degrees: int + degrees of the spherical harmonics series expansion + + Returns + ------- + ncoeff: int + number of spherical harmonics coefficients + """ + ncoeff = (2 * (degrees // 2) + 1) * ((degrees // 2) + 1) + + return ncoeff + + +@njit(cache=True) +def get_sph_harm_norm_factors(degrees): + """ + Estimate the normalization factors of the real spherical harmonics series + expansion. + + Parameters + ---------- + degrees: int + degrees of the spherical harmonics series expansion + + Returns + ------- + norm_factors: numpy.ndarray (dtype: float) + 2D array of spherical harmonics normalization factors + """ + norm_factors = np.zeros(shape=(degrees + 1, 2 * degrees + 1)) + for n in range(0, degrees + 1, 2): + for m in range(0, n + 1, 1): + norm_factors[n, m] = norm_factor(n, m) + + norm_factors = norm_factors[::2] + + return norm_factors + + +@njit(cache=True) +def norm_factor(n, m): + """ + Compute the normalization factor of the term of degree n and order m + of the real-valued spherical harmonics series expansion. + + Parameters + ---------- + n: int + degree index + + m: int + order index + + Returns + ------- + nf: float + normalization factor + """ + if m == 0: + nf = np.sqrt((2 * n + 1) / (4 * np.pi)) + else: + nf = (-1)**m * np.sqrt(2) * np.sqrt(((2 * n + 1) / (4 * np.pi) * + (factorial(n - np.abs(m)) / factorial(n + np.abs(m))))) + + return nf + + +@njit(cache=True) +def sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm): + if order == -2: + return norm[2] * 3 * sin_theta**2 * np.sin(2 * phi) + elif order == -1: + return norm[1] * 3 * sin_theta * cos_theta * np.sin(phi) + elif order == 0: + return norm[0] * 0.5 * (3 * cos_theta ** 2 - 1) + elif order == 1: + 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) + + +@njit(cache=True) +def sph_harm_degree_4(order, phi, sin_theta, cos_theta, norm): + if order == -4: + return norm[4] * 105 * sin_theta**4 * np.sin(4 * phi) + elif order == -3: + return norm[3] * 105 * sin_theta**3 * cos_theta * np.sin(3 * phi) + elif order == -2: + return norm[2] * 7.5 * sin_theta**2 * (7 * cos_theta ** 2 - 1) * np.sin(2 * phi) + elif order == -1: + return norm[1] * 2.5 * sin_theta * (7 * cos_theta ** 3 - 3 * cos_theta) * np.sin(phi) + elif order == 0: + return norm[0] * 0.125 * (35 * cos_theta ** 4 - 30 * cos_theta ** 2 + 3) + elif order == 1: + return norm[1] * 2.5 * sin_theta * (7 * cos_theta ** 3 - 3 * cos_theta) * np.cos(phi) + elif order == 2: + return norm[2] * 7.5 * sin_theta**2 * (7 * cos_theta ** 2 - 1) * np.cos(2 * phi) + elif order == 3: + 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) + + +@njit(cache=True) +def sph_harm_degree_6(order, phi, sin_theta, cos_theta, norm): + if order == -6: + return norm[6] * 10395 * sin_theta**6 * np.sin(6 * phi) + elif order == -5: + return norm[5] * 10395 * sin_theta**5 * cos_theta * np.sin(5 * phi) + elif order == -4: + return norm[4] * 472.5 * sin_theta**4 * (11 * cos_theta ** 2 - 1) * np.sin(4 * phi) + elif order == -3: + return norm[3] * 157.5 * sin_theta**3 * (11 * cos_theta ** 3 - 3 * cos_theta) * np.sin(3 * phi) + elif order == -2: + return norm[2] * 13.125 * sin_theta**2 * (33 * cos_theta ** 4 - 18 * cos_theta ** 2 + 1) * np.sin(2 * phi) + elif order == -1: + return norm[1] * 2.625 * sin_theta \ + * (33 * cos_theta**5 - 30 * cos_theta**3 + 5 * cos_theta) * np.sin(phi) + elif order == 0: + return norm[0] * 0.0625 * (231 * cos_theta ** 6 - 315 * cos_theta ** 4 + 105 * cos_theta ** 2 - 5) + elif order == 1: + return norm[1] * 2.625 * sin_theta \ + * (33 * cos_theta**5 - 30 * cos_theta**3 + 5 * cos_theta) * np.cos(phi) + elif order == 2: + return norm[2] * 13.125 * sin_theta**2 * (33 * cos_theta ** 4 - 18 * cos_theta ** 2 + 1) * np.cos(2 * phi) + elif order == 3: + return norm[3] * 157.5 * sin_theta**3 * (11 * cos_theta ** 3 - 3 * cos_theta) * np.cos(3 * phi) + elif order == 4: + return norm[4] * 472.5 * sin_theta**4 * (11 * cos_theta ** 2 - 1) * np.cos(4 * phi) + elif order == 5: + 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) + + +@njit(cache=True) +def sph_harm_degree_8(order, phi, sin_theta, cos_theta, norm): + if order == -8: + return norm[8] * 2027025 * sin_theta**8 * np.sin(8 * phi) + elif order == -7: + return norm[7] * 2027025 * sin_theta**7 * cos_theta * np.sin(7 * phi) + elif order == -6: + return norm[6] * 67567.5 * sin_theta**6 * (15 * cos_theta ** 2 - 1) * np.sin(6 * phi) + elif order == -5: + return norm[5] * 67567.5 * sin_theta**5 * (5 * cos_theta ** 3 - cos_theta) * np.sin(5 * phi) + elif order == -4: + return norm[4] * 1299.375 * sin_theta**4 * (65 * cos_theta ** 4 - 26 * cos_theta ** 2 + 1) * np.sin(4 * phi) + elif order == -3: + return norm[3] * 433.125 * sin_theta**3 \ + * (39 * cos_theta**5 - 26 * cos_theta**3 + 3 * cos_theta) * np.sin(3 * phi) + elif order == -2: + return norm[2] * 19.6875 * sin_theta**2 \ + * (143 * cos_theta**6 - 143 * cos_theta**4 + 33 * cos_theta**2 - 1) * np.sin(2 * phi) + elif order == -1: + return norm[1] * 0.5625 * sin_theta \ + * (715 * cos_theta**7 - 1001 * cos_theta**5 + 385 * cos_theta**3 - 35 * cos_theta) * np.sin(phi) + elif order == 0: + return norm[0] * 0.0078125 \ + * (6435 * cos_theta**8 - 12012 * cos_theta**6 + 6930 * cos_theta**4 - 1260 * cos_theta**2 + 35) + elif order == 1: + return norm[1] * 0.5625 * sin_theta \ + * (715 * cos_theta**7 - 1001 * cos_theta**5 + 385 * cos_theta**3 - 35 * cos_theta) * np.cos(phi) + elif order == 2: + return norm[2] * 19.6875 * sin_theta**2 \ + * (143 * cos_theta**6 - 143 * cos_theta**4 + 33 * cos_theta**2 - 1) * np.cos(2 * phi) + elif order == 3: + return norm[3] * 433.125 * sin_theta**3 \ + * (39 * cos_theta**5 - 26 * cos_theta**3 + 3 * cos_theta) * np.cos(3 * phi) + elif order == 4: + return norm[4] * 1299.375 * sin_theta**4 * (65 * cos_theta ** 4 - 26 * cos_theta ** 2 + 1) * np.cos(4 * phi) + elif order == 5: + return norm[5] * 67567.5 * sin_theta**5 * (5 * cos_theta ** 3 - cos_theta) * np.cos(5 * phi) + elif order == 6: + return norm[6] * 67567.5 * sin_theta**6 * (15 * cos_theta ** 2 - 1) * np.cos(6 * phi) + elif order == 7: + 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) + + +@njit(cache=True) +def sph_harm_degree_10(order, phi, sin_theta, cos_theta, norm): + if order == -10: + return norm[10] * 654729075 * sin_theta**10 * np.sin(10 * phi) + elif order == -9: + return norm[9] * 654729075 * sin_theta**9 * cos_theta * np.sin(9 * phi) + elif order == -8: + return norm[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta ** 2 - 1) * np.sin(8 * phi) + elif order == -7: + return norm[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta ** 3 - 3 * cos_theta) * np.sin(7 * phi) + elif order == -6: + return norm[6] * 84459.375 * sin_theta**6 \ + * (323 * cos_theta**4 - 102 * cos_theta**2 + 3) * np.sin(6 * phi) + elif order == -5: + return norm[5] * 16891.875 * sin_theta**5 \ + * (323 * cos_theta**5 - 170 * cos_theta**3 + 15 * cos_theta) * np.sin(5 * phi) + elif order == -4: + return norm[4] * 2815.3125 * sin_theta**4 \ + * (323 * cos_theta**6 - 255 * cos_theta**4 + 45 * cos_theta**2 - 1) * np.sin(4 * phi) + elif order == -3: + return norm[3] * 402.1875 * sin_theta**3 \ + * (323 * cos_theta**7 - 357 * cos_theta**5 + 105 * cos_theta**3 - 7 * cos_theta) * np.sin(3 * phi) + elif order == -2: + return norm[2] * 3.8671875 * sin_theta**2 \ + * (4199 * cos_theta**8 - 6188 * cos_theta**6 + 2730 * cos_theta**4 - 364 * cos_theta**2 + 7) \ + * np.sin(2 * phi) + elif order == -1: + return norm[1] * 0.4296875 * sin_theta \ + * (4199 * cos_theta**9 - 7956 * cos_theta**7 + 4914 * cos_theta**5 - 1092 * cos_theta**3 + 63 * cos_theta) \ + * np.sin(phi) + elif order == 0: + return norm[0] * 0.00390625 \ + * (46189 * cos_theta**10 - 109395 * cos_theta**8 + 90090 * cos_theta**6 + - 30030 * cos_theta**4 + 3465 * cos_theta**2 - 63) + elif order == 1: + return norm[1] * 0.4296875 * sin_theta \ + * (4199 * cos_theta**9 - 7956 * cos_theta**7 + 4914 * cos_theta**5 + - 1092 * cos_theta**3 + 63 * cos_theta) * np.cos(phi) + elif order == 2: + return norm[2] * 3.8671875 * sin_theta**2 \ + * (4199 * cos_theta**8 - 6188 * cos_theta**6 + 2730 * cos_theta**4 - 364 * cos_theta**2 + 7) \ + * np.cos(2 * phi) + elif order == 3: + return norm[3] * 402.1875 * sin_theta**3 \ + * (323 * cos_theta**7 - 357 * cos_theta**5 + 105 * cos_theta**3 - 7 * cos_theta) * np.cos(3 * phi) + elif order == 4: + return norm[4] * 2815.3125 * sin_theta**4 \ + * (323 * cos_theta**6 - 255 * cos_theta**4 + 45 * cos_theta**2 - 1) * np.cos(4 * phi) + elif order == 5: + return norm[5] * 16891.875 * sin_theta**5 \ + * (323 * cos_theta**5 - 170 * cos_theta**3 + 15 * cos_theta) * np.cos(5 * phi) + elif order == 6: + return norm[6] * 84459.375 * sin_theta**6 \ + * (323 * cos_theta**4 - 102 * cos_theta**2 + 3) * np.cos(6 * phi) + elif order == 7: + return norm[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta ** 3 - 3 * cos_theta) * np.cos(7 * phi) + elif order == 8: + return norm[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta ** 2 - 1) * np.cos(8 * phi) + elif order == 9: + 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) diff --git a/foa3d/utils.py b/foa3d/utils.py index 77c44ee..b2c970a 100644 --- a/foa3d/utils.py +++ b/foa3d/utils.py @@ -34,7 +34,6 @@ def create_background_mask(img, method='yen', black_bg=False): bg_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) boolean background mask """ - # select thresholding method if method == 'li': thresh = threshold_li(img) @@ -87,10 +86,9 @@ 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 + '.npy') + mmap_path = path.join(tmp_dir, name + '.mmap') if path.exists(mmap_path): unlink(mmap_path) @@ -144,7 +142,6 @@ def get_available_cores(): num_cpu: int number of available cores """ - num_cpu = environ.pop('OMP_NUM_THREADS', default=None) num_cpu = cpu_count() if num_cpu is None else int(num_cpu) @@ -165,7 +162,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] @@ -225,7 +221,6 @@ def divide_nonzero(nd_array1, nd_array2, new_val=1e-10): divided: numpy.ndarray divided array """ - divisor = np.copy(nd_array2) divisor[divisor == 0] = new_val divided = np.divide(nd_array1, divisor) @@ -256,7 +251,6 @@ def elapsed_time(start_time): secs: float seconds """ - stop_time = perf_counter() tot = stop_time - start_time @@ -303,7 +297,6 @@ def get_item_bytes(data): bts: int item size in bytes """ - # type byte size try: bts = int(np.iinfo(data.dtype).bits / 8) @@ -328,10 +321,9 @@ def get_config_label(cli_args): cfg_lbl: str Frangi filter configuration label """ - cfg_lbl = '_s' for s in cli_args.scales: - cfg_lbl += f'{s}' + 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}' return cfg_lbl @@ -428,7 +420,6 @@ def normalize_image(img, min_val=None, max_val=None, max_out_val=255.0, dtype=np norm_img: numpy.ndarray normalized image """ - # get min and max values if min_val is None: min_val = np.min(img) @@ -461,7 +452,6 @@ 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]) @@ -534,7 +524,6 @@ def transform_axes(nd_array, flipped=None, swapped=None, expand=None): nd_array: numpy.ndarray transformed data array """ - if flipped is not None: nd_array = np.flip(nd_array, axis=flipped) @@ -571,7 +560,6 @@ 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)