From 4495e962bf1b13de8de9cf4e96c25d258d127a45 Mon Sep 17 00:00:00 2001 From: msorelli Date: Wed, 18 Sep 2024 18:18:38 +0200 Subject: [PATCH] Rename variables --- foa3d/__main__.py | 11 +- foa3d/frangi.py | 4 +- foa3d/input.py | 143 +++++++------- foa3d/odf.py | 54 +++--- foa3d/pipeline.py | 414 ++++++++++++++++++++--------------------- foa3d/preprocessing.py | 53 +++--- foa3d/printing.py | 129 +++++++------ foa3d/slicing.py | 8 +- foa3d/utils.py | 45 +++-- 9 files changed, 426 insertions(+), 435 deletions(-) diff --git a/foa3d/__main__.py b/foa3d/__main__.py index dd8e785..eb30c5d 100644 --- a/foa3d/__main__.py +++ b/foa3d/__main__.py @@ -6,7 +6,7 @@ def foa3d(cli_args): - # load microscopy volume image or array of fiber orientation vectors + # load 3D microscopy image or 4D array of fiber orientation vectors img, ts_msk, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args) # get resources configuration @@ -14,23 +14,22 @@ def foa3d(cli_args): # conduct parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices if not is_fiber: - fiber_vec_img, iso_fiber_img, px_sz, img_name \ + fbr_vec_img, iso_fbr_img, px_sz, img_name \ = parallel_frangi_on_slices(img, cli_args, save_dir[0], tmp_dir, img_name, ts_msk=ts_msk, ram=ram, jobs=jobs, is_tiled=is_tiled) - else: - fiber_vec_img, iso_fiber_img, px_sz = (img, None, None) + 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(fiber_vec_img, iso_fiber_img, cli_args, px_sz, save_dir[1], tmp_dir, img_name, ram=ram) + parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir[1], tmp_dir, img_name, ram=ram) # delete temporary folder delete_tmp_folder(tmp_dir) def main(): - # start Foa3D pipeline by terminal + # start Foa3D by terminal print_pipeline_heading() foa3d(cli_args=get_cli_parser()) diff --git a/foa3d/frangi.py b/foa3d/frangi.py index d8d960b..0efd4c2 100644 --- a/foa3d/frangi.py +++ b/foa3d/frangi.py @@ -214,7 +214,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 - enhanced_img = reject_vesselness_background(vesselness, eigen2, eigen3, dark) + enhanced_img = reject_fiber_background(vesselness, eigen2, eigen3, dark) # stack eigenvalues list eigenval = np.stack(eigenval, axis=-1) @@ -332,7 +332,7 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True return enhanced_img, fiber_vec, eigenval -def reject_vesselness_background(vesselness, eigen2, eigen3, dark): +def reject_fiber_background(vesselness, eigen2, eigen3, dark): """ Reject the fiber background, exploiting the sign of the "secondary" eigenvalues λ2 and λ3. diff --git a/foa3d/input.py b/foa3d/input.py index c9f9758..858a6c7 100644 --- a/foa3d/input.py +++ b/foa3d/input.py @@ -52,7 +52,7 @@ def get_cli_parser(): 'Scientific Reports, 13, pp. 4160.\n', formatter_class=CustomFormatter) cli_parser.add_argument(dest='image_path', - help='path to input microscopy volume image or to 4D array of fiber orientation vectors\n' + help='path to input 3D microscopy image or to 4D array of fiber orientation vectors\n' '* supported formats: .tif (image), ' '.npy (image or fiber vectors), .yml (ZetaStitcher stitch file)\n' '* image axes order: (Z, Y, X)\n' @@ -65,22 +65,20 @@ def get_cli_parser(): help='Frangi background score sensitivity') cli_parser.add_argument('-s', '--scales', nargs='+', type=float, default=[1.25], help='list of Frangi filter scales [μm]') - cli_parser.add_argument('-l', '--lpf-mask', action='store_true', default=False, - help='toggle lipofuscin-based neuronal body masking') cli_parser.add_argument('-j', '--jobs', type=int, default=None, help='number of parallel threads used by the Frangi filtering stage: ' '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 filtering stage [GB]: use all if None') cli_parser.add_argument('-m', '--mmap', action='store_true', default=False, - help='create a memory-mapped array of the microscopy volume image') + help='create a memory-mapped array of the 3D microscopy image') 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 the X axis [μm]') cli_parser.add_argument('--psf-fwhm-y', type=float, default=0.692, help='PSF FWHM along the Y axis [μm]') cli_parser.add_argument('--psf-fwhm-z', type=float, default=2.612, help='PSF FWHM along the Z axis [μm]') - cli_parser.add_argument('--ch-mye', type=int, default=1, help='myelinated fibers channel') - cli_parser.add_argument('--ch-lpf', type=int, default=0, help='lipofuscin channel (soma)') + cli_parser.add_argument('--fb-ch', type=int, default=1, help='myelinated fibers channel') + cli_parser.add_argument('--bc-ch', type=int, default=0, help='neuronal bodies channel') 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, @@ -91,6 +89,8 @@ def get_cli_parser(): help='degrees of the spherical harmonics series expansion (even number between 2 and 10)') cli_parser.add_argument('-o', '--out', type=str, default=None, help='output directory') + cli_parser.add_argument('-c', '--cell-msk', action='store_true', default=False, + help='toggle optional neuronal body masking') cli_parser.add_argument('-t', '--tissue-msk', action='store_true', default=False, help='apply tissue reconstruction mask (binarized MIP)') @@ -100,26 +100,26 @@ def get_cli_parser(): return cli_args -def get_image_info(img, px_sz, mask_lpf, ch_mye, ch_axis=None, is_tiled=False): +def get_image_info(img, px_sz, msk_bc, fb_ch, ch_ax=None, is_tiled=False): """ - Get information on the input microscopy volume image. + Get information on the input 3D microscopy image. Parameters ---------- img: numpy.ndarray - microscopy volume image + 3D microscopy image px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - mask_lpf: bool - if True, mask neuronal bodies exploiting the autofluorescence - signal of lipofuscin pigments + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel - ch_mye: int + fb_ch: int myelinated fibers channel - ch_axis: int + ch_ax: int channel axis is_tiled: bool @@ -127,39 +127,39 @@ def get_image_info(img, px_sz, mask_lpf, ch_mye, ch_axis=None, is_tiled=False): Returns ------- - img_shape: numpy.ndarray (shape=(3,), dtype=int) + img_shp: numpy.ndarray (shape=(3,), dtype=int) volume image shape [px] - img_shape_um: numpy.ndarray (shape=(3,), dtype=float) + img_shp_um: numpy.ndarray (shape=(3,), dtype=float) volume image shape [μm] img_item_sz: int - array item size (in bytes) + image item size [B] - ch_mye: int + fb_ch: int myelinated fibers channel - mask_lpf: bool - if True, mask neuronal bodies exploiting the autofluorescence - signal of lipofuscin pigments + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel """ # adapt channel axis - img_shape = np.asarray(img.shape) - ndim = len(img_shape) + img_shp = np.asarray(img.shape) + ndim = len(img_shp) if ndim == 4: - ch_axis = 1 if is_tiled else -1 + ch_ax = 1 if is_tiled else -1 elif ndim == 3: - ch_mye = None - mask_lpf = False + fb_ch = None + msk_bc = False - # get info on microscopy volume image - if ch_axis is not None: - img_shape = np.delete(img_shape, ch_axis) - img_shape_um = np.multiply(img_shape, px_sz) + # get info on 3D microscopy image + if ch_ax is not None: + img_shp = np.delete(img_shp, ch_ax) + img_shp_um = np.multiply(img_shp, px_sz) img_item_sz = get_item_bytes(img) - return img_shape, img_shape_um, img_item_sz, ch_mye, mask_lpf + return img_shp, img_shp_um, img_item_sz, fb_ch, msk_bc def get_file_info(cli_args): @@ -174,26 +174,26 @@ def get_file_info(cli_args): Returns ------- img_path: str - path to the microscopy volume image + path to the 3D microscopy image img_name: str - name of the microscopy volume image + name of the 3D microscopy image img_fmt: str - format of the microscopy volume image + format of the 3D microscopy image is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher is_mmap: bool - create a memory-mapped array of the microscopy volume image, + 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) - ch_mye: int + fb_ch: int myelinated fibers channel """ @@ -212,10 +212,10 @@ def get_file_info(cli_args): is_tiled = True if img_fmt == 'yml' else False # apply tissue reconstruction mask (binarized MIP) - mip_msk = cli_args.tissue_msk - ch_mye = cli_args.ch_mye + msk_mip = cli_args.tissue_msk + fb_ch = cli_args.fb_ch - return img_path, img_name, img_fmt, is_tiled, is_mmap, mip_msk, ch_mye + return img_path, img_name, img_fmt, is_tiled, is_mmap, msk_mip, fb_ch def get_frangi_config(cli_args, img_name): @@ -228,7 +228,7 @@ def get_frangi_config(cli_args, img_name): populated namespace of command line arguments img_name: str - name of the microscopy volume image + name of the 3D microscopy image Returns ------- @@ -259,15 +259,15 @@ def get_frangi_config(cli_args, img_name): z_rng: int output z-range in [px] - ch_lpf: int + bc_ch: int neuronal bodies channel - ch_mye: int + fb_ch: int myelinated fibers channel - mask_lpf: bool - if True, mask neuronal bodies exploiting the autofluorescence - signal of lipofuscin pigments + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel hsv_vec_cmap: bool @@ -287,9 +287,9 @@ def get_frangi_config(cli_args, img_name): scales_px = scales_um / px_sz_iso[0] # image channels - ch_mye = cli_args.ch_mye - ch_lpf = cli_args.ch_lpf - mask_lpf = cli_args.lpf_mask + fb_ch = cli_args.fb_ch + bc_ch = cli_args.bc_ch + msk_bc = cli_args.cell_msk # fiber orientation colormap hsv_vec_cmap = cli_args.hsv @@ -304,7 +304,7 @@ def get_frangi_config(cli_args, img_name): out_name = '{}img{}'.format(pfx, img_name) return alpha, beta, gamma, scales_px, scales_um, smooth_sigma, \ - px_sz, px_sz_iso, z_rng, ch_lpf, ch_mye, mask_lpf, hsv_vec_cmap, out_name + px_sz, px_sz_iso, z_rng, bc_ch, fb_ch, msk_bc, hsv_vec_cmap, out_name def get_resolution(cli_args): @@ -337,7 +337,7 @@ def get_resolution(cli_args): def get_resource_config(cli_args): """ - Retrieve resource usage configuration of the Foa3D pipeline. + Retrieve resource usage configuration of the Foa3D tool. Parameters ---------- @@ -347,11 +347,11 @@ def get_resource_config(cli_args): Returns ------- max_ram: float - maximum RAM available to the Frangi filtering stage [B] + maximum RAM available to the Frangi filter stage [B] jobs: int number of parallel jobs (threads) - used by the Frangi filtering stage + used by the Frangi filter stage """ # resource parameters (convert maximum RAM to bytes) @@ -365,7 +365,7 @@ def get_resource_config(cli_args): def load_microscopy_image(cli_args): """ - Load microscopy volume image from TIFF, NumPy or ZetaStitcher .yml file. + Load 3D microscopy image from TIFF, NumPy or ZetaStitcher .yml file. Alternatively, the processing pipeline accepts as input NumPy or HDF5 files of fiber orientation vector data: in this case, the Frangi filter stage will be skipped. @@ -378,7 +378,7 @@ def load_microscopy_image(cli_args): Returns ------- img: numpy.ndarray or NumPy memory-map object - microscopy volume image or array of fiber orientation vectors + 3D microscopy image or array of fiber orientation vectors tissue_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask @@ -398,16 +398,13 @@ def load_microscopy_image(cli_args): img_name: str microscopy image filename - - cli_args: see ArgumentParser.parse_args - updated namespace of command line arguments """ # create temporary directory tmp_dir = tempfile.mkdtemp() # retrieve input file information - img_path, img_name, img_fmt, is_tiled, is_mmap, mip_msk, ch_mye = get_file_info(cli_args) + img_path, img_name, img_fmt, is_tiled, is_mmap, msk_mip, fb_ch = get_file_info(cli_args) # import fiber orientation vector data tic = perf_counter() @@ -415,10 +412,10 @@ def load_microscopy_image(cli_args): img, is_fiber = load_orient(img_path, img_name, img_fmt) tissue_msk = None - # import raw microscopy volume image + # import raw 3D microscopy image else: img, tissue_msk, is_fiber = load_raw(img_path, img_name, img_fmt, is_tiled=is_tiled, is_mmap=is_mmap, - tmp_dir=tmp_dir, mip_msk=mip_msk, ch_mye=ch_mye) + tmp_dir=tmp_dir, msk_mip=msk_mip, fb_ch=fb_ch) # print import time print_import_time(tic) @@ -439,13 +436,13 @@ def load_orient(img_path, img_name, img_fmt): Parameters ---------- img_path: str - path to the microscopy volume image + path to the 3D microscopy image img_name: str - name of the microscopy volume image + name of the 3D microscopy image img_fmt: str - format of the microscopy volume image + format of the 3D microscopy image Returns ------- @@ -477,26 +474,26 @@ def load_orient(img_path, img_name, img_fmt): return img, is_fiber -def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None, mip_msk=False, ch_mye=1): +def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None, msk_mip=False, fb_ch=1): """ - Load raw microscopy volume image. + Load 3D microscopy image. Parameters ---------- img_path: str - path to the microscopy volume image + path to the 3D microscopy image img_name: str - name of the microscopy volume image + name of the 3D microscopy image img_fmt: str - format of the microscopy volume image + format of the 3D microscopy image is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher is_mmap: bool - create a memory-mapped array of the microscopy volume image, + create a memory-mapped array of the 3D microscopy image, increasing the parallel processing performance (the image will be preliminarily loaded to RAM) @@ -506,13 +503,13 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir mip_msk: bool apply tissue reconstruction mask (binarized MIP) - ch_mye: int + fb_ch: int myelinated fibers channel Returns ------- img: numpy.ndarray or NumPy memory-map object - microscopy volume image + 3D microscopy image tissue_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask @@ -546,12 +543,12 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img[:], mmap_mode='r') # compute tissue reconstruction mask (binarized MIP) - if mip_msk: + if msk_mip: dims = len(img.shape) if dims == 3: tissue_mip = np.max(img[:], axis=0) elif dims == 4: - img_mye = img[:, ch_mye, :, :] if is_tiled else img[..., ch_mye] + img_mye = img[:, fb_ch, :, :] if is_tiled else img[..., fb_ch] tissue_mip = np.max(img_mye, axis=0) tissue_msk = create_background_mask(tissue_mip, method='li', black_bg=True) diff --git a/foa3d/odf.py b/foa3d/odf.py index 471e06d..0e034e0 100644 --- a/foa3d/odf.py +++ b/foa3d/odf.py @@ -6,14 +6,14 @@ from foa3d.utils import normalize_image, transform_axes -def compute_disarray(fiber_vec): +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 ---------- - fiber_vec: numpy.ndarray (shape=(N,3), dtype=float) + 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) @@ -22,14 +22,14 @@ def compute_disarray(fiber_vec): disarray: float32 local angular disarray """ - fiber_vec = np.delete(fiber_vec, np.all(fiber_vec == 0, axis=-1), axis=0) - disarray = (1 - np.linalg.norm(np.mean(fiber_vec, axis=0))).astype(np.float32) + 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(fiber_vec, norm): +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 @@ -37,7 +37,7 @@ def compute_fiber_angles(fiber_vec, norm): Parameters ---------- - fiber_vec: numpy.ndarray (shape=(N,3), dtype=float) + 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) @@ -53,14 +53,14 @@ def compute_fiber_angles(fiber_vec, norm): fiber polar angle [rad] """ - fiber_vec = fiber_vec[norm > 0, :] - phi = np.arctan2(fiber_vec[:, 1], fiber_vec[:, 2]) - theta = np.arccos(fiber_vec[:, 0] / norm[norm > 0]) + 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(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, vec_tensor_eigen, vxl_side, odf_norm, +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): """ Compute the spherical harmonics coefficients iterating over super-voxels @@ -68,7 +68,7 @@ def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarra Parameters ---------- - fiber_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float) + fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float) fiber orientation vectors odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) @@ -114,18 +114,18 @@ def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarra """ # iterate over ODF super-voxels - ref_vxl_size = min(vxl_side, fiber_vec.shape[0]) * vxl_side**2 - for z in range(0, fiber_vec.shape[0], vxl_side): + ref_vxl_size = min(vxl_side, fbr_vec.shape[0]) * vxl_side**2 + for z in range(0, fbr_vec.shape[0], vxl_side): zmax = z + vxl_side - for y in range(0, fiber_vec.shape[1], vxl_side): + for y in range(0, fbr_vec.shape[1], vxl_side): ymax = y + vxl_side - for x in range(0, fiber_vec.shape[2], vxl_side): + for x in range(0, fbr_vec.shape[2], vxl_side): xmax = x + vxl_side # slice orientation voxel - vec_vxl = fiber_vec[z:zmax, y:ymax, x:xmax, :] + vec_vxl = fbr_vec[z:zmax, y:ymax, x:xmax, :] zerovec = np.count_nonzero(np.all(vec_vxl == 0, axis=-1)) sli_vxl_size = np.prod(vec_vxl.shape[:-1]) @@ -259,14 +259,14 @@ def compute_real_sph_harm(degree, order, phi, sin_theta, cos_theta, norm_factors return real_sph_harm -def compute_vec_tensor_eigen(fiber_vec): +def compute_vec_tensor_eigen(fbr_vec): """ Compute the eigenvalues of the 3x3 orientation tensor obtained from a reshaped super-voxel of fiber orientation vectors. Parameters ---------- - fiber_vec: numpy.ndarray (shape=(N,3), dtype=float) + fbr_vec: numpy.ndarray (shape=(N,3), dtype=float) fiber orientation vectors (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx) @@ -276,10 +276,10 @@ def compute_vec_tensor_eigen(fiber_vec): orientation tensor eigenvalues in ascending order """ - fiber_vec = np.delete(fiber_vec, np.all(fiber_vec == 0, axis=-1), axis=0) - fiber_vec.shape = (-1, 3) + 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)) - for v in fiber_vec: + for v in fbr_vec: t += np.outer(v, v) vec_tensor_eigen = sp.linalg.eigh(t, eigvals_only=True).astype(np.float32) @@ -362,14 +362,14 @@ def fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff): return real_sph_harm -def fiber_vectors_to_sph_harm(fiber_vec, degrees, norm_factors): +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 ---------- - fiber_vec: numpy.ndarray (shape=(N,3), dtype=float) + 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) @@ -385,14 +385,14 @@ def fiber_vectors_to_sph_harm(fiber_vec, degrees, norm_factors): real-valued spherical harmonics coefficients """ - fiber_vec.shape = (-1, 3) + fbr_vec.shape = (-1, 3) ncoeff = get_sph_harm_ncoeff(degrees) - norm = np.linalg.norm(fiber_vec, axis=-1) - if np.sum(norm) < np.sqrt(fiber_vec.shape[0]): + 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(fiber_vec, norm) + phi, theta = compute_fiber_angles(fbr_vec, norm) real_sph_harm = fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff) diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py index 80ad3df..e5cdb92 100644 --- a/foa3d/pipeline.py +++ b/foa3d/pipeline.py @@ -34,34 +34,33 @@ def compute_fractional_anisotropy(eigenval): Returns ------- - frac_anis: numpy.ndarray (shape=(3,), dtype=float) + fa: numpy.ndarray (shape=(3,), dtype=float) fractional anisotropy """ - frac_anis = \ - 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))) + 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 frac_anis + return fa -def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, z_rng=(0, None), mask_lpf=False, ram=None): +def init_frangi_volumes(img_shp, slc_shp, rsz_ratio, tmp_dir, z_rng=(0, None), msk_bc=False, ram=None): """ Initialize the output datasets of the Frangi filtering stage. Parameters ---------- - img_shape: numpy.ndarray (shape=(3,), dtype=int) + img_shp: numpy.ndarray (shape=(3,), dtype=int) volume image shape [px] - slice_shape: numpy.ndarray (shape=(3,), dtype=int) + slc_shp: numpy.ndarray (shape=(3,), dtype=int) shape of the basic image slices analyzed using parallel threads [px] - resize_ratio: numpy.ndarray (shape=(3,), dtype=float) + rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) 3D image resize ratio tmp_dir: str @@ -70,38 +69,38 @@ def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, z_rng=(0, z_rng: int output z-range in [px] - mask_lpf: bool - if True, mask neuronal bodies exploiting the autofluorescence - signal of lipofuscin pigments + msk_bc: bool + if True, mask neuronal bodies + in the optionally provided image channel ram: float maximum RAM available to the Frangi filtering stage [B] Returns ------- - fiber_dset_path: str + fbr_dset_path: str path to initialized fiber orientation HDF5 dataset (not fitting the available RAM) - fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) initialized fiber orientation volume image - fiber_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) + fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) initialized orientation colormap image - frac_anis_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + fa_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized fractional anisotropy image frangi_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized Frangi-enhanced image - fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized fiber mask image - iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized fiber image (isotropic resolution) - soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized soma mask image z_sel: NumPy slice object @@ -109,53 +108,53 @@ def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, z_rng=(0, """ # shape copies - img_shape = img_shape.copy() - slice_shape = slice_shape.copy() + 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 = slice_shape[0] - img_shape[0] = z_max - z_min + 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_shape) - tot_shape = tuple(np.ceil(resize_ratio * img_shape).astype(int)) - slice_shape[0] = tot_shape[0] + img_dims = len(img_shp) + tot_shape = tuple(np.ceil(rsz_ratio * img_shp).astype(int)) + slc_shp[0] = tot_shape[0] # fiber channel arrays - iso_fiber_img, _ = init_volume(tot_shape, dtype='uint8', chunks=slice_shape, name='iso_fiber', tmp=tmp_dir, ram=ram) - frangi_img, _ = init_volume(tot_shape, dtype='uint8', chunks=slice_shape, name='frangi', tmp=tmp_dir, ram=ram) - fiber_msk, _ = init_volume(tot_shape, dtype='uint8', chunks=slice_shape, name='fiber_msk', tmp=tmp_dir, ram=ram) - frac_anis_img, _ = init_volume(tot_shape, dtype='uint8', chunks=slice_shape, name='frac_anis', tmp=tmp_dir, ram=ram) + iso_fbr_img, _ = init_volume(tot_shape, dtype='uint8', chunks=slc_shp, name='iso_fiber', tmp=tmp_dir, ram=ram) + frangi_img, _ = init_volume(tot_shape, dtype='uint8', chunks=slc_shp, name='frangi', tmp=tmp_dir, ram=ram) + fbr_msk, _ = init_volume(tot_shape, dtype='uint8', chunks=slc_shp, name='fbr_msk', tmp=tmp_dir, ram=ram) + fa_img, _ = init_volume(tot_shape, dtype='uint8', chunks=slc_shp, name='fa', tmp=tmp_dir, ram=ram) # soma channel array - if mask_lpf: - soma_msk, _ = init_volume(tot_shape, dtype='uint8', chunks=slice_shape, name='soma_msk', tmp=tmp_dir, ram=ram) + if msk_bc: + bc_msk, _ = init_volume(tot_shape, dtype='uint8', chunks=slc_shp, name='bc_msk', tmp=tmp_dir, ram=ram) else: - soma_msk = None + bc_msk = None # fiber orientation arrays vec_shape = tot_shape + (img_dims,) - vec_slice_shape = tuple(slice_shape) + (img_dims,) - fiber_vec_img, fiber_dset_path = \ + vec_slice_shape = tuple(slc_shp) + (img_dims,) + fbr_vec_img, fbr_dset_path = \ init_volume(vec_shape, dtype='float32', chunks=vec_slice_shape, name='fiber_vec', tmp=tmp_dir, ram=ram) - fiber_vec_clr, _ = \ + fbr_vec_clr, _ = \ init_volume(vec_shape, dtype='uint8', chunks=vec_slice_shape, name='fiber_cmap', tmp=tmp_dir, ram=ram) - return fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_img, \ - frangi_img, fiber_msk, iso_fiber_img, soma_msk, z_sel + return fbr_dset_path, fbr_vec_img, fbr_vec_clr, fa_img, \ + frangi_img, fbr_msk, iso_fbr_img, bc_msk, z_sel -def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, ram=None): +def init_odf_volumes(vec_img_shp, tmp_dir, odf_scale, odf_degrees=6, ram=None): """ Initialize the output datasets of the ODF analysis stage. Parameters ---------- - vec_img_shape: tuple + vec_img_shp: tuple vector volume shape [px] tmp_dir: str @@ -198,37 +197,36 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, ram=None) """ # ODI maps shape - odi_shape = tuple(np.ceil(np.divide(vec_img_shape, odf_scale)).astype(int)) + odi_shp = tuple(np.ceil(np.divide(vec_img_shp, odf_scale)).astype(int)) # create downsampled background memory map - bg_shape = tuple(np.flip(odi_shape)) - bg_mrtrix, _ = init_volume(bg_shape, dtype='uint8', chunks=tuple(bg_shape[:2]) + (1,), + bg_shp = tuple(np.flip(odi_shp)) + bg_mrtrix, _ = init_volume(bg_shp, dtype='uint8', chunks=tuple(bg_shp[:2]) + (1,), name='bg_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) # create ODF memory map - num_coeff = get_sph_harm_ncoeff(odf_degrees) - odf_shape = odi_shape + (num_coeff,) - odf, _ = init_volume(odf_shape, dtype='float32', chunks=(1, 1, 1, num_coeff), + nc = get_sph_harm_ncoeff(odf_degrees) + odf_shp = odi_shp + (nc,) + odf, _ = init_volume(odf_shp, dtype='float32', chunks=(1, 1, 1, nc), name='odf_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) # create orientation tensor memory map - vec_tensor_shape = odi_shape + (3,) - vec_tensor_eigen, _ = \ - init_volume(vec_tensor_shape, dtype='float32', chunks=(1, 1, 1, 3), - name='tensor_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) + vec_tensor_shape = odi_shp + (3,) + vec_tensor_eigen, _ = init_volume(vec_tensor_shape, dtype='float32', + chunks=(1, 1, 1, 3), name='tensor_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) # create ODI memory maps - odi_pri, _ = init_volume(odi_shape, dtype='float32', + odi_pri, _ = init_volume(odi_shp, dtype='float32', name='odi_pri_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) - odi_sec, _ = init_volume(odi_shape, dtype='float32', + odi_sec, _ = init_volume(odi_shp, dtype='float32', name='odi_sec_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) - odi_tot, _ = init_volume(odi_shape, dtype='float32', + odi_tot, _ = init_volume(odi_shp, dtype='float32', name='odi_tot_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) - odi_anis, _ = init_volume(odi_shape, dtype='float32', + odi_anis, _ = init_volume(odi_shp, dtype='float32', name='odi_anis_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) # create disarray map - disarray, _ = init_volume(odi_shape, dtype='float32', + disarray, _ = init_volume(odi_shp, dtype='float32', name='disarray{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) return odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray, vec_tensor_eigen @@ -287,11 +285,10 @@ def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+', return vol, hdf5_path -def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, scales_px, px_rsz_ratio, z_sel, - fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, iso_fiber_img, fiber_msk, soma_msk, - ts_msk=None, ch_lpf=0, ch_mye=1, alpha=0.05, beta=1, gamma=100, dark=False, msk_bc=False, - hsv_vec_cmap=False, pad_mode='reflect', is_tiled=False, fiber_thresh='li', soma_thresh='yen', - mip_thr=0.005): +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, alpha=0.05, beta=1, gamma=100, dark=False, msk_bc=False, + hsv_vec_cmap=False, pad_mode='reflect', is_tiled=False, fb_thr='li', bc_thr='yen', ts_thr=0.005): """ Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image. @@ -300,25 +297,23 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc img: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X)) fiber fluorescence volume image - rng_in: NumPy slice object + in_rng: NumPy slice object input image range (fibers) - rng_in_neu: NumPy slice object + bc_rng: NumPy slice object input image range (neurons) - rng_out: NumPy slice object + out_rng: NumPy slice object output range pad: numpy.ndarray (axis order=(Z,Y,X)) image padding array - ovlp: numpy.ndarray (shape=(3,), dtype=int) - slice overlap range - along each axis side + ovlp: int + overlapping range between slices along each axis [px] smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the low-pass Gaussian filter [px] - (applied to the XY plane) + 3D standard deviation of the smoothing Gaussian filter [px] scales_px: numpy.ndarray (dtype=int) spatial scales [px] @@ -329,34 +324,34 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc z_sel: NumPy slice object selected z-depth range - fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vector image - fiber_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) + fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) orientation colormap image - frac_anis_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + 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 volume image (fiber probability volume) - iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) isotropic fiber image - fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) fiber mask image - soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) soma mask image tissue_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask - ch_lpf: int + bc_ch: int neuronal bodies channel - ch_mye: int + fb_ch: int myelinated fibers channel alpha: float @@ -376,7 +371,7 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc 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) - mask_lpf: bool + msk_bc: bool if True, mask neuronal bodies exploiting the autofluorescence signal of lipofuscin pigments @@ -386,14 +381,14 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc is_tiled: bool must be True for tiled reconstructions aligned using ZetaStitcher - fiber_thresh: str - enhanced fiber channel thresholding method + fb_thr: str + thresholding method applied to the myelinated fibers channel - soma_thresh: str - soma channel thresholding method + bc_thr: str + thresholding method applied to the neuronal bodies channel - mip_thr: float - relative threshold on non-zero tissue MIP pixels + ts_thr: float + relative threshold on non-zero tissue pixels Returns ------- @@ -401,121 +396,119 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc """ # slice fiber image slice - fiber_slice, tissue_msk_slice = slice_image(img, rng_in, ch_mye, ts_msk=ts_msk, is_tiled=is_tiled) + fbr_slc, ts_msk_slc = slice_image(img, in_rng, fb_ch, ts_msk=ts_msk, is_tiled=is_tiled) # skip background slice - crt = np.count_nonzero(tissue_msk_slice) / np.prod(tissue_msk_slice.shape) > mip_thr \ - if tissue_msk_slice is not None else np.max(fiber_slice) != 0 + 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_fiber_slice, rsz_pad, rsz_tissue_msk_slice = \ - correct_image_anisotropy(fiber_slice, px_rsz_ratio, sigma=smooth_sigma, pad=pad, ts_msk=tissue_msk_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 rsz_pad is not None: - iso_fiber_slice = np.pad(iso_fiber_slice, rsz_pad, mode=pad_mode) + 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 rsz_tissue_msk_slice is not None: - rsz_tissue_msk_slice = np.pad(rsz_tissue_msk_slice, rsz_pad, mode='constant') + 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_slice, fiber_vec_slice, eigenval_slice = \ - frangi_filter(iso_fiber_slice, scales_px=scales_px, alpha=alpha, beta=beta, gamma=gamma, dark=dark) + 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_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice, rsz_tissue_msk_slice = \ - crop_lst([iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice, rsz_tissue_msk_slice], - rng_out, ovlp) + 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 - frac_anis_slice = compute_fractional_anisotropy(eigenval_slice) + fa_slc = compute_fractional_anisotropy(eigenval_slc) # generate fiber orientation color map - fiber_clr_slice = hsv_orient_cmap(fiber_vec_slice) if hsv_vec_cmap else rgb_orient_cmap(fiber_vec_slice) + fbr_clr_slc = hsv_orient_cmap(fbr_vec_slc) if hsv_vec_cmap else rgb_orient_cmap(fbr_vec_slc) # remove background - fiber_vec_slice, fiber_clr_slice, frac_anis_slice, fiber_msk_slice = \ - mask_background(frangi_slice, fiber_vec_slice, fiber_clr_slice, frac_anis_slice, - tissue_msk=rsz_tissue_msk_slice, method=fiber_thresh, invert=False) + fbr_vec_slc, fbr_clr_slc, fa_slc, fiber_msk_slice = \ + mask_background(frangi_slc, fbr_vec_slc, fbr_clr_slc, fa_slc, + ts_msk=ts_msk_slc_rsz, method=fb_thr, invert=False) # (optional) neuronal body masking if msk_bc: # get soma image slice - soma_slice = slice_image(img, rng_in_neu, ch_lpf, is_tiled=is_tiled) + bc_slc = slice_image(img, bc_rng, bc_ch, is_tiled=is_tiled) # resize soma slice (lateral blurring and downsampling) - iso_soma_slice, _, _ = correct_image_anisotropy(soma_slice, px_rsz_ratio) + iso_bc_slc, _, _ = correct_image_anisotropy(bc_slc, px_rsz_ratio) # crop isotropized soma slice - iso_soma_slice = crop(iso_soma_slice, rng_out) + iso_bc_slc = crop(iso_bc_slc, out_rng) # mask neuronal bodies - fiber_vec_slice, fiber_clr_slice, frac_anis_slice, soma_msk_slice = \ - mask_background(iso_soma_slice, fiber_vec_slice, fiber_clr_slice, frac_anis_slice, - method=soma_thresh, invert=True) + 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) + else: - soma_msk_slice = None + bc_msk_slc = None # fill memory-mapped output arrays - fill_frangi_volumes(fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, iso_fiber_img, fiber_msk, soma_msk, - fiber_vec_slice, fiber_clr_slice, frac_anis_slice, frangi_slice, iso_fiber_slice, - fiber_msk_slice, soma_msk_slice, rng_out, z_sel) + 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, + fiber_msk_slice, bc_msk_slc, out_rng, z_sel) -def fill_frangi_volumes(fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, iso_fiber_img, fiber_msk, soma_msk, - fiber_vec_slice, fiber_clr_slice, frac_anis_slice, frangi_slice, iso_fiber_slice, - fiber_msk_slice, soma_msk_slice, rng_out, z_sel): +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): """ - Fill the memory-mapped output arrays of the Frangi filtering stage. + Fill the memory-mapped output arrays of the Frangi filter stage. Parameters ---------- - fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vector image - fiber_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) + fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) orientation colormap image - frac_anis_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + 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) - iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) isotropic fiber image - fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) fiber mask image - soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) soma mask image - fiber_vec_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + fbr_vec_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) fiber orientation vector image slice - fiber_clr_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fbr_clr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) orientation colormap image slice - frac_anis_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) fractional anisotropy image slice - frangi_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + frangi_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) Frangi-enhanced image slice - iso_fiber_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + iso_fbr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) isotropic fiber image slice - fiber_msk_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + fbr_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) fiber mask image slice - soma_msk_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + bc_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) soma mask image slice - rng_out: tuple + out_rng: tuple 3D slice output index range z_sel: NumPy slice object @@ -527,36 +520,35 @@ def fill_frangi_volumes(fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, """ # fill memory-mapped output arrays - vec_rng_out = tuple(np.append(rng_out, slice(0, 3, 1))) - fiber_vec_img[vec_rng_out] = fiber_vec_slice[z_sel, ...] - fiber_vec_clr[vec_rng_out] = fiber_clr_slice[z_sel, ...] - iso_fiber_img[rng_out] = iso_fiber_slice[z_sel, ...].astype(np.uint8) - frac_anis_img[rng_out] = (255 * frac_anis_slice[z_sel, ...]).astype(np.uint8) - frangi_img[rng_out] = (255 * frangi_slice[z_sel, ...]).astype(np.uint8) - fiber_msk[rng_out] = (255 * (1 - fiber_msk_slice[z_sel, ...])).astype(np.uint8) + 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) + frangi_img[out_rng] = (255 * frangi_slc[z_sel, ...]).astype(np.uint8) + fa_img[out_rng] = (255 * fa_slc[z_sel, ...]).astype(np.uint8) + fbr_msk[out_rng] = (255 * (1 - fbr_msk_slc[z_sel, ...])).astype(np.uint8) # fill memory-mapped output soma mask, if available - if soma_msk is not None: - soma_msk[rng_out] = (255 * soma_msk_slice[z_sel, ...]).astype(np.uint8) + if bc_msk is not None: + bc_msk[out_rng] = (255 * bc_msk_slc[z_sel, ...]).astype(np.uint8) -def mask_background(img, fiber_vec_slice, fiber_clr_slice, - frac_anis_slice=None, method='yen', invert=False, tissue_msk=None): +def mask_background(img, fbr_vec_slc, fbr_clr_slc, fa_slc=None, method='yen', invert=False, ts_msk=None): """ - Mask orientation volume arrays. + Mask fiber orientation arrays. Parameters ---------- img: numpy.ndarray (axis order=(Z,Y,X)) fiber (or neuron) fluorescence volume image - fiber_vec_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + fbr_vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) fiber orientation vector slice - fiber_clr_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) + fbr_clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) fiber orientation colormap slice - frac_anis_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) fractional anisotropy slice method: str @@ -565,57 +557,57 @@ def mask_background(img, fiber_vec_slice, fiber_clr_slice, invert: bool mask inversion flag - tissue_msk: numpy.ndarray (dtype=bool) + ts_msk: numpy.ndarray (dtype=bool) tissue reconstruction binary mask Returns ------- - fiber_vec_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) + 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) - frac_anis_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=float) + fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float) fractional anisotropy patch (masked) - background_mask: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) + bg: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) background mask """ # generate background mask - background = create_background_mask(img, method=method) + bg = create_background_mask(img, method=method) # apply tissue reconstruction mask, when provided - if tissue_msk is not None: - background = np.logical_or(background, np.logical_not(tissue_msk)) + if ts_msk is not None: + bg = np.logical_or(bg, np.logical_not(ts_msk)) # invert mask if invert: - background = np.logical_not(background) + bg = np.logical_not(bg) # apply mask to input arrays - fiber_vec_slice[background, :] = 0 - fiber_clr_slice[background, :] = 0 - frac_anis_slice[background] = 0 + fbr_vec_slc[bg, :] = 0 + fbr_clr_slc[bg, :] = 0 + fa_slc[bg] = 0 - return fiber_vec_slice, fiber_clr_slice, frac_anis_slice, background + return fbr_vec_slc, fbr_clr_slc, fa_slc, bg -def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, img_name, odf_scale_um, - odf_norm, odf_deg=6, ram=None): +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, ram=None): """ Estimate 3D fiber ODFs from basic orientation data chunks using parallel threads. Parameters ---------- - fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vectors - iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + iso_fbr_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) isotropic fiber volume - px_size_iso: numpy.ndarray (shape=(3,), dtype=float) + px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) adjusted isotropic pixel size [μm] save_dir: str @@ -645,23 +637,23 @@ def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, i """ # derive the ODF kernel size in [px] - odf_scale = int(np.ceil(odf_scale_um / px_size_iso[0])) + 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(fiber_vec_img.shape[:-1], tmp_dir, odf_scale=odf_scale, odf_degrees=odf_deg, ram=ram) + = init_odf_volumes(fbr_vec_img.shape[:-1], tmp_dir, odf_scale=odf_scale, odf_degrees=odf_deg, ram=ram) # generate downsampled background for MRtrix3 mrview - bg_img = fiber_vec_img if iso_fiber_img is None else iso_fiber_img + 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(fiber_vec_img, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, tensor, + 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_size_iso, save_dir, img_name, odf_scale_um) + px_sz_iso, save_dir, img_name, odf_scale_um) def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk=None, @@ -722,10 +714,10 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk # get Frangi filter configuration alpha, beta, gamma, frangi_sigma, frangi_sigma_um, smooth_sigma, px_sz, px_sz_iso, \ - z_rng, ch_lpf, ch_mye, msk_bc, hsv_vec_cmap, img_name = get_frangi_config(cli_args, img_name) + z_rng, bc_ch, fb_ch, msk_bc, hsv_vec_cmap, img_name = get_frangi_config(cli_args, img_name) # get info about the input microscopy image - img_shp, img_shp_um, img_item_sz, ch_mye, msk_bc = get_image_info(img, px_sz, msk_bc, ch_mye, is_tiled=is_tiled) + img_shp, img_shp_um, img_item_sz, fb_ch, msk_bc = get_image_info(img, px_sz, msk_bc, fb_ch, is_tiled=is_tiled) # 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 = \ @@ -733,11 +725,11 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk # 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_bg=msk_bc, jobs=jobs) + generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp, msk_bc=msk_bc, jobs=jobs) # initialize output arrays - fiber_dset_path, fbr_vec_img, fbr_vec_clr, frac_anis_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, mask_lpf=msk_bc, ram=ram) + fbr_dset_path, 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, ram=ram) # print Frangi filter configuration print_frangi_info(alpha, beta, gamma, frangi_sigma_um, img_shp_um, in_slc_shp_um, tot_slc_num, @@ -749,14 +741,14 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk 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, - frac_anis_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, ts_msk=ts_msk, - ch_lpf=ch_lpf, ch_mye=ch_mye, alpha=alpha, beta=beta, gamma=gamma, + fa_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, ts_msk=ts_msk, + bc_ch=bc_ch, fb_ch=fb_ch, alpha=alpha, beta=beta, gamma=gamma, msk_bc=msk_bc, hsv_vec_cmap=hsv_vec_cmap, is_tiled=is_tiled) for i in range(tot_slc_num)) # save output arrays - save_frangi_arrays(fiber_dset_path, fbr_vec_img, fbr_vec_clr, frac_anis_img, frangi_img, fbr_msk, bc_msk, - px_sz_iso, save_dir, img_name) + save_frangi_arrays(fbr_dset_path, fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, fbr_msk, bc_msk, px_sz_iso, + save_dir, img_name) # print Frangi filtering time print_analysis_time(start_time) @@ -764,7 +756,7 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk return fbr_vec_img, iso_fbr_img, px_sz_iso, img_name -def parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_size_iso, save_dir, tmp_dir, img_name, +def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz_iso, save_dir, tmp_dir, img_name, backend='loky', ram=None, verbose=100): """ Iterate over the required spatial scales and apply the parallel ODF analysis @@ -772,16 +764,16 @@ def parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_size_iso, Parameters ---------- - fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vector image - iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + 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_size_iso: numpy.ndarray (axis order=(3,), dtype=float) + px_sz_iso: numpy.ndarray (axis order=(3,), dtype=float) adjusted isotropic pixel size [μm] save_dir: str @@ -820,14 +812,14 @@ def parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_size_iso, num_cpu = get_available_cores() # generate pixel size if not provided - if px_size_iso is None: - px_size_iso = cli_args.px_size_z * np.ones((3,)) + if px_sz_iso is None: + px_sz_iso = cli_args.px_size_z * np.ones((3,)) # parallel ODF analysis of fiber orientation vectors # over the required spatial scales n_jobs = min(num_cpu, len(cli_args.odf_res)) with Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose, max_nbytes=None) as parallel: - parallel(delayed(odf_analysis)(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, img_name, + parallel(delayed(odf_analysis)(fbr_vec_img, iso_fbr_img, px_sz_iso, save_dir, tmp_dir, img_name, odf_norm=odf_norm, odf_deg=cli_args.odf_deg, odf_scale_um=s, ram=ram) for s in cli_args.odf_res) @@ -835,36 +827,36 @@ def parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_size_iso, print_analysis_time(start_time) -def save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, - fiber_msk, soma_msk, px_size, save_dir, img_name): +def save_frangi_arrays(fbr_dset_path, fbr_vec_img, fbr_vec_clr, fa_img, frangi_img, + fbr_msk, bc_msk, px_sz, save_dir, img_name): """ Save the output arrays of the Frangi filter stage to TIF files. Parameters ---------- - fiber_dset_path: str + fbr_dset_path: str path to initialized fiber orientation HDF5 dataset (not fitting the available RAM) - fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) + fbr_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) fiber orientation vector image - fiber_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) + fbr_vec_clr: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) orientation colormap image - frac_anis_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + 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 volume image (fiber probability) - fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + fbr_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) fiber mask image - soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + bc_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) neuron mask image - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size (Z,Y,X) [μm] save_dir: str @@ -879,30 +871,30 @@ def save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_ """ # move large fiber orientation dataset to saving directory - if fiber_dset_path is not None: - shutil.move(fiber_dset_path, path.join(save_dir, 'fiber_vec_{0}.h5'.format(img_name))) + if fbr_dset_path is not None: + shutil.move(fbr_dset_path, path.join(save_dir, 'fiber_vec_{0}.h5'.format(img_name))) # or save orientation vectors to TIFF else: - save_array('fiber_vec_{0}'.format(img_name), save_dir, np.moveaxis(fiber_vec_img, -1, 1), px_size) + save_array('fiber_vec_{0}'.format(img_name), save_dir, np.moveaxis(fbr_vec_img, -1, 1), px_sz) # save orientation color map to TIFF - save_array('fiber_cmap_{0}'.format(img_name), save_dir, fiber_vec_clr, px_size) + save_array('fiber_cmap_{0}'.format(img_name), save_dir, fbr_vec_clr, px_sz) # save fractional anisotropy map to TIFF - save_array('frac_anis_{0}'.format(img_name), save_dir, frac_anis_img, px_size) + save_array('frac_anis_{0}'.format(img_name), save_dir, fa_img, px_sz) # save Frangi-enhanced fiber volume to TIFF - save_array('frangi_{0}'.format(img_name), save_dir, frangi_img, px_size) + save_array('frangi_{0}'.format(img_name), save_dir, frangi_img, px_sz) # save masked fiber volume to TIFF - save_array('fiber_msk_{0}'.format(img_name), save_dir, fiber_msk, px_size) + save_array('fiber_msk_{0}'.format(img_name), save_dir, fbr_msk, px_sz) # save masked soma volume to TIFF - if soma_msk is not None: - save_array('soma_msk_{0}'.format(img_name), save_dir, soma_msk, px_size) + if bc_msk is not None: + save_array('soma_msk_{0}'.format(img_name), save_dir, bc_msk, px_sz) -def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, disarray, px_size, save_dir, img_name, odf_scale_um): +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 @@ -931,7 +923,7 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, disarray, px_s disarray: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32) - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size (Z,Y,X) [μm] save_dir: str @@ -950,8 +942,8 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, disarray, px_s # ODF analysis volumes to Nifti files (adjusted view for MRtrix3) save_array('bg_mrtrixview_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, bg, fmt='nii') save_array('odf_mrtrixview_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odf, fmt='nii') - save_array('odi_pri_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_pri, px_size, odi=True) - save_array('odi_sec_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_sec, px_size, odi=True) - save_array('odi_tot_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_tot, px_size, odi=True) - save_array('odi_anis_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_anis, px_size, odi=True) - save_array('disarray_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, disarray, px_size, odi=True) + save_array('odi_pri_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_pri, px_sz, odi=True) + save_array('odi_sec_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_sec, px_sz, odi=True) + save_array('odi_tot_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_tot, px_sz, odi=True) + save_array('odi_anis_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_anis, px_sz, odi=True) + save_array('disarray_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, disarray, px_sz, odi=True) diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py index 30d610b..aa7e5df 100644 --- a/foa3d/preprocessing.py +++ b/foa3d/preprocessing.py @@ -49,10 +49,10 @@ def config_anisotropy_correction(px_sz, psf_fwhm): # adjust PSF anisotropy via lateral Gaussian blurring if cndt_1: - # estimate the PSF variance from input FWHM values [μm^2] + # estimate the PSF variance from input FWHM values [μm²] psf_var = np.square(fwhm_to_sigma(psf_fwhm)) - # estimate the in-plane filter variance [μm^2] + # estimate the in-plane filter variance [μm²] gauss_var = np.max(psf_var) - psf_var # ...and the corresponding standard deviation [px] @@ -72,23 +72,22 @@ def config_anisotropy_correction(px_sz, psf_fwhm): return smooth_sigma, px_sz_iso -def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mode='reflect', +def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, pad_mode='reflect', anti_alias=True, trunc=4, ts_msk=None): """ - Smooth the input volume image along the X and Y axes so that the lateral - and longitudinal sizes of the optical system's PSF become equal. - Downsample data in the XY plane in order to uniform the 3D pixel size. + Smooth and downsample the raw microscopy image so as to uniform the lateral sizes + of the optical system's PSF and the original voxel size. Parameters ---------- img: numpy.ndarray (axis order=(Z,Y,X)) - microscopy volume image + 3D microscopy image rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) 3D resize ratio sigma: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the low-pass Gaussian filter [px] + 3D standard deviation of the smoothing Gaussian filter [px] (resolution anisotropy correction) pad: numpy.ndarray (shape=(3,2), dtype=int) @@ -98,57 +97,57 @@ def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mo image padding mode adopted for the smoothing Gaussian filter anti_alias: bool - if True, apply an antialiasing filter when downsampling the XY plane + if True, apply an antialiasing filter when downsampling trunc: int truncate the Gaussian kernel at this many standard deviations tissue_msk: numpy.ndarray (dtype=bool) - tissue reconstruction binary mask + tissue binary mask Returns ------- iso_img: numpy.ndarray (axis order=(Z,Y,X)) - isotropic microscopy volume image + isotropic microscopy image - rsz_pad: numpy.ndarray (shape=(3,2), dtype=int) + pad_rsz: numpy.ndarray (shape=(3,2), dtype=int) resized image padding array - rsz_tissue_msk: numpy.ndarray (dtype=bool) - resized tissue reconstruction binary mask + ts_msk_rsz: numpy.ndarray (dtype=bool) + resized tissue binary mask """ # no resizing if np.all(rsz_ratio == 1): # resize tissue mask, when available if ts_msk is not None: - tissue_msk = ts_msk[np.newaxis, ...] - tissue_msk = np.repeat(tissue_msk, img.shape[0], axis=0) + ts_msk = ts_msk[np.newaxis, ...] + ts_msk = np.repeat(ts_msk, img.shape[0], axis=0) return img, pad, ts_msk # lateral blurring else: if sigma is not None: - img = gaussian_filter(img, sigma=sigma, mode=smooth_pad_mode, truncate=trunc, output=np.float32) + img = gaussian_filter(img, sigma=sigma, mode=pad_mode, truncate=trunc, output=np.float32) # downsampling - iso_shape = np.ceil(np.multiply(np.asarray(img.shape), rsz_ratio)).astype(int) - iso_img = np.zeros(shape=iso_shape, dtype=img.dtype) - for z in range(iso_shape[0]): + iso_shp = np.ceil(np.multiply(np.asarray(img.shape), rsz_ratio)).astype(int) + iso_img = np.zeros(shape=iso_shp, dtype=img.dtype) + for z in range(iso_shp[0]): iso_img[z, ...] = \ - resize(img[z, ...], output_shape=tuple(iso_shape[1:]), anti_aliasing=anti_alias, preserve_range=True) + resize(img[z, ...], output_shape=tuple(iso_shp[1:]), anti_aliasing=anti_alias, preserve_range=True) # resize padding array accordingly - rsz_pad = np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad)).astype(int) \ + pad_rsz = np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad)).astype(int) \ if pad is not None else None # resize tissue mask, when available if ts_msk is not None: - rsz_tissue_msk = resize(ts_msk, output_shape=tuple(iso_shape[1:]), preserve_range=True) - rsz_tissue_msk = rsz_tissue_msk[np.newaxis, ...] - rsz_tissue_msk = np.repeat(rsz_tissue_msk, iso_shape[0], axis=0) + ts_msk_rsz = resize(ts_msk, output_shape=tuple(iso_shp[1:]), preserve_range=True) + ts_msk_rsz = ts_msk_rsz[np.newaxis, ...] + ts_msk_rsz = np.repeat(ts_msk_rsz, iso_shp[0], axis=0) else: - rsz_tissue_msk = None + ts_msk_rsz = None - return iso_img, rsz_pad, rsz_tissue_msk + return iso_img, pad_rsz, ts_msk_rsz diff --git a/foa3d/printing.py b/foa3d/printing.py index e3f33db..871c32a 100644 --- a/foa3d/printing.py +++ b/foa3d/printing.py @@ -39,10 +39,10 @@ def color_text(r, g, b, text): return clr_text -def print_frangi_info(alpha, beta, gamma, scales_um, image_shape_um, in_slice_shape_um, tot_slice_num, - px_size, image_item_size, lpf_soma_mask): +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): """ - Print Frangi filter heading. + Print Frangi filter stage heading. Parameters ---------- @@ -58,23 +58,24 @@ def print_frangi_info(alpha, beta, gamma, scales_um, image_shape_um, in_slice_sh scales_um: list (dtype=float) analyzed spatial scales [μm] - image_shape_um: numpy.ndarray (shape=(3,), dtype=float) + img_shp_um: numpy.ndarray (shape=(3,), dtype=float) volume image shape [μm] - in_slice_shape_um: numpy.ndarray (shape=(3,), dtype=float) + in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=float) shape of the analyzed image slices [μm] - tot_slice_num: int + tot_slc_num: int total number of analyzed image slices - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - image_item_size: int + img_item_sz: int image item size (in bytes) - lpf_soma_mask: bool - neuronal body masking flag + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel Returns ------- @@ -93,15 +94,15 @@ def print_frangi_info(alpha, beta, gamma, scales_um, image_shape_um, in_slice_sh print("Enhanced diameters [\u03BCm]: {}\n".format(4 * scales_um)) # print iterative analysis information - print_slicing_info(image_shape_um, in_slice_shape_um, tot_slice_num, px_size, image_item_size) + print_slicing_info(img_shp_um, in_slc_shp_um, tot_slc_num, px_sz, img_item_sz) # print neuron masking info - print_soma_masking(lpf_soma_mask) + print_soma_masking(msk_bc) def print_analysis_time(start_time): """ - Print volume image analysis time. + Print image processing time. Parameters ---------- @@ -112,18 +113,19 @@ def print_analysis_time(start_time): ------- None """ - _, mins, secs = elapsed_time(start_time) - print("\nMicroscopy image analyzed in: {0} min {1:3.1f} s\n".format(mins, secs)) + _, hrs, mins, secs = elapsed_time(start_time) + print("\nMicroscopy image analyzed in: {0} min {1:3.1f} s\n".format(mins, secs)) if hrs == 0 else \ + print("\nMicroscopy image analyzed in: {0} hrs {1} min {2:3.1f} s\n".format(hrs, mins, secs)) def print_blur(smooth_sigma_um, psf_fwhm): """ - Print gaussian lowpass filter standard deviation. + Print the standard deviation of the smoothing Gaussian filter. Parameters ---------- smooth_sigma_um: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the low-pass Gaussian filter [μm] + 3D standard deviation of the smoothing Gaussian filter [μm] (resolution anisotropy correction) psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) @@ -140,7 +142,7 @@ def print_blur(smooth_sigma_um, psf_fwhm): def print_import_time(start_time): """ - Print volume image import time. + Print image import time. Parameters ---------- @@ -151,7 +153,7 @@ def print_import_time(start_time): ------- None """ - _, mins, secs = elapsed_time(start_time) + _, _, mins, secs = elapsed_time(start_time) print("Image loaded in: {0} min {1:3.1f} s".format(mins, secs)) @@ -178,7 +180,7 @@ def print_odf_info(odf_scales_um, odf_degrees): def print_pipeline_heading(): """ - Print Foa3D pipeline heading. + Print Foa3D tool heading. Returns ------- @@ -199,13 +201,13 @@ def print_prepro_heading(): print("\n Z Y X") -def print_native_res(px_size, psf_fwhm): +def print_native_res(px_sz, psf_fwhm): """ - Print pixel and optical resolution of the microscopy system. + Print the native pixel size and optical resolution of the microscopy system. Parameters ---------- - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) @@ -215,13 +217,13 @@ def print_native_res(px_size, psf_fwhm): ------- None """ - print("Pixel size [μm]: ({0:.3f}, {1:.3f}, {2:.3f})".format(px_size[0], px_size[1], px_size[2])) + print("Pixel size [μm]: ({0:.3f}, {1:.3f}, {2:.3f})".format(px_sz[0], px_sz[1], px_sz[2])) print("PSF FWHM [μm]: ({0:.3f}, {1:.3f}, {2:.3f})".format(psf_fwhm[0], psf_fwhm[1], psf_fwhm[2])) def print_new_res(px_sz_iso): """ - Print adjusted isotropic spatial resolution. + Print the adjusted isotropic spatial resolution. Parameters ---------- @@ -235,26 +237,25 @@ def print_new_res(px_sz_iso): print("Adjusted pixel size [μm]: ({0:.3f}, {1:.3f}, {2:.3f})\n".format(px_sz_iso[0], px_sz_iso[1], px_sz_iso[2])) -def print_slicing_info(image_shape_um, slice_shape_um, tot_slice_num, px_size, image_item_size): +def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz): """ - Print information on the slicing of the basic image sub-volumes - iteratively processed by the Foa3D pipeline. + Print information on the slicing of the basic image sub-volumes processed by the Foa3D tool. Parameters ---------- - image_shape_um: numpy.ndarray (shape=(3,), dtype=float) - volume image shape [μm] + img_shp_um: numpy.ndarray (shape=(3,), dtype=float) + 3D microscopy image [μm] - slice_shape_um: numpy.ndarray (shape=(3,), dtype=float) + slc_shp_um: numpy.ndarray (shape=(3,), dtype=float) shape of the analyzed image slices [μm] - tot_slice_num: int + tot_slc_num: int total number of analyzed image slices - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - image_item_size: int + img_item_sz: int image item size (in bytes) Returns @@ -263,49 +264,45 @@ def print_slicing_info(image_shape_um, slice_shape_um, tot_slice_num, px_size, i """ # adjust slice shape - if np.any(image_shape_um < slice_shape_um): - slice_shape_um = image_shape_um + if np.any(img_shp_um < slc_shp_um): + slc_shp_um = img_shp_um # get image memory size - image_size = image_item_size * np.prod(np.divide(image_shape_um, px_size)) + img_sz = img_item_sz * np.prod(np.divide(img_shp_um, px_sz)) # get slice memory size - max_slice_size = image_item_size * np.prod(np.divide(slice_shape_um, px_size)) + max_slc_sz = img_item_sz * np.prod(np.divide(slc_shp_um, px_sz)) # print info print("\n Z Y X") - print("Total image shape [μm]: ({0:.1f}, {1:.1f}, {2:.1f})" - .format(image_shape_um[0], image_shape_um[1], image_shape_um[2])) - print("Total image size [MB]: {0}\n" - .format(np.ceil(image_size / 1024**2).astype(int))) - print("Image slice shape [μm]: ({0:.1f}, {1:.1f}, {2:.1f})" - .format(slice_shape_um[0], slice_shape_um[1], slice_shape_um[2])) - print("Image slice size [MB]: {0}" - .format(np.ceil(max_slice_size / 1024**2).astype(int))) - print("Image slice number: {0}\n" - .format(tot_slice_num)) - - -def print_soma_masking(lpf_soma_mask): + print("Total image shape [μm]: ({0:.1f}, {1:.1f}, {2:.1f})".format(img_shp_um[0], img_shp_um[1], img_shp_um[2])) + print("Total image size [MB]: {0}\n".format(np.ceil(img_sz / 1024**2).astype(int))) + print("Image slice shape [μm]: ({0:.1f}, {1:.1f}, {2:.1f})".format(slc_shp_um[0], slc_shp_um[1], slc_shp_um[2])) + print("Image slice size [MB]: {0}".format(np.ceil(max_slc_sz / 1024**2).astype(int))) + print("Image slice number: {0}\n".format(tot_slc_num)) + + +def print_soma_masking(msk_bc): """ - Print info on lipofuscin-based neuronal body masking. + Print information on the optional masking of neuronal bodies. Parameters ---------- - lpf_soma_mask: bool - neuronal body masking flag + msk_bc: bool + if True, mask neuronal bodies within + the optionally provided channel Returns ------- None """ prt = 'Soma mask: ' - print('{}active\n'.format(prt)) if lpf_soma_mask else print('{}not active\n'.format(prt)) + print('{}active\n'.format(prt)) if msk_bc else print('{}not active\n'.format(prt)) -def print_image_shape(cli_args, img, is_tiled, channel_ax=None): +def print_image_shape(cli_args, img, is_tiled, ch_ax=None): """ - Print volume image shape. + Print 3D microscopy image shape. Parameters ---------- @@ -318,7 +315,7 @@ def print_image_shape(cli_args, img, is_tiled, channel_ax=None): is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher - channel_ax: int + ch_ax: int channel axis (if ndim == 4) Returns @@ -327,18 +324,18 @@ def print_image_shape(cli_args, img, is_tiled, channel_ax=None): """ # get pixel size - px_size_z = cli_args.px_size_z - px_size_xy = cli_args.px_size_xy + px_sz_z = cli_args.px_size_z + px_sz_xy = cli_args.px_size_xy # adapt axis order - img_shape = img.shape - if len(img_shape) == 4: - channel_ax = 1 if is_tiled else -1 + img_shp = img.shape + if len(img_shp) == 4: + ch_ax = 1 if is_tiled else -1 # get image shape - if channel_ax is not None: - img_shape = np.delete(img_shape, channel_ax) + if ch_ax is not None: + img_shp = np.delete(img_shp, ch_ax) print("\n Z Y X") print("Image shape [μm]: ({0:.1f}, {1:.1f}, {2:.1f})" - .format(img_shape[0] * px_size_z, img_shape[1] * px_size_xy, img_shape[2] * px_size_xy)) + .format(img_shp[0] * px_sz_z, img_shp[1] * px_sz_xy, img_shp[2] * px_sz_xy)) diff --git a/foa3d/slicing.py b/foa3d/slicing.py index 248fe84..8d6aad5 100644 --- a/foa3d/slicing.py +++ b/foa3d/slicing.py @@ -387,7 +387,7 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi return batch_sz, in_slc_shp, in_slc_shp_um, px_rsz_ratio, ovlp, ovlp_rsz -def crop(img, rng, ovlp, flip=()): +def crop(img, rng, ovlp=np.zeros((3,)), flip=()): """ Shrink image slice at total volume boundaries, for overall shape consistency. @@ -459,7 +459,7 @@ def crop_lst(img_lst, rng, ovlp=None, flip=()): return img_lst -def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bg=False, jobs=None): +def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bc=False, jobs=None): """ Generate image slice ranges for the Frangi filtering stage. @@ -481,7 +481,7 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms ovlp: int overlapping range between image slices along each axis [px] - msk_bg: bool + msk_bc: bool if True, mask neuronal bodies in the optionally provided image channel @@ -540,7 +540,7 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms out_rng_lst.append(out_rng) # (optional) neuronal body channel - if msk_bg: + if msk_bc: bc_rng, _ = compute_slice_range((z, y, x), in_slc_shp, img_shp, slc_per_dim) bc_rng_lst.append(bc_rng) else: diff --git a/foa3d/utils.py b/foa3d/utils.py index c5440f6..a763026 100644 --- a/foa3d/utils.py +++ b/foa3d/utils.py @@ -38,8 +38,8 @@ def create_background_mask(img, method='yen', black_bg=False): # select thresholding method if method == 'li': - initial_li_guess = np.mean(img[img != 0]) - thresh = threshold_li(img, initial_guess=initial_li_guess) + init_li = np.mean(img[img != 0]) + thresh = threshold_li(img, initial_guess=init_li) elif method == 'niblack': thresh = threshold_niblack(img, window_size=15, k=0.2) elif method == 'sauvola': @@ -171,7 +171,7 @@ def get_item_size(dtype): Returns ------- - item_size: int + item_sz: int item size in bytes """ @@ -182,17 +182,17 @@ def get_item_size(dtype): lst_4 = ['uint64', 'int64', 'float64', np.float64] if dtype in lst_1: - item_size = 1 + item_sz = 1 elif dtype in lst_2: - item_size = 2 + item_sz = 2 elif dtype in lst_3: - item_size = 4 + item_sz = 4 elif dtype in lst_4: - item_size = 8 + item_sz = 8 else: raise ValueError("Unsupported data type!") - return item_size + return item_sz def delete_tmp_folder(tmp_dir): @@ -214,7 +214,7 @@ def delete_tmp_folder(tmp_dir): pass -def divide_nonzero(nd_array1, nd_array2, new_value=1e-10): +def divide_nonzero(nd_array1, nd_array2, new_val=1e-10): """ Divide two arrays handling zero denominator values. @@ -226,7 +226,7 @@ def divide_nonzero(nd_array1, nd_array2, new_value=1e-10): nd_array2: numpy.ndarray divisor array - new_value: float + new_val: float substituted value Returns @@ -235,9 +235,9 @@ def divide_nonzero(nd_array1, nd_array2, new_value=1e-10): divided array """ - denominator = np.copy(nd_array2) - denominator[denominator == 0] = new_value - divided = np.divide(nd_array1, denominator) + divisor = np.copy(nd_array2) + divisor[divisor == 0] = new_val + divided = np.divide(nd_array1, divisor) return divided @@ -253,10 +253,13 @@ def elapsed_time(start_time): Returns ------- - total: float + tot: float total time [s] - mins: float + hrs: int + hours + + mins: int minutes secs: float @@ -264,11 +267,15 @@ def elapsed_time(start_time): """ stop_time = perf_counter() - total = stop_time - start_time - mins = total // 60 - secs = total % 60 + tot = stop_time - start_time + + secs = tot % 86400 + hrs = int(secs // 3600) + secs %= 3600 + mins = int(secs // 60) + secs %= 60 - return total, mins, secs + return tot, hrs, mins, secs def fwhm_to_sigma(fwhm):