diff --git a/foa3d/__main__.py b/foa3d/__main__.py
index 698cc32..dd8e785 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, tissue_msk, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args)
+    img, ts_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,9 @@ 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, tissue_msk=tissue_msk,
-                                        ram=ram, jobs=jobs, is_tiled=is_tiled)
+            = 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)
+
     else:
         fiber_vec_img, iso_fiber_img, px_sz = (img, None, None)
 
diff --git a/foa3d/frangi.py b/foa3d/frangi.py
index f620780..d8d960b 100644
--- a/foa3d/frangi.py
+++ b/foa3d/frangi.py
@@ -17,7 +17,7 @@ def analyze_hessian_eigen(img, sigma, trunc=4):
     Parameters
     ----------
     img: numpy.ndarray (axis order=(Z,Y,X))
-        microscopy volume image
+        3D microscopy image
 
     sigma: int
         spatial scale [px]
diff --git a/foa3d/input.py b/foa3d/input.py
index 4974a59..c9f9758 100644
--- a/foa3d/input.py
+++ b/foa3d/input.py
@@ -193,7 +193,7 @@ def get_file_info(cli_args):
     mip_msk: bool
         apply tissue reconstruction mask (binarized MIP)
 
-    ch_mye: int  
+    ch_mye: int
         myelinated fibers channel
     """
 
@@ -248,8 +248,7 @@ def get_frangi_config(cli_args, img_name):
         Frangi filter scales [μm]
 
     smooth_sigma: numpy.ndarray (shape=(3,), dtype=int)
-        3D standard deviation of the low-pass Gaussian filter [px]
-        (applied to the XY plane)
+        3D standard deviation of the smoothing Gaussian filter [px]
 
     px_sz: numpy.ndarray (shape=(3,), dtype=float)
         pixel size [μm]
@@ -507,7 +506,7 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir
     mip_msk: bool
         apply tissue reconstruction mask (binarized MIP)
 
-    ch_mye: int  
+    ch_mye: int
         myelinated fibers channel
 
     Returns
@@ -524,7 +523,7 @@ def load_raw(img_path, img_name, img_fmt, is_tiled=False, is_mmap=False, tmp_dir
     """
 
     # print heading
-    print(color_text(0, 191, 255, "\nMicroscopy Volume Image Import\n"))
+    print(color_text(0, 191, 255, "\nMicroscopy Image Import\n"))
 
     # load microscopy tiled reconstruction (aligned using ZetaStitcher)
     if is_tiled:
diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py
index 9e2db01..748e0c1 100644
--- a/foa3d/pipeline.py
+++ b/foa3d/pipeline.py
@@ -14,8 +14,8 @@
 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, crop_slice_lst, slice_channel)
+from foa3d.slicing import (config_frangi_batch, generate_slice_lists,
+                           crop, crop_lst, slice_image)
 from foa3d.utils import (create_background_mask, create_hdf5_dset,
                          create_memory_map, divide_nonzero,
                          get_available_cores, get_item_size, hsv_orient_cmap,
@@ -289,7 +289,7 @@ 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,
-                   tissue_msk=None, ch_lpf=0, ch_mye=1, alpha=0.05, beta=1, gamma=100, dark=False, mask_lpf=False,
+                   ts_msk=None, ch_lpf=0, ch_mye=1, alpha=0.05, beta=1, gamma=100, dark=False, msk_bc=False,
                    hsv_vec_cmap=False, pad_mode='reflect', is_tiled=False, fiber_thresh='li', soma_thresh='yen',
                    mip_thr=0.005):
     """
@@ -401,7 +401,7 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
     """
 
     # slice fiber image slice
-    fiber_slice, tissue_msk_slice = slice_channel(img, rng_in, channel=ch_mye, tissue_msk=tissue_msk, is_tiled=is_tiled)
+    fiber_slice, tissue_msk_slice = slice_image(img, rng_in, ch_mye, ts_msk=ts_msk, is_tiled=is_tiled)
 
     # skip background slice
     crt = np.count_nonzero(tissue_msk_slice) / np.prod(tissue_msk_slice.shape) > mip_thr \
@@ -410,8 +410,7 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
 
         # preprocess fiber slice
         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)
+            correct_image_anisotropy(fiber_slice, px_rsz_ratio, sigma=smooth_sigma, pad=pad, ts_msk=tissue_msk_slice)
 
         # pad fiber slice if required
         if rsz_pad is not None:
@@ -427,8 +426,8 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
 
         # crop resulting slices
         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)
+            crop_lst([iso_fiber_slice, frangi_slice, fiber_vec_slice, eigenval_slice, rsz_tissue_msk_slice],
+                             rng_out, ovlp)
 
         # generate fractional anisotropy image
         frac_anis_slice = compute_fractional_anisotropy(eigenval_slice)
