diff --git a/foa3d/__main__.py b/foa3d/__main__.py
index ff736b9..a39242c 100644
--- a/foa3d/__main__.py
+++ b/foa3d/__main__.py
@@ -1,7 +1,6 @@
 from foa3d.input import get_cli_parser, load_microscopy_image
-from foa3d.pipeline import parallel_odf_over_scales, parallel_frangi_over_slices
+from foa3d.pipeline import parallel_frangi_over_slices, parallel_odf_over_scales
 from foa3d.printing import print_pipeline_heading
-from foa3d.utils import delete_tmp_folder
 
 
 def foa3d(cli_args):
@@ -15,9 +14,6 @@ def foa3d(cli_args):
     # generate 3D fiber ODF maps over the spatial scales of interest using concurrent workers
     parallel_odf_over_scales(cli_args, save_dirs, out_img['vec'], out_img['iso'], out_img['px_sz'], in_img['name'])
 
-    # delete temporary folder
-    delete_tmp_folder(save_dirs['tmp'], (in_img, out_img))
-
 
 def main():
     # start Foa3D by terminal
diff --git a/foa3d/frangi.py b/foa3d/frangi.py
index 39d9bcc..0dfa719 100644
--- a/foa3d/frangi.py
+++ b/foa3d/frangi.py
@@ -5,8 +5,8 @@
 from joblib import Parallel, delayed
 from scipy import ndimage as ndi
 
-from foa3d.utils import (create_background_mask, create_memory_map,
-                         divide_nonzero, hsv_orient_cmap, rgb_orient_cmap)
+from foa3d.utils import (create_background_mask, divide_nonzero,
+                         hsv_orient_cmap, rgb_orient_cmap)
 
 
 def analyze_hessian_eigen(img, sigma, trunc=4):
@@ -333,7 +333,7 @@ def compute_parall_scaled_orientation(scales_px, img, alpha=0.001, beta=1, gamma
         Hessian eigenvalues sorted by absolute value (ascending order)
     """
     n_scales = len(scales_px)
-    with Parallel(n_jobs=n_scales, prefer='threads', require='sharedmem') as parallel:
+    with Parallel(n_jobs=n_scales, prefer='threads') as parallel:
         par_res = parallel(
                     delayed(compute_scaled_orientation)(
                         scales_px[i], img=img,
@@ -357,12 +357,58 @@ def compute_parall_scaled_orientation(scales_px, img, alpha=0.001, beta=1, gamma
     return frangi, fbr_vec, eigenval
 
 
-def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False):
+def init_frangi_arrays(in_img, cli_args, cfg, slc_shp, rsz_ratio):
     """
-    Initialize the output datasets of the Frangi filtering stage.
+    Initialize the output datasets of the Frangi filter stage.
 
     Parameters
     ----------
+    in_img: dict
+        input image dictionary
+
+            data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X))
+                3D microscopy image
+
+            ts_msk: numpy.ndarray (dtype=bool)
+                tissue reconstruction binary mask
+
+            ch_ax: int
+                RGB image channel axis (either 1, 3, or None for grayscale images)
+
+            fb_ch: int
+                neuronal fibers channel
+
+            bc_ch: int
+                brain cell soma channel
+
+            msk_bc: bool
+                if True, mask neuronal bodies within
+                the optionally provided channel
+
+            psf_fwhm: numpy.ndarray (shape=(3,), dtype=float)
+                3D FWHM of the PSF [μm]
+
+            px_sz: numpy.ndarray (shape=(3,), dtype=float)
+                pixel size [μm]
+
+            name: str
+                name of the 3D microscopy image
+
+            is_vec: bool
+                vector field flag
+
+            shape: numpy.ndarray (shape=(3,), dtype=int)
+                total image shape
+
+            shape_um: numpy.ndarray (shape=(3,), dtype=float)
+                total image shape [μm]
+
+            item_sz: int
+                image item size [B]
+
+    cli_args: see ArgumentParser.parse_args
+        populated namespace of command line arguments
+
     cfg: dict
         Frangi filter configuration
 
@@ -406,9 +452,6 @@ def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False):
             exp_all: bool
                 export all images
 
-    img_shp: numpy.ndarray (shape=(3,), dtype=int)
-        3D image shape [px]
-
     slc_shp: numpy.ndarray (shape=(3,), dtype=int)
         shape of the basic image slices
         analyzed using parallel threads [px]
@@ -416,68 +459,62 @@ def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False):
     rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
         3D image resize ratio
 
-    tmp_dir: str
-        temporary file directory
-
-    msk_bc: bool
-        if True, mask neuronal bodies
-        in the optionally provided image channel
-
     Returns
     -------
     out_img: dict
         output image dictionary
 
-        vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
-            initialized fiber orientation 3D image
+            vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
+                initialized fiber orientation 3D image
 
-        clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8)
-            initialized orientation colormap image
+            clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8)
+                initialized orientation colormap image
 
-        fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            initialized fractional anisotropy image
+            fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
+                initialized fractional anisotropy image
 
-        frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            initialized Frangi-enhanced image
+            frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
+                initialized Frangi-enhanced image
 
-        iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            initialized fiber image (isotropic resolution)
+            iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
+                initialized fiber image (isotropic resolution)
 
-        fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            initialized fiber mask image
+            fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
+                initialized fiber mask image
 
-        bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            initialized soma mask image
+            bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
+                initialized soma mask image
 
-    z_sel: NumPy slice object
-        selected z-depth range
+            z_out: NumPy slice object
+                selected z-depth range
     """
     # shape copies
-    img_shp = img_shp.copy()
-    slc_shp = slc_shp.copy()
+    img_shp = in_img['shape'].copy()
+    slc_shp = np.ceil(np.multiply(rsz_ratio, slc_shp)).astype(int)
 
-    # adapt output z-axis shape if required
-    z_min, z_max = cfg['z_rng']
+    # adapt output z-axis range if required
+    z_min = int(np.floor(cli_args.z_min / np.max(in_img['px_sz'])))
+    z_max = int(np.ceil(cli_args.z_max / np.max(in_img['px_sz']))) if cli_args.z_max is not None else cli_args.z_max
     if z_min != 0 or z_max is not None:
         if z_max is None:
             z_max = slc_shp[0]
         img_shp[0] = z_max - z_min
-    z_sel = slice(z_min, z_max, 1)
+    z_out = slice(z_min, z_max, 1)
 
     # output shape
     img_dims = len(img_shp)
     tot_shp = tuple(np.ceil(rsz_ratio * img_shp).astype(int))
 
     # fiber channel arrays
-    iso_fbr_img = create_memory_map(tot_shp, dtype='uint8', name='iso_fiber', tmp_dir=tmp_dir)
+    iso_fbr_img = np.zeros(tot_shp, dtype='uint8')
     if cfg['exp_all']:
-        frangi_img = create_memory_map(tot_shp, dtype='uint8', name='frangi', tmp_dir=tmp_dir)
-        fbr_msk = create_memory_map(tot_shp, dtype='uint8', name='fbr_msk', tmp_dir=tmp_dir)
-        fa_img = create_memory_map(tot_shp, dtype='float32', name='fa', tmp_dir=tmp_dir)
+        frangi_img = np.zeros(tot_shp, dtype='uint8')
+        fbr_msk = np.zeros(tot_shp, dtype='uint8')
+        fa_img = np.zeros(tot_shp, dtype='float32')
 
         # soma channel array
-        if msk_bc:
-            bc_msk = create_memory_map(tot_shp, dtype='uint8', name='bc_msk', tmp_dir=tmp_dir)
+        if in_img['msk_bc']:
+            bc_msk = np.zeros(tot_shp, dtype='uint8')
         else:
             bc_msk = None
     else:
@@ -485,20 +522,14 @@ def init_frangi_arrays(cfg, img_shp, slc_shp, rsz_ratio, tmp_dir, msk_bc=False):
 
     # fiber orientation arrays
     vec_shape = tot_shp + (img_dims,)
-    fbr_vec_img = create_memory_map(vec_shape, dtype='float32', name='fiber_vec', tmp_dir=tmp_dir)
-    fbr_vec_clr = create_memory_map(vec_shape, dtype='uint8', name='fiber_cmap', tmp_dir=tmp_dir)
+    fbr_vec_img = np.zeros(vec_shape, dtype='float32')
+    fbr_vec_clr = np.zeros(vec_shape, dtype='uint8')
 
     # fill output image dictionary
-    out_img = dict()
-    out_img['vec'] = fbr_vec_img
-    out_img['clr'] = fbr_vec_clr
-    out_img['fa'] = fa_img
-    out_img['frangi'] = frangi_img
-    out_img['iso'] = iso_fbr_img
-    out_img['fbr_msk'] = fbr_msk
-    out_img['bc_msk'] = bc_msk
+    out_img = {'vec': fbr_vec_img, 'clr': fbr_vec_clr, 'fa': fa_img, 'frangi': frangi_img,
+               'iso': iso_fbr_img, 'fbr_msk': fbr_msk, 'bc_msk': bc_msk, 'px_sz': cfg['px_sz'], 'z_out': z_out}
 
-    return out_img, z_sel
+    return out_img
 
 
 def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=False, hsv=False):
@@ -527,69 +558,91 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=Fals
         (i.e., negative contrast polarity)
 
     hsv: bool
+        generate an HSV colormap of 3D fiber orientations
 
     Returns
     -------
-    orient_slc: dict
-        slice orientation dictionary
+    out_slc: dict
+        slice output data dictionary
 
-        vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
-            3D fiber orientation field
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+                3D fiber orientation field
 
-        clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
-            orientation colormap image
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap image
 
-        fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-            fractional anisotropy image
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy image
 
-    frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
-        Frangi's vesselness likelihood image
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
+                Frangi's vesselness likelihood image
     """
