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,12 @@ 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)
+
+    ch_mye: int  
+        myelinated fibers channel
     """
 
     # get microscopy image path and name
@@ -202,7 +211,11 @@ 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
+    ch_mye = cli_args.ch_mye
+
+    return img_path, img_name, img_fmt, is_tiled, is_mmap, mip_msk, ch_mye
 
 
 def get_frangi_config(cli_args, img_name):
@@ -368,6 +381,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 +408,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, ch_mye = 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, ch_mye=ch_mye)
 
     # print import time
     print_import_time(tic)
@@ -412,7 +430,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 +478,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, ch_mye=1):
     """
     Load raw microscopy volume image.
 
@@ -486,11 +504,20 @@ 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)
+
+    ch_mye: int  
+        myelinated fibers channel
+
     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 +546,21 @@ 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:
+            img_mye = img[:, ch_mye, :, :] if is_tiled else img[..., ch_mye]
+            tissue_mip = np.max(img_mye, axis=0)
+
+        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/output.py

@@ -101,7 +101,7 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', odi=False):
         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
+        bigtiff = True if nd_array.itemsize * nd_array.size >= 4294967296 else False
 
         # save array to TIFF file
         out_name = '{}.{}'.format(fname, fmt)

foa3d/pipeline.py

@@ -275,7 +275,7 @@ def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+',
     if ram is None:
         ram = psutil.virtual_memory()[1]
     item_sz = get_item_size(dtype)
-    vol_sz = item_sz * np.prod(shape)
+    vol_sz = item_sz * np.prod(shape).astype(np.float64)
 
     # create memory-mapped array or HDF5 file depending on available memory resources
     if vol_sz >= ram:
@@ -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.005):
     """
     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
 
@@ -369,7 +373,7 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
         (i.e., negative contrast polarity)
 
     hsv_vec_cmap: bool
-        if True, generate a HSV orientation color map based on XY-plane orientation angles
+        if True, generate an HSV orientation color map based on XY-plane orientation angles
         (instead of an RGB map using the cartesian components of the estimated vectors)
 
     mask_lpf: bool
@@ -388,32 +392,43 @@ 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 = \
-            correct_image_anisotropy(fiber_slice, px_rsz_ratio, sigma=smooth_sigma, pad=pad)
+        iso_fiber_slice, rsz_pad, rsz_tissue_msk_slice = \
+            correct_image_anisotropy(fiber_slice, px_rsz_ratio,
+                                     sigma=smooth_sigma, pad=pad, tissue_msk=tissue_msk_slice)
 
         # pad fiber slice if required
         if rsz_pad is not None:
             iso_fiber_slice = np.pad(iso_fiber_slice, rsz_pad, mode=pad_mode)
 
+            # pad tissue mask if available
+            if rsz_tissue_msk_slice is not None:
+                rsz_tissue_msk_slice = np.pad(rsz_tissue_msk_slice, rsz_pad, mode='constant')
+
         # 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, frangi_slice, fiber_vec_slice, eigenval_slice = \
-            crop_slice_lst([iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice], rng_out, ovlp=ovlp)
+        iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice, rsz_tissue_msk_slice = \
+            crop_slice_lst([iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice, rsz_tissue_msk_slice],
+                           rng_out, ovlp=ovlp)
 
         # generate fractional anisotropy image
         frac_anis_slice = compute_fractional_anisotropy(eigenval_slice)
@@ -423,7 +438,8 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
 
         # 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)
+            mask_background(frangi_slice, fiber_vec_slice, fiber_clr_slice,
+                            tissue_msk=rsz_tissue_msk_slice, method=fiber_thresh, invert=False)
 
         # (optional) neuronal body masking
         if mask_lpf:
@@ -432,7 +448,7 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
             soma_slice = slice_channel(img, rng_in_neu, channel=ch_lpf, is_tiled=is_tiled)
 
             # resize soma slice (lateral blurring and downsampling)
-            iso_soma_slice, _ = correct_image_anisotropy(soma_slice, px_rsz_ratio)
+            iso_soma_slice, _, _ = correct_image_anisotropy(soma_slice, px_rsz_ratio)
 
             # crop isotropized soma slice
             iso_soma_slice = crop_slice(iso_soma_slice, rng_out)
@@ -525,7 +541,8 @@ def fill_frangi_volumes(fiber_vec_img, fiber_vec_clr, frac_anis_img, frangi_img,
         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):
+def mask_background(img, fiber_vec_slice, fiber_clr_slice,
+                    frac_anis_slice=None, method='yen', invert=False, tissue_msk=None):
     """
     Mask orientation volume arrays.
 
@@ -549,6 +566,9 @@ def mask_background(img, fiber_vec_slice, fiber_clr_slice, frac_anis_slice=None,
     invert: bool
         mask inversion flag
 
+    tissue_msk: numpy.ndarray (dtype=bool)
+        tissue reconstruction binary mask
+
     Returns
     -------
     fiber_vec_slice: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
@@ -567,6 +587,10 @@ def mask_background(img, fiber_vec_slice, fiber_clr_slice, frac_anis_slice=None,
     # generate background mask
     background = create_background_mask(img, method=method)
 
+    # apply tissue reconstruction mask, when provided
+    if tissue_msk is not None:
+        background = np.logical_or(background, np.logical_not(tissue_msk))
+
     # invert mask
     if invert:
         background = np.logical_not(background)
@@ -646,7 +670,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 +692,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 +771,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