@@ -437,21 +436,21 @@ def fiber_analysis(img, rng_in, rng_in_neu, rng_out, pad, ovlp, smooth_sigma, sc
         fiber_clr_slice = hsv_orient_cmap(fiber_vec_slice) if hsv_vec_cmap else rgb_orient_cmap(fiber_vec_slice)
 
         # remove background
-        fiber_vec_slice, fiber_clr_slice, fiber_msk_slice = \
-            mask_background(frangi_slice, fiber_vec_slice, fiber_clr_slice,
+        fiber_vec_slice, fiber_clr_slice, frac_anis_slice, fiber_msk_slice = \
+            mask_background(frangi_slice, fiber_vec_slice, fiber_clr_slice, frac_anis_slice,
                             tissue_msk=rsz_tissue_msk_slice, method=fiber_thresh, invert=False)
 
         # (optional) neuronal body masking
-        if mask_lpf:
+        if msk_bc:
 
             # get soma image slice
-            soma_slice = slice_channel(img, rng_in_neu, channel=ch_lpf, is_tiled=is_tiled)
+            soma_slice = slice_image(img, rng_in_neu, ch_lpf, is_tiled=is_tiled)
 
             # resize soma slice (lateral blurring and downsampling)
             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)
+            iso_soma_slice = crop(iso_soma_slice, rng_out)
 
             # mask neuronal bodies
             fiber_vec_slice, fiber_clr_slice, frac_anis_slice, soma_msk_slice = \
@@ -598,14 +597,9 @@ def mask_background(img, fiber_vec_slice, fiber_clr_slice,
     # apply mask to input arrays
     fiber_vec_slice[background, :] = 0
     fiber_clr_slice[background, :] = 0
+    frac_anis_slice[background] = 0
 
-    # (optional) mask fractional anisotropy
-    if frac_anis_slice is not None:
-        frac_anis_slice[background] = 0
-        return fiber_vec_slice, fiber_clr_slice, frac_anis_slice, background
-
-    else:
-        return fiber_vec_slice, fiber_clr_slice, background
+    return fiber_vec_slice, fiber_clr_slice, frac_anis_slice, background
 
 
 def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, img_name, odf_scale_um,
@@ -670,8 +664,8 @@ 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, tissue_msk=None,
-                              ram=None, jobs=4, backend='threading', is_tiled=False, verbose=100):
+def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk=None,
+                              ram=None, jobs=4, backend='threading', is_tiled=False, verbose=10):
     """
     Apply 3D Frangi filtering to basic TPFM image slices using parallel threads.
 
@@ -713,77 +707,61 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, tissue
 
     Returns
     -------
-    fiber_vec_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32)
+    fbr_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 volume image (fiber probability volume)
-
-    iso_fiber_img: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
+    iso_fbr_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
+    px_sz_iso: int
+        isotropic pixel size [μm]
 
-    soma_msk: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
-        soma mask image
+    img_name: str
+        name of the input microscopy 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)
+    alpha, beta, gamma, frangi_sigma, frangi_sigma_um, smooth_sigma, px_sz, px_sz_iso, \
+        z_rng, ch_lpf, ch_mye, msk_bc, 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_mye, mask_lpf = \
-        get_image_info(img, px_size, mask_lpf, ch_mye, is_tiled=is_tiled)
+    # get info about the input microscopy image
+    img_shp, img_shp_um, img_item_sz, ch_mye, msk_bc = get_image_info(img, px_sz, msk_bc, 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, ram=ram, jobs=jobs)
+    # configure the batch of basic image slices analyzed using parallel threads
+    batch_sz, in_slc_shp, in_slc_shp_um, px_rsz_ratio, ovlp, ovlp_rsz = \
+        config_frangi_batch(px_sz, px_sz_iso, img_shp, img_item_sz, smooth_sigma, frangi_sigma, ram=ram, jobs=jobs)
 
-    # get info on the processed image slices
-    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_um, mask_lpf, batch_size, slice_size=max_slice_size)
+    # get info about the processed image slices
+    in_rng_lst, in_pad_lst, out_rng_lst, bc_rng_lst, out_slc_shp, tot_slc_num, batch_sz = \
+        generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp, msk_bg=msk_bc, jobs=jobs)
 
     # initialize output arrays
-    fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_img, \
-        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_rng=z_rng, mask_lpf=mask_lpf, ram=ram)
+    fiber_dset_path, fbr_vec_img, fbr_vec_clr, frac_anis_img, frangi_img, fbr_msk, iso_fbr_img, bc_msk, z_sel = \
+        init_frangi_volumes(img_shp, out_slc_shp, px_rsz_ratio, tmp_dir, z_rng=z_rng, mask_lpf=msk_bc, 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, mask_lpf)
+    print_frangi_info(alpha, beta, gamma, frangi_sigma_um, img_shp_um, in_slc_shp_um, tot_slc_num,
+                      px_sz, img_item_sz, msk_bc)
 
-    # parallel Frangi filter-based fiber orientation analysis of microscopy sub-volumes
+    # parallel Frangi filter-based fiber orientation analysis of microscopy image slices
     start_time = perf_counter()
-    with Parallel(n_jobs=batch_size, backend=backend, verbose=verbose, max_nbytes=None) as parallel:
+    with Parallel(n_jobs=batch_sz, backend=backend, verbose=verbose, max_nbytes=None) as parallel:
         parallel(
-            delayed(fiber_analysis)(img,
-                                    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, tissue_msk=tissue_msk,
+            delayed(fiber_analysis)(img, in_rng_lst[i], bc_rng_lst[i], out_rng_lst[i], in_pad_lst[i], ovlp_rsz,
+                                    smooth_sigma, frangi_sigma, px_rsz_ratio, z_sel, fbr_vec_img, fbr_vec_clr,
+                                    frac_anis_img, frangi_img, iso_fbr_img, fbr_msk, bc_msk, ts_msk=ts_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))
+                                    msk_bc=msk_bc, hsv_vec_cmap=hsv_vec_cmap, is_tiled=is_tiled)
+            for i in range(tot_slc_num))
 
-    # save Frangi output arrays
-    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)
+    # save output arrays
+    save_frangi_arrays(fiber_dset_path, fbr_vec_img, fbr_vec_clr, frac_anis_img, frangi_img, fbr_msk, bc_msk,
+                       px_sz_iso, save_dir, img_name)
 
     # print Frangi filtering time
     print_analysis_time(start_time)
 
-    return fiber_vec_img, iso_fiber_img, px_size_iso, img_name
+    return fbr_vec_img, iso_fbr_img, px_sz_iso, img_name
 
 
 def parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_size_iso, save_dir, tmp_dir, img_name,
@@ -893,7 +871,7 @@ def save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_
         saving directory string path
 
     img_name: str
-        name of the input microscopy volume image
+        name of the input microscopy image
 
     Returns
     -------
diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py
index c59ea8d..fb6ff22 100644
--- a/foa3d/preprocessing.py
+++ b/foa3d/preprocessing.py
@@ -63,8 +63,7 @@ def config_anisotropy_correction(px_sz, psf_fwhm):
             print_blur(smooth_sigma_um, psf_fwhm)
 
         # print pixel resize info
-        if cndt_2:
-            print_new_res(px_sz_iso)
+        print_new_res(px_sz_iso) if cndt_2 else print()
 
     # skip line
     else:
@@ -74,7 +73,7 @@ def config_anisotropy_correction(px_sz, psf_fwhm):
 
 
 def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mode='reflect',
-                             anti_alias=True, trunc=4, tissue_msk=None):
+                             anti_alias=True, trunc=4, ts_msk=None):
     """
     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.