-    # single-scale vesselness analysis
-    if len(scales_px) == 1:
-        frangi, fbr_vec, eigenval = \
-            compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark)
-
-    # parallel scaled vesselness analysis
-    else:
-        frangi, fbr_vec, eigenval = \
-            compute_parall_scaled_orientation(scales_px, img, alpha=alpha, beta=beta, gamma=gamma, dark=dark)
-
-    # generate fractional anisotropy image
+    # single-scale or parallel multi-scale vesselness analysis
+    ns = len(scales_px)
+    frangi = np.zeros((ns,) + img.shape, dtype='float32')
+    eigenvec = np.zeros((ns,) + img.shape + (3,), dtype='float32')
+    eigenval = np.zeros((ns,) + img.shape + (3,), dtype='float32')
+    for s in range(ns):
+        frangi[s], eigenvec[s], eigenval[s] = \
+            compute_scaled_orientation(scales_px[s], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark)
+
+    # get maximum response across the requested scales
+    best_idx = np.expand_dims(np.argmax(frangi, axis=0), axis=0)
+    frangi = np.take_along_axis(frangi, best_idx, axis=0).squeeze(axis=0)
+
+    # select fiber orientation vectors (and the associated eigenvalues) among different scales
+    best_idx = np.expand_dims(best_idx, axis=-1)
+    eigenval = np.take_along_axis(eigenval, best_idx, axis=0).squeeze(axis=0)
+    fbr_vec = np.take_along_axis(eigenvec, best_idx, axis=0).squeeze(axis=0)
+
+    # compute fractional anisotropy image and fiber orientation color map
     fa = compute_fractional_anisotropy(eigenval)
-
-    # generate fiber orientation color map
     fbr_clr = hsv_orient_cmap(fbr_vec) if hsv else rgb_orient_cmap(fbr_vec)
 
-    # fill orientation dictionary
-    orient_slc = dict()
-    orient_slc['vec_slc'] = fbr_vec
-    orient_slc['clr_slc'] = fbr_clr
-    orient_slc['fa_slc'] = fa
+    # fill slice output dictionary
+    out_slc = {'vec': fbr_vec, 'clr': fbr_clr, 'fa': fa, 'frangi': frangi}
 
-    return orient_slc, frangi
+    return out_slc
 
 
-def mask_background(img, orient, method='yen', invert=False, ts_msk=None):
+def mask_background(out_slc, ref_img=None, method='yen', invert=False, ornt_keys=('vec', 'clr', 'fa')):
     """
-    Mask fiber orientation arrays.
+    Mask fiber orientation data arrays.
 
     Parameters
     ----------
-    img: numpy.ndarray (axis order=(Z,Y,X))
-        fiber (or neuron) fluorescence 3D image
+    out_slc: dict
+        slice output dictionary
+
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+                3D fiber orientation field
+
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap
+
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy
+
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image slice (fiber probability image)
+
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image slice
 
-    orient: dict
-        slice orientation dictionary
+            ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                tissue mask slice
 
-            vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
-                fiber orientation vector slice
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask slice
 
-            clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
-                orientation colormap slice
+            bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool)
+                brain cell mask slice
 
-            fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
-                fractional anisotropy slice
+            rng: NumPy slice object
+                output range
+
+    ref_img: numpy.ndarray (axis order=(Z,Y,X)
+        reference image used for thresholding
 
     method: str
         thresholding method (refer to skimage.filters)
@@ -597,45 +650,38 @@ def mask_background(img, orient, method='yen', invert=False, ts_msk=None):
     invert: bool
         mask inversion flag
 
-    ts_msk: numpy.ndarray (dtype=bool)
-        tissue reconstruction binary mask
+    ornt_keys: tuple
+        fiber orientation data keys
 
     Returns
     -------
-    orient: dict
-        (masked) slice orientation dictionary
-
-            vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
-                fiber orientation vector slice (masked)
-
-            clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
-                orientation colormap slice (masked)
-
-            fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
-                fractional anisotropy slice (masked)
+    out_slc: dict
 
     bg: numpy.ndarray (axis order=(Z,Y,X), dtype=bool)
         background mask
     """
     # generate background mask
-    bg = create_background_mask(img, method=method)
+    if ref_img is None:
+        ref_img = out_slc['frangi']
+    bg = create_background_mask(ref_img, method=method)
 
     # apply tissue reconstruction mask, when provided
-    if ts_msk is not None:
-        bg = np.logical_or(bg, np.logical_not(ts_msk))
+    if out_slc['ts_msk'] is not None:
+        bg = np.logical_or(bg, np.logical_not(out_slc['ts_msk']))
 
     # invert mask
     if invert:
         bg = np.logical_not(bg)
 
     # apply mask to input orientation data dictionary
-    for key in orient.keys():
-        if orient[key].ndim == 3:
-            orient[key][bg] = 0
-        else:
-            orient[key][bg, :] = 0
+    for key in out_slc.keys():
+        if key in ornt_keys:
+            if out_slc[key].ndim == 3:
+                out_slc[key][bg] = 0
+            else:
+                out_slc[key][bg, :] = 0
 
-    return orient, bg
+    return out_slc, bg
 
 
 def reject_vesselness_background(vesselness, eigen2, eigen3, dark):
@@ -707,83 +753,94 @@ def sort_eigen(eigenval, eigenvec, axis=-1):
     return sorted_eigenval, sorted_eigenvec
 
 
-def write_frangi_arrays(out_rng, z_sel, iso_slc, frangi_slc, fbr_msk_slc, bc_msk_slc, vec_slc, clr_slc, fa_slc,
-                        vec, clr, fa, frangi, iso, fbr_msk, bc_msk):
+def write_frangi_arrays(out_img, out_slc):
     """
     Fill the memory-mapped output arrays of the Frangi filter stage.
 
     Parameters
     ----------
-    out_rng: tuple
-        3D slice output index range
+    out_img: dict
+        output image dictionary
+
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32)
+                fiber orientation vector field
+
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap image
 
-    z_sel: NumPy slice object
-        selected z-depth range
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy image
 
-    iso_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        isotropic fiber image slice
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image (fiber probability image)
 
-    frangi_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32)
-        Frangi-enhanced image slice
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image
 
-    fbr_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        fiber mask image slice
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask image
 
-    bc_msk_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        soma mask image slice
+            bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                soma mask image
+
+            px_sz: numpy.ndarray (shape=(3,), dtype=float)
+                output pixel size [μm]
 
-    vec_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32)
-        fiber orientation vector image slice
+    out_slc: dict
+        slice output dictionary
 
-    clr_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        orientation colormap image slice
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+                3D fiber orientation field
 
-    fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=float32)
-        fractional anisotropy image slice
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap
 
-    vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
-        fiber orientation vector field
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy
 
-    clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8)
-        orientation colormap image
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image slice (fiber probability image)
 
-    fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-        fractional anisotropy image
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image slice
 
-    frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-        Frangi-enhanced image (fiber probability image)
+            ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                tissue mask slice
 
-    iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-        isotropic fiber image
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask slice
 
-    fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-        fiber mask image
+            bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool)
+                brain cell mask slice
 
-    bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-        soma mask image
+            rng: NumPy slice object
+                output range
 
     Returns
     -------
     None
     """
-    # fill memory-mapped output arrays
-    vec_rng_out = tuple(np.append(out_rng, slice(0, 3, 1)))
-    vec[vec_rng_out] = vec_slc[z_sel, ...]
-    clr[vec_rng_out] = clr_slc[z_sel, ...]
-    iso[out_rng] = iso_slc[z_sel, ...].astype(np.uint8)
-
-    # optional output images: Frangi filter response
-    if frangi is not None:
-        frangi[out_rng] = (255 * frangi_slc[z_sel, ...]).astype(np.uint8)
-
-    # optional output images: fractional anisotropy
-    if fa is not None:
-        fa[out_rng] = fa_slc[z_sel, ...]
-
-    # optional output images: fiber mask
-    if fbr_msk is not None:
-        fbr_msk[out_rng] = (255 * (1 - fbr_msk_slc[z_sel, ...])).astype(np.uint8)
-
-    # optional output images: neuronal soma mask
-    if bc_msk is not None:
-        bc_msk[out_rng] = (255 * bc_msk_slc[z_sel, ...]).astype(np.uint8)
+    # write valid slice output
+    if out_slc is not None:
+        out_rng = out_slc['rng']
+        z_out = out_img['z_out']
+        vec_rng_out = tuple(np.append(out_rng, slice(0, 3, 1)))
+        out_img['vec'][vec_rng_out] = out_slc['vec'][z_out, ...]
+        out_img['clr'][vec_rng_out] = out_slc['clr'][z_out, ...]
+        out_img['iso'][out_rng] = out_slc['iso'][z_out, ...].astype(np.uint8)
+
+        # optional output images: fractional anisotropy
+        if out_img['fa'] is not None:
+            out_img['fa'][out_rng] = out_slc['fa'][z_out, ...]
+
+        # optional output images: Frangi filter response
+        if out_img['frangi'] is not None:
+            out_img['frangi'][out_rng] = (255 * out_slc['frangi'][z_out, ...]).astype(np.uint8)
+
+        # optional output images: fiber mask
+        if out_img['fbr_msk'] is not None:
+            out_img['fbr_msk'][out_rng] = (255 * (1 - out_slc['fbr_msk'][z_out, ...])).astype(np.uint8)
+
+        # optional output images: neuronal soma mask
+        if out_img['bc_msk'] is not None:
+            out_img['bc_msk'][out_rng] = (255 * out_slc['bc_msk'][z_out, ...]).astype(np.uint8)
diff --git a/foa3d/input.py b/foa3d/input.py
index a3f61e9..726565f 100644
--- a/foa3d/input.py
+++ b/foa3d/input.py
@@ -15,8 +15,8 @@
 from foa3d.preprocessing import config_anisotropy_correction
 from foa3d.printing import (color_text, print_flsh, print_image_info,
                             print_import_time)
-from foa3d.utils import (create_background_mask, create_memory_map,
-                         detect_ch_axis, get_item_bytes, get_config_label)
+from foa3d.utils import (create_background_mask, detect_ch_axis, get_item_bytes,
+                         get_config_label)
 
 
 class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawTextHelpFormatter):
