diff --git a/foa3d/__main__.py b/foa3d/__main__.py index 563aea3..eb892ad 100644 --- a/foa3d/__main__.py +++ b/foa3d/__main__.py @@ -1,5 +1,5 @@ -from foa3d.input import get_cli_parser, get_pipeline_config, load_microscopy_image -from foa3d.pipeline import (parallel_odf_on_scales, parallel_frangi_on_slices) +from foa3d.input import get_cli_parser, get_resource_config, load_microscopy_image +from foa3d.pipeline import (parallel_odf_at_scales, parallel_frangi_on_slices) from foa3d.printing import print_pipeline_heading from foa3d.utils import delete_tmp_folder @@ -7,29 +7,22 @@ def foa3d(cli_args): # load microscopy volume image or array of fiber orientation vectors - img, mosaic, skip_frangi, cli_args, save_subdirs, tmp_dir, img_name = load_microscopy_image(cli_args) + img, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args) - # get the fiber orientation analysis pipeline configuration - alpha, beta, gamma, scales_um, smooth_sigma, px_size, px_size_iso, \ - odf_scales_um, odf_degrees, z_min, z_max, ch_neuron, ch_fiber, \ - lpf_soma_mask, max_ram_mb, jobs, img_name = get_pipeline_config(cli_args, skip_frangi, img_name) + # get resources configuration + ram, jobs = get_resource_config(cli_args) # conduct parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices - if not skip_frangi: - fiber_vec_img, iso_fiber_img \ - = parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_subdirs[0], tmp_dir, img_name, - alpha=alpha, beta=beta, gamma=gamma, frangi_sigma_um=scales_um, - z_min=z_min, z_max=z_max, lpf_soma_mask=lpf_soma_mask, - ch_neuron=ch_neuron, ch_fiber=ch_fiber, mosaic=mosaic, - max_ram_mb=max_ram_mb, jobs=jobs) - - # estimate 3D fiber ODF maps over the spatial scales of interest using concurrent workers - if odf_scales_um: - if skip_frangi: - fiber_vec_img = img - iso_fiber_img = None - parallel_odf_on_scales(fiber_vec_img, iso_fiber_img, px_size_iso, save_subdirs[1], tmp_dir, img_name, - odf_scales_um=odf_scales_um, odf_degrees=odf_degrees, max_ram_mb=max_ram_mb) + if not is_fiber: + fiber_vec_img, iso_fiber_img, px_sz, img_name \ + = parallel_frangi_on_slices(img, cli_args, save_dir[0], tmp_dir, img_name, ram=ram, jobs=jobs, + is_tiled=is_tiled) + else: + fiber_vec_img, iso_fiber_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) # delete temporary folder delete_tmp_folder(tmp_dir) diff --git a/foa3d/frangi.py b/foa3d/frangi.py index 9ca7780..f620780 100644 --- a/foa3d/frangi.py +++ b/foa3d/frangi.py @@ -141,9 +141,6 @@ def compute_scaled_hessian(img, sigma=1, trunc=4): Hessian matrix of image second derivatives """ - # get number of dimensions - ndim = img.ndim - # scale selection scaled_img = ndi.gaussian_filter(img, sigma=sigma, output=np.float32, truncate=trunc) @@ -152,15 +149,15 @@ def compute_scaled_hessian(img, sigma=1, trunc=4): # compute the Hessian matrix elements hessian_elements = [np.gradient(gradient_list[ax0], axis=ax1) - for ax0, ax1 in combinations_with_replacement(range(ndim), 2)] + for ax0, ax1 in combinations_with_replacement(range(img.ndim), 2)] # scale the elements of the Hessian matrix corr_factor = sigma ** 2 hessian_elements = [corr_factor * element for element in hessian_elements] # create the Hessian matrix from its basic elements - hessian = np.zeros((ndim, ndim) + scaled_img.shape, dtype=scaled_img.dtype) - for index, (ax0, ax1) in enumerate(combinations_with_replacement(range(ndim), 2)): + hessian = np.zeros((img.ndim, img.ndim) + scaled_img.shape, dtype=scaled_img.dtype) + for index, (ax0, ax1) in enumerate(combinations_with_replacement(range(img.ndim), 2)): element = hessian_elements[index] hessian[ax0, ax1, ...] = element if ax0 != ax1: @@ -264,30 +261,6 @@ def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma): return vesselness -def convert_frangi_scales(scales_um, px_size): - """ - Compute the Frangi filter scales in pixel. - - Parameters - ---------- - scales_um: list (dtype=float) - Frangi filter scales [μm] - - px_size: int - isotropic pixel size [μm] - - Returns - ------- - scales_px: numpy.ndarray (dtype=int) - Frangi filter scales [px] - """ - - scales_um = np.asarray(scales_um) - scales_px = scales_um / px_size - - return scales_px - - def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True): """ Apply 3D Frangi filter to microscopy volume image. @@ -329,7 +302,7 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True n_scales = len(scales_px) if n_scales == 1: enhanced_img, fiber_vec, eigenval \ - = compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark) + = compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark) # parallel scaled vesselness analysis else: @@ -381,17 +354,13 @@ def reject_vesselness_background(vesselness, eigen2, eigen3, dark): Returns ------- - vesselness: numpy.ndarray (axis order=(Z,Y,X)) + vesselness: numpy.ndarray (axis order=(Z,Y,X), dtype=float) masked Frangi's vesselness likelihood image """ - if dark: - vesselness[eigen2 < 0] = 0 - vesselness[eigen3 < 0] = 0 - else: - vesselness[eigen2 > 0] = 0 - vesselness[eigen3 > 0] = 0 - vesselness[np.isnan(vesselness)] = 0 + 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 return vesselness diff --git a/foa3d/input.py b/foa3d/input.py index 14b345c..8af4622 100644 --- a/foa3d/input.py +++ b/foa3d/input.py @@ -16,7 +16,7 @@ from foa3d.output import create_save_dirs from foa3d.preprocessing import config_anisotropy_correction from foa3d.printing import (color_text, print_image_shape, print_import_time, - print_resolution) + print_native_res) from foa3d.utils import create_memory_map, get_item_bytes, get_output_prefix @@ -64,8 +64,8 @@ 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('-n', '--neuron-mask', action='store_true', default=False, - help='lipofuscin-based neuronal body masking') + 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') @@ -78,10 +78,12 @@ def get_cli_parser(): 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-fiber', type=int, default=1, help='myelinated fibers channel') - cli_parser.add_argument('--ch-neuron', type=int, default=0, help='neuronal soma channel') + 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('--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='toggle HSV colormap for 3D fiber orientations') cli_parser.add_argument('-o', '--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, @@ -93,7 +95,7 @@ def get_cli_parser(): return cli_args -def get_image_info(img, px_size, ch_fiber, mosaic=False, ch_axis=None): +def get_image_info(img, px_sz, mask_lpf, ch_mye, ch_axis=None, is_tiled=False): """ Get information on the input microscopy volume image. @@ -102,18 +104,22 @@ def get_image_info(img, px_size, ch_fiber, mosaic=False, ch_axis=None): img: numpy.ndarray microscopy volume image - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - ch_fiber: int - myelinated fibers channel + mask_lpf: bool + if True, mask neuronal bodies exploiting the autofluorescence + signal of lipofuscin pigments - mosaic: bool - True for tiled reconstructions aligned using ZetaStitcher + ch_mye: int + myelinated fibers channel ch_axis: int channel axis + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher + Returns ------- img_shape: numpy.ndarray (shape=(3,), dtype=int) @@ -122,28 +128,36 @@ def get_image_info(img, px_size, ch_fiber, mosaic=False, ch_axis=None): img_shape_um: numpy.ndarray (shape=(3,), dtype=float) volume image shape [μm] - img_item_size: int + img_item_sz: int array item size (in bytes) + + ch_mye: int + myelinated fibers channel + + mask_lpf: bool + if True, mask neuronal bodies exploiting the autofluorescence + signal of lipofuscin pigments """ # adapt channel axis img_shape = np.asarray(img.shape) ndim = len(img_shape) if ndim == 4: - ch_axis = 1 if mosaic else -1 + ch_axis = 1 if is_tiled else -1 elif ndim == 3: - ch_fiber = None + ch_mye = None + mask_lpf = 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_size) - img_item_size = get_item_bytes(img) + img_shape_um = np.multiply(img_shape, px_sz) + img_item_sz = get_item_bytes(img) - return img_shape, img_shape_um, img_item_size, ch_fiber + return img_shape, img_shape_um, img_item_sz, ch_mye, mask_lpf -def get_image_file(cli_args): +def get_file_info(cli_args): """ Get microscopy image file path and format. @@ -163,17 +177,17 @@ def get_image_file(cli_args): img_fmt: str format of the microscopy volume image - mosaic: bool + is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher - in_mmap: bool + is_mmap: bool create a memory-mapped array of the microscopy volume image, increasing the parallel processing performance (the image will be preliminarily loaded to RAM) """ # get microscopy image path and name - in_mmap = cli_args.mmap + is_mmap = cli_args.mmap img_path = cli_args.image_path img_name = path.basename(img_path) split_name = img_name.split('.') @@ -183,26 +197,23 @@ def get_image_file(cli_args): raise ValueError('Format must be specified for input volume images!') else: img_fmt = split_name[-1] - img_name = img_name.replace('.' + split_name[-1], '') - mosaic = True if img_fmt == 'yml' else False + img_name = img_name.replace('.{}'.format(img_fmt), '') + is_tiled = True if img_fmt == 'yml' else False - return img_path, img_name, img_fmt, mosaic, in_mmap + return img_path, img_name, img_fmt, is_tiled, is_mmap -def get_pipeline_config(cli_args, vector, img_name): +def get_frangi_config(cli_args, img_name): """ - Retrieve the Foa3D pipeline configuration. + Get Frangi filter configuration. Parameters ---------- cli_args: see ArgumentParser.parse_args populated namespace of command line arguments - vector: bool - True for fiber orientation vector data - img_name: str - name of the input volume image + name of the microscopy volume image Returns ------- @@ -215,91 +226,74 @@ def get_pipeline_config(cli_args, vector, img_name): 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_sigma: numpy.ndarray (shape=(3,), dtype=int) 3D standard deviation of the low-pass Gaussian filter [px] - (resolution anisotropy correction) + (applied to the XY plane) - px_size: numpy.ndarray (shape=(3,), dtype=float) + px_sz: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - px_size_iso: numpy.ndarray (shape=(3,), dtype=float) - adjusted isotropic pixel size [μm] - - odf_scales_um: list (dtype: float) - list of fiber ODF resolution values (super-voxel sides [μm]) + px_sz_iso: int + isotropic pixel size [μm] - odf_degrees: int - degrees of the spherical harmonics series expansion + z_rng: int + output z-range in [px] - z_min: int - minimum output z-depth [px] - - z_max: int - maximum output z-depth [px] - - ch_neuron: int + ch_lpf: int neuronal bodies channel - ch_fiber: int + ch_mye: int myelinated fibers channel - lpf_soma_mask: bool - neuronal body masking + mask_lpf: bool + if True, mask neuronal bodies exploiting the autofluorescence + signal of lipofuscin pigments - max_ram_mb: float - maximum RAM available to the Frangi filtering stage [MB] + hsv_vec_cmap: bool - jobs: int - number of parallel jobs (threads) - used by the Frangi filtering stage - - img_name: str - microscopy image filename + out_name: str + output file name """ + # microscopy image pixel and PSF size + px_sz, psf_fwhm = get_resolution(cli_args) + + # preprocessing configuration (in-plane smoothing) + smooth_sigma, px_sz_iso = config_anisotropy_correction(px_sz, psf_fwhm) + # Frangi filter parameters - alpha = cli_args.alpha - beta = cli_args.beta - gamma = cli_args.gamma - scales_um = cli_args.scales - if type(scales_um) is not list: - scales_um = [scales_um] - - # pipeline parameters - lpf_soma_mask = cli_args.neuron_mask - ch_neuron = cli_args.ch_neuron - ch_fiber = cli_args.ch_fiber - max_ram = cli_args.ram - max_ram_mb = None if max_ram is None else max_ram * 1000 - jobs = cli_args.jobs + 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] - # ODF parameters - odf_scales_um = cli_args.odf_res - odf_degrees = cli_args.odf_deg + # image channels + ch_mye = cli_args.ch_mye + ch_lpf = cli_args.ch_lpf + mask_lpf = cli_args.lpf_mask - # microscopy image pixel size and PSF FWHM - px_size, psf_fwhm = get_resolution(cli_args, vector) + # fiber orientation colormap + hsv_vec_cmap = cli_args.hsv # forced output z-range - z_min = cli_args.z_min - z_max = cli_args.z_max - z_min = int(np.floor(z_min / px_size[0])) - if z_max is not None: - z_max = int(np.ceil(z_max / px_size[0])) + 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) - # preprocessing configuration - smooth_sigma, px_size_iso = config_anisotropy_correction(px_size, psf_fwhm, vector) + # add Frangi filter configuration prefix to output filenames + pfx = get_output_prefix(scales_um, alpha, beta, gamma) + out_name = '{}img{}'.format(pfx, img_name) - # add pipeline configuration prefix to output filenames - if not vector: - pfx = get_output_prefix(scales_um, alpha, beta, gamma) - img_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 - return alpha, beta, gamma, scales_um, smooth_sigma, px_size, px_size_iso, odf_scales_um, odf_degrees, \ - z_min, z_max, ch_neuron, ch_fiber, lpf_soma_mask, max_ram_mb, jobs, img_name - -def get_resolution(cli_args, vector): +def get_resolution(cli_args): """ Retrieve microscopy resolution information from command line arguments. @@ -308,42 +302,57 @@ def get_resolution(cli_args, vector): cli_args: see ArgumentParser.parse_args populated namespace of command line arguments - vector: bool - True for fiber orientation vector data - Returns ------- - 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) 3D PSF FWHM [μm] """ - # pixel size - px_size_z = cli_args.px_size_z - px_size_xy = cli_args.px_size_xy + # 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]) - # psf size - psf_fwhm_z = cli_args.psf_fwhm_z - psf_fwhm_y = cli_args.psf_fwhm_y - psf_fwhm_x = cli_args.psf_fwhm_x + # print resolution info + print_native_res(px_sz, psf_fwhm) - # create arrays - px_size = np.array([px_size_z, px_size_xy, px_size_xy]) - psf_fwhm = np.array([psf_fwhm_z, psf_fwhm_y, psf_fwhm_x]) + return px_sz, psf_fwhm - # print resolution info - if not vector: - print_resolution(px_size, psf_fwhm) - return px_size, psf_fwhm +def get_resource_config(cli_args): + """ + Retrieve resource usage configuration of the Foa3D pipeline. + + Parameters + ---------- + cli_args: see ArgumentParser.parse_args + populated namespace of command line arguments + + Returns + ------- + max_ram: float + maximum RAM available to the Frangi filtering stage [B] + + jobs: int + number of parallel jobs (threads) + used by the Frangi filtering stage + """ + + # resource parameters (convert maximum RAM to bytes) + jobs = cli_args.jobs + max_ram = cli_args.ram + if max_ram is not None: + max_ram *= 1024**3 + + return max_ram, jobs def load_microscopy_image(cli_args): """ Load microscopy volume image from TIFF, NumPy or ZetaStitcher .yml file. - Alternatively, the processing pipeline accepts as input .npy or HDF5 + 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. @@ -357,17 +366,14 @@ def load_microscopy_image(cli_args): img: numpy.ndarray or NumPy memory-map object microscopy volume image or array of fiber orientation vectors - mosaic: bool + is_tiled: bool True for tiled microscopy reconstructions aligned using ZetaStitcher - skip_frangi: bool + is_fiber: bool True when pre-estimated fiber orientation vectors are directly provided to the pipeline - cli_args: see ArgumentParser.parse_args - updated namespace of command line arguments - - save_subdirs: list (dtype=str) + save_dir: list (dtype=str) saving subdirectory string paths tmp_dir: str @@ -375,75 +381,143 @@ def load_microscopy_image(cli_args): img_name: str microscopy image filename + + cli_args: see ArgumentParser.parse_args + updated namespace of command line arguments """ - # initialize variables - skip_frangi = False - skip_odf = cli_args.odf_res is None + # create temporary directory tmp_dir = tempfile.mkdtemp() - # retrieve volume path and name - img_path, img_name, img_fmt, mosaic, in_mmap = get_image_file(cli_args) + # retrieve input file information + img_path, img_name, img_fmt, is_tiled, is_mmap = get_file_info(cli_args) - # import start time + # import fiber orientation vector data tic = perf_counter() - - # fiber orientation vector data if img_fmt == 'npy' or img_fmt == 'h5': + img, is_fiber = load_orient(img_path, img_name, img_fmt) - # print heading - print(color_text(0, 191, 255, "\nFiber Orientation Data Import\n")) + # import raw microscopy volume image + else: + img, is_fiber = load_raw(img_path, img_name, img_fmt, is_tiled=is_tiled, is_mmap=is_mmap, tmp_dir=tmp_dir) - # load fiber orientations - if img_fmt == 'npy': - img = np.load(img_path, mmap_mode='r') - else: - img_file = h5py.File(img_path, 'r') - img = img_file.get(img_file.keys()[0]) + # print import time + print_import_time(tic) + + # print volume image shape + print_image_shape(cli_args, img, is_tiled) if not is_fiber else print() + + # create saving directory + save_dir = create_save_dirs(img_path, img_name, cli_args, is_fiber=is_fiber) + + return img, is_tiled, is_fiber, save_dir, tmp_dir, img_name - # check dimensions - if img.ndim != 4: - raise ValueError('Invalid fiber orientation dataset (ndim != 4)') - else: - skip_frangi = True - print("Loading {} orientation dataset...\n".format(img_name)) - # microscopy volume image +def load_orient(img_path, img_name, img_fmt): + """ + Load array of 3D fiber orientations. + + Parameters + ---------- + img_path: str + path to the microscopy volume image + + img_name: str + name of the microscopy volume image + + img_fmt: str + format of the microscopy volume image + + Returns + ------- + img: numpy.ndarray + 3D fiber orientation vectors + + is_fiber: bool + True when pre-estimated fiber orientation vectors + are directly provided to the pipeline + """ + + # print heading + print(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') else: - # print heading - print(color_text(0, 191, 255, "\nMicroscopy Volume Image Import\n")) + img_file = h5py.File(img_path, 'r') + img = img_file.get(img_file.keys()[0]) - # load microscopy tiled reconstruction (aligned using ZetaStitcher) - if mosaic: - print("Loading {} tiled reconstruction...\n".format(img_name)) - img = VirtualFusedVolume(img_path) + # check array shape + if img.ndim != 4: + raise ValueError('Invalid 3D fiber orientation dataset (ndim != 4)!') + else: + is_fiber = True + print("Loading {} orientation dataset...\n".format(img_name)) - # load microscopy z-stack - else: - print("Loading {} z-stack...\n".format(img_name)) - img_fmt = img_fmt.lower() - if img_fmt == 'npy': - img = np.load(img_path) - elif img_fmt == 'tif' or img_fmt == 'tiff': - img = tiff.imread(img_path) - else: - raise ValueError('Unsupported image format!') - - # grey channel fiber image - if len(img.shape) == 3: - cli_args.neuron_mask = False - - # create image memory map - if in_mmap: - img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp=tmp_dir, arr=img[:], mmap_mode='r') + return img, is_fiber - # print import time - print_import_time(tic) - # print volume image shape - print_image_shape(cli_args, img, mosaic) if not skip_frangi else print() +def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None): + """ + Load raw microscopy volume image. - # create saving directory - save_subdirs = create_save_dirs(img_path, img_name, skip_frangi=skip_frangi, skip_odf=skip_odf) + Parameters + ---------- + img_path: str + path to the microscopy volume image + + img_name: str + name of the microscopy volume image + + img_fmt: str + format of the microscopy volume image + + is_tiled: bool + True for tiled reconstructions aligned using ZetaStitcher + + is_mmap: bool + create a memory-mapped array of the microscopy volume image, + increasing the parallel processing performance + (the image will be preliminarily loaded to RAM) + + tmp_dir: str + temporary file directory + + Returns + ------- + img: numpy.ndarray or NumPy memory-map object + microscopy volume image + + is_fiber: bool + True when pre-estimated fiber orientation vectors + are directly provided to the pipeline + """ + + # print heading + print(color_text(0, 191, 255, "\nMicroscopy Volume Image Import\n")) + + # load microscopy tiled reconstruction (aligned using ZetaStitcher) + if is_tiled: + print("Loading {} tiled reconstruction...\n".format(img_name)) + img = VirtualFusedVolume(img_path) + + # load microscopy z-stack + else: + print("Loading {} z-stack...\n".format(img_name)) + img_fmt = img_fmt.lower() + if img_fmt == 'npy': + img = np.load(img_path) + elif 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') + + # raw image + is_fiber = False - return img, mosaic, skip_frangi, cli_args, save_subdirs, tmp_dir, img_name + return img, is_fiber diff --git a/foa3d/odf.py b/foa3d/odf.py index ff9f385..ad5db01 100644 --- a/foa3d/odf.py +++ b/foa3d/odf.py @@ -38,82 +38,153 @@ def compute_fiber_angles(fiber_vec, norm): return phi, theta -@njit(cache=True) -def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis): +def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, vec_tensor_eigen, vxl_side, odf_norm, + odf_deg=6, vxl_thr=0.5, vec_thr=-1): """ - Compute orientation dispersion parameters. + Compute the spherical harmonics coefficients iterating over super-voxels + of fiber orientation vectors. Parameters ---------- - vec_tensor_eigen: numpy.ndarray (axis orde=(Z,Y,X,C), dtype=float32) - orientation tensor eigenvalues - computed from an ODF super-voxel + fiber_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float) + fiber orientation vectors - Returns - ------- - odi_pri: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) + 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_sec: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) secondary orientation dispersion index - odi_tot: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) total orientation dispersion index - odi_anis: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) + odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) 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 + + vxl_side: int + side of the ODF super-voxel [px] + + odf_norm: numpy.ndarray (dtype: float) + 2D array of spherical harmonics normalization factors + + odf_deg: int + degrees of the spherical harmonics series expansion + + vxl_thr: float + minimum relative threshold on the sliced voxel volume + + vec_thr: float + minimum relative threshold on non-zero orientation vectors + + Returns + ------- + odf: numpy.ndarray (axis order=(X,Y,Z,C), dtype=float32) + volumetric map of real-valued spherical harmonics coefficients """ - k1 = np.abs(vec_tensor_eigen[..., 2]) - k2 = np.abs(vec_tensor_eigen[..., 1]) - k3 = np.abs(vec_tensor_eigen[..., 0]) - diff = np.abs(vec_tensor_eigen[..., 1] - vec_tensor_eigen[..., 0]) + # 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): + zmax = z + vxl_side - # primary dispersion (1.2732395447351628 = 4/π) - odi_pri = 2 - 1.2732395447351628 * np.arctan2(k1, k2) + for y in range(0, fiber_vec.shape[1], vxl_side): + ymax = y + vxl_side - # secondary dispersion - odi_sec = 2 - 1.2732395447351628 * np.arctan2(k1, k3) + for x in range(0, fiber_vec.shape[2], vxl_side): + xmax = x + vxl_side - # total dispersion - odi_tot = 2 - 1.2732395447351628 * np.arctan2(k1, np.sqrt(np.abs(np.multiply(k2, k3)))) + # slice orientation voxel + vec_vxl = fiber_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]) - # dispersion anisotropy - odi_anis = 2 - 1.2732395447351628 * np.arctan2(k1, diff) + # skip boundary voxels and voxels without enough non-zero orientation vectors + if sli_vxl_size / ref_vxl_size > vxl_thr and 1 - zerovec / sli_vxl_size > vec_thr: + + # flatten orientation supervoxel + vec_arr = vec_vxl.ravel() + + # compute ODF and orientation tensor eigenvalues + zv, yv, xv = z // vxl_side, y // vxl_side, x // vxl_side + odf[zv, yv, xv, :] = fiber_vectors_to_sph_harm(vec_arr, odf_deg, odf_norm) + vec_tensor_eigen[zv, yv, xv, :] = compute_vec_tensor_eigen(vec_arr) + + # compute slice-wise dispersion and anisotropy parameters + compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis) - return odi_pri, odi_sec, odi_tot, odi_anis + # keep a black background + mask_orientation_dispersion([odi_pri, odi_sec, odi_tot, odi_anis], vec_tensor_eigen) + + # transform axes + odf = transform_axes(odf, swapped=(0, 2), flipped=(1, 2)) + + return odf -def mask_orientation_dispersion(odi_tpl, vec_tensor_eigen): +@njit(cache=True) +def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis): """ - Keep a black background where no fiber orientation vector - is available. + Compute orientation dispersion parameters. Parameters ---------- - odi_tpl: tuple - tuple including the ODI maps estimated - in compute_orientation_dispersion() - - vec_tensor_eigen: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) + 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_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_tpl: tuple - updated tuple including the masked ODI maps + 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 """ - bg_msk = np.all(vec_tensor_eigen == 0, axis=-1) + k1 = np.abs(vec_tensor_eigen[..., 2]) + k2 = np.abs(vec_tensor_eigen[..., 1]) + k3 = np.abs(vec_tensor_eigen[..., 0]) + diff = np.abs(vec_tensor_eigen[..., 1] - vec_tensor_eigen[..., 0]) - odi_lst = list(odi_tpl) - for odi in odi_lst: - odi[bg_msk] = 0 + # primary dispersion (1.2732395447351628 = 4/π) + odi_pri[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, k2))) \ + .astype(np.uint8) + + # secondary dispersion + odi_sec[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, k3))) \ + .astype(np.uint8) - odi_tpl = tuple(odi_lst) + # total dispersion + odi_tot[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, np.sqrt(np.abs(np.multiply(k2, k3)))))) \ + .astype(np.uint8) - return odi_tpl + # dispersion anisotropy + odi_anis[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, diff))) \ + .astype(np.uint8) @njit(cache=True) @@ -194,114 +265,6 @@ def compute_vec_tensor_eigen(fiber_vec): return vec_tensor_eigen -def estimate_odf_coeff(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, vec_tensor_eigen, vxl_side, odf_norm, - odf_degrees=6, vxl_thr=0.5, vec_thr=-1): - """ - Estimate the spherical harmonics coefficients iterating over super-voxels - of fiber orientation vectors. - - Parameters - ---------- - fiber_vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) - fiber orientation vectors - - odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32) - initialized array of ODF spherical harmonics coefficients - - odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of primary orientation dispersion parameters - - odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of secondary orientation dispersion parameters - - odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of total orientation dispersion parameters - - odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8) - initialized array of orientation dispersion anisotropy parameters - - vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8) - initialized array of orientation tensor eigenvalues - - vxl_side: int - side of the ODF super-voxel [px] - - odf_norm: numpy.ndarray (dtype: float) - 2D array of spherical harmonics normalization factors - - odf_degrees: int - degrees of the spherical harmonics series expansion - - vxl_thr: float - minimum relative threshold on the sliced voxel volume - - vec_thr: float - minimum relative threshold on non-zero orientation vectors - - Returns - ------- - odf: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32) - volumetric map of real-valued spherical harmonics coefficients - - odi_pri: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - primary orientation dispersion index - - odi_sec: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - secondary orientation dispersion index - - odi_tot: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - total orientation dispersion index - - odi_anis: numpy.ndarray (axis order=(Z,Y,X), dtype=float32) - orientation dispersion index anisotropy - """ - - # compute spherical harmonics normalization factors (once) - odf_norm = get_sph_harm_norm_factors(odf_degrees) - - # iterate over ODF super-voxels - fiber_vec_shape = np.array(fiber_vec.shape) - ref_vxl_size = min(vxl_side, fiber_vec_shape[0]) * vxl_side**2 - for z in range(0, fiber_vec_shape[0], vxl_side): - zmax = z + vxl_side - - for y in range(0, fiber_vec_shape[1], vxl_side): - ymax = y + vxl_side - - for x in range(0, fiber_vec_shape[2], vxl_side): - xmax = x + vxl_side - - # slice vector voxel (skip boundary voxels and voxels without enough non-zero orientation vectors) - vec_vxl = fiber_vec[z:zmax, y:ymax, x:xmax, :] - zero_vecs = np.count_nonzero(np.all(vec_vxl == 0, axis=-1)) - sli_vxl_size = np.prod(vec_vxl.shape[:-1]) - if sli_vxl_size / ref_vxl_size > vxl_thr and 1 - zero_vecs / sli_vxl_size > vec_thr: - vec_arr = vec_vxl.ravel() - odf[z // vxl_side, y // vxl_side, x // vxl_side, :] \ - = fiber_vectors_to_sph_harm(vec_arr, odf_degrees, odf_norm) - vec_tensor_eigen[z // vxl_side, y // vxl_side, x // vxl_side, :] \ - = compute_vec_tensor_eigen(vec_arr) - - # compute dispersion and anisotropy parameters slice-wise - odi_pri, odi_sec, odi_tot, odi_anis \ - = compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis) - - # keep a black background - odi_pri, odi_sec, odi_tot, odi_anis \ - = mask_orientation_dispersion((odi_pri, odi_sec, odi_tot, odi_anis), vec_tensor_eigen) - - # transform axes - odf = transform_axes(odf, swapped=(0, 2), flipped=(1, 2)) - - # adjust data type - odi_pri = (255.0 * odi_pri).astype(np.uint8) - odi_sec = (255.0 * odi_sec).astype(np.uint8) - odi_tot = (255.0 * odi_tot).astype(np.uint8) - odi_anis = (255.0 * odi_anis).astype(np.uint8) - - return odf, odi_pri, odi_sec, odi_tot, odi_anis - - factorial_lut = np.array([ 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, @@ -366,8 +329,8 @@ def fiber_angles_to_sph_harm(phi, theta, degrees, norm_factors, ncoeff): real_sph_harm = np.zeros(ncoeff) i = 0 - for n in np.arange(0, degrees + 1, 2): - for m in np.arange(-n, n + 1, 1): + for n in range(0, degrees + 1, 2): + for m in range(-n, n + 1, 1): for j, (p, t) in enumerate(zip(phi, theta)): real_sph_harm[i] += compute_real_sph_harm(n, m, p, np.sin(t), np.cos(t), norm_factors) i += 1 @@ -442,22 +405,20 @@ def generate_odf_background(bg_img, bg_mrtrix_mmap, vxl_side): new_shape = bg_mrtrix_mmap.shape[:-1] # normalize - dims = bg_img.ndim - if dims == 3: + if bg_img.ndim == 3: bg_img = normalize_image(bg_img) # loop over z-slices, and resize them for z in range(0, bg_img.shape[0], vxl_side): - if dims == 3: + if bg_img.ndim == 3: tmp_slice = np.mean(bg_img[z:z + vxl_side, ...], axis=0) - elif dims == 4: + 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 = np.swapaxes(tmp_slice, 0, 1).astype(np.uint8) - bg_mrtrix_mmap[..., z // vxl_side] = resize(tmp_slice, output_shape=new_shape, - anti_aliasing=True, preserve_range=True) - return bg_mrtrix_mmap + bg_mrtrix_mmap[..., z // vxl_side] = \ + resize(tmp_slice, output_shape=new_shape, anti_aliasing=True, preserve_range=True) @njit(cache=True) @@ -499,8 +460,8 @@ def get_sph_harm_norm_factors(degrees): """ norm_factors = np.zeros(shape=(degrees + 1, 2 * degrees + 1)) - for n in np.arange(0, degrees + 1, 2): - for m in np.arange(0, n + 1, 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] @@ -508,6 +469,32 @@ def get_sph_harm_norm_factors(degrees): return norm_factors +def mask_orientation_dispersion(odi_lst, vec_tensor_eigen): + """ + Keep a black background where no fiber orientation vector + is available. + + Parameters + ---------- + odi_lst: list + list including the ODI maps estimated + in compute_orientation_dispersion() + + 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 + """ + + bg_msk = np.all(vec_tensor_eigen == 0, axis=-1) + + for odi in odi_lst: + odi[bg_msk] = 0 + + @njit(cache=True) def norm_factor(n, m): """ @@ -538,177 +525,177 @@ def norm_factor(n, m): @njit(cache=True) -def sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm_factor): +def sph_harm_degree_2(order, phi, sin_theta, cos_theta, norm): if order == -2: - return norm_factor[2] * 3 * sin_theta**2 * np.sin(2 * phi) + return norm[2] * 3 * sin_theta**2 * np.sin(2 * phi) elif order == -1: - return norm_factor[1] * 3 * sin_theta * cos_theta * np.sin(phi) + return norm[1] * 3 * sin_theta * cos_theta * np.sin(phi) elif order == 0: - return norm_factor[0] * 0.5 * (3 * cos_theta**2 - 1) + return norm[0] * 0.5 * (3 * cos_theta ** 2 - 1) elif order == 1: - return norm_factor[1] * 3 * sin_theta * cos_theta * np.cos(phi) + return norm[1] * 3 * sin_theta * cos_theta * np.cos(phi) elif order == 2: - return norm_factor[2] * 3 * sin_theta**2 * np.cos(2 * phi) + 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_factor): +def sph_harm_degree_4(order, phi, sin_theta, cos_theta, norm): if order == -4: - return norm_factor[4] * 105 * sin_theta**4 * np.sin(4 * phi) + return norm[4] * 105 * sin_theta**4 * np.sin(4 * phi) elif order == -3: - return norm_factor[3] * 105 * sin_theta**3 * cos_theta * np.sin(3 * phi) + return norm[3] * 105 * sin_theta**3 * cos_theta * np.sin(3 * phi) elif order == -2: - return norm_factor[2] * 7.5 * sin_theta**2 * (7 * cos_theta**2 - 1) * np.sin(2 * phi) + return norm[2] * 7.5 * sin_theta**2 * (7 * cos_theta ** 2 - 1) * np.sin(2 * phi) elif order == -1: - return norm_factor[1] * 2.5 * sin_theta * (7 * cos_theta**3 - 3 * cos_theta) * np.sin(phi) + return norm[1] * 2.5 * sin_theta * (7 * cos_theta ** 3 - 3 * cos_theta) * np.sin(phi) elif order == 0: - return norm_factor[0] * 0.125 * (35 * cos_theta**4 - 30 * cos_theta**2 + 3) + return norm[0] * 0.125 * (35 * cos_theta ** 4 - 30 * cos_theta ** 2 + 3) elif order == 1: - return norm_factor[1] * 2.5 * sin_theta * (7 * cos_theta**3 - 3 * cos_theta) * np.cos(phi) + return norm[1] * 2.5 * sin_theta * (7 * cos_theta ** 3 - 3 * cos_theta) * np.cos(phi) elif order == 2: - return norm_factor[2] * 7.5 * sin_theta**2 * (7 * cos_theta**2 - 1) * np.cos(2 * phi) + return norm[2] * 7.5 * sin_theta**2 * (7 * cos_theta ** 2 - 1) * np.cos(2 * phi) elif order == 3: - return norm_factor[3] * 105 * sin_theta**3 * cos_theta * np.cos(3 * phi) + return norm[3] * 105 * sin_theta**3 * cos_theta * np.cos(3 * phi) elif order == 4: - return norm_factor[4] * 105 * sin_theta**4 * np.cos(4 * phi) + 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_factor): +def sph_harm_degree_6(order, phi, sin_theta, cos_theta, norm): if order == -6: - return norm_factor[6] * 10395 * sin_theta**6 * np.sin(6 * phi) + return norm[6] * 10395 * sin_theta**6 * np.sin(6 * phi) elif order == -5: - return norm_factor[5] * 10395 * sin_theta**5 * cos_theta * np.sin(5 * phi) + return norm[5] * 10395 * sin_theta**5 * cos_theta * np.sin(5 * phi) elif order == -4: - return norm_factor[4] * 472.5 * sin_theta**4 * (11 * cos_theta**2 - 1) * np.sin(4 * phi) + return norm[4] * 472.5 * sin_theta**4 * (11 * cos_theta ** 2 - 1) * np.sin(4 * phi) elif order == -3: - return norm_factor[3] * 157.5 * sin_theta**3 * (11 * cos_theta**3 - 3 * cos_theta) * np.sin(3 * phi) + return norm[3] * 157.5 * sin_theta**3 * (11 * cos_theta ** 3 - 3 * cos_theta) * np.sin(3 * phi) elif order == -2: - return norm_factor[2] * 13.125 * sin_theta**2 * (33 * cos_theta**4 - 18 * cos_theta**2 + 1) * np.sin(2 * phi) + 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_factor[1] * 2.625 * sin_theta \ + 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_factor[0] * 0.0625 * (231 * cos_theta**6 - 315 * cos_theta**4 + 105 * cos_theta**2 - 5) + return norm[0] * 0.0625 * (231 * cos_theta ** 6 - 315 * cos_theta ** 4 + 105 * cos_theta ** 2 - 5) elif order == 1: - return norm_factor[1] * 2.625 * sin_theta \ + 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_factor[2] * 13.125 * sin_theta**2 * (33 * cos_theta**4 - 18 * cos_theta**2 + 1) * np.cos(2 * phi) + 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_factor[3] * 157.5 * sin_theta**3 * (11 * cos_theta**3 - 3 * cos_theta) * np.cos(3 * phi) + return norm[3] * 157.5 * sin_theta**3 * (11 * cos_theta ** 3 - 3 * cos_theta) * np.cos(3 * phi) elif order == 4: - return norm_factor[4] * 472.5 * sin_theta**4 * (11 * cos_theta**2 - 1) * np.cos(4 * phi) + return norm[4] * 472.5 * sin_theta**4 * (11 * cos_theta ** 2 - 1) * np.cos(4 * phi) elif order == 5: - return norm_factor[5] * 10395 * sin_theta**5 * cos_theta * np.cos(5 * phi) + return norm[5] * 10395 * sin_theta**5 * cos_theta * np.cos(5 * phi) elif order == 6: - return norm_factor[6] * 10395 * sin_theta**6 * np.cos(6 * phi) + 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_factor): +def sph_harm_degree_8(order, phi, sin_theta, cos_theta, norm): if order == -8: - return norm_factor[8] * 2027025 * sin_theta**8 * np.sin(8 * phi) + return norm[8] * 2027025 * sin_theta**8 * np.sin(8 * phi) elif order == -7: - return norm_factor[7] * 2027025 * sin_theta**7 * cos_theta * np.sin(7 * phi) + return norm[7] * 2027025 * sin_theta**7 * cos_theta * np.sin(7 * phi) elif order == -6: - return norm_factor[6] * 67567.5 * sin_theta**6 * (15 * cos_theta**2 - 1) * np.sin(6 * phi) + return norm[6] * 67567.5 * sin_theta**6 * (15 * cos_theta ** 2 - 1) * np.sin(6 * phi) elif order == -5: - return norm_factor[5] * 67567.5 * sin_theta**5 * (5 * cos_theta**3 - cos_theta) * np.sin(5 * phi) + return norm[5] * 67567.5 * sin_theta**5 * (5 * cos_theta ** 3 - cos_theta) * np.sin(5 * phi) elif order == -4: - return norm_factor[4] * 1299.375 * sin_theta**4 * (65 * cos_theta**4 - 26 * cos_theta**2 + 1) * np.sin(4 * phi) + 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_factor[3] * 433.125 * sin_theta**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_factor[2] * 19.6875 * sin_theta**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_factor[1] * 0.5625 * sin_theta \ + 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_factor[0] * 0.0078125 \ + 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_factor[1] * 0.5625 * sin_theta \ + 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_factor[2] * 19.6875 * sin_theta**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_factor[3] * 433.125 * sin_theta**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_factor[4] * 1299.375 * sin_theta**4 * (65 * cos_theta**4 - 26 * cos_theta**2 + 1) * np.cos(4 * phi) + 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_factor[5] * 67567.5 * sin_theta**5 * (5 * cos_theta**3 - cos_theta) * np.cos(5 * phi) + return norm[5] * 67567.5 * sin_theta**5 * (5 * cos_theta ** 3 - cos_theta) * np.cos(5 * phi) elif order == 6: - return norm_factor[6] * 67567.5 * sin_theta**6 * (15 * cos_theta**2 - 1) * np.cos(6 * phi) + return norm[6] * 67567.5 * sin_theta**6 * (15 * cos_theta ** 2 - 1) * np.cos(6 * phi) elif order == 7: - return norm_factor[7] * 2027025 * sin_theta**7 * cos_theta * np.cos(7 * phi) + return norm[7] * 2027025 * sin_theta**7 * cos_theta * np.cos(7 * phi) elif order == 8: - return norm_factor[8] * 2027025 * sin_theta**8 * np.cos(8 * phi) + 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_factor): +def sph_harm_degree_10(order, phi, sin_theta, cos_theta, norm): if order == -10: - return norm_factor[10] * 654729075 * sin_theta**10 * np.sin(10 * phi) + return norm[10] * 654729075 * sin_theta**10 * np.sin(10 * phi) elif order == -9: - return norm_factor[9] * 654729075 * sin_theta**9 * cos_theta * np.sin(9 * phi) + return norm[9] * 654729075 * sin_theta**9 * cos_theta * np.sin(9 * phi) elif order == -8: - return norm_factor[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta**2 - 1) * np.sin(8 * phi) + return norm[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta ** 2 - 1) * np.sin(8 * phi) elif order == -7: - return norm_factor[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta**3 - 3 * cos_theta) * np.sin(7 * phi) + return norm[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta ** 3 - 3 * cos_theta) * np.sin(7 * phi) elif order == -6: - return norm_factor[6] * 84459.375 * sin_theta**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_factor[5] * 16891.875 * sin_theta**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_factor[4] * 2815.3125 * sin_theta**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_factor[3] * 402.1875 * sin_theta**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_factor[2] * 3.8671875 * sin_theta**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_factor[1] * 0.4296875 * sin_theta \ + 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_factor[0] * 0.00390625 \ + 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_factor[1] * 0.4296875 * sin_theta \ + 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_factor[2] * 3.8671875 * sin_theta**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_factor[3] * 402.1875 * sin_theta**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_factor[4] * 2815.3125 * sin_theta**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_factor[5] * 16891.875 * sin_theta**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_factor[6] * 84459.375 * sin_theta**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_factor[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta**3 - 3 * cos_theta) * np.cos(7 * phi) + return norm[7] * 5743237.5 * sin_theta**7 * (19 * cos_theta ** 3 - 3 * cos_theta) * np.cos(7 * phi) elif order == 8: - return norm_factor[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta**2 - 1) * np.cos(8 * phi) + return norm[8] * 17229712.5 * sin_theta**8 * (19 * cos_theta ** 2 - 1) * np.cos(8 * phi) elif order == 9: - return norm_factor[9] * 654729075 * sin_theta**9 * cos_theta * np.cos(9 * phi) + return norm[9] * 654729075 * sin_theta**9 * cos_theta * np.cos(9 * phi) elif order == 10: - return norm_factor[10] * 654729075 * sin_theta**10 * np.cos(10 * phi) + return norm[10] * 654729075 * sin_theta**10 * np.cos(10 * phi) diff --git a/foa3d/output.py b/foa3d/output.py index 0e5942b..345b470 100644 --- a/foa3d/output.py +++ b/foa3d/output.py @@ -6,7 +6,7 @@ import tifffile as tiff -def create_save_dirs(img_path, img_name, skip_frangi=False, skip_odf=False): +def create_save_dirs(img_path, img_name, cli_args, is_fiber=False): """ Create saving directory. @@ -18,14 +18,13 @@ def create_save_dirs(img_path, img_name, skip_frangi=False, skip_odf=False): img_name: str name of the input volume image - skip_frangi: bool + cli_args: see ArgumentParser.parse_args + updated namespace of command line arguments + + is_fiber: bool True when fiber orientation vectors are provided as input to the pipeline - skip_odf: bool - True when no ODF analysis is required following - the Frangi filtering stage - Returns ------- save_subdirs: list (dtype=str) @@ -39,31 +38,31 @@ def create_save_dirs(img_path, img_name, skip_frangi=False, skip_odf=False): base_path = path.dirname(img_path) # create saving directory - save_dir = path.join(base_path, time_stamp + '_' + img_name) - save_subdirs = [] - if not path.isdir(save_dir): - mkdir(save_dir) + base_dir = path.join(base_path, time_stamp + '_' + img_name) + save_dir = list() + if not path.isdir(base_dir): + mkdir(base_dir) # create Frangi filter output subdirectory - if not skip_frangi: - frangi_dir = path.join(save_dir, 'frangi') + if not is_fiber: + frangi_dir = path.join(base_dir, 'frangi') mkdir(frangi_dir) - save_subdirs.append(frangi_dir) + save_dir.append(frangi_dir) else: - save_subdirs.append(None) + save_dir.append(None) # create ODF analysis output subdirectory - if not skip_odf: - odf_dir = path.join(save_dir, 'odf') + if cli_args.odf_res is not None: + odf_dir = path.join(base_dir, 'odf') mkdir(odf_dir) - save_subdirs.append(odf_dir) + save_dir.append(odf_dir) else: - save_subdirs.append(None) + save_dir.append(None) - return save_subdirs + return save_dir -def save_array(fname, save_dir, nd_array, px_size=None, format='tif', odi=False): +def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', odi=False): """ Save array to file. @@ -78,10 +77,10 @@ def save_array(fname, save_dir, nd_array, px_size=None, format='tif', odi=False) nd_array: NumPy memory-map object or HDF5 dataset data - px_size: tuple + px_sz: tuple pixel size (Z,Y,X) [um] - format: str + fmt: str output format odi: bool @@ -93,31 +92,32 @@ def save_array(fname, save_dir, nd_array, px_size=None, format='tif', odi=False) """ # check output format - format = format.lower() - if format == 'tif' or format == 'tiff': + fmt = fmt.lower() + if fmt == 'tif' or fmt == 'tiff': # retrieve image pixel size - px_size_z, px_size_y, px_size_x = px_size + px_sz_z, px_sz_y, px_sz_x = px_sz # adjust bigtiff optional argument bigtiff = True if nd_array.itemsize * np.prod(nd_array.shape) >= 4294967296 else False # save array to TIFF file + out_name = '{}.{}'.format(fname, fmt) if odi: - tiff.imwrite(path.join(save_dir, fname + '.' + format), nd_array, imagej=True, bigtiff=bigtiff, - resolution=(1 / px_size_x, 1 / px_size_y), - metadata={'axes': 'ZYX', 'spacing': px_size_z, 'unit': 'um'}, compression='zlib') + tiff.imwrite(path.join(save_dir, out_name), nd_array, imagej=True, bigtiff=bigtiff, + resolution=(1 / px_sz_x, 1 / px_sz_y), + metadata={'axes': 'ZYX', 'spacing': px_sz_z, 'unit': 'um'}, compression='zlib') else: - tiff.imwrite(path.join(save_dir, fname + '.' + format), nd_array, imagej=True, bigtiff=bigtiff, - resolution=(1 / px_size_x, 1 / px_size_y), - metadata={'spacing': px_size_z, 'unit': 'um'}, compression='zlib') + tiff.imwrite(path.join(save_dir, out_name), nd_array, imagej=True, bigtiff=bigtiff, + resolution=(1 / px_sz_x, 1 / px_sz_y), + metadata={'spacing': px_sz_z, 'unit': 'um'}, compression='zlib') # save array to NumPy file - elif format == 'npy': + elif fmt == 'npy': np.save(path.join(save_dir, fname + '.npy'), nd_array) # save array to NIfTI file - elif format == 'nii': + elif fmt == 'nii': nd_array = nib.Nifti1Image(nd_array, np.eye(4)) nd_array.to_filename(path.join(save_dir, fname + '.nii')) diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py index 6563901..42c7ec1 100644 --- a/foa3d/pipeline.py +++ b/foa3d/pipeline.py @@ -6,20 +6,20 @@ import psutil from joblib import Parallel, delayed -from foa3d.frangi import convert_frangi_scales, frangi_filter -from foa3d.input import get_image_info -from foa3d.odf import (estimate_odf_coeff, generate_odf_background, +from foa3d.frangi import frangi_filter +from foa3d.input import get_image_info, get_frangi_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_analysis_time, print_frangi_info, print_odf_info) from foa3d.slicing import (config_frangi_batch, config_frangi_slicing, - crop_slice, slice_channel) + crop_slice, crop_slice_lst, slice_channel) from foa3d.utils import (create_background_mask, create_hdf5_dset, create_memory_map, divide_nonzero, - get_available_cores, get_item_size, orient_colormap, - vector_colormap) + get_available_cores, get_item_size, hsv_orient_cmap, + rgb_orient_cmap) def compute_fractional_anisotropy(eigenval): @@ -40,16 +40,15 @@ def compute_fractional_anisotropy(eigenval): frac_anis = \ np.sqrt(0.5 * divide_nonzero( - (eigenval[..., 0] - eigenval[..., 1]) ** 2 + - (eigenval[..., 0] - eigenval[..., 2]) ** 2 + - (eigenval[..., 1] - eigenval[..., 2]) ** 2, + 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 -def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, - z_min=0, z_max=None, lpf_soma_mask=False, max_ram_mb=None): +def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, z_rng=(0, None), mask_lpf=False, ram=None): """ Initialize the output datasets of the Frangi filtering stage. @@ -68,17 +67,15 @@ def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, tmp_dir: str temporary file directory - z_min: int - minimum output z-depth in [px] + z_rng: int + output z-range in [px] - z_max: int - maximum output z-depth in [px] - - lpf_soma_mask: bool - neuronal body masking flag + mask_lpf: bool + if True, mask neuronal bodies exploiting the autofluorescence + signal of lipofuscin pigments - max_ram_mb: float - maximum RAM available [MB] + ram: float + maximum RAM available to the Frangi filtering stage [B] Returns ------- @@ -104,8 +101,8 @@ def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized fiber image (isotropic resolution) - neuron_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - initialized neuron mask image + soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + initialized soma mask image z_sel: NumPy slice object selected z-depth range @@ -116,6 +113,7 @@ def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, slice_shape = slice_shape.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] @@ -124,45 +122,40 @@ def init_frangi_volumes(img_shape, slice_shape, resize_ratio, tmp_dir, # output shape img_dims = len(img_shape) - dset_shape = tuple(np.ceil(resize_ratio * img_shape).astype(int)) - slice_shape[0] = dset_shape[0] + tot_shape = tuple(np.ceil(resize_ratio * img_shape).astype(int)) + slice_shape[0] = tot_shape[0] # fiber channel arrays - iso_fiber_img, _ = init_volume(dset_shape, dtype='uint8', chunks=slice_shape, - name='iso_fiber', tmp=tmp_dir, max_ram_mb=max_ram_mb) - frangi_img, _ = init_volume(dset_shape, dtype='uint8', chunks=slice_shape, - name='frangi', tmp=tmp_dir, max_ram_mb=max_ram_mb) - fiber_mask, _ = init_volume(dset_shape, dtype='uint8', chunks=slice_shape, - name='fiber_msk', tmp=tmp_dir, max_ram_mb=max_ram_mb) - frac_anis_img, _ = init_volume(dset_shape, dtype='uint8', chunks=slice_shape, - name='frac_anis', tmp=tmp_dir, max_ram_mb=max_ram_mb) - - # neuron channel array - if lpf_soma_mask: - neuron_msk, _ = init_volume(dset_shape, dtype='uint8', chunks=slice_shape, - name='neuron_msk', tmp=tmp_dir, max_ram_mb=max_ram_mb) + 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) + + # 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) else: - neuron_msk = None + soma_msk = None # fiber orientation arrays - vec_dset_shape = dset_shape + (img_dims,) + vec_shape = tot_shape + (img_dims,) vec_slice_shape = tuple(slice_shape) + (img_dims,) - fiber_vec_img, fiber_dset_path = init_volume(vec_dset_shape, dtype='float32', chunks=vec_slice_shape, - name='fiber_vec', tmp=tmp_dir, max_ram_mb=max_ram_mb) - fiber_vec_clr, _ = init_volume(vec_dset_shape, dtype='uint8', chunks=vec_slice_shape, - name='fiber_cmap', tmp=tmp_dir, max_ram_mb=max_ram_mb) + fiber_vec_img, fiber_dset_path = \ + init_volume(vec_shape, dtype='float32', chunks=vec_slice_shape, name='fiber_vec', tmp=tmp_dir, ram=ram) + fiber_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_mask, iso_fiber_img, neuron_msk, z_sel + frangi_img, fiber_msk, iso_fiber_img, soma_msk, z_sel -def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, max_ram_mb=None): +def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, ram=None): """ Initialize the output datasets of the ODF analysis stage. Parameters ---------- - vec_img_shape: numpy.ndarray (shape=(3,), dtype=int) + vec_img_shape: tuple vector volume shape [px] tmp_dir: str @@ -174,8 +167,8 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, max_ram_m odf_degrees: int degrees of the spherical harmonics series expansion - max_ram_mb: float - maximum RAM available [MB] + ram: float + maximum RAM available to the ODF generation stage [B] Returns ------- @@ -197,7 +190,7 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, max_ram_m odi_anis: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) initialized array of orientation dispersion anisotropy parameters - vec_tensor_eigen: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=uint8) + vec_tensor_eigen: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) initialized array of fiber orientation tensor eigenvalues """ @@ -207,34 +200,34 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, max_ram_m # 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,), - name='bg_tmp{0}'.format(odf_scale), tmp=tmp_dir, max_ram_mb=max_ram_mb) + 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), - name='odf_tmp{0}'.format(odf_scale), tmp=tmp_dir, max_ram_mb=max_ram_mb) + 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, max_ram_mb=max_ram_mb) + name='tensor_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) # create ODI memory maps odi_pri, _ = init_volume(odi_shape, dtype='uint8', - name='odi_pri_tmp{0}'.format(odf_scale), tmp=tmp_dir, max_ram_mb=max_ram_mb) + name='odi_pri_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) odi_sec, _ = init_volume(odi_shape, dtype='uint8', - name='odi_sec_tmp{0}'.format(odf_scale), tmp=tmp_dir, max_ram_mb=max_ram_mb) + name='odi_sec_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) odi_tot, _ = init_volume(odi_shape, dtype='uint8', - name='odi_tot_tmp{0}'.format(odf_scale), tmp=tmp_dir, max_ram_mb=max_ram_mb) + name='odi_tot_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) odi_anis, _ = init_volume(odi_shape, dtype='uint8', - name='odi_anis_tmp{0}'.format(odf_scale), tmp=tmp_dir, max_ram_mb=max_ram_mb) + name='odi_anis_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram) return odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, vec_tensor_eigen -def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+', max_ram_mb=None): +def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+', ram=None): """ Initialize output volume as an empty HDF5 dataset or a memory-mapped array depending on available RAM. @@ -259,37 +252,38 @@ def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+', mmap_mode: str file opening mode (memory-mapped object only) - max_ram_mb: float - maximum RAM available [MB] + ram: float + maximum RAM available to the Frangi filtering stage [B] Returns ------- vol: NumPy memory-map object or HDF5 dataset initialized data volume - file_path: str + hdf5_path: str path to the HDF5 file """ # get maximum RAM and initialized array memory size - max_ram = psutil.virtual_memory()[1] if max_ram_mb is None else max_ram_mb * 1024**2 + if ram is None: + ram = psutil.virtual_memory()[1] item_sz = get_item_size(dtype) vol_sz = item_sz * np.prod(shape) # create memory-mapped array or HDF5 file depending on available memory resources - if vol_sz >= max_ram: + if vol_sz >= ram: vol, hdf5_path = create_hdf5_dset(shape, dtype, chunks=chunks, name=name, tmp=tmp) else: - vol = create_memory_map(shape, dtype, name=name, tmp=tmp, mmap_mode=mmap_mode) + vol = create_memory_map(shape, dtype, name=name, tmp_dir=tmp, mmap_mode=mmap_mode) hdf5_path = None return vol, hdf5_path -def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad_mat, 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, neuron_msk, - ch_neuron=0, ch_fiber=1, alpha=0.05, beta=1, gamma=100, - dark=False, orient_cmap=False, lpf_soma_mask=False, mosaic=False): +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, + ch_lpf=0, ch_mye=1, alpha=0.05, beta=1, gamma=100, dark=False, mask_lpf=False, hsv_vec_cmap=False, + pad_mode='reflect', is_tiled=False, fiber_thresh='li', soma_thresh='yen'): """ Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image. @@ -307,8 +301,12 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad_mat, smooth_sigma, scal rng_out: NumPy slice object output range - pad_mat: numpy.ndarray (axis order=(Z,Y,X)) - 3D image padding 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 smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) 3D standard deviation of the low-pass Gaussian filter [px] @@ -341,13 +339,13 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad_mat, smooth_sigma, scal fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) fiber mask image - neuron_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - neuron mask image + soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + soma mask image - ch_neuron: int + ch_lpf: int neuronal bodies channel - ch_fiber: int + ch_mye: int myelinated fibers channel alpha: float @@ -363,87 +361,163 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad_mat, smooth_sigma, scal if True, enhance black 3D tubular structures (i.e., negative contrast polarity) - orient_cmap: bool - if True, generate color maps based on XY-plane orientation angles - (instead of using the cartesian components of the estimated vectors) + hsv_vec_cmap: bool + if True, generate a HSV orientation color map based on XY-plane orientation angles + (instead of an RGB map using the cartesian components of the estimated vectors) - lpf_soma_mask: bool + mask_lpf: bool if True, mask neuronal bodies exploiting the autofluorescence signal of lipofuscin pigments - mosaic: bool + pad_mode: str + image padding mode + + is_tiled: bool must be True for tiled reconstructions aligned using ZetaStitcher + fiber_thresh: str + enhanced fiber channel thresholding method + + soma_thresh: str + soma channel thresholding method + Returns ------- None """ # slice fiber image slice - fiber_slice = slice_channel(img, rng_in, channel=ch_fiber, mosaic=mosaic) + fiber_slice = slice_channel(img, rng_in, channel=ch_mye, is_tiled=is_tiled) # skip background slice if np.max(fiber_slice) != 0: # preprocess fiber slice - iso_fiber_slice, rsz_pad_mat = \ - correct_image_anisotropy(fiber_slice, px_rsz_ratio, sigma=smooth_sigma, pad_mat=pad_mat) + iso_fiber_slice, rsz_pad = \ + correct_image_anisotropy(fiber_slice, px_rsz_ratio, sigma=smooth_sigma, pad=pad) + + # pad fiber slice if required + iso_fiber_slice = np.pad(iso_fiber_slice, rsz_pad, mode=pad_mode) # 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) # crop resulting slices - iso_fiber_slice = crop_slice(iso_fiber_slice, rng_out, rsz_pad_mat) - frangi_slice = crop_slice(frangi_slice, rng_out, rsz_pad_mat) - fiber_vec_slice = crop_slice(fiber_vec_slice, rng_out, rsz_pad_mat) - eigenval_slice = crop_slice(eigenval_slice, rng_out, rsz_pad_mat) + iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice = \ + crop_slice_lst([iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice], rng_out, ovlp=ovlp) # generate fractional anisotropy image frac_anis_slice = compute_fractional_anisotropy(eigenval_slice) - # generate RGB orientation color map - if orient_cmap: - orientcol_slice = orient_colormap(fiber_vec_slice) - else: - orientcol_slice = vector_colormap(fiber_vec_slice) + # 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) - # mask background - fiber_vec_slice, orientcol_slice, fiber_mask_slice = \ - mask_background(frangi_slice, fiber_vec_slice, orientcol_slice, thresh_method='li', invert=False) + # remove background + fiber_vec_slice, fiber_clr_slice, fiber_msk_slice = \ + mask_background(frangi_slice, fiber_vec_slice, fiber_clr_slice, method=fiber_thresh, invert=False) # (optional) neuronal body masking - neuron_mask_slice = None - if lpf_soma_mask: + if mask_lpf: - # get neuron image slice - neuron_slice = slice_channel(img, rng_in_neu, channel=ch_neuron, mosaic=mosaic) + # get soma image slice + soma_slice = slice_channel(img, rng_in_neu, channel=ch_lpf, is_tiled=is_tiled) - # resize neuron slice (lateral blurring and downsampling) - iso_neuron_slice, _ = correct_image_anisotropy(neuron_slice, px_rsz_ratio) + # resize soma slice (lateral blurring and downsampling) + iso_soma_slice, _ = correct_image_anisotropy(soma_slice, px_rsz_ratio) - # crop isotropized neuron slice - iso_neuron_slice = crop_slice(iso_neuron_slice, rng_out) + # crop isotropized soma slice + iso_soma_slice = crop_slice(iso_soma_slice, rng_out) # mask neuronal bodies - fiber_vec_slice, orientcol_slice, frac_anis_slice, neuron_mask_slice = \ - mask_background(iso_neuron_slice, fiber_vec_slice, orientcol_slice, frac_anis_slice, - thresh_method='yen', invert=True) - - # fill memory-mapped output neuron mask - neuron_msk[rng_out] = (255 * neuron_mask_slice[z_sel, ...]).astype(np.uint8) + 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) + else: + soma_msk_slice = None # 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] = orientcol_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_mask_slice[z_sel, ...])).astype(np.uint8) + 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 mask_background(img, fiber_vec_slice, orientcol_slice, frac_anis_slice=None, thresh_method='yen', invert=False): +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): + """ + Fill the memory-mapped output arrays of the Frangi filtering stage. + + Parameters + ---------- + fiber_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) + orientation colormap image + + frac_anis_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) + isotropic fiber image + + fiber_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) + soma mask image + + fiber_vec_slice: 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) + orientation colormap image slice + + frac_anis_slice: 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-enhanced image slice + + iso_fiber_slice: 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) + fiber mask image slice + + soma_msk_slice: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8) + soma mask image slice + + rng_out: tuple + 3D slice output index range + + z_sel: NumPy slice object + selected z-depth range + + Returns + ------- + None + """ + + # 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) + + # 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) + + +def mask_background(img, fiber_vec_slice, fiber_clr_slice, frac_anis_slice=None, method='yen', invert=False): """ Mask orientation volume arrays. @@ -455,13 +529,13 @@ def mask_background(img, fiber_vec_slice, orientcol_slice, frac_anis_slice=None, fiber_vec_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) fiber orientation vector slice - orientcol_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) - orientation colormap slice + fiber_clr_slice: 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) fractional anisotropy slice - thresh_method: str + method: str thresholding method (refer to skimage.filters) invert: bool @@ -483,7 +557,7 @@ def mask_background(img, fiber_vec_slice, orientcol_slice, frac_anis_slice=None, """ # generate background mask - background = create_background_mask(img, thresh_method=thresh_method) + background = create_background_mask(img, method=method) # invert mask if invert: @@ -491,19 +565,19 @@ def mask_background(img, fiber_vec_slice, orientcol_slice, frac_anis_slice=None, # apply mask to input arrays fiber_vec_slice[background, :] = 0 - orientcol_slice[background, :] = 0 + fiber_clr_slice[background, :] = 0 # (optional) mask fractional anisotropy if frac_anis_slice is not None: frac_anis_slice[background] = 0 - return fiber_vec_slice, orientcol_slice, frac_anis_slice, background + return fiber_vec_slice, fiber_clr_slice, frac_anis_slice, background else: - return fiber_vec_slice, orientcol_slice, background + return fiber_vec_slice, fiber_clr_slice, background -def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, img_name, odf_scale_um, odf_norm, - odf_degrees=6, max_ram_mb=None): +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): """ Estimate 3D fiber ODFs from basic orientation data chunks using parallel threads. @@ -533,47 +607,38 @@ def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, i odf_norm: numpy.ndarray (dtype: float) 2D array of spherical harmonics normalization factors - odf_degrees: int + odf_deg: int degrees of the spherical harmonics series expansion - max_ram_mb: float - maximum RAM available [MB] + ram: float + maximum RAM available to the ODF analysis stage [B] Returns ------- None """ - # get info on the input volume of orientation vectors - vec_img_shape = np.asarray(fiber_vec_img.shape) - # derive the ODF kernel size in [px] odf_scale = int(np.ceil(odf_scale_um / px_size_iso[0])) # initialize ODF analysis output volumes odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, tensor \ - = init_odf_volumes(vec_img_shape[:-1], tmp_dir, - odf_scale=odf_scale, odf_degrees=odf_degrees, max_ram_mb=max_ram_mb) + = init_odf_volumes(fiber_vec_img.shape[:-1], tmp_dir, odf_scale=odf_scale, odf_degrees=odf_deg, ram=ram) # generate downsampled background for MRtrix3 mrview - if iso_fiber_img is None: - bg_mrtrix = generate_odf_background(fiber_vec_img, bg_mrtrix, vxl_side=odf_scale) - else: - bg_mrtrix = generate_odf_background(iso_fiber_img, bg_mrtrix, vxl_side=odf_scale) + bg_img = fiber_vec_img if iso_fiber_img is None else iso_fiber_img + generate_odf_background(bg_img, bg_mrtrix, vxl_side=odf_scale) # compute ODF coefficients - odf, odi_pri, odi_sec, odi_tot, odi_anis = \ - estimate_odf_coeff(fiber_vec_img, odf, odi_pri, odi_sec, odi_tot, odi_anis, tensor, - vxl_side=odf_scale, odf_norm=odf_norm, odf_degrees=odf_degrees) + odf = compute_odf_map(fiber_vec_img, odf, odi_pri, odi_sec, odi_tot, odi_anis, 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, px_size_iso, save_dir, img_name, odf_scale_um) -def parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_dir, tmp_dir, img_name, frangi_sigma_um, - ch_neuron=0, ch_fiber=1, alpha=0.05, beta=1, gamma=100, dark=False, z_min=0, z_max=None, - orient_cmap=False, lpf_soma_mask=False, mosaic=False, - max_ram_mb=None, jobs=4, backend='threading'): +def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, + ram=None, jobs=4, backend='threading', is_tiled=False, verbose=100): """ Apply 3D Frangi filtering to basic TPFM image slices using parallel threads. @@ -582,15 +647,8 @@ def parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_dir, img: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X)) microscopy volume image - px_size: numpy.ndarray (shape=(3,), dtype=float) - pixel size [μm] - - px_size_iso: numpy.ndarray (shape=(3,), dtype=float) - adjusted isotropic pixel size [μm] - - smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) - 3D standard deviation of the low-pass Gaussian filter [px] - (applied to the XY plane) + cli_args: see ArgumentParser.parse_args + populated namespace of command line arguments save_dir: str saving directory string path @@ -599,49 +657,10 @@ def parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_dir, temporary file directory img_name: str - name of the input volume image - - frangi_sigma_um: list (dtype=float) - Frangi filter scales in [μm] - - ch_neuron: int - neuronal bodies channel + name of the microscopy volume image - ch_fiber: int - myelinated fibers channel - - alpha: float - plate-like score sensitivity - - beta: float - blob-like score sensitivity - - gamma: float - background score sensitivity - - dark: bool - if True, enhance black 3D tubular structures - (i.e., negative contrast polarity) - - z_min: int - minimum output z-depth in [px] - - z_max: int - maximum output z-depth in [px] - - orient_cmap: bool - if True, generate color maps based on XY-plane orientation angles - (instead of using the cartesian components of the estimated vectors) - - lpf_soma_mask: bool - if True, mask neuronal bodies exploiting the autofluorescence - signal of lipofuscin pigments - - mosaic: bool - must be True for tiled reconstructions aligned using ZetaStitcher - - max_ram_mb: float - maximum RAM available to the Frangi filtering stage [MB] + ram: float + maximum RAM available to the Frangi filtering stage [B] jobs: int number of parallel jobs (threads) @@ -650,6 +669,12 @@ def parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_dir, backend: str backend module employed by joblib.Parallel + is_tiled: bool + must be True for tiled reconstructions aligned using ZetaStitcher + + verbose: int + joblib verbosity level + Returns ------- fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32) @@ -670,62 +695,63 @@ def parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_dir, fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) fiber mask image - neuron_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) - neuron mask image + soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + soma mask image """ + # get Frangi filter configuration + alpha, beta, gamma, frangi_sigma_px, frangi_sigma_um, smooth_sigma, px_size, px_size_iso, \ + z_rng, ch_lpf, ch_mye, mask_lpf, hsv_vec_cmap, img_name = get_frangi_config(cli_args, img_name) + # get info on the input volume image - img_shape, img_shape_um, img_item_size, ch_fiber = \ - get_image_info(img, px_size, ch_fiber, mosaic=mosaic) + img_shape, img_shape_um, img_item_size, ch_mye, mask_lpf = \ + get_image_info(img, px_size, mask_lpf, ch_mye, is_tiled=is_tiled) # configure batch of basic image slices analyzed in parallel batch_size, max_slice_size = \ - config_frangi_batch(frangi_sigma_um, max_ram_mb=max_ram_mb, jobs=jobs) - - # convert the spatial scales of the Frangi filter to pixel - frangi_sigma_px = convert_frangi_scales(frangi_sigma_um, px_size_iso[0]) + config_frangi_batch(frangi_sigma_um, ram=ram, jobs=jobs) # get info on the processed image slices - rng_in_lst, rng_in_neu_lst, rng_out_lst, pad_mat_lst, \ - in_slice_shape_um, out_slice_shape, px_rsz_ratio, tot_slice_num, batch_size = \ + rng_in_lst, rng_in_lpf_lst, rng_out_lst, pad_lst, \ + in_slice_shape_um, out_slice_shape, px_rsz_ratio, rsz_ovlp, tot_slice_num, batch_size = \ config_frangi_slicing(img_shape, img_item_size, px_size, px_size_iso, - smooth_sigma, frangi_sigma_px, lpf_soma_mask, batch_size, slice_size=max_slice_size) + smooth_sigma, frangi_sigma_um, mask_lpf, batch_size, slice_size=max_slice_size) # initialize output arrays fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_img, \ - frangi_img, fiber_msk, iso_fiber_img, neuron_msk, z_sel = \ + frangi_img, fiber_msk, iso_fiber_img, soma_msk, z_sel = \ init_frangi_volumes(img_shape, out_slice_shape, px_rsz_ratio, tmp_dir, - z_min=z_min, z_max=z_max, lpf_soma_mask=lpf_soma_mask, max_ram_mb=max_ram_mb) + z_rng=z_rng, mask_lpf=mask_lpf, ram=ram) # print Frangi filter configuration print_frangi_info(alpha, beta, gamma, frangi_sigma_um, img_shape_um, in_slice_shape_um, tot_slice_num, - px_size, img_item_size, lpf_soma_mask) + px_size, img_item_size, mask_lpf) # parallel Frangi filter-based fiber orientation analysis of microscopy sub-volumes start_time = perf_counter() - with Parallel(n_jobs=batch_size, backend=backend, verbose=100, max_nbytes=None) as parallel: + with Parallel(n_jobs=batch_size, backend=backend, verbose=verbose, max_nbytes=None) as parallel: parallel( delayed(fiber_analysis)(img, - rng_in_lst[i], rng_in_neu_lst[i], rng_out_lst[i], pad_mat_lst[i], + rng_in_lst[i], rng_in_lpf_lst[i], rng_out_lst[i], pad_lst[i], rsz_ovlp, smooth_sigma, frangi_sigma_px, px_rsz_ratio, z_sel, fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, - iso_fiber_img, fiber_msk, neuron_msk, - ch_neuron=ch_neuron, ch_fiber=ch_fiber, alpha=alpha, beta=beta, gamma=gamma, - dark=dark, orient_cmap=orient_cmap, lpf_soma_mask=lpf_soma_mask, mosaic=mosaic) + iso_fiber_img, fiber_msk, soma_msk, ch_lpf=ch_lpf, ch_mye=ch_mye, + alpha=alpha, beta=beta, gamma=gamma, mask_lpf=mask_lpf, + hsv_vec_cmap=hsv_vec_cmap, is_tiled=is_tiled) for i in range(tot_slice_num)) # save Frangi output arrays - save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, fiber_msk, neuron_msk, + save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img, fiber_msk, soma_msk, px_size_iso, save_dir, img_name) # print Frangi filtering time print_analysis_time(start_time) - return fiber_vec_img, iso_fiber_img + return fiber_vec_img, iso_fiber_img, px_size_iso, img_name -def parallel_odf_on_scales(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, img_name, - odf_scales_um, odf_degrees=6, backend='threading', max_ram_mb=None): +def parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_size_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 implemented in parallel_odf_on_slices(). @@ -738,6 +764,9 @@ def parallel_odf_on_scales(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, iso_fiber_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) adjusted isotropic pixel size [μm] @@ -750,17 +779,14 @@ def parallel_odf_on_scales(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, img_name: str name of the input volume image - odf_scales_um: list (dtype: float) - list of fiber ODF resolution values (super-voxel sides [μm]) - - odf_degrees: int - degrees of the spherical harmonics series expansion - backend: str backend module employed by joblib.Parallel - max_ram_mb: float - maximum RAM available [MB] + ram: float + maximum RAM available to the ODF analysis stage [B] + + verbose: int + joblib verbosity level Returns ------- @@ -771,29 +797,32 @@ def parallel_odf_on_scales(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, start_time = perf_counter() # print ODF analysis heading - print_odf_info(odf_scales_um, odf_degrees) + print_odf_info(cli_args.odf_res, cli_args.odf_deg) # compute spherical harmonics normalization factors (once for all scales) - norm_factors = get_sph_harm_norm_factors(odf_degrees) + odf_norm = get_sph_harm_norm_factors(cli_args.odf_deg) - # number of logical cores + # get number of logical cores 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,)) + # parallel ODF analysis of fiber orientation vectors # over the required spatial scales - n_jobs = min(num_cpu, len(odf_scales_um)) - with Parallel(n_jobs=n_jobs, backend=backend, verbose=100, max_nbytes=None) as parallel: + 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, - odf_norm=norm_factors, odf_degrees=odf_degrees, odf_scale_um=s, - max_ram_mb=max_ram_mb) - for s in odf_scales_um) + odf_norm=odf_norm, odf_deg=cli_args.odf_deg, odf_scale_um=s, ram=ram) + for s in cli_args.odf_res) # print ODF analysis time 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, neuron_msk, px_size, save_dir, img_name): + fiber_msk, soma_msk, px_size, save_dir, img_name): """ Save the output arrays of the Frangi filter stage to TIF files. @@ -818,7 +847,7 @@ def save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_ fiber_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) fiber mask image - neuron_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + soma_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) @@ -840,23 +869,23 @@ def save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_ shutil.move(fiber_dset_path, path.join(save_dir, 'fiber_vec_{0}.h5'.format(img_name))) # or save orientation vectors to NumPy file else: - save_array('fiber_vec_{0}'.format(img_name), save_dir, fiber_vec_img, format='npy') + save_array('fiber_vec_{0}'.format(img_name), save_dir, fiber_vec_img, fmt='npy') - # save orientation color map to TIF + # save orientation color map to TIFF save_array('fiber_cmap_{0}'.format(img_name), save_dir, fiber_vec_clr, px_size) - # save fractional anisotropy map to TIF + # save fractional anisotropy map to TIFF save_array('frac_anis_{0}'.format(img_name), save_dir, frac_anis_img, px_size) - # save Frangi-enhanced fiber volume to TIF + # save Frangi-enhanced fiber volume to TIFF save_array('frangi_{0}'.format(img_name), save_dir, frangi_img, px_size) - # save masked fiber volume to TIF + # save masked fiber volume to TIFF save_array('fiber_msk_{0}'.format(img_name), save_dir, fiber_msk, px_size) - # save masked neuron volume to TIF - if neuron_msk is not None: - save_array('neuron_msk_{0}'.format(img_name), save_dir, neuron_msk, px_size) + # 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) def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_dir, img_name, odf_scale_um): @@ -868,22 +897,22 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_ Parameters ---------- - odf_img: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z,C), dtype=float32) + odf: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z,C), dtype=float32) ODF spherical harmonics coefficients - bg_mrtrix_img: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z), dtype=uint8) + bg: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z), dtype=uint8) background for ODF visualization in MRtrix3 - odi_pri_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + odi_pri: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) primary orientation dispersion parameter - odi_sec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + odi_sec: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) secondary orientation dispersion parameter - odi_tot_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + odi_tot: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) total orientation dispersion parameter - odi_anis_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) + odi_anis: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) orientation dispersion anisotropy parameter px_size: numpy.ndarray (shape=(3,), dtype=float) @@ -903,8 +932,8 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_ None """ # 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, format='nii') - save_array('odf_mrtrixview_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odf, format='nii') + 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) diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py index e4ed482..9589bee 100644 --- a/foa3d/preprocessing.py +++ b/foa3d/preprocessing.py @@ -2,11 +2,11 @@ from scipy.ndimage import gaussian_filter from skimage.transform import resize -from foa3d.printing import print_prepro_heading +from foa3d.printing import print_blur, print_new_res, print_prepro_heading from foa3d.utils import fwhm_to_sigma -def config_anisotropy_correction(px_size, psf_fwhm, vector): +def config_anisotropy_correction(px_sz, psf_fwhm): """ Scanning and light-sheet fluorescence microscopes provide 3D data characterized by a lower resolution along the optical axis @@ -18,53 +18,43 @@ def config_anisotropy_correction(px_size, psf_fwhm, vector): 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) 3D FWHM of the PSF [μm] - vector: bool - Returns ------- smooth_sigma: numpy.ndarray (shape=(3,), dtype=int) 3D standard deviation of the low-pass Gaussian filter [px] (resolution anisotropy correction) - px_size_iso: numpy.ndarray (shape=(3,), dtype=float) + px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) new isotropic pixel size [μm] """ # print preprocessing heading - if not vector: - print_prepro_heading() + print_prepro_heading() # set the isotropic pixel resolution equal to the z-sampling step - px_size_iso = px_size[0] * np.ones(shape=(3,)) + px_sz_iso = px_sz[0] * np.ones(shape=(3,)) # adjust PSF anisotropy via lateral Gaussian blurring if not np.all(psf_fwhm == psf_fwhm[0]): - # estimate the PSF variance from input FWHM values [μm**2] + # estimate the PSF variance from input FWHM values [μm^2] psf_var = np.square(fwhm_to_sigma(psf_fwhm)) - # estimate the in-plane filter variance [μm**2] - gauss_var_x = psf_var[0] - psf_var[2] - gauss_var_y = psf_var[0] - psf_var[1] + # estimate the in-plane filter variance [μm^2] + gauss_var = np.array([0, psf_var[0] - psf_var[1], psf_var[0] - psf_var[2]]) # ...and the corresponding standard deviation [px] - gauss_sigma_x = np.sqrt(gauss_var_x) / px_size[2] - gauss_sigma_y = np.sqrt(gauss_var_y) / px_size[1] - gauss_sigma_z = 0 - smooth_sigma = np.array([gauss_sigma_z, gauss_sigma_y, gauss_sigma_x]) + smooth_sigma = np.divide(np.sqrt(gauss_var), px_sz) # print preprocessing info - gauss_sigma_um = np.multiply(smooth_sigma, px_size) - if not vector: - print("\n Z Y X") - print("Gaussian blur \u03C3 [μm]: ({0:.3f}, {1:.3f}, {2:.3f})" - .format(gauss_sigma_um[0], gauss_sigma_um[1], gauss_sigma_um[2]), end='\r') + smooth_sigma_um = np.multiply(smooth_sigma, px_sz) + print_blur(smooth_sigma_um) # (no blurring) else: @@ -72,17 +62,12 @@ def config_anisotropy_correction(px_size, psf_fwhm, vector): smooth_sigma = None # print pixel resize info - if not vector: - print("\nOriginal pixel size [μm]: ({0:.3f}, {1:.3f}, {2:.3f})" - .format(px_size[0], px_size[1], px_size[2])) - print("Adjusted pixel size [μm]: ({0:.3f}, {1:.3f}, {2:.3f})\n" - .format(px_size_iso[0], px_size_iso[1], px_size_iso[2])) + print_new_res(px_sz_iso, psf_fwhm) - return smooth_sigma, px_size_iso + return smooth_sigma, px_sz_iso -def correct_image_anisotropy(img, rsz_ratio, - sigma=None, pad_mat=None, pad_mode='reflect', anti_aliasing=True, trunc=4): +def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mode='reflect', anti_alias=True, trunc=4): """ 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. @@ -100,13 +85,13 @@ def correct_image_anisotropy(img, rsz_ratio, 3D standard deviation of the low-pass Gaussian filter [px] (resolution anisotropy correction) - pad_mat: numpy.ndarray (shape=(3,2), dtype=int) - padding range array + pad: numpy.ndarray (shape=(3,2), dtype=int) + image padding array to be resized - pad_mode: str - data padding mode adopted for the Gaussian filter + smooth_pad_mode: str + image padding mode adopted for the smoothing Gaussian filter - anti_aliasing: bool + anti_alias: bool if True, apply an anti-aliasing filter when downsampling the XY plane trunc: int @@ -117,8 +102,8 @@ def correct_image_anisotropy(img, rsz_ratio, iso_img: numpy.ndarray (axis order=(Z,Y,X)) isotropic microscopy volume image - rsz_pad_mat: numpy.ndarray (shape=(3,2), dtype=int) - resized padding range array + rsz_pad: numpy.ndarray (shape=(3,2), dtype=int) + resized image padding array """ # no resizing if np.all(rsz_ratio == 1): @@ -127,19 +112,19 @@ def correct_image_anisotropy(img, rsz_ratio, # lateral blurring else: if sigma is not None: - img = gaussian_filter(img, sigma=sigma, mode=pad_mode, truncate=trunc, output=np.float32) + img = gaussian_filter(img, sigma=sigma, mode=smooth_pad_mode, truncate=trunc, output=np.float32) - # lateral downsampling + # 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_img[z, ...] = resize(img[z, ...], output_shape=tuple(iso_shape[1:]), - anti_aliasing=anti_aliasing, preserve_range=True) + iso_img[z, ...] = \ + resize(img[z, ...], output_shape=tuple(iso_shape[1:]), anti_aliasing=anti_alias, preserve_range=True) - # resize padding matrix - if pad_mat is not None: - rsz_pad_mat = (np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad_mat))).astype(int) - return iso_img, rsz_pad_mat + # resize padding array accordingly + if pad is not None: + rsz_pad = np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad)).astype(int) + return iso_img, rsz_pad else: return iso_img, None diff --git a/foa3d/printing.py b/foa3d/printing.py index 357da0a..fa23e1d 100644 --- a/foa3d/printing.py +++ b/foa3d/printing.py @@ -89,8 +89,8 @@ def print_frangi_info(alpha, beta, gamma, scales_um, image_shape_um, in_slice_sh print(u"\n\u03B1: {0:.3f}\n".format(alpha) + u"\u03B2: {0:.3f}\n".format(beta) + u"\u03B3: {0}\n".format(gamma)) - print("Scales [\u03BCm]: {}".format(scales_um)) - print("Diameters [\u03BCm]: {}".format(4 * scales_um)) + print("Enhanced scales [\u03BCm]: {}".format(scales_um)) + 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) @@ -116,6 +116,25 @@ def print_analysis_time(start_time): print("\nVolume image analyzed in: {0} min {1:3.1f} s\n".format(mins, secs)) +def print_blur(smooth_sigma_um): + """ + Print gaussian lowpass filter standard deviation. + + Parameters + ---------- + smooth_sigma_um: numpy.ndarray (shape=(3,), dtype=int) + 3D standard deviation of the low-pass Gaussian filter [μm] + (resolution anisotropy correction) + + Returns + ------- + None + """ + print("\n Z Y X") + print("Gaussian blur \u03C3 [μm]: ({0:.3f}, {1:.3f}, {2:.3f})" + .format(smooth_sigma_um[0], smooth_sigma_um[1], smooth_sigma_um[2]), end='\r') + + def print_import_time(start_time): """ Print volume image import time. @@ -158,10 +177,6 @@ def print_pipeline_heading(): """ Print Foa3D pipeline heading. - Parameters - ---------- - None - Returns ------- None @@ -173,18 +188,14 @@ def print_prepro_heading(): """ Print preprocessing heading. - Parameters - ---------- - None - Returns ------- None """ - print(color_text(0, 191, 255, "\n\nMicroscopy Volume Image Preprocessing"), end='\r') + print(color_text(0, 191, 255, "\n\nMicroscopy Volume Image Preprocessing")) -def print_resolution(px_size, psf_fwhm): +def print_native_res(px_size, psf_fwhm): """ Print pixel and optical resolution of the microscopy system. @@ -204,6 +215,26 @@ def print_resolution(px_size, psf_fwhm): 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, psf_fwhm): + """ + Print adjusted isotropic spatial resolution. + + Parameters + ---------- + px_sz_iso: numpy.ndarray (shape=(3,), dtype=float) + new isotropic pixel size [μm] + + psf_fwhm: numpy.ndarray (shape=(3,), dtype=float) + 3D FWHM of the PSF [μm] + + Returns + ------- + None + """ + print("\nAdjusted pixel size [μm]: ({0:.3f}, {1:.3f}, {2:.3f})".format(px_sz_iso[0], px_sz_iso[1], px_sz_iso[2])) + print("Adjusted PSF FWHM [μm]: ({0:.3f}, {0:.3f}, {0:.3f})\n".format(psf_fwhm[0])) + + def print_slicing_info(image_shape_um, slice_shape_um, tot_slice_num, px_size, image_item_size): """ Print information on the slicing of the basic image sub-volumes @@ -247,11 +278,11 @@ def print_slicing_info(image_shape_um, slice_shape_um, tot_slice_num, px_size, i .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("Basic slice shape [μm]: ({0:.1f}, {1:.1f}, {2:.1f})" + 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("Basic slice size [MB]: {0}" + print("Image slice size [MB]: {0}" .format(np.ceil(max_slice_size / 1024**2).astype(int))) - print("Basic slice number: {0}\n" + print("Image slice number: {0}\n" .format(tot_slice_num)) @@ -268,13 +299,11 @@ def print_soma_masking(lpf_soma_mask): ------- None """ - if lpf_soma_mask: - print("Lipofuscin-based soma masking: ON\n") - else: - print("Lipofuscin-based soma masking: OFF\n") + prt = 'Lipofuscin soma mask: ' + print('{}active\n'.format(prt)) if lpf_soma_mask else print('{}not active\n'.format(prt)) -def print_image_shape(cli_args, img, mosaic, channel_ax=None): +def print_image_shape(cli_args, img, is_tiled, channel_ax=None): """ Print volume image shape. @@ -286,7 +315,7 @@ def print_image_shape(cli_args, img, mosaic, channel_ax=None): img: numpy.ndarray (shape=(Z,Y,X)) microscopy volume image - mosaic: bool + is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher channel_ax: int @@ -304,7 +333,7 @@ def print_image_shape(cli_args, img, mosaic, channel_ax=None): # adapt axis order img_shape = img.shape if len(img_shape) == 4: - channel_ax = 1 if mosaic else -1 + channel_ax = 1 if is_tiled else -1 # get image shape if channel_ax is not None: diff --git a/foa3d/slicing.py b/foa3d/slicing.py index ec0b53b..0c67018 100644 --- a/foa3d/slicing.py +++ b/foa3d/slicing.py @@ -6,7 +6,7 @@ from foa3d.utils import get_available_cores -def compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip=False, ovlp_rng=0): +def compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip=False, ovlp=0): """ Adjust image slice coordinates at boundaries. @@ -30,8 +30,8 @@ def compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip= flip: bool if True, flip axis coordinates - ovlp_rng: int - optional slice sampling range + ovlp: int + optional slicing range extension along each axis (on each side) Returns @@ -43,7 +43,7 @@ def compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip= adjusted stop index ovlp: numpy.ndarray (shape=(2,), dtype=int) - lower and upper overlapped boundaries + lower and upper padded boundaries along axis """ @@ -58,18 +58,18 @@ def compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip= stop = img_shape[ax] - start_tmp # adjust start and stop coordinates - start, stop, ovlp = trim_axis_range(img_shape, start, stop, ax=ax, ovlp_rng=ovlp_rng) + start, stop, pad = adjust_axis_range(img_shape, start, stop, ax=ax, ovlp=ovlp) # handle image shape residuals at boundaries if ax_iter[ax] == slice_per_dim[ax] - 1: if np.remainder(img_shape[ax], slice_shape[ax]) > 0: stop = img_shape[ax] - ovlp[1] = 0 + pad[1] = ovlp - return start, stop, ovlp + return start, stop, pad -def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp_rng=0, flip=False): +def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp=0, flip=False): """ Compute basic slice coordinates from microscopy volume image. @@ -88,9 +88,9 @@ def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp_rng slice_per_dim: numpy.ndarray (shape=(3,), dtype=int) number of image slices along each axis - ovlp_rng: int - optional slice sampling range - extension along each axis (on each side) + ovlp: int + slice overlap range + along each axis side flip: bool if True, flip axes coordinates @@ -100,8 +100,8 @@ def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp_rng rng: tuple 3D slice index ranges - ovlp: np.ndarray (shape=(3,2), dtype=int) - overlapped boundaries + pad: np.ndarray (shape=(3,2), dtype=int) + padded boundaries """ # adjust original image patch coordinates @@ -109,11 +109,11 @@ def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp_rng dims = len(ax_iter) start = np.zeros((dims,), dtype=int) stop = np.zeros((dims,), dtype=int) - ovlp = np.zeros((dims, 2), dtype=int) + pad = np.zeros((dims, 2), dtype=int) slc = tuple() for ax in range(dims): - start[ax], stop[ax], ovlp[ax] = \ - compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip=flip, ovlp_rng=ovlp_rng) + start[ax], stop[ax], pad[ax] = \ + compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip=flip, ovlp=ovlp) slc += (slice(start[ax], stop[ax], 1),) # generate tuple of slice index ranges @@ -122,12 +122,12 @@ def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp_rng # handle invalid slices for r in rng: if r.start is None: - return None, ovlp + return None, pad - return rng, ovlp + return rng, pad -def compute_slice_shape(img_shape, item_size, px_size=None, max_size=100, ovlp_rng=0): +def compute_slice_shape(img_shape, item_size, px_size=None, max_size=1e6, ovlp=0): """ Compute basic image chunk shape depending on its maximum size (in bytes). @@ -142,12 +142,13 @@ def compute_slice_shape(img_shape, item_size, px_size=None, max_size=100, ovlp_r px_size: numpy.ndarray (shape=(3,), dtype=float) pixel size [μm] - max_slice_size: float + max_size: float maximum memory size (in bytes) of the basic image slices analyzed using parallel threads - ovlp_rng: int - slice overlap range along each axis + ovlp: int + slice overlap range + along each axis side Returns ------- @@ -161,11 +162,12 @@ def compute_slice_shape(img_shape, item_size, px_size=None, max_size=100, ovlp_r (if px_size is provided) """ if len(img_shape) == 4: - item_size = item_size * 3 + item_size *= 3 - slice_depth = img_shape[0] - slice_side = np.round(1024 * np.sqrt((max_size / (slice_depth * item_size))) - 2 * ovlp_rng) - slice_shape = np.array([slice_depth, slice_side, slice_side]).astype(int) + tot_ovlp = 2 * ovlp + slice_depth = img_shape[0] + tot_ovlp + slice_side = np.round(np.sqrt((max_size / (slice_depth * item_size)))) + slice_shape = np.array([slice_depth, slice_side, slice_side]).astype(int) - tot_ovlp slice_shape = np.min(np.stack((img_shape[:3], slice_shape)), axis=0) if px_size is not None: @@ -176,7 +178,7 @@ def compute_slice_shape(img_shape, item_size, px_size=None, max_size=100, ovlp_r return slice_shape -def compute_overlap_range(smooth_sigma, frangi_sigma, truncate=4): +def compute_overlap_range(smooth_sigma, frangi_sigma, px_rsz_ratio, truncate=4): """ Compute lateral slice extension range for coping with smoothing-related boundary artifacts. @@ -190,24 +192,37 @@ def compute_overlap_range(smooth_sigma, frangi_sigma, truncate=4): frangi_sigma: numpy.ndarray (dtype=int) analyzed spatial scales [px] + px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) + 3D image resize ratio + truncate: int neglect the Gaussian smoothing kernel after this many standard deviations Returns ------- - ext_rng: int - slice extended index range (on each side) + ovlp: int + slice overlap range + along each axis side + + rsz_ovlp: numpy.ndarray (shape=(3,), dtype=int) + resized image slice overlap range (if px_rsz_ratio is not None) """ max_sigma = np.max(np.concatenate((smooth_sigma, frangi_sigma))) - ext_rng = int(np.ceil(2 * truncate * max_sigma) // 2) if smooth_sigma is not None else 0 + ovlp = int(np.ceil(2 * truncate * max_sigma) // 2) if smooth_sigma is not None else 0 - return ext_rng + if px_rsz_ratio is not None: + rsz_ovlp = np.multiply(ovlp * np.ones((3,)), px_rsz_ratio).astype(int) + + return ovlp, rsz_ovlp + + else: + return ovlp def config_frangi_batch(frangi_scales, mem_growth_factor=149.7, mem_fudge_factor=1.0, - min_slice_size_mb=-1, jobs=None, max_ram_mb=None): + min_slice_size=-1, jobs=None, ram=None): """ Compute size and number of the batches of basic microscopy image slices analyzed in parallel. @@ -224,28 +239,28 @@ def config_frangi_batch(frangi_scales, mem_growth_factor=149.7, mem_fudge_factor mem_fudge_factor: float memory fudge factor - min_slice_size_mb: float - minimum slice size in [MB] + min_slice_size: float + minimum slice size in [B] jobs: int number of parallel jobs (threads) used by the Frangi filtering stage - max_ram_mb: float - maximum RAM available to the Frangi filtering stage [MB] + ram: float + maximum RAM available to the Frangi filtering stage [B] Returns ------- slice_batch_size: int slice batch size - slice_size_mb: float - memory size (in megabytes) of the basic image slices + slice_size: float + memory size (in bytes) of the basic image slices fed to the Frangi filter """ # maximum RAM not provided: use all - if max_ram_mb is None: - max_ram_mb = psutil.virtual_memory()[1] / 1e6 + if ram is None: + ram = psutil.virtual_memory()[1] # number of logical cores num_cpu = get_available_cores() @@ -261,15 +276,15 @@ def config_frangi_batch(frangi_scales, mem_growth_factor=149.7, mem_fudge_factor slice_batch_size = 1 # get image slice size - slice_size_mb = get_slice_size(max_ram_mb, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales) - while slice_size_mb < min_slice_size_mb: + slice_size = get_slice_size(ram, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales) + while slice_size < min_slice_size: slice_batch_size -= 1 - slice_size_mb = get_slice_size(max_ram_mb, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales) + slice_size = get_slice_size(ram, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales) - return slice_batch_size, slice_size_mb + return slice_batch_size, slice_size -def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sigma, frangi_sigma, lpf_soma_mask, +def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sigma, frangi_sigma, mask_lpf, batch_size, slice_size): """ Image slicing configuration for the parallel Frangi filtering of basic chunks @@ -296,6 +311,10 @@ def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sig frangi_sigma: numpy.ndarray (or int) Frangi filter scales [px] + mask_lpf: bool + if True, mask neuronal bodies exploiting the autofluorescence + signal of lipofuscin pigments + batch_size: int slice batch size @@ -314,8 +333,8 @@ def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sig out_rng_lst: list list of output slice index ranges - in_ovlp_lst: list - list of slice overlap arrays + in_pad_lst: list + list of slice padding arrays in_slice_shape_um: numpy.ndarray (shape=(3,), dtype=float) shape of the basic image slices analyzed iteratively [μm] @@ -326,6 +345,9 @@ def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sig px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float) 3D image resize ratio + rsz_ovlp: numpy.ndarray (shape=(3,), dtype=int) + resized image slice overlap range + tot_slice_num: int total number of analyzed image slices @@ -333,17 +355,19 @@ def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sig adjusted slice batch size """ - # compute input patch padding range (border artifacts suppression) - ovlp_rng = compute_overlap_range(smooth_sigma, frangi_sigma) + # get pixel resize ratio + px_rsz_ratio = np.divide(px_size, px_size_iso) + + # compute input slice overlap (border artifacts suppression) + ovlp, rsz_ovlp = compute_overlap_range(smooth_sigma, frangi_sigma, px_rsz_ratio) - # shape of processed TPFM slices + # shape of processed microscopy image slices in_slice_shape, in_slice_shape_um = compute_slice_shape(img_shape, item_size, - px_size=px_size, max_size=slice_size, ovlp_rng=ovlp_rng) + px_size=px_size, max_size=slice_size, ovlp=ovlp) # adjust shapes according to the anisotropy correction - px_rsz_ratio = np.divide(px_size, px_size_iso) - out_slice_shape = np.ceil(px_rsz_ratio * in_slice_shape).astype(int) - out_image_shape = np.ceil(px_rsz_ratio * img_shape).astype(int) + out_slice_shape = np.ceil(np.multiply(px_rsz_ratio, in_slice_shape)).astype(int) + out_image_shape = np.ceil(np.multiply(px_rsz_ratio, img_shape)).astype(int) # iteratively define input/output slice 3D ranges slice_per_dim = np.ceil(np.divide(img_shape, in_slice_shape)).astype(int) @@ -351,26 +375,26 @@ def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sig # initialize empty range lists out_rng_lst = list() - in_ovlp_lst = list() + in_pad_lst = list() in_rng_lst = list() in_rng_lst_neu = list() for z, y, x in product(range(slice_per_dim[0]), range(slice_per_dim[1]), range(slice_per_dim[2])): # index ranges of analyzed fiber image slices (with padding) - in_rng, ovlp = compute_slice_range((z, y, x), in_slice_shape, img_shape, slice_per_dim, ovlp_rng=ovlp_rng) + in_rng, pad = compute_slice_range((z, y, x), in_slice_shape, img_shape, slice_per_dim, ovlp=ovlp) # get output slice ranges # for valid input range instances if in_rng is not None: in_rng_lst.append(in_rng) - in_ovlp_lst.append(ovlp) + in_pad_lst.append(pad) # output index ranges out_rng, _ = compute_slice_range((z, y, x), out_slice_shape, out_image_shape, slice_per_dim) out_rng_lst.append(out_rng) # optional neuron masking - if lpf_soma_mask: + if mask_lpf: in_rng_neu, _ = compute_slice_range((z, y, x), in_slice_shape, img_shape, slice_per_dim) in_rng_lst_neu.append(in_rng_neu) else: @@ -382,13 +406,13 @@ def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sig if batch_size > tot_slice_num: batch_size = tot_slice_num - return in_rng_lst, in_rng_lst_neu, out_rng_lst, in_ovlp_lst, \ - in_slice_shape_um, out_slice_shape, px_rsz_ratio, tot_slice_num, batch_size + return in_rng_lst, in_rng_lst_neu, out_rng_lst, in_pad_lst, \ + in_slice_shape_um, out_slice_shape, px_rsz_ratio, rsz_ovlp, tot_slice_num, batch_size def crop_slice(img_slice, rng, ovlp=None, flipped=()): """ - Shrink image slice at volume boundaries, for overall shape consistency. + Shrink image slice at total volume boundaries, for overall shape consistency. Parameters ---------- @@ -398,8 +422,8 @@ def crop_slice(img_slice, rng, ovlp=None, flipped=()): rng: tuple 3D slice index range - ovlp: np.ndarray (shape=(3,2), dtype=int) - overlapped boundaries + ovlp: numpy.ndarray (shape=(3,), dtype=int) + slice overlap range (on each axis side) flipped: tuple flipped axes @@ -412,9 +436,9 @@ def crop_slice(img_slice, rng, ovlp=None, flipped=()): # delete overlapping slice boundaries if ovlp is not None: - img_slice = img_slice[ovlp[0, 0]:img_slice.shape[0] - ovlp[0, 1], - ovlp[1, 0]:img_slice.shape[1] - ovlp[1, 1], - ovlp[2, 0]:img_slice.shape[2] - ovlp[2, 1]] + img_slice = img_slice[ovlp[0]:img_slice.shape[0] - ovlp[0], + ovlp[1]:img_slice.shape[1] - ovlp[1], + ovlp[2]:img_slice.shape[2] - ovlp[2]] # check slice shape and output index ranges out_slice_shape = img_slice.shape @@ -433,6 +457,37 @@ def crop_slice(img_slice, rng, ovlp=None, flipped=()): return cropped +def crop_slice_lst(img_slice_lst, rng, ovlp=None, flipped=()): + """ + Shrink list of image slices at total volume boundaries, for overall shape consistency. + + Parameters + ---------- + img_slice_lst: list + list of image slices generated by + the Frangi filtering stage to be cropped + + rng: tuple + 3D slice index range + + ovlp: numpy.ndarray (shape=(3,), dtype=int) + slice overlap range (on each axis side) + + flipped: tuple + flipped axes + + Returns + ------- + img_slice_lst: list + list of cropped image slices + """ + + for s, img_slice in enumerate(img_slice_lst): + img_slice_lst[s] = crop_slice(img_slice, rng, ovlp=ovlp, flipped=flipped) + + return img_slice_lst + + def get_slice_size(max_ram, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales=1): """ Compute the size of the basic microscopy image slices fed to the Frangi filtering stage. @@ -440,7 +495,7 @@ def get_slice_size(max_ram, mem_growth_factor, mem_fudge_factor, slice_batch_siz Parameters ---------- max_ram: float - available RAM + available RAM [B] mem_growth_factor: float empirical memory growth factor @@ -458,7 +513,7 @@ def get_slice_size(max_ram, mem_growth_factor, mem_fudge_factor, slice_batch_siz Returns ------- slice_size: float - memory size (in megabytes) of the basic image slices + memory size (in bytes) of the basic image slices fed to the pipeline stage """ slice_size = max_ram / (slice_batch_size * mem_growth_factor * mem_fudge_factor * num_scales) @@ -466,7 +521,7 @@ def get_slice_size(max_ram, mem_growth_factor, mem_fudge_factor, slice_batch_siz return slice_size -def slice_channel(img, rng, channel, mosaic=False): +def slice_channel(img, rng, channel, is_tiled=False): """ Slice desired channel from input image volume. @@ -481,28 +536,28 @@ def slice_channel(img, rng, channel, mosaic=False): channel: int image channel axis - mosaic: bool + is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher Returns ------- - image_slice: numpy.ndarray + img_slice: numpy.ndarray sliced image patch """ z_rng, r_rng, c_rng = rng if channel is None: - image_slice = img[z_rng, r_rng, c_rng] + img_slice = img[z_rng, r_rng, c_rng] else: - if mosaic: - image_slice = img[z_rng, channel, r_rng, c_rng] + if is_tiled: + img_slice = img[z_rng, channel, r_rng, c_rng] else: - image_slice = img[z_rng, r_rng, c_rng, channel] + img_slice = img[z_rng, r_rng, c_rng, channel] - return image_slice + return img_slice -def trim_axis_range(img_shape, start, stop, ax, ovlp_rng=0): +def adjust_axis_range(img_shape, start, stop, ax, ovlp=0): """ Trim slice axis range at total image boundaries. @@ -520,9 +575,9 @@ def trim_axis_range(img_shape, start, stop, ax, ovlp_rng=0): ax: int image axis - ovlp_rng: int - optional slice sampling range - extension along each axis (on each side) + ovlp: int + slice overlap range + along each axis side Returns ------- @@ -532,30 +587,23 @@ def trim_axis_range(img_shape, start, stop, ax, ovlp_rng=0): stop: int adjusted slice stop coordinate - ovlp: numpy.ndarray (shape=(2,), dtype=int) - lower and upper overlap ranges + pad: numpy.ndarray (shape=(2,), dtype=int) + lower and upper padding ranges """ - # initialize slice overlap array - ovlp = np.zeros((2,), dtype=int) + # initialize slice padding array + pad = np.zeros((2,), dtype=int) - # skip z-axis - if ax != 0: + # adjust start coordinate + start -= ovlp + if start < 0: + pad[0] = -start + start = 0 - # adjust start coordinate - start -= ovlp_rng - if start < 0: - ovlp[0] = ovlp_rng + start - start = 0 - else: - ovlp[0] = ovlp_rng - - # adjust stop coordinate - stop += ovlp_rng - if stop > img_shape[ax]: - ovlp[1] = img_shape[ax] - (stop - ovlp_rng) - stop = img_shape[ax] - else: - ovlp[1] = ovlp_rng + # adjust stop coordinate + stop += ovlp + if stop > img_shape[ax]: + pad[1] = stop - img_shape[ax] + stop = img_shape[ax] - return start, stop, ovlp + return start, stop, pad diff --git a/foa3d/utils.py b/foa3d/utils.py index 4754a82..7182e7e 100644 --- a/foa3d/utils.py +++ b/foa3d/utils.py @@ -15,43 +15,43 @@ threshold_yen) -def create_background_mask(img, thresh_method='yen'): +def create_background_mask(img, method='yen'): """ Compute background mask. Parameters ---------- - img: numpy.ndarray (shape=(Z,Y,X)) + img: numpy.ndarray (axis order=(Z,Y,X)) microscopy volume image - thresh_method: str + method: str image thresholding method Returns ------- - background_mask: numpy.ndarray (shape=(Z,Y,X), dtype=bool) + bg_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) boolean background mask """ # select thresholding method - if thresh_method == 'li': + if method == 'li': initial_li_guess = np.mean(img[img != 0]) thresh = threshold_li(img, initial_guess=initial_li_guess) - elif thresh_method == 'niblack': + elif method == 'niblack': thresh = threshold_niblack(img, window_size=15, k=0.2) - elif thresh_method == 'sauvola': + elif method == 'sauvola': thresh = threshold_sauvola(img, window_size=15, k=0.2, r=None) - elif thresh_method == 'triangle': + elif method == 'triangle': thresh = threshold_triangle(img, nbins=256) - elif thresh_method == 'yen': + elif method == 'yen': thresh = threshold_yen(img, nbins=256) else: - raise ValueError(" Unsupported thresholding method!!!") + raise ValueError("Unsupported thresholding method!!!") # compute mask - background_mask = img < thresh + bg_msk = img < thresh - return background_mask + return bg_msk def create_hdf5_dset(dset_shape, dtype, chunks=True, name='tmp', tmp=None): @@ -69,6 +69,12 @@ def create_hdf5_dset(dset_shape, dtype, chunks=True, name='tmp', tmp=None): chunks: tuple (dtype: int) or bool shape of the chunked storage layout (default: auto chunking) + name: str + filename + + tmp: str + path to existing temporary saving directory + Returns ------- dset: @@ -77,18 +83,17 @@ def create_hdf5_dset(dset_shape, dtype, chunks=True, name='tmp', tmp=None): file_path: str path to the HDF5 file """ - if tmp is None: tmp = tempfile.mkdtemp() - file_path = path.join(tmp, name + '.h5') + file_path = path.join(tmp, '{}.h5'.format(name)) file = File(file_path, 'w') - dset = file.create_dataset(None, dset_shape, chunks=chunks, dtype=dtype) + dset = file.create_dataset(None, dset_shape, chunks=tuple(chunks), dtype=dtype) return dset, file_path -def create_memory_map(shape, dtype, name='tmp', tmp=None, arr=None, mmap_mode='r+'): +def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mode='r+'): """ Create a memory-map to an array stored in a binary file on disk. @@ -103,7 +108,7 @@ def create_memory_map(shape, dtype, name='tmp', tmp=None, arr=None, mmap_mode='r name: str optional temporary filename - tmp: str + tmp_dir: str temporary file directory arr: numpy.ndarray @@ -118,14 +123,16 @@ def create_memory_map(shape, dtype, name='tmp', tmp=None, arr=None, mmap_mode='r memory-mapped array """ - if tmp is None: - tmp = tempfile.mkdtemp() - mmap_path = path.join(tmp, name + '.mmap') + if tmp_dir is None: + tmp_dir = tempfile.mkdtemp() + mmap_path = path.join(tmp_dir, name + '.mmap') if path.exists(mmap_path): unlink(mmap_path) + if arr is None: arr = np.zeros(tuple(shape), dtype=dtype) + _ = dump(arr, mmap_path) mmap = load(mmap_path, mmap_mode=mmap_mode) del arr @@ -138,10 +145,6 @@ def get_available_cores(): """ Return the number of available logical cores. - Parameters - ---------- - None - Returns ------- num_cpu: int @@ -149,10 +152,7 @@ def get_available_cores(): """ num_cpu = environ.pop('OMP_NUM_THREADS', default=None) - if num_cpu is None: - num_cpu = cpu_count() - else: - num_cpu = int(num_cpu) + num_cpu = cpu_count() if num_cpu is None else int(num_cpu) return num_cpu @@ -299,20 +299,17 @@ def get_item_bytes(data): Returns ------- - bytes: int + bts: int item size in bytes """ - # get data type - data_type = data.dtype - # type byte size try: - bytes = int(np.iinfo(data_type).bits / 8) + bts = int(np.iinfo(data.dtype).bits / 8) except ValueError: - bytes = int(np.finfo(data_type).bits / 8) + bts = int(np.finfo(data.dtype).bits / 8) - return bytes + return bts def get_output_prefix(scales_um, alpha, beta, gamma): @@ -342,8 +339,9 @@ def get_output_prefix(scales_um, alpha, beta, gamma): pfx = 'sc' for s in scales_um: - pfx = pfx + str(s) + '_' - pfx = 'a' + str(alpha) + '_b' + str(beta) + '_g' + str(gamma) + '_' + pfx + pfx += '{}_'.format(s) + + pfx = 'a{}_b{}_g{}_{}'.format(alpha, beta, gamma, pfx) return pfx @@ -377,12 +375,12 @@ def normalize_angle(angle, lower=0.0, upper=360.0, dtype=None): if lower >= upper """ # convert to array if needed - isList = False + is_list = False if np.isscalar(angle): angle = np.array(angle) elif isinstance(angle, list): angle = np.array(angle) - isList = True + is_list = True # check limits if lower >= upper: @@ -407,7 +405,7 @@ def normalize_angle(angle, lower=0.0, upper=360.0, dtype=None): norm_angle = norm_angle.astype(dtype) # convert back to list - if isList: + if is_list: norm_angle = list(norm_angle) return norm_angle @@ -438,7 +436,7 @@ def normalize_image(img, max_out_val=255.0, dtype=np.uint8): min_val = np.min(img) max_val = np.max(img) - # normalization + # normalize if max_val != 0: if max_val != min_val: norm_img = (((img - min_val) / (max_val - min_val)) * max_out_val).astype(dtype) @@ -450,7 +448,7 @@ def normalize_image(img, max_out_val=255.0, dtype=np.uint8): return norm_img -def orient_colormap(vec_img): +def hsv_orient_cmap(vec_img): """ Compute HSV colormap of vector orientations from 3D vector field. @@ -465,12 +463,8 @@ def orient_colormap(vec_img): orientation color map """ - # get input array shape - vec_img_shape = vec_img.shape - - # select planar components - vy = vec_img[..., 1] - vx = vec_img[..., 2] + # extract planar components + vy, vx = (vec_img[..., 1], vec_img[..., 2]) # compute the in-plane versor length vxy_abs = np.sqrt(np.square(vx) + np.square(vy)) @@ -481,16 +475,13 @@ def orient_colormap(vec_img): vxy_ang = np.divide(vxy_ang, np.pi) # initialize colormap - rgb_map = np.zeros(shape=tuple(list(vec_img_shape[:-1]) + [3]), dtype=np.uint8) - for z in range(vec_img_shape[0]): + rgb_map = np.zeros(shape=tuple(list(vec_img.shape[:-1]) + [3]), dtype=np.uint8) + for z in range(vec_img.shape[0]): - # generate colormap slice by slice - h = vxy_ang[z] - s = vxy_abs[z] - v = s - hsv_map = np.stack((h, s, v), axis=-1) + # generate HSV colormap slice by slice + hsv_map = np.stack((vxy_ang[z], vxy_abs[z], vxy_abs[z]), axis=-1) - # conversion to 8-bit precision + # convert to RGB map with 8-bit precision rgb_map[z] = (255.0 * hsv_to_rgb(hsv_map)).astype(np.uint8) return rgb_map @@ -558,7 +549,7 @@ def transform_axes(nd_array, flipped=None, swapped=None, expand=None): return nd_array -def vector_colormap(vec_img): +def rgb_orient_cmap(vec_img, minimum=0, stretch=1, q=8): """ Compute RGB colormap of orientation vector components from 3D vector field. @@ -567,26 +558,30 @@ def vector_colormap(vec_img): vec_img: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float) n-dimensional array of orientation vectors + minimum: int + intensity that should be mapped to black (a scalar or array for R, G, B) + + stretch: int + linear stretch of the image + + q: int + asinh softening parameter + Returns ------- rgb_map: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8) orientation color map """ - # get input array shape - vec_img_shape = vec_img.shape - # take absolute value vec_img = np.abs(vec_img) # initialize colormap - rgb_map = np.zeros(shape=vec_img_shape, dtype=np.uint8) - for z in range(vec_img_shape[0]): + rgb_map = np.zeros(shape=vec_img.shape, dtype=np.uint8) + for z in range(vec_img.shape[0]): # generate colormap slice by slice - img_r = vec_img[z, :, :, 2] - img_g = vec_img[z, :, :, 1] - img_b = vec_img[z, :, :, 0] - rgb_map[z] = make_lupton_rgb(img_r, img_g, img_b, minimum=0, stretch=1, Q=8) + img_r, img_g, img_b = (vec_img[z, :, :, 2], vec_img[z, :, :, 1], vec_img[z, :, :, 0]) + rgb_map[z] = make_lupton_rgb(img_r, img_g, img_b, minimum=minimum, stretch=stretch, Q=q) return rgb_map diff --git a/requirements.in b/requirements.in index 29171e3..36a724e 100644 --- a/requirements.in +++ b/requirements.in @@ -7,3 +7,4 @@ nibabel numba psutil scikit-image +zetastitcher==0.7.0.post1 diff --git a/requirements.txt b/requirements.txt index b5b708c..7e8828f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,6 @@ # -# This file is autogenerated by pip-compile with python 3.8 -# To update, run: +# This file is autogenerated by pip-compile with Python 3.9 +# by the following command: # # pip-compile requirements.in # @@ -10,20 +10,39 @@ alive-progress==2.4.1 # via -r requirements.in astropy==5.1 # via -r requirements.in +cachetools==5.3.1 + # via zetastitcher +clarabel==0.6.0 + # via cvxpy +coloredlogs==15.0.1 + # via zetastitcher contourpy==1.0.5 # via matplotlib +cvxpy==1.4.1 + # via zetastitcher cycler==0.11.0 # via matplotlib +daqp==0.5.1 + # via qpsolvers +ecos==2.0.12 + # via + # cvxpy + # qpsolvers fonttools==4.37.4 # via matplotlib grapheme==0.6.0 # via alive-progress h5py==3.7.0 # via -r requirements.in +humanfriendly==10.0 + # via coloredlogs +humanize==4.8.0 + # via zetastitcher imageio==2.22.1 - # via scikit-image -importlib-metadata==5.0.0 - # via numba + # via + # pims + # scikit-image + # zetastitcher joblib==1.2.0 # via -r requirements.in kiwisolver==1.4.4 @@ -33,7 +52,9 @@ llvmlite==0.39.1 matplotlib==3.6.0 # via -r requirements.in networkx==2.8.7 - # via scikit-image + # via + # scikit-image + # zetastitcher nibabel==4.0.2 # via -r requirements.in numba==0.56.2 @@ -41,52 +62,106 @@ numba==0.56.2 numpy==1.23.3 # via # astropy + # clarabel # contourpy + # cvxpy + # ecos # h5py # imageio # matplotlib # nibabel # numba + # opencv-python + # osqp + # pandas + # pims # pyerfa # pywavelets + # qpsolvers # scikit-image # scipy + # scs # tifffile + # zetastitcher +opencv-python==4.8.1.78 + # via zetastitcher +osqp==0.6.3 + # via + # cvxpy + # qpsolvers packaging==21.3 # via # astropy # matplotlib # nibabel # scikit-image +pandas==2.1.1 + # via zetastitcher pillow==9.2.0 # via # imageio # matplotlib # scikit-image +pims==0.6.1 + # via zetastitcher psutil==5.9.4 - # via -r requirements.in + # via + # -r requirements.in + # zetastitcher +pybind11==2.11.1 + # via cvxpy pyerfa==2.0.0.1 # via astropy pyparsing==3.0.9 # via # matplotlib # packaging +pyreadline3==3.4.1 + # via humanfriendly python-dateutil==2.8.2 - # via matplotlib + # via + # matplotlib + # pandas +pytz==2023.3.post1 + # via pandas pywavelets==1.4.1 # via scikit-image pyyaml==6.0 - # via astropy + # via + # astropy + # zetastitcher +qdldl==0.1.7.post0 + # via osqp +qpsolvers==4.0.0 + # via zetastitcher scikit-image==0.19.3 # via -r requirements.in scipy==1.9.1 - # via scikit-image + # via + # clarabel + # cvxpy + # ecos + # osqp + # qpsolvers + # scikit-image + # scs + # zetastitcher +scs==3.2.3 + # via + # cvxpy + # qpsolvers six==1.16.0 # via python-dateutil +slicerator==1.1.0 + # via pims tifffile==2022.8.12 - # via scikit-image -zipp==3.8.1 - # via importlib-metadata + # via + # scikit-image + # zetastitcher +tzdata==2023.3 + # via pandas +zetastitcher==0.7.0.post1 + # via -r requirements.in # The following packages are considered to be unsafe in a requirements file: # setuptools