@@ -122,11 +121,11 @@ def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mo
     if np.all(rsz_ratio == 1):
 
         # resize tissue mask, when available
-        if tissue_msk is not None:
-            tissue_msk = tissue_msk[np.newaxis, ...]
+        if ts_msk is not None:
+            tissue_msk = ts_msk[np.newaxis, ...]
             tissue_msk = np.repeat(tissue_msk, img.shape[0], axis=0)
 
-        return img, pad, tissue_msk
+        return img, pad, ts_msk
 
     # lateral blurring
     else:
@@ -145,8 +144,8 @@ def correct_image_anisotropy(img, rsz_ratio, sigma=None, pad=None, smooth_pad_mo
             if pad is not None else None
 
         # resize tissue mask, when available
-        if tissue_msk is not None:
-            rsz_tissue_msk = resize(tissue_msk, output_shape=tuple(iso_shape[1:]), preserve_range=True)
+        if ts_msk is not None:
+            rsz_tissue_msk = resize(ts_msk, output_shape=tuple(iso_shape[1:]), preserve_range=True)
             rsz_tissue_msk = rsz_tissue_msk[np.newaxis, ...]
             rsz_tissue_msk = np.repeat(rsz_tissue_msk, iso_shape[0], axis=0)
         else:
diff --git a/foa3d/printing.py b/foa3d/printing.py
index 0ab081c..ad14613 100644
--- a/foa3d/printing.py
+++ b/foa3d/printing.py
@@ -113,7 +113,7 @@ def print_analysis_time(start_time):
     None
     """
     _, mins, secs = elapsed_time(start_time)
-    print("\nVolume image analyzed in: {0} min {1:3.1f} s\n".format(mins, secs))
+    print("\nMicroscopy image analyzed in: {0} min {1:3.1f} s\n".format(mins, secs))
 
 
 def print_blur(smooth_sigma_um, psf_fwhm):
@@ -135,7 +135,7 @@ def print_blur(smooth_sigma_um, psf_fwhm):
     """
     print("Gaussian blur \u03C3      [μm]: ({0:.3f}, {1:.3f}, {2:.3f})"
           .format(smooth_sigma_um[0], smooth_sigma_um[1], smooth_sigma_um[2]))
-    print("Adjusted PSF FWHM    [μm]: ({0:.3f}, {0:.3f}, {0:.3f})\n".format(psf_fwhm[0]))
+    print("Adjusted PSF FWHM    [μm]: ({0:.3f}, {0:.3f}, {0:.3f})".format(psf_fwhm[0]))
 
 
 def print_import_time(start_time):
@@ -152,7 +152,7 @@ def print_import_time(start_time):
     None
     """
     _, mins, secs = elapsed_time(start_time)
-    print("Volume image loaded in: {0} min {1:3.1f} s".format(mins, secs))
+    print("Image loaded in: {0} min {1:3.1f} s".format(mins, secs))
 
 
 def print_odf_info(odf_scales_um, odf_degrees):
@@ -195,7 +195,7 @@ def print_prepro_heading():
     -------
     None
     """
-    print(color_text(0, 191, 255, "\n\nMicroscopy Volume Image Preprocessing"))
+    print(color_text(0, 191, 255, "\n\nMicroscopy Image Preprocessing"))
     print("\n                              Z      Y      X")
 
 
@@ -299,7 +299,7 @@ def print_soma_masking(lpf_soma_mask):
     -------
     None
     """
-    prt = 'Lipofuscin soma mask: '
+    prt = 'Soma mask: '
     print('{}active\n'.format(prt)) if lpf_soma_mask else print('{}not active\n'.format(prt))
 
 
diff --git a/foa3d/slicing.py b/foa3d/slicing.py
index 967dc07..2bb3dd8 100644
--- a/foa3d/slicing.py
+++ b/foa3d/slicing.py
@@ -1,184 +1,181 @@
-from itertools import product
-
 import numpy as np
 import psutil
 
+from itertools import product
+from joblib import Parallel, delayed
+from numba import jit, njit
+
 from foa3d.utils import get_available_cores
 
 
-def compute_axis_range(ax, ax_iter, slice_shape, img_shape, slice_per_dim, flip=False, ovlp=0):
+@njit(cache=True)
+def adjust_axis_range(ax_iter, img_shp, slc_per_ax, start, stop, ovlp=0):
     """
-    Adjust image slice coordinates at boundaries.
+    Trim slice axis range at image boundaries
+    and adjust padded ranges accordingly.
 
     Parameters
     ----------
-    ax: int
-        axis index
-
     ax_iter: int
-        iteration counter along axis (see pipeline.iterate_frangi_on_slices)
+        iteration counter along axis
 
-    slice_shape: numpy.ndarray (shape=(3,), dtype=int)
-        shape of the basic image slices analyzed iteratively [px]
+    img_shp: int
+        total image shape along axis [px]
 
-    img_shape: numpy.ndarray (shape=(3,), dtype=int)
-        total image shape [px]
+    slc_per_ax: int
+        slices along axis
 
-    slice_per_dim: numpy.ndarray (shape=(3,), dtype=int)
-        number of image slices along each axis
+    start: int
+        start coordinate [px]
 
-    flip: bool
-        if True, flip axis coordinates
+    stop: int
+        stop coordinate [px]
 
     ovlp: int
-        optional slicing range
-        extension along each axis (on each side)
+        overlapping range between slices along each axis [px]
 
     Returns
     -------
     start: int
-        adjusted start index
+        adjusted start coordinate [px]
 
     stop: int
-        adjusted stop index
+        adjusted stop coordinate [px]
 
-    ovlp: numpy.ndarray (shape=(2,), dtype=int)
-        lower and upper padded boundaries
-        along axis
+    pad: numpy.ndarray (shape=(2,), dtype=int)
+        lower and upper padding ranges [px]
     """
 
-    # compute start and stop coordinates
-    start = ax_iter[ax] * slice_shape[ax]
-    stop = start + slice_shape[ax]
+    # initialize slice padding array
+    pad = np.zeros((2,), dtype=np.int64)
 
-    # flip coordinates if required
-    if flip:
-        start_tmp = start
-        start = img_shape[ax] - stop
-        stop = img_shape[ax] - start_tmp
+    # adjust start coordinate
+    start -= ovlp
+    if start < 0:
+        pad[0] = -start
+        start = 0
 
-    # adjust start and stop coordinates
-    start, stop, pad = adjust_axis_range(img_shape, start, stop, ax=ax, ovlp=ovlp)
+    # adjust stop coordinate
+    stop += ovlp
+    if stop > img_shp:
+        pad[1] = stop - img_shp
+        stop = img_shp
 
-    # 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]
-            pad[1] = ovlp
+    # handle image residuals at boundaries
+    if ax_iter == slc_per_ax - 1:
+        pad[1] = ovlp
+        stop = img_shp
 
     return start, stop, pad
 
 
-def compute_slice_range(ax_iter, slice_shape, img_shape, slice_per_dim, ovlp=0, flip=False):
+@njit(cache=True)
+def compute_axis_range(ax, ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip=False):
     """
-    Compute basic slice coordinates from microscopy volume image.
+    Adjust image slice coordinates at boundaries.
 
     Parameters
     ----------
-    ax_iter: tuple
-        iteration counters along axes
+    ax: int
+        image axis
 
-    slice_shape: numpy.ndarray (shape=(3,), dtype=int)
+    ax_iter: int
+        iteration counter along axis (see pipeline.iterate_frangi_on_slices)
+
+    slc_shp: numpy.ndarray (shape=(3,), dtype=int)
         shape of the basic image slices
         analyzed using parallel threads [px]
 
-    img_shape: numpy.ndarray (shape=(3,), dtype=int)
+    img_shp: numpy.ndarray (shape=(3,), dtype=int)
         total image shape [px]
 
-    slice_per_dim: numpy.ndarray (shape=(3,), dtype=int)
-        number of image slices along each axis
+    slc_per_dim: numpy.ndarray (shape=(3,), dtype=int)
+        image slices along each axis
 
     ovlp: int
-        slice overlap range
-        along each axis side
+        overlapping range between slices along each axis [px]
 
     flip: bool
-        if True, flip axes coordinates
+        flip axes
 
     Returns
     -------
-    rng: tuple
-        3D slice index ranges
+    start: int
+        start index [px]
 
-    pad: np.ndarray (shape=(3,2), dtype=int)
-        padded boundaries
+    stop: int
+        stop index [px]
+
+    pad: numpy.ndarray (shape=(2,), dtype=int)
+        lower and upper padding ranges [px]
     """
 
-    # adjust original image patch coordinates
-    # and generate padding range matrix
-    dims = len(ax_iter)
-    start = np.zeros((dims,), dtype=int)
-    stop = np.zeros((dims,), dtype=int)
-    pad = np.zeros((dims, 2), dtype=int)
-    slc = tuple()
-    for ax in range(dims):
-        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),)
+    # compute start and stop coordinates
+    start = ax_iter[ax] * slc_shp[ax]
+    stop = start + slc_shp[ax]
 
-    # generate tuple of slice index ranges
-    rng = np.index_exp[slc]
+    # flip coordinates if required
+    if flip:
+        start_tmp = start
+        start = img_shp[ax] - stop
+        stop = img_shp[ax] - start_tmp
 
-    # handle invalid slices
-    for r in rng:
-        if r.start is None:
-            return None, pad
+    # adjust start and stop coordinates
+    start, stop, pad = adjust_axis_range(ax_iter[ax], img_shp[ax], slc_per_dim[ax], start, stop, ovlp=ovlp)
 
-    return rng, pad
+    return start, stop, pad
 
 
-def compute_slice_shape(img_shape, item_size, px_size=None, max_size=1e6, ovlp=0):
+def compute_slice_range(ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip=False):
     """
-    Compute basic image chunk shape depending on its maximum size (in bytes).
+    Compute basic slice coordinates from microscopy volume image.
 
     Parameters
     ----------
-    img_shape: numpy.ndarray (shape=(3,), dtype=int)
-        total image shape [px]
+    ax_iter: tuple
+        iteration counters along each axis
 
-    item_size: int
-        image item size (in bytes)
+    slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the basic image slices
+        analyzed using parallel threads [px]
 
-    px_size: numpy.ndarray (shape=(3,), dtype=float)
-        pixel size [μm]
+    img_shp: numpy.ndarray (shape=(3,), dtype=int)
+        total image shape [px]
 
-    max_size: float
-        maximum memory size (in bytes) of the basic image slices
-        analyzed using parallel threads
+    slc_per_dim: numpy.ndarray (shape=(3,), dtype=int)
+        image slices along each axis
 
     ovlp: int
-        slice overlap range
-        along each axis side
+        overlapping range between slices along each axis [px]
+
+    flip: bool
+        flip axes
 
     Returns
     -------
-    slice_shape: numpy.ndarray (shape=(3,), dtype=int)
-        shape of the basic image slices
-        analyzed using parallel threads [px]
+    rng: np.ndarray
+        3D slice index ranges [px]
 
-    slice_shape_um: numpy.ndarray (shape=(3,), dtype=float)
-        shape of the basic image slices
-        analyzed using parallel threads [μm]
-        (if px_size is provided)
+    pad: np.ndarray (shape=(3,2), dtype=int)
+        slice padding range [px]
     """
-    if len(img_shape) == 4:
-        item_size *= 3
 
-    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)
+    # generate axis range and padding array
+    dim = len(ax_iter)
+    slc = tuple()
+    pad = np.zeros((dim, 2), dtype=np.int64)
+    for ax in range(dim):
+        start, stop, pad[ax] = compute_axis_range(ax, ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=ovlp, flip=flip)
+        slc += (slice(start, stop, 1),)
 
-    if px_size is not None:
-        slice_shape_um = np.multiply(slice_shape, px_size)
-        return slice_shape, slice_shape_um
+    # generate range array
+    rng = np.index_exp[slc]
 
-    else:
-        return slice_shape
+    return rng, pad
 
 
-def compute_overlap_range(smooth_sigma, frangi_sigma, px_rsz_ratio, truncate=2):
+@njit(cache=True)
+def compute_overlap(smooth_sigma, frangi_sigma, px_rsz_ratio=None, trunc=2):
     """
     Compute lateral slice extension range
     for coping with smoothing-related boundary artifacts.
@@ -186,8 +183,7 @@ def compute_overlap_range(smooth_sigma, frangi_sigma, px_rsz_ratio, truncate=2):
     Parameters
     ----------
     smooth_sigma: numpy.ndarray (shape=(3,), dtype=int)
-        3D standard deviation of the low-pass Gaussian filter [px]
-        (applied to the XY plane)
+        3D standard deviation of the smoothing Gaussian filter [px]
 
     frangi_sigma: numpy.ndarray (dtype=int)
         analyzed spatial scales [px]
@@ -195,422 +191,472 @@ def compute_overlap_range(smooth_sigma, frangi_sigma, px_rsz_ratio, truncate=2):
     px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
         3D image resize ratio
 
-    truncate: int
+    trunc: int
         neglect the Gaussian smoothing kernel
         after this many standard deviations
 
     Returns
     -------
     ovlp: int
-        slice overlap range
-        along each axis side
+        overlapping range between image slices along each axis [px]
 
-    rsz_ovlp: numpy.ndarray (shape=(3,), dtype=int)
-        resized image slice overlap range (if px_rsz_ratio is not None)
+    ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int)
+        resized overlapping range between image slices
+        (if px_rsz_ratio is not None) [px]
     """
     if smooth_sigma is not None:
         frangi_sigma = np.concatenate((smooth_sigma, frangi_sigma))
 
-    ovlp = int(np.ceil(2 * truncate * np.max(frangi_sigma)) // 2)
+    ovlp = int(np.ceil(2 * trunc * np.max(frangi_sigma)) // 2)
+    ovlp_rsz = np.multiply(ovlp * np.ones((3,)), px_rsz_ratio).astype(np.int64) if px_rsz_ratio is not None else None
 
-    if px_rsz_ratio is not None:
-        rsz_ovlp = np.multiply(ovlp * np.ones((3,)), px_rsz_ratio).astype(int)
+    return ovlp, ovlp_rsz
 
-        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=-1, jobs=None, ram=None):
+def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0):
     """