@@ -389,15 +389,15 @@ def get_frangi_config(cli_args, in_img):
 
     """
     # initialize Frangi filter configuration dictionary
-    frangi_cfg = dict()
+    frangi_cfg = {}
 
     # preprocessing configuration (in-plane smoothing)
-    frangi_cfg['smooth_sd'], px_sz_iso = config_anisotropy_correction(in_img['px_sz'], in_img['psf_fwhm'])
+    frangi_cfg['smooth_sd'], frangi_cfg['px_sz'] = config_anisotropy_correction(in_img['px_sz'], in_img['psf_fwhm'])
 
     # Frangi filter parameters
     frangi_cfg['alpha'], frangi_cfg['beta'], frangi_cfg['gamma'] = (cli_args.alpha, cli_args.beta, cli_args.gamma)
     frangi_cfg['scales_um'] = np.array(cli_args.scales)
-    frangi_cfg['scales_px'] = frangi_cfg['scales_um'] / px_sz_iso[0]
+    frangi_cfg['scales_px'] = frangi_cfg['scales_um'] / np.max(frangi_cfg['px_sz'])
 
     # Frangi filter response threshold
     frangi_cfg['fb_thr'] = cli_args.fb_thr
@@ -413,7 +413,7 @@ def get_frangi_config(cli_args, in_img):
     # export all flag
     frangi_cfg['exp_all'] = cli_args.exp_all
 
-    return frangi_cfg, px_sz_iso
+    return frangi_cfg
 
 
 def get_resolution(cli_args):
@@ -434,8 +434,8 @@ def get_resolution(cli_args):
         3D PSF FWHM [μm]
     """
     # create pixel and psf size arrays
-    px_sz = np.array([cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy])
-    psf_fwhm = np.array([cli_args.psf_fwhm_z, cli_args.psf_fwhm_y, cli_args.psf_fwhm_x])
+    px_sz = (cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy)
+    psf_fwhm = (cli_args.psf_fwhm_z, cli_args.psf_fwhm_y, cli_args.psf_fwhm_x)
 
     return px_sz, psf_fwhm
 
@@ -513,30 +513,21 @@ def load_microscopy_image(cli_args):
     # import fiber orientation vector data
     tic = perf_counter()
     if is_vec:
-        in_img = load_orient(img_path, img_name, img_fmt, tmp_dir=save_dirs['tmp'])
+        in_img = load_orient(img_path, img_fmt)
 
     # import raw 3D microscopy image
     else:
-        in_img = load_raw(img_path, img_name, img_fmt, psf_fwhm,
-                          is_tiled=is_tiled, tmp_dir=save_dirs['tmp'],
-                          msk_mip=msk_mip, msk_bc=msk_bc, fb_ch=fb_ch, bc_ch=bc_ch)
+        in_img = load_raw(img_path, img_fmt, psf_fwhm,
+                          is_tiled=is_tiled, msk_mip=msk_mip, msk_bc=msk_bc, fb_ch=fb_ch, bc_ch=bc_ch)
 
-    # add image pixel size to input data dictionary
-    in_img['px_sz'] = px_sz
-
-    # add image name to input data dictionary
-    in_img['name'] = img_name
-
-    # add vector field flag
-    in_img['is_vec'] = is_vec
+    # populate input image dictionary
+    in_img.update({'px_sz': px_sz, 'name': img_name, 'is_vec': is_vec})
 
     # get microscopy image size information
     in_img = get_image_size(in_img)
 
-    # print import time
+    # print input data information
     print_import_time(tic)
-
-    # print volume image shape and resolution
     if not is_vec:
         print_image_info(in_img)
     else:
@@ -545,7 +536,7 @@ def load_microscopy_image(cli_args):
     return in_img, save_dirs
 
 
-def load_orient(img_path, img_name, img_fmt, tmp_dir=None):
+def load_orient(img_path, img_fmt):
     """
     Load array of 3D fiber orientations.
 
@@ -554,21 +545,15 @@ def load_orient(img_path, img_name, img_fmt, tmp_dir=None):
     img_path: str
         path to the 3D microscopy image
 
-    img_name: str
-        name of the 3D microscopy image
-
     img_fmt: str
         format of the 3D microscopy image
 
-    tmp_dir: str
-        temporary file directory
-
     Returns
     -------
     in_img: dict
         input image dictionary
 
-            data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X))
+            data: numpy.ndarray (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X))
                 3D microscopy image
 
             ts_msk: numpy.ndarray (dtype=bool)
@@ -589,10 +574,6 @@ def load_orient(img_path, img_name, img_fmt, tmp_dir=None):
         if ch_ax != 3:
             vec_img = np.moveaxis(vec_img, ch_ax, -1)
 
-        # memory-map the input TIFF image
-        vec_img = create_memory_map(vec_img.shape, dtype=vec_img.dtype, name=img_name,
-                                    tmp_dir=tmp_dir, arr=vec_img, mmap_mode='r')
-
     # check array shape
     if vec_img.ndim != 4:
         raise ValueError('Invalid 3D fiber orientation dataset (ndim != 4)!')
@@ -600,16 +581,13 @@ def load_orient(img_path, img_name, img_fmt, tmp_dir=None):
         print_flsh(f"Loading {img_path} orientation vector field...\n")
 
         # populate input image dictionary
-        in_img = dict()
-        in_img['data'] = vec_img
-        in_img['ts_msk'] = None
-        in_img['ch_ax'] = None
+        in_img = {'data': vec_img, 'ts_msk': None, 'ch_ax': None}
 
         return in_img
 
 
-def load_raw(img_path, img_name, img_fmt, psf_fwhm,
-             is_tiled=False, tmp_dir=None, msk_mip=False, msk_bc=False, fb_ch=1, bc_ch=0):
+def load_raw(img_path, img_fmt, psf_fwhm,
+             is_tiled=False, msk_mip=False, msk_bc=False, fb_ch=1, bc_ch=0):
     """
     Load 3D microscopy image.
 
@@ -618,9 +596,6 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm,
     img_path: str
         path to the 3D microscopy image
 
-    img_name: str
-        name of the 3D microscopy image
-
     img_fmt: str
         format of the 3D microscopy image
 
@@ -630,9 +605,6 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm,
     is_tiled: bool
         True for tiled reconstructions aligned using ZetaStitcher
 
-    tmp_dir: str
-        temporary file directory
-
     msk_mip: bool
         apply tissue reconstruction mask (binarized MIP)
 
@@ -651,7 +623,7 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm,
     in_img: dict
         input image dictionary
 
-            data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X))
+            data: numpy.ndarray (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X))
                 3D microscopy image
 
             ts_msk: numpy.ndarray (dtype=bool)
@@ -690,20 +662,14 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm,
         else:
             raise ValueError('Unsupported image format!')
 
-        # create image memory map
-        img = create_memory_map(img.shape, dtype=img.dtype, name=img_name, tmp_dir=tmp_dir, arr=img, mmap_mode='r')
-
     # detect channel axis (RGB images)
     ch_ax = detect_ch_axis(img)
 
     # generate tissue reconstruction mask
     if msk_mip:
         dims = len(img.shape)
-
-        # grayscale image
         if dims == 3:
             img_fbr = img
-        # RGB image
         elif dims == 4:
             img_fbr = img[:, fb_ch, :, :] if ch_ax == 1 else img[..., fb_ch]
         else:
@@ -720,13 +686,7 @@ def load_raw(img_path, img_name, img_fmt, psf_fwhm,
         ts_msk = None
 
     # populate input image dictionary
-    in_img = dict()
-    in_img['data'] = img
-    in_img['ts_msk'] = ts_msk
-    in_img['ch_ax'] = ch_ax
-    in_img['fb_ch'] = fb_ch
-    in_img['bc_ch'] = bc_ch
-    in_img['msk_bc'] = msk_bc
-    in_img['psf_fwhm'] = psf_fwhm
+    in_img = {'data': img, 'ts_msk': ts_msk, 'ch_ax': ch_ax, 'fb_ch': fb_ch,
+              'bc_ch': bc_ch,'msk_bc': msk_bc, 'psf_fwhm': psf_fwhm}
 
     return in_img
diff --git a/foa3d/odf.py b/foa3d/odf.py
index 89a42fa..16cf071 100644
--- a/foa3d/odf.py
+++ b/foa3d/odf.py
@@ -4,7 +4,7 @@
 from skimage.transform import resize
 
 from foa3d.spharm import fiber_vectors_to_sph_harm, get_sph_harm_ncoeff
-from foa3d.utils import create_memory_map, normalize_image, transform_axes
+from foa3d.utils import normalize_image, transform_axes
 
 
 def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tensor_eigen, vxl_side, odf_norm,
@@ -230,7 +230,7 @@ def generate_odf_background(bg_mrtrix, fbr_vec, vxl_side, iso_fbr=None):
             resize(tmp_slice, output_shape=new_shape, anti_aliasing=True, preserve_range=True)
 
 
-def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False):
+def init_odf_arrays(vec_img_shp, odf_scale, odf_deg=6, exp_all=False):
     """
     Initialize the output datasets of the ODF analysis stage.
 
@@ -285,36 +285,31 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False):
 
     # create downsampled background memory map
     bg_shp = tuple(np.flip(odi_shp))
-    bg_mrtrix = create_memory_map(bg_shp, dtype='uint8', name=f'bg_tmp{odf_scale}', tmp_dir=tmp_dir)
+    bg_mrtrix = np.zeros(bg_shp, dtype='uint8')
 
     # create ODF memory map
     nc = get_sph_harm_ncoeff(odf_deg)
     odf_shp = odi_shp + (nc,)
-    odf = create_memory_map(odf_shp, dtype='float32', name=f'odf_tmp{odf_scale}', tmp_dir=tmp_dir)
+    odf = np.zeros(odf_shp, dtype='float32')
 
     # create orientation tensor memory map
     vec_tensor_shape = odi_shp + (3,)
-    vec_tensor_eigen = create_memory_map(vec_tensor_shape, dtype='float32',
-                                         name=f'tensor_tmp{odf_scale}', tmp_dir=tmp_dir)
+    vec_tensor_eigen = np.zeros(vec_tensor_shape, dtype='float32')
 
     # fiber density memory map
-    fbr_dnst = create_memory_map(odi_shp, dtype='float32', name=f'fbr_dnst_tmp{odf_scale}', tmp_dir=tmp_dir)
+    fbr_dnst = np.zeros(odi_shp, dtype='float32')
 
     # create ODI memory maps
-    odi_tot = create_memory_map(odi_shp, dtype='float32', name=f'odi_tot_tmp{odf_scale}', tmp_dir=tmp_dir)
+    odi_tot = np.zeros(odi_shp, dtype='float32')
     if exp_all:
-        odi_pri = create_memory_map(odi_shp, dtype='float32', name=f'odi_pri_tmp{odf_scale}', tmp_dir=tmp_dir)
-        odi_sec = create_memory_map(odi_shp, dtype='float32', name=f'odi_sec_tmp{odf_scale}', tmp_dir=tmp_dir)
-        odi_anis = create_memory_map(odi_shp, dtype='float32', name=f'odi_anis_tmp{odf_scale}', tmp_dir=tmp_dir)
+        odi_pri = np.zeros(odi_shp, dtype='float32')
+        odi_sec = np.zeros(odi_shp, dtype='float32')
+        odi_anis = np.zeros(odi_shp, dtype='float32')
     else:
         odi_pri, odi_sec, odi_anis = (None, None, None)
 
     # fill output image dictionary
-    odi = dict()
-    odi['odi_pri'] = odi_pri
-    odi['odi_sec'] = odi_sec
-    odi['odi_tot'] = odi_tot
-    odi['odi_anis'] = odi_anis
+    odi = {'odi_pri': odi_pri, 'odi_sec': odi_sec, 'odi_tot': odi_tot, 'odi_anis': odi_anis}
 
     return odf, odi, fbr_dnst, bg_mrtrix, vec_tensor_eigen
 
diff --git a/foa3d/output.py b/foa3d/output.py
index c58467b..eaf2ba0 100644
--- a/foa3d/output.py
+++ b/foa3d/output.py
@@ -1,5 +1,4 @@
 import psutil
-import tempfile
 
 from datetime import datetime
 from os import makedirs, mkdir, path
@@ -35,7 +34,7 @@ def create_save_dirs(img_path, img_name, cli_args, is_vec=False):
     -------
     save_dirs: dict
         saving directories
-        ('frangi': Frangi filter, 'odf': ODF analysis, 'tmp': temporary files)
+        ('frangi': Frangi filter, 'odf': ODF analysis)
     """
     # get current time
     time_stamp = datetime.now().strftime("%Y%m%d-%H%M%S")
