Minor fixes

foa3d/__main__.py

@@ -7,16 +7,18 @@
 def foa3d(cli_args):
 
     # load 3D microscopy image or 4D array of fiber orientation vectors
-    img, ts_msk, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args)
+    img, ts_msk, is_tiled, is_fovec, 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:
+    if not is_fovec:
         fbr_vec_img, iso_fbr_img, px_sz, img_name \
             = parallel_frangi_on_slices(img, cli_args, save_dir[0], tmp_dir, img_name,
                                         ts_msk=ts_msk, ram=ram, jobs=jobs, is_tiled=is_tiled)
+
+    # skip Frangi filter stage if orientation vectors were directly provided as input
     else:
         fbr_vec_img, iso_fbr_img, px_sz = (img, None, None)
 

foa3d/frangi.py

@@ -107,7 +107,6 @@ def compute_frangi_features(eigen1, eigen2, eigen3, gamma):
         background score sensitivity
         (automatically computed if not provided as input)
     """
-
     ra = divide_nonzero(np.abs(eigen2), np.abs(eigen3))
     rb = divide_nonzero(np.abs(eigen1), np.sqrt(np.abs(np.multiply(eigen2, eigen3))))
     s = compute_structureness(eigen1, eigen2, eigen3)

foa3d/input.py

@@ -52,47 +52,56 @@ def get_cli_parser():
                     'Scientific Reports, 13, pp. 4160.\n',
         formatter_class=CustomFormatter)
     cli_parser.add_argument(dest='image_path',
-                            help='path to input 3D microscopy image or to 4D array of fiber orientation vectors\n'
-                                 '* supported formats: .tif (image), '
-                                 '.npy (image or fiber vectors), .yml (ZetaStitcher stitch file)\n'
-                                 '* image  axes order: (Z, Y, X)\n'
-                                 '* vector axes order: (Z, Y, X, C)')
+                            help='path to input 3D microscopy image or 4D array of fiber orientation vectors\n'
+                                 '* supported formats:\n'
+                                 '  - .tif .tiff (microscopy image or fiber orientation vectors)\n'
+                                 '  - .yml (ZetaStitcher\'s stitch file of tiled microscopy reconstruction)\n'
+                                 '  - .npy (fiber orientation vectors)\n'
+                                 '* image axes order:\n'
+                                 '  - grayscale image:      (Z, Y, X)\n'
+                                 '  - RGB image:            (Z, Y, X, C)\n'
+                                 '  - RGB tiled image:      (Z, C, Y, X)\n'                                 
+                                 '  - NumPy vector image:   (Z, Y, X, C)\n'
+                                 '  - HDF5  vector image:   (Z, Y, X, C)\n'
+                                 '  - TIFF  vector image:   (Z, C, Y, X)\n')
     cli_parser.add_argument('-a', '--alpha', type=float, default=0.001,
-                            help='Frangi plate-like object sensitivity')
+                            help='Frangi\'s plate-like object sensitivity')
     cli_parser.add_argument('-b', '--beta', type=float, default=1.0,
-                            help='Frangi blob-like object sensitivity')
+                            help='Frangi\'s blob-like object sensitivity')
     cli_parser.add_argument('-g', '--gamma', type=float, default=None,
-                            help='Frangi background score sensitivity')
+                            help='Frangi\'s 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('-j', '--jobs', type=int, default=None,
-                            help='number of parallel threads used by the Frangi filtering stage: '
+                            help='number of parallel threads used by the Frangi filter stage: '
                                  'use one thread per logical core if None')
     cli_parser.add_argument('-r', '--ram', type=float, default=None,
-                            help='maximum RAM available to the Frangi filtering stage [GB]: use all if None')
+                            help='maximum RAM available to the Frangi filter stage [GB]: use all if None')
     cli_parser.add_argument('-m', '--mmap', action='store_true', default=False,
-                            help='create a memory-mapped array of the 3D microscopy image')
+                            help='create a memory-mapped array of input microscopy or fiber orientation images')
     cli_parser.add_argument('--px-size-xy', type=float, default=0.878, help='lateral pixel size [μm]')
     cli_parser.add_argument('--px-size-z', type=float, default=1.0, help='longitudinal pixel size [μm]')
-    cli_parser.add_argument('--psf-fwhm-x', type=float, default=0.692, help='PSF FWHM along the X axis [μm]')
-    cli_parser.add_argument('--psf-fwhm-y', type=float, default=0.692, help='PSF FWHM along the Y axis [μm]')
-    cli_parser.add_argument('--psf-fwhm-z', type=float, default=2.612, help='PSF FWHM along the Z axis [μm]')
-    cli_parser.add_argument('--fb-ch', type=int, default=1, help='myelinated fibers channel')
+    cli_parser.add_argument('--psf-fwhm-x', type=float, default=0.692, help='PSF FWHM along horizontal x-axis [μm]')
+    cli_parser.add_argument('--psf-fwhm-y', type=float, default=0.692, help='PSF FWHM along vertical x-axis [μm]')
+    cli_parser.add_argument('--psf-fwhm-z', type=float, default=2.612, help='PSF FWHM along depth z-axis [μm]')
+    cli_parser.add_argument('--fb-ch', type=int, default=1, help='neuronal fibers channel')
     cli_parser.add_argument('--bc-ch', type=int, default=0, help='neuronal bodies channel')
     cli_parser.add_argument('--z-min', type=float, default=0, help='forced minimum output z-depth [μm]')
     cli_parser.add_argument('--z-max', type=float, default=None, help='forced maximum output z-depth [μm]')
     cli_parser.add_argument('--hsv', action='store_true', default=False,
-                            help='toggle HSV colormap for 3D fiber orientations')
+                            help='generate HSV colormap for 3D fiber orientations')
     cli_parser.add_argument('--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,
                             help='degrees of the spherical harmonics series expansion (even number between 2 and 10)')
     cli_parser.add_argument('-o', '--out', type=str, default=None,
                             help='output directory')
     cli_parser.add_argument('-c', '--cell-msk', action='store_true', default=False,
-                            help='toggle optional neuronal body masking')
+                            help='apply neuronal body mask (the optional channel of neuronal bodies must be available)')
     cli_parser.add_argument('-t', '--tissue-msk', action='store_true', default=False,
-                            help='apply tissue reconstruction mask (binarized MIP)')
+                            help='apply tissue background mask')
+    cli_parser.add_argument('-v', '--vec', action='store_true', default=False,
+                            help='fiber orientation vector image')
 
     # parse arguments
     cli_args = cli_parser.parse_args()
@@ -185,6 +194,10 @@ def get_file_info(cli_args):
     is_tiled: bool
         True for tiled reconstructions aligned using ZetaStitcher
 
+    is_fovec: bool
+        True when pre-estimated fiber orientation vectors
+        are directly provided to the pipeline
+
     is_mmap: bool
         create a memory-mapped array of the 3D microscopy image,
         increasing the parallel processing performance
@@ -210,12 +223,13 @@ def get_file_info(cli_args):
         img_fmt = split_name[-1]
         img_name = img_name.replace('.{}'.format(img_fmt), '')
         is_tiled = True if img_fmt == 'yml' else False
+        is_fovec = cli_args.vec
 
     # apply tissue reconstruction mask (binarized MIP)
     msk_mip = cli_args.tissue_msk
     fb_ch = cli_args.fb_ch
 
-    return img_path, img_name, img_fmt, is_tiled, is_mmap, msk_mip, fb_ch
+    return img_path, img_name, img_fmt, is_tiled, is_fovec, is_mmap, msk_mip, fb_ch
 
 
 def get_frangi_config(cli_args, img_name):
@@ -386,7 +400,7 @@ def load_microscopy_image(cli_args):
     is_tiled: bool
         True for tiled microscopy reconstructions aligned using ZetaStitcher
 
-    is_fiber: bool
+    is_fovec: bool
         True when pre-estimated fiber orientation vectors
         are directly provided to the pipeline
 
@@ -404,32 +418,32 @@ def load_microscopy_image(cli_args):
     tmp_dir = tempfile.mkdtemp()
 
     # retrieve input file information
-    img_path, img_name, img_fmt, is_tiled, is_mmap, msk_mip, fb_ch = get_file_info(cli_args)
+    img_path, img_name, img_fmt, is_tiled, is_fovec, is_mmap, msk_mip, fb_ch = 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)
+    if is_fovec:
+        img = load_orient(img_path, img_name, img_fmt, is_mmap=is_mmap, tmp_dir=tmp_dir)
         tissue_msk = None
 
     # import raw 3D microscopy image
     else:
-        img, tissue_msk, is_fiber = load_raw(img_path, img_name, img_fmt, is_tiled=is_tiled, is_mmap=is_mmap,
-                                             tmp_dir=tmp_dir, msk_mip=msk_mip, fb_ch=fb_ch)
+        img, tissue_msk = load_raw(img_path, img_name, img_fmt,
+                                   is_tiled=is_tiled, is_mmap=is_mmap, tmp_dir=tmp_dir, msk_mip=msk_mip, fb_ch=fb_ch)
 
     # 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()
+    print_image_shape(cli_args, img, is_tiled) if not is_fovec else print()
 
     # create saving directory
-    save_dir = create_save_dirs(img_path, img_name, cli_args, is_fiber=is_fiber)
+    save_dir = create_save_dirs(img_path, img_name, cli_args, is_fovec=is_fovec)
 
-    return img, tissue_msk, is_tiled, is_fiber, save_dir, tmp_dir, img_name
+    return img, tissue_msk, is_tiled, is_fovec, save_dir, tmp_dir, img_name
 
 
-def load_orient(img_path, img_name, img_fmt):
+def load_orient(img_path, img_name, img_fmt, is_mmap=False, tmp_dir=None):
     """
     Load array of 3D fiber orientations.
 