-    Compute size and number of the batches of basic microscopy image slices
-    analyzed in parallel.
+    Compute basic image chunk shape depending on its maximum size (in bytes).
 
     Parameters
     ----------
-    frangi_scales: list (dtype=float)
-        analyzed spatial scales in [μm]
-
-    mem_growth_factor: float
-        empirical memory growth factor
-        of the Frangi filtering stage
+    img_shp: numpy.ndarray (shape=(3,), dtype=int)
+        total image shape [px]
 
-    mem_fudge_factor: float
-        memory fudge factor
+    item_sz: int
+        image item size [B]
 
-    min_slice_size: float
-        minimum slice size in [B]
+    px_sz: numpy.ndarray (shape=(3,), dtype=float)
+        pixel size [μm]
 
-    jobs: int
-        number of parallel jobs (threads)
-        used by the Frangi filtering stage
+    slc_sz: float
+        maximum memory size of the basic image slices
+        analyzed using parallel threads [B]
 
-    ram: float
-        maximum RAM available to the Frangi filtering stage [B]
+    ovlp: int
+        overlapping range between image slices along each axis [px]
 
     Returns
     -------
-    slice_batch_size: int
-        slice batch size
+    slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the basic image slices
+        analyzed using parallel threads [px]
 
-    slice_size: float
-        memory size (in bytes) of the basic image slices
-        fed to the Frangi filter
+    slc_shp_um: numpy.ndarray (shape=(3,), dtype=float)
+        shape of the basic image slices
+        analyzed using parallel threads [μm]
+        (if px_size is provided)
     """
-    # maximum RAM not provided: use all
-    if ram is None:
-        ram = psutil.virtual_memory()[1]
+    if len(img_shp) == 4:
+        item_sz *= 3
 
-    # number of logical cores
-    num_cpu = get_available_cores()
-    if jobs is None:
-        jobs = num_cpu
+    side = np.round((slc_sz / item_sz)**(1/3))
+    slc_shp = (side * np.ones((3,)) - 2 * ovlp).astype(np.int64)
+    slc_shp = np.min(np.stack((img_shp[:3], slc_shp)), axis=0)
+    slc_shp_um = np.multiply(slc_shp, px_sz) if px_sz is not None else None
 
-    # number of spatial scales
-    num_scales = len(frangi_scales)
+    return slc_shp, slc_shp_um
 
-    # initialize slice batch size
-    slice_batch_size = np.min([jobs // num_scales, num_cpu]).astype(int)
-    if slice_batch_size == 0:
-        slice_batch_size = 1
 
-    # get image slice size
-    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 = get_slice_size(ram, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales)
+def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1):
+    """
+    Compute the size of the basic microscopy image slices fed to the Frangi filtering stage.
 