@@ -69,10 +68,6 @@ def create_save_dirs(img_path, img_name, cli_args, is_vec=False):
     else:
         save_dirs['odf'] = None
 
-    # create temporary directory
-    tmp_dir = tempfile.mkdtemp(dir=base_out_dir)
-    save_dirs['tmp'] = tmp_dir
-
     return save_dirs
 
 
@@ -113,7 +108,7 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None):
 
     # check output format
     fmt = fmt.lower()
-    if fmt == 'tif' or fmt == 'tiff':
+    if fmt in ('tif', 'tiff'):
 
         # retrieve image pixel size
         px_sz_z, px_sz_y, px_sz_x = px_sz
@@ -144,7 +139,7 @@ def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', ram=None):
         raise ValueError("Unsupported data format!!!")
 
 
-def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None):
+def save_frangi_arrays(save_dir, img_name, out_img, ram=None):
     """
     Save the output arrays of the Frangi filter stage to TIF files.
 
@@ -178,8 +173,8 @@ def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None):
         bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
             neuron mask image
 
-    px_sz: numpy.ndarray (shape=(3,), dtype=float)
-        pixel size (Z,Y,X) [μm]
+        px_sz: numpy.ndarray (shape=(3,), dtype=float)
+            pixel size (Z,Y,X) [μm]
 
     ram: float
         maximum RAM available
@@ -190,8 +185,8 @@ def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None):
     """
     # loop over output image dictionary fields and save to TIFF
     for img_key in out_img.keys():
-        if out_img[img_key] is not None:
-            save_array(f'{img_key}_{img_name}', save_dir, out_img[img_key], px_sz, ram=ram)
+        if isinstance(out_img[img_key], np.ndarray) and out_img[img_key] is not None:
+            save_array(f'{img_key}_{img_name}', save_dir, out_img[img_key], out_img['px_sz'], ram=ram)
 
     # print output directory
     print_flsh(f"\nFrangi filter arrays saved to: {save_dir}\n")
@@ -247,28 +242,21 @@ def save_odf_arrays(save_dir, img_name, odf_scale_um, px_sz, odf, bg, fbr_dnst,
     sbfx = f'{odf_scale_um}_{img_name}'
     save_array(f'bg_mrtrixview_sv{sbfx}', save_dir, bg, fmt='nii')
     save_array(f'odf_mrtrixview_sv{sbfx}', save_dir, odf, fmt='nii')
-    del bg
-    del odf
 
     # save fiber density
     save_array(f'fbr_dnst_sv{sbfx}', save_dir, fbr_dnst, px_sz)
-    del fbr_dnst
 
     # save total orientation dispersion
     save_array(f'odi_tot_sv{sbfx}', save_dir, odi_tot, px_sz)
-    del odi_tot
 
     # save primary orientation dispersion
     if odi_pri is not None:
         save_array(f'odi_pri_sv{sbfx}', save_dir, odi_pri, px_sz)
-        del odi_pri
 
     # save secondary orientation dispersion
     if odi_sec is not None:
         save_array(f'odi_sec_sv{sbfx}', save_dir, odi_sec, px_sz)
-        del odi_sec
 
     # save orientation dispersion anisotropy
     if odi_anis is not None:
         save_array(f'odi_anis_sv{sbfx}', save_dir, odi_anis, px_sz)
-        del odi_anis
diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py
index f79fed9..b6d2f51 100644
--- a/foa3d/pipeline.py
+++ b/foa3d/pipeline.py
@@ -1,7 +1,6 @@
-from time import perf_counter
-
 import numpy as np
 from joblib import Parallel, delayed
+from time import perf_counter
 
 from foa3d.frangi import (frangi_filter, init_frangi_arrays, mask_background,
                           write_frangi_arrays)
@@ -10,10 +9,10 @@
                        init_odf_arrays)
 from foa3d.output import save_frangi_arrays, save_odf_arrays
 from foa3d.preprocessing import correct_anisotropy
-from foa3d.printing import (print_flsh, print_frangi_info, print_odf_info,
-                            print_frangi_progress)
-from foa3d.slicing import (check_background, config_frangi_batch,
-                           generate_slice_lists, crop, crop_lst, slice_image)
+from foa3d.printing import (print_flsh, print_frangi_config,
+                            print_frangi_progress, print_odf_info)
+from foa3d.slicing import (check_background, get_slice_config,
+                           slice_generator, crop, crop_img_dict)
 from foa3d.spharm import get_sph_harm_norm_factors
 from foa3d.utils import get_available_cores
 
@@ -78,157 +77,82 @@ def parallel_frangi_over_slices(cli_args, save_dirs, in_img):
     out_img: dict
         output image dictionary
 
-        fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
-            fiber orientation vector field
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float32)
+                fiber orientation vector field
 
-        fbr_clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8)
-            orientation colormap image
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap image
 
-        fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            fractional anisotropy image
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy image
 
-        frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            Frangi-enhanced image (fiber probability image)
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image (fiber probability image)
 
-        iso_fbr_img: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            isotropic fiber image
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image
 
-        fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            fiber mask image
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask image
 
-        bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-            soma mask image
+            bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                soma mask image
 
-        px_sz: numpy.ndarray (shape=(3,), dtype=float)
-            output pixel size [μm]
+            px_sz: numpy.ndarray (shape=(3,), dtype=float)
+                output pixel size [μm]
     """
-    # skip Frangi filter stage if orientation vectors were directly provided as input
+    # skip the Frangi filter stage if an orientation vector field was provided as input
     if in_img['is_vec']:
-        out_img = dict()
-        out_img['fbr_vec'] = in_img['data']
-        out_img['iso_fbr'] = None
-        out_img['px_sz'] = None
+        out_img = {'vec': in_img['data'], 'iso': None, 'px_sz': None}
 
     else:
-
-        # get resources configuration
+        # get processing configuration
         ram, jobs = get_resource_config(cli_args)
+        frangi_cfg = get_frangi_config(cli_args, in_img)
+        batch_sz, slc_shp, rsz_ratio, ovlp, tot_slc = get_slice_config(in_img, frangi_cfg, ram=ram, jobs=jobs)
 
-        # get Frangi filter configuration
-        frangi_cfg, px_sz_iso = get_frangi_config(cli_args, in_img)
-
-        # configure the batch of basic image slices analyzed using parallel threads
-        batch_sz, slc_shp, slc_shp_um, rsz_ratio, ovlp, ovlp_rsz = \
-            config_frangi_batch(in_img, frangi_cfg, px_sz_iso, ram=ram, jobs=jobs)
-
-        # generate image range lists
-        slc_rng, out_slc_shp, tot_slc, batch_sz = \
-            generate_slice_lists(slc_shp, in_img['shape'], batch_sz, rsz_ratio, ovlp, in_img['msk_bc'])
-
-        # initialize output arrays
-        out_img, z_sel = init_frangi_arrays(frangi_cfg, in_img['shape'], out_slc_shp, rsz_ratio, save_dirs['tmp'],
-                                            msk_bc=in_img['msk_bc'])
-
-        # print Frangi filter configuration
-        print_frangi_info(in_img, frangi_cfg, slc_shp_um, tot_slc)
-
-        # parallel Frangi filter-based fiber orientation analysis of microscopy image slices
-        print_flsh(f"[Parallel(n_jobs={batch_sz})]: Using backend ThreadingBackend with {batch_sz} concurrent workers.")
+        # conduct a Frangi-filter-based analysis of fiber orientations using concurrent workers
+        slc_gen = slice_generator(in_img, slc_shp, rsz_ratio, ovlp, in_img['msk_bc'])
+        out_img = init_frangi_arrays(in_img, cli_args, frangi_cfg, slc_shp, rsz_ratio)
         t_start = perf_counter()
-        with Parallel(n_jobs=batch_sz, prefer='threads', require='sharedmem') as parallel:
-            parallel(
-                delayed(frangi_analysis)(s, in_img, out_img, frangi_cfg, ovlp_rsz, rsz_ratio, z_sel, **slc_rng,
-                                         print_info=(t_start, batch_sz, tot_slc))
-                for s in range(tot_slc))
+        print_frangi_config(in_img, frangi_cfg, slc_shp, batch_sz)
+        with Parallel(n_jobs=batch_sz, prefer='threads', return_as="generator_unordered") as parallel:
+            out_gen = parallel(delayed(frangi_analysis)(s, frangi_cfg, ovlp, rsz_ratio, (t_start, batch_sz, tot_slc))
+                               for s in slc_gen)
 
-        # save output arrays
-        save_frangi_arrays(save_dirs['frangi'], in_img['name'], out_img, px_sz_iso, ram=ram)
+            # write results to output arrays once available
+            for out_slc in out_gen:
+                write_frangi_arrays(out_img, out_slc)
 
-        # add isotropic pixel size
-        out_img['px_sz'] = px_sz_iso
+        # save output arrays to file
+        save_frangi_arrays(save_dirs['frangi'], in_img['name'], out_img, ram=ram)
 
     return out_img
 
 
-def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_rng, out_rng, in_pad,
-                    print_info=None, pad_mode='reflect', bc_thr='yen', ts_thr=0.0001, verbose=10):
+def frangi_analysis(img_slc_data, cfg, ovlp, rsz_ratio, print_info):
     """
     Conduct a Frangi-based fiber orientation analysis on basic slices selected from the whole microscopy volume image.
 
     Parameters
     ----------
