From 8a894f44bcb0dc44323ecc4e87178f11f723e971 Mon Sep 17 00:00:00 2001 From: msorelli Date: Fri, 28 Jun 2024 11:16:14 +0200 Subject: [PATCH] Add new background rejection mechanism --- foa3d/__main__.py | 6 +++--- foa3d/input.py | 49 +++++++++++++++++++++++++++++++++++------- foa3d/pipeline.py | 29 ++++++++++++++++++------- foa3d/preprocessing.py | 8 +++---- foa3d/slicing.py | 12 +++++++++-- foa3d/utils.py | 7 ++++-- 6 files changed, 83 insertions(+), 28 deletions(-) diff --git a/foa3d/__main__.py b/foa3d/__main__.py index eb892ad..698cc32 100644 --- a/foa3d/__main__.py +++ b/foa3d/__main__.py @@ -7,7 +7,7 @@ def foa3d(cli_args): # load microscopy volume image or array of fiber orientation vectors - img, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args) + img, tissue_msk, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args) # get resources configuration ram, jobs = get_resource_config(cli_args) @@ -15,8 +15,8 @@ def foa3d(cli_args): # conduct parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices if not is_fiber: fiber_vec_img, iso_fiber_img, px_sz, img_name \ - = parallel_frangi_on_slices(img, cli_args, save_dir[0], tmp_dir, img_name, ram=ram, jobs=jobs, - is_tiled=is_tiled) + = parallel_frangi_on_slices(img, cli_args, save_dir[0], tmp_dir, img_name, tissue_msk=tissue_msk, + ram=ram, jobs=jobs, is_tiled=is_tiled) else: fiber_vec_img, iso_fiber_img, px_sz = (img, None, None) diff --git a/foa3d/input.py b/foa3d/input.py index 78c1e8e..461adea 100644 --- a/foa3d/input.py +++ b/foa3d/input.py @@ -17,7 +17,8 @@ from foa3d.preprocessing import config_anisotropy_correction from foa3d.printing import (color_text, print_image_shape, print_import_time, print_native_res) -from foa3d.utils import create_memory_map, get_item_bytes, get_output_prefix +from foa3d.utils import (create_background_mask, create_memory_map, + get_item_bytes, get_output_prefix) class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawTextHelpFormatter): @@ -90,7 +91,9 @@ def get_cli_parser(): help='degrees of the spherical harmonics series expansion (even number between 2 and 10)') cli_parser.add_argument('-o', '--out', type=str, default=None, help='output directory') - + cli_parser.add_argument('-t', '--tissue-msk', action='store_true', default=False, + help='apply tissue reconstruction mask (binarized MIP)') + # parse arguments cli_args = cli_parser.parse_args() @@ -186,6 +189,9 @@ def get_file_info(cli_args): create a memory-mapped array of the microscopy volume image, increasing the parallel processing performance (the image will be preliminarily loaded to RAM) + + mip_msk: bool + apply tissue reconstruction mask (binarized MIP) """ # get microscopy image path and name @@ -202,7 +208,10 @@ def get_file_info(cli_args): img_name = img_name.replace('.{}'.format(img_fmt), '') is_tiled = True if img_fmt == 'yml' else False - return img_path, img_name, img_fmt, is_tiled, is_mmap + # apply tissue reconstruction mask (binarized MIP) + mip_msk = cli_args.tissue_msk + + return img_path, img_name, img_fmt, is_tiled, is_mmap, mip_msk def get_frangi_config(cli_args, img_name): @@ -368,6 +377,9 @@ def load_microscopy_image(cli_args): img: numpy.ndarray or NumPy memory-map object microscopy volume image or array of fiber orientation vectors + tissue_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + is_tiled: bool True for tiled microscopy reconstructions aligned using ZetaStitcher @@ -392,16 +404,18 @@ def load_microscopy_image(cli_args): tmp_dir = tempfile.mkdtemp() # retrieve input file information - img_path, img_name, img_fmt, is_tiled, is_mmap = get_file_info(cli_args) + img_path, img_name, img_fmt, is_tiled, is_mmap, mip_msk = get_file_info(cli_args) # import fiber orientation vector data tic = perf_counter() if img_fmt == 'npy' or img_fmt == 'h5': img, is_fiber = load_orient(img_path, img_name, img_fmt) + tissue_msk = None # 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) + img, tissue_msk, is_fiber = load_raw(img_path, img_name, img_fmt, + is_tiled=is_tiled, is_mmap=is_mmap, tmp_dir=tmp_dir, mip_msk=mip_msk) # print import time print_import_time(tic) @@ -412,7 +426,7 @@ def load_microscopy_image(cli_args): # 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 + return img, tissue_msk, is_tiled, is_fiber, save_dir, tmp_dir, img_name def load_orient(img_path, img_name, img_fmt): @@ -460,7 +474,7 @@ def load_orient(img_path, img_name, img_fmt): return img, is_fiber -def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None): +def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None, mip_msk=False): """ Load raw microscopy volume image. @@ -486,11 +500,17 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir tmp_dir: str temporary file directory + mip_msk: bool + apply tissue reconstruction mask (binarized MIP) + Returns ------- img: numpy.ndarray or NumPy memory-map object microscopy volume image + tissue_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + is_fiber: bool True when pre-estimated fiber orientation vectors are directly provided to the pipeline @@ -519,7 +539,20 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir if is_mmap: img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img[:], mmap_mode='r') + # compute tissue reconstruction mask (binarized MIP) + if mip_msk: + dims = len(img.shape) + if dims == 3: + tissue_mip = np.max(img[:], axis=0) + elif dims == 4: + tissue_mip = np.max(img[:], axis=(0, -1)) + + tissue_msk = create_background_mask(tissue_mip, method='li', black_bg=True) + + else: + tissue_msk = None + # raw image is_fiber = False - return img, is_fiber + return img, tissue_msk, is_fiber diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py index 4cb432d..6b50a24 100644 --- a/foa3d/pipeline.py +++ b/foa3d/pipeline.py @@ -289,8 +289,9 @@ def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+', 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'): + tissue_msk=None, 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', + mip_thr=0.995): """ Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image. @@ -349,6 +350,9 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8) soma mask image + tissue_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + ch_lpf: int neuronal bodies channel @@ -388,16 +392,22 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc soma_thresh: str soma channel thresholding method + mip_thr: float + relative threshold on non-zero tissue MIP pixels + Returns ------- None """ # slice fiber image slice - fiber_slice = slice_channel(img, rng_in, channel=ch_mye, is_tiled=is_tiled) + fiber_slice, tissue_msk_slice = slice_channel(img, rng_in, channel=ch_mye, tissue_msk=tissue_msk, is_tiled=is_tiled) # skip background slice - if np.max(fiber_slice) != 0: + crt = np.count_nonzero(tissue_msk_slice) / np.prod(tissue_msk_slice.shape) > mip_thr \ + if tissue_msk_slice is not None else np.max(fiber_slice) != 0 + + if crt: # preprocess fiber slice iso_fiber_slice, rsz_pad = \ @@ -646,7 +656,7 @@ def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, i px_size_iso, save_dir, img_name, odf_scale_um) -def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, +def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, tissue_msk=None, ram=None, jobs=4, backend='threading', is_tiled=False, verbose=100): """ Apply 3D Frangi filtering to basic TPFM image slices using parallel threads. @@ -668,6 +678,9 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, img_name: str name of the microscopy volume image + tissue_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + ram: float maximum RAM available to the Frangi filtering stage [B] @@ -744,9 +757,9 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, 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, 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) + iso_fiber_img, fiber_msk, soma_msk, tissue_msk=tissue_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 diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py index dacf28b..f9451c2 100644 --- a/foa3d/preprocessing.py +++ b/foa3d/preprocessing.py @@ -128,9 +128,7 @@ def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mo resize(img[z, ...], output_shape=tuple(iso_shape[1:]), anti_aliasing=anti_alias, preserve_range=True) # 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 + rsz_pad = np.floor(np.multiply(np.array([rsz_ratio, rsz_ratio]).transpose(), pad)).astype(int) \ + if pad is not None else None - else: - return iso_img, None + return iso_img, rsz_pad diff --git a/foa3d/slicing.py b/foa3d/slicing.py index 360ed8f..9bb10dc 100644 --- a/foa3d/slicing.py +++ b/foa3d/slicing.py @@ -522,7 +522,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, is_tiled=False): +def slice_channel(img, rng, channel, tissue_msk=None, is_tiled=False): """ Slice desired channel from input image volume. @@ -537,6 +537,9 @@ def slice_channel(img, rng, channel, is_tiled=False): channel: int image channel axis + tissue_msk: numpy.ndarray (dtype=bool) + tissue reconstruction binary mask + is_tiled: bool True for tiled reconstructions aligned using ZetaStitcher @@ -544,6 +547,9 @@ def slice_channel(img, rng, channel, is_tiled=False): ------- img_slice: numpy.ndarray sliced image patch + + tissue_msk_slice: numpy.ndarray (dtype=bool) + sliced tissue reconstruction binary mask """ z_rng, r_rng, c_rng = rng @@ -555,7 +561,9 @@ def slice_channel(img, rng, channel, is_tiled=False): else: img_slice = img[z_rng, r_rng, c_rng, channel] - return img_slice + tissue_msk_slice = tissue_msk[r_rng, c_rng] if tissue_msk is not None else None + + return img_slice, tissue_msk_slice def adjust_axis_range(img_shape, start, stop, ax, ovlp=0): diff --git a/foa3d/utils.py b/foa3d/utils.py index c462035..c5440f6 100644 --- a/foa3d/utils.py +++ b/foa3d/utils.py @@ -15,7 +15,7 @@ threshold_yen) -def create_background_mask(img, method='yen'): +def create_background_mask(img, method='yen', black_bg=False): """ Compute background mask. @@ -27,6 +27,9 @@ def create_background_mask(img, method='yen'): method: str image thresholding method + black_bg: bool + generate foreground mask + Returns ------- bg_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool) @@ -49,7 +52,7 @@ def create_background_mask(img, method='yen'): raise ValueError("Unsupported thresholding method!!!") # compute mask - bg_msk = img < thresh + bg_msk = img >= thresh if black_bg else img < thresh return bg_msk