-    return slice_batch_size, slice_size
+    Parameters
+    ----------
+    max_ram: float
+        available RAM [B]
 
+    mem_growth: float
+        empirical memory growth factor
+
+    mem_fudge: float
+        memory fudge factor
+
+    batch_sz: int
+        slice batch size
 
-def config_frangi_slicing(img_shape, item_size, px_size, px_size_iso, smooth_sigma, frangi_sigma, mask_lpf,
-                          batch_size, slice_size):
+    ns: int
+        number of spatial scales
+
+    Returns
+    -------
+    slc_sz: float
+        memory size of the basic image slices
+        analyzed using parallel threads [B]
+    """
+    slc_sz = max_ram / (batch_sz * mem_growth * mem_fudge * ns)
+
+    return slc_sz
+
+
+def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi_sigma,
+                        mem_growth=149.7, mem_fudge=1.0, jobs=None, ram=None, shp_thr=7):
     """
-    Image slicing configuration for the parallel Frangi filtering of basic chunks
-    of the input microscopy volume using concurrent threads.
+    Compute size and number of the batches of basic microscopy image slices
+    analyzed in parallel.
 
     Parameters
     ----------
-    img_shape: numpy.ndarray (shape=(3,), dtype=int)
-        total image shape [px]
+    px_sz: numpy.ndarray (shape=(3,), dtype=float)
+        pixel size [μm]
 