-    s: int
-        image slice index
+    img_slc_data: tuple
+        input image slice data
 
-    in_img: dict
-        input image dictionary (extended)
-
-            data: numpy.ndarray or NumPy memory-map object (axis order=(Z,Y,X) or (Z,Y,X,C) or (Z,C,Y,X))
-                3D microscopy image
+            fbr_slc: numpy.ndarray (axis order=(Z,Y,X))
+                fiber image slice
 
             ts_msk: numpy.ndarray (dtype=bool)
-                tissue reconstruction binary mask
-
-            ch_ax: int
-                RGB image channel axis (either 1, 3, or None for grayscale images)
-
-            fb_ch: int
-                neuronal fibers channel
-
-            bc_ch: int
-                brain cell soma channel
-
-            msk_bc: bool
-                if True, mask neuronal bodies within
-                the optionally provided channel
-
-            psf_fwhm: numpy.ndarray (shape=(3,), dtype=float)
-                3D FWHM of the PSF [μm]
-
-            px_sz: numpy.ndarray (shape=(3,), dtype=float)
-                pixel size [μm]
+                tissue mask slice
 
-            name: str
-                name of the 3D microscopy image
-
-            is_vec: bool
-                vector field flag
-
-            shape: numpy.ndarray (shape=(3,), dtype=int)
-                total image shape
+            bc_slc: numpy.ndarray (axis order=(Z,Y,X))
+                brain cell image slice
 
-            shape_um: numpy.ndarray (shape=(3,), dtype=float)
-                total image shape [μm]
+            pad_rng: numpy.ndarray (axis order=(Z,Y,X))
+                image padding range
 
-            item_sz: int
-                image item size [B]
-
-    out_img: dict
-        output image dictionary
-
-            vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
-                fiber orientation vector field
-
-            clr: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=uint8)
-                orientation colormap image
-
-            fa: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-                fractional anisotropy image
-
-            frangi: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-                Frangi-enhanced image (fiber probability image)
-
-            iso: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-                isotropic fiber image
-
-            fbr_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-                fiber mask image
-
-            bc_msk: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
-                soma mask image
-
-            px_sz: numpy.ndarray (shape=(3,), dtype=float)
-                output pixel size [μm]
+            out_rng: NumPy slice object
+                output range
 
     cfg: dict
         Frangi filter configuration
@@ -254,9 +178,6 @@ def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_
             px_sz: numpy.ndarray (shape=(3,), dtype=float)
                 pixel size [μm]
 
-            z_rng: int
-                output z-range in [px]
-
             bc_ch: int
                 neuronal bodies channel
 
@@ -279,72 +200,60 @@ def frangi_analysis(s, in_img, out_img, cfg, ovlp, rsz_ratio, z_sel, in_rng, bc_
     rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
         3D image resize ratio
 
-    z_sel: NumPy slice object
-        selected z-depth range
-
-    in_rng: NumPy slice object
-        input image range (fibers)
+    Returns
+    -------
+    out_slc: dict
+        slice output dictionary
 
-    bc_rng: NumPy slice object
-        input image range (neurons)
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+                3D fiber orientation field
 
-    out_rng: NumPy slice object
-        output range
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap
 
-    in_pad: numpy.ndarray (axis order=(Z,Y,X))
-        image padding range
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy
 
-    print_info: tuple
-        optional printed information
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image slice (fiber probability image)
 
-    pad_mode: str
-        image padding mode
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image slice
 
-    bc_thr: str
-        thresholding method applied to the neuronal bodies channel
+            ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                tissue mask slice
 
-    ts_thr: float
-        relative threshold on non-zero tissue pixels
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask slice
 
-    verbose: int
-        verbosity level (print progress every N=verbose image slices)
+            bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool)
+                brain cell mask slice
 
-    Returns
-    -------
-    None
+            rng: NumPy slice object
+                output range
     """
-    # select fiber image slice
-    fbr_slc, ts_msk_slc = slice_image(in_img['data'], in_rng[s], in_img['ch_ax'], ch=in_img['fb_ch'],
-                                      ts_msk=in_img['ts_msk'])
+    # unpack image slice variables
+    fbr_slc, ts_msk_slc, bc_slc, pad_rng, out_rng = img_slc_data
 
     # process non-background image slice
-    is_valid = check_background(fbr_slc, ts_msk=ts_msk_slc, ts_thr=ts_thr)
+    is_valid = check_background(fbr_slc, ts_msk=ts_msk_slc)
     if is_valid:
-        orient_slc, frangi_slc, iso_slc, fbr_msk_slc = \
-            analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng[s], in_pad[s], pad_mode=pad_mode, ts_msk=ts_msk_slc)
-
-        # (optional) neuronal body masking
-        if in_img['msk_bc']:
-
-            # get soma image slice
-            bc_slc, _ = slice_image(in_img['data'], bc_rng[s], in_img['ch_ax'], ch=in_img['bc_ch'])
+        out_slc = analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad_rng, ts_msk=ts_msk_slc)
 
-            # suppress soma contribution
-            orient_slc, bc_msk_slc = reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng[s], bc_thr)
+        # (optional) neuronal cell body masking
+        if bc_slc is not None:
+            out_slc = reject_brain_cells(bc_slc, out_slc, rsz_ratio, out_rng)
 
-        # soma mask not available
-        else:
-            bc_msk_slc = None
-
-        # fill memory-mapped output arrays
-        write_frangi_arrays(out_rng[s], z_sel, iso_slc, frangi_slc, fbr_msk_slc, bc_msk_slc, **orient_slc, **out_img)
+        out_slc['rng'] = out_rng
+    else:
+        out_slc = None
+        
+    print_frangi_progress(print_info, is_valid)
 
-    # print progress
-    if print_info is not None:
-        print_frangi_progress(print_info, is_valid, verbose=verbose)
+    return out_slc
 
 
-def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflect', ts_msk=None):
+def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad_rng, ts_msk=None):
     """
     Analyze 3D fiber orientations exploiting a Frangi-filter-based
     unsupervised enhancement of tubular structures.
@@ -378,9 +287,6 @@ def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflec
             px_sz: numpy.ndarray (shape=(3,), dtype=float)
                 pixel size [μm]
 
-            z_rng: int
-                output z-range in [px]
-
             bc_ch: int
                 neuronal bodies channel
 
@@ -406,67 +312,67 @@ def analyze_fibers(fbr_slc, cfg, ovlp, rsz_ratio, out_rng, pad, pad_mode='reflec
     out_rng: NumPy slice object
         output range
 
-    pad: numpy.ndarray (axis order=(Z,Y,X))
+    pad_rng: numpy.ndarray (axis order=(Z,Y,X))
         image padding range
 
-    pad_mode: str
-        image padding mode
-
     ts_msk: numpy.ndarray (dtype=bool)
         tissue reconstruction binary mask
 
     Returns
     -------
-    orient_slc: dict
-        slice orientation dictionary
+    out_slc: dict
+        slice output dictionary
+
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+                3D fiber orientation field
 
-        vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
-            3D fiber orientation field
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap
 
-        clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
-            orientation colormap image
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy
 
-        fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-            fractional anisotropy image
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image slice (fiber probability image)
 
-    frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        Frangi-enhanced image slice (fiber probability image)
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image slice
 
-    iso_fbr_img: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        isotropic fiber image slice
+            ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                tissue mask slice
 
-    fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-        fiber mask image slice
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask slice
     """
     # preprocess fiber slice (lateral blurring and downsampling)
     iso_fbr_slc, pad_rsz, ts_msk_slc_rsz = \
-        correct_anisotropy(fbr_slc, rsz_ratio, sigma=cfg['smooth_sd'], pad=pad, ts_msk=ts_msk)
+        correct_anisotropy(fbr_slc, rsz_ratio, sigma=cfg['smooth_sd'], pad=pad_rng, ts_msk=ts_msk)
 
     # pad fiber slice if required
     if pad_rsz is not None:
-        iso_fbr_slc = np.pad(iso_fbr_slc, pad_rsz, mode=pad_mode)
+        iso_fbr_slc = np.pad(iso_fbr_slc, pad_rsz, mode='reflect')
 
         # pad tissue mask if available
         if ts_msk_slc_rsz is not None:
             ts_msk_slc_rsz = np.pad(ts_msk_slc_rsz, pad_rsz, mode='constant')
 
     # 3D Frangi filter
-    orient_slc, frangi_slc = \
-        frangi_filter(iso_fbr_slc, scales_px=cfg['scales_px'],
-                      alpha=cfg['alpha'], beta=cfg['beta'], gamma=cfg['gamma'], hsv=cfg['hsv_cmap'])
+    out_slc = frangi_filter(iso_fbr_slc, scales_px=cfg['scales_px'],
+                            alpha=cfg['alpha'], beta=cfg['beta'], gamma=cfg['gamma'], hsv=cfg['hsv_cmap'])
 
     # crop resulting slices
-    orient_slc, frangi_slc, iso_fbr_slc, ts_msk_slc_rsz = \
-        crop_lst([orient_slc, frangi_slc, iso_fbr_slc, ts_msk_slc_rsz], out_rng, ovlp)
+    out_slc.update({'iso': iso_fbr_slc, 'ts_msk': ts_msk_slc_rsz})
+    ovlp_rsz = np.multiply(ovlp * np.ones((3,)), rsz_ratio).astype(int)
+    out_slc = crop_img_dict(out_slc, out_rng, ovlp_rsz)
 
     # remove Frangi filter background
-    orient_slc, fbr_msk_slc = \
-        mask_background(frangi_slc, orient_slc, ts_msk=ts_msk_slc_rsz, method=cfg['fb_thr'], invert=False)
+    out_slc, msk = mask_background(out_slc, method=cfg['fb_thr'], invert=False)
+    out_slc['fbr_msk'] = msk
 
-    return orient_slc, frangi_slc, iso_fbr_slc, fbr_msk_slc
+    return out_slc
 
 
-def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'):
+def reject_brain_cells(bc_slc, out_slc, rsz_ratio, out_rng, bc_thr='yen'):
     """
     Suppress soma contribution using the optional image channel
     of neuronal bodies.
@@ -476,17 +382,29 @@ def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'):
     bc_slc: numpy.ndarray (axis order=(Z,Y,X))
         brain cell soma image slice
 
-    orient_slc: dict
-        slice orientation dictionary
+    out_slc: dict
+        slice output dictionary
 
-            vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
                 3D fiber orientation field
 
-            clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
-                orientation colormap image
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap
 
-            fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-                fractional anisotropy image
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy
+
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image slice (fiber probability image)
+
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image slice
+
+            ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                tissue mask slice
+
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask slice
 
     rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
         3D image resize ratio
@@ -495,24 +413,36 @@ def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'):
         output range
 
     bc_thr: str
-        thresholding method applied to the neuronal bodies channel
+        thresholding method applied to the channel of brain cell bodies
 
     Returns
     -------
-    orient_slc: dict
-        (masked) slice orientation dictionary
+    out_slc: dict
+        (masked) slice output dictionary
+
+            vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+                3D fiber orientation field
+
+            clr: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
+                orientation colormap
+
+            fa: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fractional anisotropy
+
+            frangi: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                Frangi-enhanced image slice (fiber probability image)
 
-        vec_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
-            (masked) 3D fiber orientation field
+            iso: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                isotropic fiber image slice
 
-        clr_slc: numpy.ndarray (axis order=(Z,Y,X,C), dtype=uint8)
-            (masked) orientation colormap image
+            ts_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                tissue mask slice
 
-        fa_slc: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
-            (masked) fractional anisotropy image
+            fbr_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=uint8)
+                fiber mask slice
 