@@ -444,14 +458,18 @@ def load_orient(img_path, img_name, img_fmt):
     img_fmt: str
         format of the 3D microscopy image
 
+    is_mmap: bool
+        create a memory-mapped array of the 3D microscopy image,
+        increasing the parallel processing performance
+        (the image will be preliminarily loaded to RAM)
+
+    tmp_dir: str
+        temporary file directory
+
     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
@@ -460,18 +478,24 @@ def load_orient(img_path, img_name, img_fmt):
     # load fiber orientations
     if img_fmt == 'npy':
         img = np.load(img_path, mmap_mode='r')
-    else:
+        if is_mmap:
+            img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r')
+    elif img_fmt == 'h5':
         img_file = h5py.File(img_path, 'r')
         img = img_file.get(img_file.keys()[0])
+    elif img_fmt == 'tif' or img_fmt == 'tiff':
+        img = tiff.imread(img_path)
+        img = np.moveaxis(img, 1, -1)
+        if is_mmap:
+            img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r')
 
     # 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))
 
-        return img, is_fiber
+        return img
 
 
 def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir=None, msk_mip=False, fb_ch=1):
@@ -511,12 +535,8 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir
     img: numpy.ndarray or NumPy memory-map object
         3D microscopy image
 
-    tissue_msk: numpy.ndarray (dtype=bool)
+    ts_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
     """
 
     # print heading
@@ -538,25 +558,26 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir
         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')
+        # 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')
 
-    # compute tissue reconstruction mask (binarized MIP)
+    # generate tissue reconstruction mask
     if msk_mip:
         dims = len(img.shape)
         if dims == 3:
-            tissue_mip = np.max(img[:], axis=0)
+            img_fbr = img
         elif dims == 4:
-            img_mye = img[:, fb_ch, :, :] if is_tiled else img[..., fb_ch]
-            tissue_mip = np.max(img_mye, axis=0)
+            img_fbr = img[:, fb_ch, :, :] if is_tiled else img[..., fb_ch]
 
-        tissue_msk = create_background_mask(tissue_mip, method='li', black_bg=True)
+        # compute MIP (naive for loop to minimize the required RAM)
+        ts_mip = np.zeros(img_fbr.shape[1:], dtype=img_fbr.dtype)
+        for z in range(img_fbr.shape[0]):
+            stk = np.stack((ts_mip, img_fbr[z]))
+            ts_mip = np.max(stk, axis=0)
 
+        ts_msk = create_background_mask(ts_mip, method='li', black_bg=True)
     else:
-        tissue_msk = None
-
-    # raw image
-    is_fiber = False
+        ts_msk = None
 
-    return img, tissue_msk, is_fiber
+    return img, ts_msk

foa3d/output.py

@@ -6,7 +6,7 @@
 import tifffile as tiff
 
 
-def create_save_dirs(img_path, img_name, cli_args, is_fiber=False):
+def create_save_dirs(img_path, img_name, cli_args, is_fovec=False):
     """
     Create saving directory.
 