-    item_size: int
-        image item size (in bytes)
+    px_sz_iso: int
+        isotropic pixel size [μm]
 
-    px_size: numpy.ndarray (shape=(3,), dtype=float)
-        pixel size [μm]
+    img_shp: numpy.ndarray (shape=(3,), dtype=int)
+        total image shape [px]
 
-    px_size_iso: numpy.ndarray (shape=(3,), dtype=float)
-        adjusted isotropic pixel size [μm]
+    item_sz: int
+        image item size [B]
 
     smooth_sigma: numpy.ndarray (shape=(3,), dtype=int)
-        3D standard deviation of the low-pass Gaussian filter [px]
-        (resolution anisotropy correction)
+        3D standard deviation of the smoothing Gaussian filter [px]
 
-    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
+    frangi_sigma: list (dtype=float)
+        analyzed spatial scales in [μm]
 
-    batch_size: int
-        slice batch size
+    mem_growth: float
+        empirical memory growth factor
+        of the Frangi filtering stage
 
-    slice_size: float
-        maximum memory size (in megabytes) of the basic image slices
-        analyzed iteratively
+    mem_fudge: float
+        memory fudge factor
 
-    Returns
-    -------
-    in_rng_lst: list
-        list of input slice index ranges
+    jobs: int
+        number of parallel jobs (threads)
+        used by the Frangi filtering stage
 
-    in_rng_lst_neu: list
-        list of soma channel slice index ranges
+    ram: float
+        maximum RAM available to the Frangi filtering stage [B]
 
-    out_rng_lst: list
-        list of output slice index ranges
+    shp_thr: int
+        minimum slice side [px]
 
-    in_pad_lst: list
-        list of slice padding arrays
+    Returns
+    -------
+    batch_sz: int
+        slice batch size
 
-    in_slice_shape_um: numpy.ndarray (shape=(3,), dtype=float)
-        shape of the basic image slices analyzed iteratively [μm]
+    in_slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the basic image slices
+        analyzed using parallel threads [px]
 
-    out_slice_shape: numpy.ndarray (shape=(3,), dtype=int)
-        shape of the processed image slices [px]
+    in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the basic image slices
+        analyzed using parallel threads [μm]
 
     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
+    ovlp: int
+        overlapping range between image slices along each axis [px]
 
-    batch_size: int
-        adjusted slice batch size
+    ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int)
+        resized overlapping range between image slices [px]
     """
+    # maximum RAM not provided: use all
+    if ram is None:
+        ram = psutil.virtual_memory()[1]
 
-    # 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 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=ovlp)
+    # number of logical cores
+    num_cpu = get_available_cores()
+    if jobs is None:
+        jobs = num_cpu
 
-    # adjust shapes according to the anisotropy correction
-    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)
+    # number of spatial scales
+    ns = len(frangi_sigma)
 
-    # iteratively define input/output slice 3D ranges
-    slice_per_dim = np.ceil(np.divide(img_shape, in_slice_shape)).astype(int)
-    tot_slice_num = int(np.prod(slice_per_dim))
+    # initialize slice batch size
+    batch_sz = np.min([jobs // ns, num_cpu]).astype(int) + 1
 
-    # initialize empty range lists
-    out_rng_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, 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_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 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:
-                in_rng_lst_neu.append(None)
+    # get pixel resize ratio
+    px_rsz_ratio = np.divide(px_sz, px_sz_iso)
+
+    # compute slice overlap (boundary artifacts suppression)
+    ovlp, ovlp_rsz = compute_overlap(smooth_sigma, frangi_sigma, px_rsz_ratio=px_rsz_ratio)
+
+    # compute the shape of basic microscopy image slices
+    in_slc_shp = np.array([-1])
+    while np.any(in_slc_shp < shp_thr):
+        batch_sz -= 1
+        if batch_sz == 0:
+            raise ValueError(
+                "Basic image slices do not fit the available RAM: decrease spatial scales and/or maximum scale value!")
         else:
-            tot_slice_num -= 1
-
-    # adjust slice batch size
-    if batch_size > tot_slice_num:
-        batch_size = tot_slice_num
+            slc_sz = compute_slice_size(ram, mem_growth, mem_fudge, batch_sz, ns)
+            in_slc_shp, in_slc_shp_um = compute_slice_shape(img_shp, item_sz, px_sz=px_sz, slc_sz=slc_sz, ovlp=ovlp)
 
-    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
+    return batch_sz, in_slc_shp, in_slc_shp_um, px_rsz_ratio, ovlp, ovlp_rsz
 
 
-def crop_slice(img_slice, rng, ovlp=None, flipped=()):
+@njit(cache=True)
+def crop(img, rng, ovlp, flip=()):
     """
     Shrink image slice at total volume boundaries, for overall shape consistency.
 
     Parameters
     ----------
-    img_slice: numpy.ndarray (axis order: (Z,Y,X))
-        image slice
+    img: numpy.ndarray (axis order=(Z,Y,X))
+        3D microscopy image
 
     rng: tuple
         3D slice index range
 
-    ovlp: numpy.ndarray (shape=(3,), dtype=int)
-        slice overlap range (on each axis side)
+    ovlp: numpy.ndarray (axis order=(Z,Y,X), dtype=int)
+        overlapping range between image slices [px]
 