-    bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool)
-        brain cell mask
+            bc_msk: numpy.ndarray (axis order=(Z,Y,X), dtype=bool)
+                brain cell mask
     """
     # resize soma slice (lateral downsampling)
     iso_bc, _, _ = correct_anisotropy(bc_slc, rsz_ratio)
@@ -521,12 +451,13 @@ def reject_brain_cells(bc_slc, orient_slc, rsz_ratio, out_rng, bc_thr='yen'):
     iso_bc = crop(iso_bc, out_rng)
 
     # mask neuronal bodies
-    orient_slc, bc_msk = mask_background(iso_bc, orient_slc, method=bc_thr, invert=True)
+    out_slc, msk = mask_background(out_slc, ref_img=iso_bc, method=bc_thr, invert=True)
+    out_slc['bc_msk'] = msk
 
-    return orient_slc, bc_msk
+    return out_slc
 
 
-def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_name, backend='loky', verbose=100):
+def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_name):
     """
     Iterate over the required spatial scales and apply the parallel ODF analysis
     implemented in parallel_odf_on_slices().
@@ -552,12 +483,6 @@ def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_n
     img_name: str
         name of the input volume image
 
-    backend: str
-        backend module employed by joblib.Parallel
-
-    verbose: int
-        joblib verbosity level
-
     Returns
     -------
     None
@@ -565,26 +490,21 @@ def parallel_odf_over_scales(cli_args, save_dirs, fbr_vec, iso_fbr, px_sz, img_n
     # print ODF analysis heading
     print_odf_info(cli_args.odf_res, cli_args.odf_deg)
 
-    # export all images (optional)
-    exp_all = cli_args.exp_all
-
     # compute spherical harmonics normalization factors (once for all scales)
     odf_norm = get_sph_harm_norm_factors(cli_args.odf_deg)
 
-    # get number of logical cores
-    num_cpu = get_available_cores()
-
     # get pixel size from CLI arguments if not provided
     if px_sz is None:
-        px_sz = np.array([cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy])
+        px_sz = (cli_args.px_size_z, cli_args.px_size_xy, cli_args.px_size_xy)
 
     # parallel ODF analysis of fiber orientation vectors
     # over spatial scales of interest
     n_scales = len(cli_args.odf_res)
+    num_cpu = get_available_cores()
     batch_sz = min(num_cpu, n_scales)
-    with Parallel(n_jobs=batch_sz, backend=backend, verbose=verbose) as parallel:
-        parallel(delayed(odf_analysis)(fbr_vec, iso_fbr, px_sz, save_dirs, img_name,
-                                       odf_norm=odf_norm, odf_deg=cli_args.odf_deg, odf_scale_um=s, exp_all=exp_all)
+    with Parallel(n_jobs=batch_sz, verbose=10, prefer='threads') as parallel:
+        parallel(delayed(odf_analysis)(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_norm=odf_norm,
+                                       odf_deg=cli_args.odf_deg, odf_scale_um=s, exp_all=cli_args.exp_all)
                  for s in cli_args.odf_res)
 
     # print output directory
@@ -634,9 +554,9 @@ def odf_analysis(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_scale_um, odf
 
     # initialize ODF analysis output arrays
     odf, odi, dnst, bg_mrtrix, vec_tensor_eigen = \
-        init_odf_arrays(fbr_vec.shape[:-1], save_dirs['tmp'], odf_scale=odf_scale, odf_deg=odf_deg, exp_all=exp_all)
+        init_odf_arrays(fbr_vec.shape[:-1], odf_scale=odf_scale, odf_deg=odf_deg, exp_all=exp_all)
 
-    # generate downsampled background for MRtrix3 mrview
+    # generate down-sampled background for visualization in MRtrix3
     generate_odf_background(bg_mrtrix, fbr_vec, vxl_side=odf_scale, iso_fbr=iso_fbr)
 
     # compute ODF coefficients
diff --git a/foa3d/preprocessing.py b/foa3d/preprocessing.py
index 4f65851..d6148b7 100644
--- a/foa3d/preprocessing.py
+++ b/foa3d/preprocessing.py
@@ -35,7 +35,7 @@ def config_anisotropy_correction(px_sz, psf_fwhm):
         new isotropic pixel size [μm]
     """
     # set the isotropic pixel resolution equal to the z-sampling step
-    px_sz_iso = np.max(px_sz) * np.ones(shape=(3,))
+    px_sz_iso = tuple(np.max(px_sz) * np.ones(shape=(3,)))
 
     # detect preprocessing requirement
     cndt_1 = not np.all(psf_fwhm == psf_fwhm[0])
diff --git a/foa3d/printing.py b/foa3d/printing.py
index 5777e41..72290ce 100644
--- a/foa3d/printing.py
+++ b/foa3d/printing.py
@@ -76,7 +76,7 @@ def print_blur(sigma_um, psf_fwhm):
                f"Adjusted PSF FWHM    [μm]: ({psf_sz:.3f}, {psf_sz:.3f}, {psf_sz:.3f})")
 
 
-def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num):
+def print_frangi_config(in_img, cfg, slc_shp, batch_sz):
     """
     Print Frangi filter stage heading.
 
@@ -125,7 +125,7 @@ def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num):
             item_sz: int
                 image item size [B]
 
-    frangi_cfg: dict
+    cfg: dict
         Frangi filter configuration
 
             alpha: float
@@ -168,37 +168,36 @@ def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num):
             exp_all: bool
                 export all images
 
-    in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=float)
-        shape of the analyzed image slices [μm]
+    slc_shp: numpy.ndarray (shape=(3,), dtype=int)
+        shape of the analyzed image slices [px]
 
-    tot_slc_num: int
-        total number of analyzed image slices
+    batch_sz: int
+        slice batch size
 
     Returns
     -------
     None
     """
     # print Frangi filter sensitivity and scales
-    alpha = frangi_cfg['alpha']
-    beta = frangi_cfg['beta']
-    gamma = frangi_cfg['gamma']
-
-    scales_um = np.asarray(frangi_cfg['scales_um'])
+    alpha = cfg['alpha']
+    beta = cfg['beta']
+    gamma = cfg['gamma']
+    scales_um = np.asarray(cfg['scales_um'])
     if gamma is None:
         gamma = 'auto'
 
-    print_flsh(color_text(0, 191, 255, "\n\n\n3D Frangi Filter\n") + "\nSensitivity\n" +
+    print_flsh(color_text(0, 191, 255, "\n3D Frangi Filter\n") + "\nSensitivity\n" +
                f"• plate-like \u03B1: {alpha:.1e}\n• blob-like  \u03B2: {beta:.1e}\n• background \u03B3: {gamma}\n\n" +
                f"Enhanced scales      [μm]: {scales_um}\nEnhanced diameters   [μm]: {4 * scales_um}\n")
 
-    # print iterative analysis information
-    print_slicing_info(in_img['shape_um'], in_slc_shp_um, tot_slc_num, in_img['px_sz'], in_img['item_sz'])
-
-    # print neuron masking info
+    # print parallel processing information
+    slc_shp_um = np.multiply(cfg['px_sz'], slc_shp)
+    print_slicing_info(in_img['shape_um'], slc_shp_um, in_img['px_sz'], in_img['item_sz'])
     print_soma_masking(in_img['msk_bc'])
+    print_flsh(f"[Parallel(n_jobs={batch_sz})]: Using backend ThreadingBackend with {batch_sz} concurrent workers.")
 
 
-def print_frangi_progress(info, is_valid, verbose=1):
+def print_frangi_progress(info, is_valid, verbose=5):
     """
     Print Frangi filter progress.
 
@@ -338,7 +337,7 @@ def print_prepro_heading():
                "\n                              Z      Y      X")
 
 