@@ -21,7 +21,7 @@ def create_save_dirs(img_path, img_name, cli_args, is_fiber=False):
     cli_args: see ArgumentParser.parse_args
         updated namespace of command line arguments
 
-    is_fiber: bool
+    is_fovec: bool
         True when fiber orientation vectors are provided as input
         to the pipeline
 
@@ -46,7 +46,7 @@ def create_save_dirs(img_path, img_name, cli_args, is_fiber=False):
         makedirs(base_out_dir)
 
     # create Frangi filter output subdirectory
-    if not is_fiber:
+    if not is_fovec:
         frangi_dir = path.join(base_out_dir, 'frangi')
         makedirs(frangi_dir)
         save_dir_lst.append(frangi_dir)

foa3d/pipeline.py

@@ -756,7 +756,7 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk
     return fbr_vec_img, iso_fbr_img, px_sz_iso, img_name
 
 
-def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz_iso, save_dir, tmp_dir, img_name,
+def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz, save_dir, tmp_dir, img_name,
                            backend='loky', ram=None, verbose=100):
     """
     Iterate over the required spatial scales and apply the parallel ODF analysis
@@ -773,8 +773,8 @@ def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz_iso, save_d
     cli_args: see ArgumentParser.parse_args
         populated namespace of command line arguments
 
-    px_sz_iso: numpy.ndarray (axis order=(3,), dtype=float)
-        adjusted isotropic pixel size [μm]
+    px_sz: numpy.ndarray (axis order=(3,), dtype=float)
+        pixel size [μm]
 
     save_dir: str
         saving directory string path
@@ -812,14 +812,14 @@ def parallel_odf_at_scales(fbr_vec_img, iso_fbr_img, cli_args, px_sz_iso, save_d
     num_cpu = get_available_cores()
 
     # generate pixel size if not provided
-    if px_sz_iso is None:
-        px_sz_iso = cli_args.px_size_z * np.ones((3,))
+    if px_sz is None:
+        px_sz = np.array([cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy])
 
     # parallel ODF analysis of fiber orientation vectors
     # over the required spatial scales
     n_jobs = min(num_cpu, len(cli_args.odf_res))
     with Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose, max_nbytes=None) as parallel:
-        parallel(delayed(odf_analysis)(fbr_vec_img, iso_fbr_img, px_sz_iso, save_dir, tmp_dir, img_name,
+        parallel(delayed(odf_analysis)(fbr_vec_img, iso_fbr_img, px_sz, save_dir, tmp_dir, img_name,
                                        odf_norm=odf_norm, odf_deg=cli_args.odf_deg, odf_scale_um=s, ram=ram)
                  for s in cli_args.odf_res)
 

foa3d/slicing.py

@@ -367,7 +367,7 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi
     ns = len(frangi_sigma_um)
 
     # initialize slice batch size
-    batch_sz = np.min([jobs, num_cpu]).astype(int) + 1
+    batch_sz = np.min([jobs // ns, num_cpu]).astype(int) + 1
 
     # get pixel resize ratio
     px_rsz_ratio = np.divide(px_sz, px_sz_iso)

foa3d/utils.py

@@ -38,8 +38,7 @@ def create_background_mask(img, method='yen', black_bg=False):
 
     # select thresholding method
     if method == 'li':
-        init_li = np.mean(img[img != 0])
-        thresh = threshold_li(img, initial_guess=init_li)
+        thresh = threshold_li(img)
     elif method == 'niblack':
         thresh = threshold_niblack(img, window_size=15, k=0.2)
     elif method == 'sauvola':