-    flipped: tuple
+    flip: tuple
         flipped axes
 
     Returns
     -------
-    cropped: numpy.ndarray
-        cropped image slice
+    img_out: numpy.ndarray
+        cropped microscopy image
     """
 
-    # delete overlapping slice boundaries
+    # delete overlapping boundaries
     if ovlp is not None:
-        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]]
+        img = img[ovlp[0]:img.shape[0] - ovlp[0], ovlp[1]:img.shape[1] - ovlp[1], ovlp[2]:img.shape[2] - ovlp[2]]
 
-    # check slice shape and output index ranges
-    out_slice_shape = img_slice.shape
-    crop_rng = np.zeros(shape=(3,), dtype=int)
+    # check image shape and output index ranges
+    crop_rng = np.zeros(shape=(3,), dtype=np.int64)
     for s in range(3):
-        crop_rng[s] = out_slice_shape[s] - np.arange(rng[s].start, rng[s].stop, rng[s].step).size
+        crop_rng[s] = img.shape[s] - np.arange(rng[s].start, rng[s].stop, rng[s].step).size
 
-    # crop slice if required
-    cropped = \
-        img_slice[crop_rng[0] or None:, ...] if 0 in flipped else img_slice[:-crop_rng[0] or None, ...]
-    cropped = \
-        cropped[:, crop_rng[1] or None:, ...] if 1 in flipped else cropped[:, :-crop_rng[1] or None, ...]
-    cropped = \
-        cropped[:, :, crop_rng[2] or None:, ...] if 2 in flipped else cropped[:, :, :-crop_rng[2] or None, ...]
+    # crop image if required
+    img_out = img[crop_rng[0] or None:, ...] if 0 in flip else img[:-crop_rng[0] or None, ...]
+    img_out = img_out[:, crop_rng[1] or None:, ...] if 1 in flip else img_out[:, :-crop_rng[1] or None, ...]
+    img_out = img_out[:, :, crop_rng[2] or None:, ...] if 2 in flip else img_out[:, :, :-crop_rng[2] or None, ...]
 
-    return cropped
+    return img_out
 
 
-def crop_slice_lst(img_slice_lst, rng, ovlp=None, flipped=()):
+def crop_lst(img_lst, rng, ovlp=None, flip=()):
     """
     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
+    img_lst: list
+        list of images to be cropped
 
     rng: tuple
         3D slice index range
 
-    ovlp: numpy.ndarray (shape=(3,), dtype=int)
-        slice overlap range (on each axis side)
+    ovlp: int
+        overlapping range between image slices along each axis [px]
 
-    flipped: tuple
+    flip: tuple
         flipped axes
 
     Returns
     -------
-    img_slice_lst: list
+    img_lst: list
         list of cropped image slices
     """
 
-    for s, img_slice in enumerate(img_slice_lst):
-        if img_slice is not None:
-            img_slice_lst[s] = crop_slice(img_slice, rng, ovlp=ovlp, flipped=flipped)
+    for s, img in enumerate(img_lst):
+        if img is not None:
+            img_lst[s] = crop(img, rng, ovlp=ovlp, flip=flip)
 
-    return img_slice_lst
+    return img_lst
 
 
-def get_slice_size(max_ram, mem_growth_factor, mem_fudge_factor, slice_batch_size, num_scales=1):
+def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bg=False, jobs=None):
     """
-    Compute the size of the basic microscopy image slices fed to the Frangi filtering stage.
+    Generate image slice ranges for the Frangi filtering stage.
 
     Parameters
     ----------
-    max_ram: float
-        available RAM [B]
-
-    mem_growth_factor: float
-        empirical memory growth factor
-        of the pipeline stage (Frangi filter or ODF estimation)
+    in_slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the basic image slices
+        analyzed using parallel threads [px]
 
-    mem_fudge_factor: float
-        memory fudge factor
+    img_shp: numpy.ndarray (shape=(3,), dtype=int)
+        total image shape [px]
 
-    slice_batch_size: int
+    batch_sz: int
         slice batch size
 
-    num_scales: int
-        number of spatial scales analyzed in parallel
+    px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
+        3D image resize ratio
+
+    ovlp: int
+        overlapping range between image slices along each axis [px]
+
+    msk_bg: bool
+        if True, mask neuronal bodies
+        in the optionally provided image channel
+
+    jobs: int
+        number of parallel jobs (threads)
 
     Returns
     -------
-    slice_size: float
-        memory size (in bytes) of the basic image slices
-        fed to the pipeline stage
+    in_rng_lst: list
+        list of input slice index ranges
+
+    in_pad_lst: list
+        list of slice padding arrays
+
+    out_rng_lst: list
+        list of output slice index ranges
+
+    bc_rng_lst: list
+        (optional) list of neuronal body slice index ranges
+
+    out_slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the processed image slices [px]
+
+    tot_slc_num: int
+        total number of analyzed image slices
+
+    batch_sz: int
+        adjusted slice batch size
     """
-    slice_size = max_ram / (slice_batch_size * mem_growth_factor * mem_fudge_factor * num_scales)
 
-    return slice_size
+    # adjust output shapes
+    out_slc_shp = np.ceil(np.multiply(px_rsz_ratio, in_slc_shp)).astype(int)
+    out_img_shp = np.ceil(np.multiply(px_rsz_ratio, img_shp)).astype(int)
+
+    # number of image slices
+    slc_per_dim = np.floor(np.divide(img_shp, in_slc_shp)).astype(int)
+    tot_slc_num = int(np.prod(slc_per_dim))
+
+    # number of logical cores
+    if jobs is None:
+        jobs = get_available_cores()
+
+    # initialize empty range lists
+    in_pad_lst = list()
+    in_rng_lst = list()
+    bc_rng_lst = list()
+    out_rng_lst = list()
+    with Parallel(n_jobs=jobs, backend='threading', verbose=0) as parallel:
+        parallel(
+            delayed(generate_slice_range)(z, y, x, in_slc_shp, out_slc_shp, img_shp, out_img_shp, slc_per_dim,
+                                          in_rng_lst, in_pad_lst, out_rng_lst, bc_rng_lst, ovlp=ovlp, msk_bg=msk_bg)
+            for z, y, x in product(range(slc_per_dim[0]), range(slc_per_dim[1]), range(slc_per_dim[2])))
+
+    # adjust slice batch size
+    if batch_sz > tot_slc_num:
+        batch_sz = tot_slc_num
 
+    return in_rng_lst, in_pad_lst, out_rng_lst, bc_rng_lst, out_slc_shp, tot_slc_num, batch_sz
 
-def slice_channel(img, rng, channel, tissue_msk=None, is_tiled=False):
+
+def generate_slice_range(z, y, x, in_slc_shp, out_slc_shp, img_shp, out_img_shp, slc_per_dim,
+                         in_rng_lst, in_pad_lst, out_rng_lst, bc_rng_lst, ovlp=0, msk_bg=False):
     """