-def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz):
+def print_slicing_info(img_shp_um, slc_shp_um, px_sz, img_item_sz):
     """
     Print information on the slicing of the basic image sub-volumes processed by the Foa3D tool.
 
@@ -378,8 +377,7 @@ def print_slicing_info(img_shp_um, slc_shp_um, tot_slc_num, px_sz, img_item_sz):
     print_flsh(f"Total image shape    [μm]: ({img_shp_um[0]:.1f}, {img_shp_um[1]:.1f}, {img_shp_um[2]:.1f})\n" +
                f"Total image size     [MB]: {np.ceil(img_sz / 1024**2).astype(int)}\n\n" +
                f"Image slice shape    [μm]: ({slc_shp_um[0]:.1f}, {slc_shp_um[1]:.1f}, {slc_shp_um[2]:.1f})\n" +
-               f"Image slice size     [MB]: {np.ceil(max_slc_sz / 1024**2).astype(int)}\n" +
-               f"Image slice number:        {tot_slc_num}\n")
+               f"Image slice size     [MB]: {np.ceil(max_slc_sz / 1024**2).astype(int)}\n")
 
 
 def print_soma_masking(msk_bc):
diff --git a/foa3d/slicing.py b/foa3d/slicing.py
index 32f5516..246ebda 100644
--- a/foa3d/slicing.py
+++ b/foa3d/slicing.py
@@ -4,7 +4,6 @@
 import psutil
 
 from numba import njit
-from foa3d.printing import print_flsh
 from foa3d.utils import get_available_cores
 
 
@@ -187,7 +186,7 @@ def compute_slice_range(ax_iter, slc_shp, img_shp, slc_per_dim, ovlp=0, flip=Fal
     return rng, pad
 
 
-def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2):
+def compute_overlap(smooth_sigma, frangi_sigma_um, rsz_ratio=None, trunc=2):
     """
     Compute lateral slice extension range
     for coping with smoothing-related boundary artifacts.
@@ -200,7 +199,7 @@ def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2):
     frangi_sigma_um: numpy.ndarray (dtype=float)
         analyzed spatial scales [μm]
 
-    px_rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
+    rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
         3D image resize ratio
 
     trunc: int
@@ -211,22 +210,17 @@ def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2):
     -------
     ovlp: int
         overlapping range between image slices along each axis [px]
-
-    ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int)
-        resized overlapping range between image slices
-        (if px_rsz_ratio is not None) [px]
     """
-    sigma = np.max(frangi_sigma_um) / px_rsz_ratio
+    sigma = np.max(frangi_sigma_um) / rsz_ratio
     if smooth_sigma is not None:
         sigma = np.concatenate((smooth_sigma, sigma))
 
     ovlp = int(np.ceil(2 * trunc * np.max(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
 
-    return ovlp, ovlp_rsz
+    return ovlp
 
 
-def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0):
+def compute_slice_shape(img_shp, item_sz, slc_sz=1e6, ovlp=0):
     """
     Compute basic image chunk shape depending on its maximum size (in bytes).
 
@@ -238,9 +232,6 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0):
     item_sz: int
         image item size [B]
 
-    px_sz: numpy.ndarray (shape=(3,), dtype=float)
-        pixel size [μm]
-
     slc_sz: float
         maximum memory size of the basic image slices
         analyzed using parallel threads [B]
@@ -253,11 +244,6 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0):
     slc_shp: numpy.ndarray (shape=(3,), dtype=int)
         shape of the basic image slices
         analyzed using parallel threads [px]
-
-    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)
     """
     if len(img_shp) == 4:
         item_sz *= 3
@@ -267,12 +253,11 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0):
     side_z = min(img_shp[0] + tot_ovlp, side_xyz)
     side_xy = np.floor(np.sqrt(np.abs(slc_sz / side_z / item_sz))).astype(int)
     slc_shp = np.array([side_z, side_xy, side_xy]) - tot_ovlp
-    slc_shp_um = np.multiply(slc_shp, px_sz) if px_sz is not None else None
 
-    return slc_shp, slc_shp_um
+    return slc_shp
 
 
-def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1):
+def compute_slice_size(max_ram, mem_growth, batch_sz):
     """
     Compute the size of the basic microscopy image slices fed to the Frangi filtering stage.
 
@@ -284,27 +269,21 @@ def compute_slice_size(max_ram, mem_growth, mem_fudge, batch_sz, ns=1):
     mem_growth: float
         empirical memory growth factor
 
-    mem_fudge: float
-        memory fudge factor
-
     batch_sz: int
         slice batch 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)
+    slc_sz = max_ram / (batch_sz * mem_growth)
 
     return slc_sz
 
 
-def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fudge=1.0, jobs=None, ram=None, shp_thr=7):
+def get_slice_config(in_img, frangi_cfg, mem_growth=149.7, jobs=None, ram=None, shp_thr=7):
     """
     Compute size and number of the batches of basic microscopy image slices
     analyzed in parallel.
@@ -312,7 +291,6 @@ def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fud
     Parameters
     ----------
 
-
     px_sz_iso: int
         isotropic pixel size [μm]
 
@@ -342,18 +320,14 @@ def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fud
         shape of the basic image slices
         analyzed using parallel threads [px]
 
-    in_slc_shp_um: numpy.ndarray (shape=(3,), dtype=int)
-        shape of the basic image slices
-        analyzed using parallel threads [μm]
-
     rsz_ratio: numpy.ndarray (shape=(3,), dtype=float)
         3D image resize ratio
 
     ovlp: int
         overlapping range between image slices along each axis [px]
 
-    ovlp_rsz: numpy.ndarray (shape=(3,), dtype=int)
-        resized overlapping range between image slices [px]
+    tot_slc: int
+        total number of image slices
     """
     # maximum RAM not provided: use all
     if ram is None:
@@ -362,22 +336,14 @@ def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fud
     # number of logical cores
     num_cpu = get_available_cores()
     if jobs is None:
