Add new background rejection mechanism

foa3d/__main__.py

@@ -7,16 +7,16 @@
 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)
 
     # 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)
 

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

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

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

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,13 +537,19 @@ 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
 
     Returns
     -------
     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):

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