-    Slice desired channel from input image volume.
+    Generate image slice range.
 
     Parameters
     ----------
-    img: numpy.ndarray
-        microscopy volume image
+    z: int
+        slice depth
 
-    rng: tuple (dtype=int)
-        3D index ranges
+    y: int
+        slice vertical position
 
-    channel: int
-        image channel axis
+    x: int
+        slice horizontal position
 
-    tissue_msk: numpy.ndarray (dtype=bool)
-        tissue reconstruction binary mask
+    in_slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the basic image slices
+        analyzed using parallel threads [px]
 
-    is_tiled: bool
-        True for tiled reconstructions aligned using ZetaStitcher
+    out_slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the image slices
+        returned by the Frangi filter stage [px]
+
+    img_shp: numpy.ndarray (shape=(3,), dtype=int)
+        total image shape [px]
+
+    out_img_shp: numpy.ndarray (shape=(3,), dtype=int)
+        output image shape [px]
+
+    slc_per_dim: numpy.ndarray (shape=(3,), dtype=int)
+        image slices along each axis
+
+    in_rng_lst: list
+        list of input slice index ranges
+
+    in_pad_lst: list
+        list of slice padding arrays
+
+    out_rng_lst: list
+        list of output slice index ranges
+
+    bc_rng_lst: list
+        (optional) list of neuronal body slice index ranges
+
+    ovlp: int
+        overlapping range between slices along each axis [px]
+
+    msk_bg: bool
+        if True, mask neuronal bodies
+        in the optionally provided image channel
 
     Returns
     -------
-    img_slice: numpy.ndarray
-        sliced image patch
-    
-    tissue_msk_slice: numpy.ndarray (dtype=bool)
-        sliced tissue reconstruction binary mask
+    None
     """
-    z_rng, r_rng, c_rng = rng
 
-    if channel is None:
-        img_slice = img[z_rng, r_rng, c_rng]
-    else:
-        img_slice = img[z_rng, channel, r_rng, c_rng] if is_tiled else img[z_rng, r_rng, c_rng, channel]
+    # index ranges of analyzed image slices (with padding)
+    in_rng, pad = compute_slice_range((z, y, x), in_slc_shp, img_shp, slc_per_dim, ovlp=ovlp)
+    in_rng_lst.append(in_rng)
+    in_pad_lst.append(pad)
 
-    tissue_msk_slice = tissue_msk[r_rng, c_rng] if tissue_msk is not None else None
+    # output index ranges
+    out_rng, _ = compute_slice_range((z, y, x), out_slc_shp, out_img_shp, slc_per_dim)
+    out_rng_lst.append(out_rng)
 
-    return img_slice, tissue_msk_slice
+    # (optional) neuronal body channel
+    if msk_bg:
+        bc_rng, _ = compute_slice_range((z, y, x), in_slc_shp, img_shp, slc_per_dim)
+        bc_rng_lst.append(bc_rng)
+    else:
+        bc_rng_lst.append(None)
 
 
-def adjust_axis_range(img_shape, start, stop, ax, ovlp=0):
+#@jit(cache=True, forceobj=True)
+def slice_image(img, rng, ch, ts_msk=None, is_tiled=False):
     """
-    Trim slice axis range at total image boundaries.
+    Slice desired channel from input image volume.
 
     Parameters
     ----------
-    img_shape: tuple or np.ndarray
-        image shape
+    img: numpy.ndarray
+        3D microscopy image
 
-    start: int
-        slice start coordinate
+    rng: tuple (dtype=int)
+        3D index ranges
 
-    stop: int
-        slice stop coordinate
+    ch: int
+        channel
 
-    ax: int
-        image axis
+    ts_msk: numpy.ndarray (dtype=bool)
+        tissue binary mask
 
-    ovlp: int
-        slice overlap range
-        along each axis side
+    is_tiled: bool
+        True for tiled reconstructions
+        aligned using ZetaStitcher
 
     Returns
     -------
-    start: int
-        adjusted slice start coordinate
-
-    stop: int
-        adjusted slice stop coordinate
+    img_slc: numpy.ndarray
+        image slice
 
-    pad: numpy.ndarray (shape=(2,), dtype=int)
-        lower and upper padding ranges
+    ts_msk_slc: numpy.ndarray (dtype=bool)
+        tissue mask slice
     """
+    z_rng, r_rng, c_rng = rng
 
-    # initialize slice padding array
-    pad = np.zeros((2,), dtype=int)
-
-    # adjust start coordinate
-    start -= ovlp
-    if start < 0:
-        pad[0] = -start
-        start = 0
+    if ch is None:
+        img_slc = img[z_rng, r_rng, c_rng]
+    else:
+        img_slc = img[z_rng, ch, r_rng, c_rng] if is_tiled else img[z_rng, r_rng, c_rng, ch]
 
-    # adjust stop coordinate
-    stop += ovlp
-    if stop > img_shape[ax]:
-        pad[1] = stop - img_shape[ax]
-        stop = img_shape[ax]
+    ts_msk_slc = ts_msk[r_rng, c_rng] if ts_msk is not None else None
 
-    return start, stop, pad
+    return img_slc, ts_msk_slc