-        jobs = num_cpu
-
-    # number of spatial scales
-    ns = len(frangi_cfg['scales_um'])
-
-    # initialize slice batch size
-    batch_sz = np.min([jobs // ns, num_cpu]).astype(int) + 1
-
-    # get pixel resize ratio
-    rsz_ratio = np.divide(in_img['px_sz'], px_sz_iso)
+        jobs = num_cpu   
 
     # compute slice overlap (boundary artifacts suppression)
-    ovlp, ovlp_rsz = compute_overlap(frangi_cfg['smooth_sd'], frangi_cfg['scales_um'],
-                                     px_rsz_ratio=rsz_ratio)
+    rsz_ratio = np.divide(in_img['px_sz'], frangi_cfg['px_sz'])
+    ovlp = compute_overlap(frangi_cfg['smooth_sd'], frangi_cfg['scales_um'], rsz_ratio=rsz_ratio)
 
     # compute the shape of basic microscopy image slices
+    batch_sz = np.min([jobs, num_cpu]).astype(int) + 1
     in_slc_shp = np.array([-1])
     while np.any(in_slc_shp < shp_thr):
         batch_sz -= 1
@@ -385,11 +351,14 @@ def config_frangi_batch(in_img, frangi_cfg, px_sz_iso, mem_growth=149.7, mem_fud
             raise ValueError(
                 "Basic image slices do not fit the available RAM: decrease spatial scales and/or maximum scale value!")
         else:
-            slc_sz = compute_slice_size(ram, mem_growth, mem_fudge, batch_sz, ns)
-            in_slc_shp, in_slc_shp_um = compute_slice_shape(in_img['shape'], in_img['item_sz'],
-                                                            px_sz=in_img['px_sz'], slc_sz=slc_sz, ovlp=ovlp)
+            slc_sz = compute_slice_size(ram, mem_growth, batch_sz)
+            in_slc_shp = compute_slice_shape(in_img['shape'], in_img['item_sz'], slc_sz=slc_sz, ovlp=ovlp)
+
+    # adjust batch size
+    tot_slc = np.prod(np.floor(np.divide(in_img['shape'], in_slc_shp)).astype(int))
+    batch_sz = min(batch_sz, tot_slc)
 
-    return batch_sz, in_slc_shp, in_slc_shp_um, rsz_ratio, ovlp, ovlp_rsz
+    return batch_sz, in_slc_shp, rsz_ratio, ovlp, tot_slc
 
 
 def crop(img, rng, ovlp=None, flip=()):
@@ -432,14 +401,14 @@ def crop(img, rng, ovlp=None, flip=()):
     return img_out
 
 
-def crop_lst(img_lst, rng, ovlp=None, flip=()):
+def crop_img_dict(img_dict, rng, ovlp=None, flip=()):
     """
     Shrink list of image slices at total volume boundaries, for overall shape consistency.
 
     Parameters
     ----------
-    img_lst: list
-        list of images to be cropped
+    img_dict: list
+        dictionary of image data to be cropped
 
     rng: tuple
         3D slice index range
@@ -452,21 +421,17 @@ def crop_lst(img_lst, rng, ovlp=None, flip=()):
 
     Returns
     -------
-    img_lst: list
-        list of cropped image slices
+    img_dict: list
+        dictionary of cropped image slices
     """
-    for s, img in enumerate(img_lst):
-        if img is not None:
-            if isinstance(img, dict):
-                for key in img.keys():
-                    img_lst[s][key] = crop(img[key], rng, ovlp=ovlp, flip=flip)
-            else:
-                img_lst[s] = crop(img, rng, ovlp=ovlp, flip=flip)
+    for key in img_dict.keys():
+        if img_dict[key] is not None:
+            img_dict[key] = crop(img_dict[key], rng, ovlp=ovlp, flip=flip)
 
-    return img_lst
+    return img_dict
 
 
-def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, msk_bc=False):
+def slice_generator(in_img, in_slc_shp, px_rsz_ratio, ovlp=0, msk_bc=False):
     """
     Generate image slice ranges for the Frangi filtering stage.
 
@@ -518,48 +483,32 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms
     """
     # 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)
+    out_img_shp = np.ceil(np.multiply(px_rsz_ratio, in_img['shape'])).astype(int)
 
     # number of image slices along each axis
-    slc_per_dim = np.floor(np.divide(img_shp, in_slc_shp)).astype(int)
-
-    # initialize empty range lists
-    # fill slice range dictionary
-    slc_rng = dict()
-    slc_rng['in_rng'] = list()
-    slc_rng['in_pad'] = list()
-    slc_rng['out_rng'] = list()
-    slc_rng['bc_rng'] = list()
-    tot_slc = np.prod(slc_per_dim)
-    for i, zyx in enumerate(product(range(slc_per_dim[0]), range(slc_per_dim[1]), range(slc_per_dim[2]))):
+    slc_per_dim = np.floor(np.divide(in_img['shape'], in_slc_shp)).astype(int)
+
+    # generate image slices
+    for zyx in product(range(slc_per_dim[0]), range(slc_per_dim[1]), range(slc_per_dim[2])):
 
         # index ranges of analyzed image slices (with padding)
-        in_rng, pad = compute_slice_range(zyx, in_slc_shp, img_shp, slc_per_dim, ovlp=ovlp)
-        slc_rng['in_rng'].append(in_rng)
-        slc_rng['in_pad'].append(pad)
+        in_rng, pad = compute_slice_range(zyx, in_slc_shp, in_img['shape'], slc_per_dim, ovlp=ovlp)
 
         # output index ranges
         out_rng, _ = compute_slice_range(zyx, out_slc_shp, out_img_shp, slc_per_dim)
-        slc_rng['out_rng'].append(out_rng)
 
         # (optional) neuronal body channel
         if msk_bc:
-            bc_rng, _ = compute_slice_range(zyx, in_slc_shp, img_shp, slc_per_dim)
-            slc_rng['bc_rng'].append(bc_rng)
+            bc_rng, _ = compute_slice_range(zyx, in_slc_shp, in_img['shape'], slc_per_dim)
+            bc_slc, _ = slice_image(in_img['data'], bc_rng, in_img['ch_ax'], ch=in_img['bc_ch'])
         else:
-            slc_rng['bc_rng'].append(None)
-
-        # show progress
-        print_flsh(f"Generating image slice ranges:\t{100 * (i+1) / tot_slc:.1f}%", end='\r')
-
-    # total number of slices
-    tot_slc = len(slc_rng['in_rng'])
+            bc_slc = None
 
-    # adjust slice batch size
-    if batch_sz > tot_slc:
-        batch_sz = tot_slc
+        # slice input image reconstruction here
+        fbr_slc, ts_msk_slc = slice_image(in_img['data'], in_rng, in_img['ch_ax'], ch=in_img['fb_ch'],
+                                          ts_msk=in_img['ts_msk'])
 
-    return slc_rng, out_slc_shp, tot_slc, batch_sz
+        yield fbr_slc, ts_msk_slc, bc_slc, pad, out_rng
 
 
 def slice_image(img, rng, ch_ax, ch, ts_msk=None):
diff --git a/foa3d/utils.py b/foa3d/utils.py
index 0695187..296e7f5 100644
--- a/foa3d/utils.py
+++ b/foa3d/utils.py
@@ -1,8 +1,5 @@
-import gc
-import tempfile
 from multiprocessing import cpu_count
-from os import environ, path, unlink
-from shutil import rmtree
+from os import environ
 from time import perf_counter
 
 import numpy as np
@@ -57,55 +54,6 @@ def create_background_mask(img, method='yen', black_bg=False):
     return bg_msk
 
 
-def create_memory_map(shape, dtype, name='tmp', tmp_dir=None, arr=None, mmap_mode='w+'):
-    """
-    Create a memory-map to an array stored in a binary file on disk.
-
-    Parameters
-    ----------
-    shape: tuple
-        shape of the stored array
-
-    dtype:
-        data-type used to interpret the file contents
-
-    name: str
-        optional temporary filename
-
-    tmp_dir: str
-        temporary file directory
-
-    arr: numpy.ndarray
-        array to be mapped
-
-    mmap_mode: str
-        file opening mode
-
-    Returns
-    -------
-    mmap: NumPy memory map
-        memory-mapped array
-    """
-    if tmp_dir is None:
-        tmp_dir = tempfile.mkdtemp()
-    mmap_path = path.join(tmp_dir, name + '.mmap')
-
-    if path.exists(mmap_path):
-        unlink(mmap_path)   
-
-    if arr is None:
-        _ = open(mmap_path, mode='w+')
-        mmap = np.memmap(mmap_path, dtype=dtype, mode=mmap_mode, shape=shape)
-        mmap[:] = 0
-    else:
-        _ = dump(arr, mmap_path)
-        mmap = load(mmap_path, mmap_mode=mmap_mode)
-        del arr
-        _ = gc.collect()
-
-    return mmap
-
-
 def detect_ch_axis(img):
     """
     Detect image channel axis.
@@ -182,33 +130,6 @@ def get_item_size(dtype):
     return item_sz
 
 
-def delete_tmp_folder(tmp_dir, tmp_data):
-    """
-    Delete temporary folder.
-
-    Parameters
-    ----------
-    tmp_dir: str
-        path to temporary folder to be removed
-    
-    tmp_data: tuple
-        temporary data dictionaries
-
-    Returns
-    -------
-    None
-    """
-    try:
-        for data in tmp_data:
-            for k in list(data.keys()):
-                del data[k]
-
-        rmtree(tmp_dir)
-
-    except OSError:
-        pass
-
-
 def divide_nonzero(nd_array1, nd_array2, new_val=1e-10):
     """
     Divide two arrays handling zero denominator values.
diff --git a/requirements.in b/requirements.in
index f89bf76..fe1f595 100644
--- a/requirements.in
+++ b/requirements.in
@@ -1,5 +1,5 @@
 astropy
-joblib
+joblib==1.4.2
 matplotlib
 nibabel
 numba
diff --git a/requirements.txt b/requirements.txt
index cd43c64..4f1022d 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,68 +1,67 @@
 #
-# This file is autogenerated by pip-compile with Python 3.9
+# This file is autogenerated by pip-compile with Python 3.8
 # by the following command:
 #
 #    pip-compile requirements.in
 #
-about-time==3.1.1
-    # via alive-progress
-astropy==5.1
+astropy==5.2.2
     # via -r requirements.in
-cachetools==5.3.1
+cachetools==5.5.0
     # via zetastitcher
-clarabel==0.6.0
+clarabel==0.9.0
     # via cvxpy
 coloredlogs==15.0.1
     # via zetastitcher
-contourpy==1.0.5
+contourpy==1.1.1
     # via matplotlib
-cvxpy==1.4.1
+cvxpy==1.5.2
     # via zetastitcher
-cycler==0.11.0
+cycler==0.12.1
     # via matplotlib
-daqp==0.5.1
-    # via qpsolvers
-ecos==2.0.12
-    # via
-    #   cvxpy
-    #   qpsolvers
-fonttools==4.37.4
+ecos==2.0.14
+    # via cvxpy
+fonttools==4.55.0
     # via matplotlib
-grapheme==0.6.0
-    # via alive-progress
 humanfriendly==10.0
     # via coloredlogs
-humanize==4.8.0
+humanize==4.10.0
     # via zetastitcher
-imageio==2.22.1
+imageio==2.35.1
     # via
     #   pims
     #   scikit-image
     #   zetastitcher
-joblib==1.2.0
+importlib-metadata==8.5.0
+    # via numba
+importlib-resources==6.4.5
+    # via
+    #   matplotlib
+    #   nibabel
+joblib==1.4.2
     # via -r requirements.in
-kiwisolver==1.4.4
+kiwisolver==1.4.7
     # via matplotlib
-llvmlite==0.39.1
+lazy-loader==0.4
+    # via scikit-image
+llvmlite==0.41.1
     # via numba
-matplotlib==3.6.0
+matplotlib==3.7.5
     # via -r requirements.in
-networkx==2.8.7
+networkx==3.1
     # via
     #   scikit-image
     #   zetastitcher
-nibabel==4.0.2
+nibabel==5.2.1
     # via -r requirements.in
-numba==0.56.2
+numba==0.58.1
     # via -r requirements.in
-numpy==1.23.3
+numpy==1.24.4
     # via
     #   astropy
     #   clarabel
     #   contourpy
     #   cvxpy
     #   ecos
-    #   h5py
     #   imageio
     #   matplotlib
     #   nibabel
@@ -73,91 +72,89 @@ numpy==1.23.3
     #   pims
     #   pyerfa
     #   pywavelets
+    #   qdldl
     #   qpsolvers
     #   scikit-image
     #   scipy
     #   scs
     #   tifffile
     #   zetastitcher
-opencv-python==4.8.1.78
+opencv-python==4.10.0.84
     # via zetastitcher
-osqp==0.6.3
-    # via
-    #   cvxpy
-    #   qpsolvers
-packaging==21.3
+osqp==0.6.7.post3
+    # via cvxpy
+packaging==24.2
     # via
     #   astropy
+    #   lazy-loader
     #   matplotlib
     #   nibabel
+    #   pims
     #   scikit-image
-pandas==2.1.1
+pandas==2.0.3
     # via zetastitcher
-pillow==9.2.0
+pillow==10.4.0
     # via
     #   imageio
     #   matplotlib
     #   scikit-image
-pims==0.6.1
+pims==0.7
     # via zetastitcher
-psutil==5.9.4
+psutil==6.1.0
     # via
     #   -r requirements.in
     #   zetastitcher
-pybind11==2.11.1
-    # via cvxpy
-pyerfa==2.0.0.1
+pyerfa==2.0.0.3
     # via astropy
-pyparsing==3.0.9
-    # via
-    #   matplotlib
-    #   packaging
-pyreadline3==3.4.1
+pyparsing==3.1.4
+    # via matplotlib
+pyreadline3==3.5.4
     # via humanfriendly
-python-dateutil==2.8.2
+python-dateutil==2.9.0.post0
     # via
     #   matplotlib
     #   pandas
-pytz==2023.3.post1
+pytz==2024.2
     # via pandas
 pywavelets==1.4.1
     # via scikit-image
-pyyaml==6.0
+pyyaml==6.0.2
     # via
     #   astropy
     #   zetastitcher
-qdldl==0.1.7.post0
+qdldl==0.1.7.post4
     # via osqp
-qpsolvers==4.0.0
+qpsolvers==4.4.0
     # via zetastitcher
-scikit-image==0.19.3
+scikit-image==0.21.0
     # via -r requirements.in
-scipy==1.9.1
+scipy==1.10.1
     # via
     #   clarabel
     #   cvxpy
     #   ecos
     #   osqp
+    #   qdldl
     #   qpsolvers
     #   scikit-image
     #   scs
     #   zetastitcher
-scs==3.2.3
-    # via
-    #   cvxpy
-    #   qpsolvers
+scs==3.2.7
+    # via cvxpy
 six==1.16.0
     # via python-dateutil
 slicerator==1.1.0
     # via pims
-tifffile==2022.8.12
+tifffile==2023.7.10
     # via
+    #   pims
     #   scikit-image
     #   zetastitcher
-tzdata==2023.3
+tzdata==2024.2
     # via pandas
 zetastitcher==0.7.0.post1
     # via -r requirements.in
-
-# The following packages are considered to be unsafe in a requirements file:
-# setuptools
+zipp==3.20.2
+    # via
+    #   importlib-metadata
+    #   importlib-resources