Minor restorations

foa3d/__main__.py

@@ -1,35 +1,28 @@
-from foa3d.input import get_cli_parser, get_pipeline_config, load_microscopy_image
-from foa3d.pipeline import (parallel_odf_on_scales, parallel_frangi_on_slices)
+from foa3d.input import get_cli_parser, get_resource_config, load_microscopy_image
+from foa3d.pipeline import (parallel_odf_at_scales, parallel_frangi_on_slices)
 from foa3d.printing import print_pipeline_heading
 from foa3d.utils import delete_tmp_folder
 
 
 def foa3d(cli_args):
 
     # load microscopy volume image or array of fiber orientation vectors
-    img, mosaic, skip_frangi, cli_args, save_subdirs, tmp_dir, img_name = load_microscopy_image(cli_args)
+    img, is_tiled, is_fiber, save_dir, tmp_dir, img_name = load_microscopy_image(cli_args)
 
-    # get the fiber orientation analysis pipeline configuration
-    alpha, beta, gamma, scales_um, smooth_sigma, px_size, px_size_iso, \
-        odf_scales_um, odf_degrees, z_min, z_max, ch_neuron, ch_fiber, \
-        lpf_soma_mask, max_ram_mb, jobs, img_name = get_pipeline_config(cli_args, skip_frangi, img_name)
+    # get resources configuration
+    ram, jobs = get_resource_config(cli_args)
 
     # conduct parallel 3D Frangi-based fiber orientation analysis on batches of basic image slices
-    if not skip_frangi:
-        fiber_vec_img, iso_fiber_img \
-            = parallel_frangi_on_slices(img, px_size, px_size_iso, smooth_sigma, save_subdirs[0], tmp_dir, img_name,
-                                        alpha=alpha, beta=beta, gamma=gamma, frangi_sigma_um=scales_um,
-                                        z_min=z_min, z_max=z_max, lpf_soma_mask=lpf_soma_mask,
-                                        ch_neuron=ch_neuron, ch_fiber=ch_fiber, mosaic=mosaic,
-                                        max_ram_mb=max_ram_mb, jobs=jobs)
-
-    # estimate 3D fiber ODF maps over the spatial scales of interest using concurrent workers
-    if odf_scales_um:
-        if skip_frangi:
-            fiber_vec_img = img
-            iso_fiber_img = None
-        parallel_odf_on_scales(fiber_vec_img, iso_fiber_img, px_size_iso, save_subdirs[1], tmp_dir, img_name,
-                               odf_scales_um=odf_scales_um, odf_degrees=odf_degrees, max_ram_mb=max_ram_mb)
+    if not is_fiber:
+        fiber_vec_img, iso_fiber_img, px_sz, img_name \
+            = parallel_frangi_on_slices(img, cli_args, save_dir[0], tmp_dir, img_name, ram=ram, jobs=jobs,
+                                        is_tiled=is_tiled)
+    else:
+        fiber_vec_img, iso_fiber_img, px_sz = (img, None, None)
+
+    # estimate 3D fiber ODF maps at the spatial scales of interest using concurrent workers
+    if cli_args.odf_res:
+        parallel_odf_at_scales(fiber_vec_img, iso_fiber_img, cli_args, px_sz, save_dir[1], tmp_dir, img_name, ram=ram)
 
     # delete temporary folder
     delete_tmp_folder(tmp_dir)

foa3d/frangi.py

@@ -141,9 +141,6 @@ def compute_scaled_hessian(img, sigma=1, trunc=4):
         Hessian matrix of image second derivatives
     """
 
-    # get number of dimensions
-    ndim = img.ndim
-
     # scale selection
     scaled_img = ndi.gaussian_filter(img, sigma=sigma, output=np.float32, truncate=trunc)
 
@@ -152,15 +149,15 @@ def compute_scaled_hessian(img, sigma=1, trunc=4):
 
     # compute the Hessian matrix elements
     hessian_elements = [np.gradient(gradient_list[ax0], axis=ax1)
-                        for ax0, ax1 in combinations_with_replacement(range(ndim), 2)]
+                        for ax0, ax1 in combinations_with_replacement(range(img.ndim), 2)]
 
     # scale the elements of the Hessian matrix
     corr_factor = sigma ** 2
     hessian_elements = [corr_factor * element for element in hessian_elements]
 
     # create the Hessian matrix from its basic elements
-    hessian = np.zeros((ndim, ndim) + scaled_img.shape, dtype=scaled_img.dtype)
-    for index, (ax0, ax1) in enumerate(combinations_with_replacement(range(ndim), 2)):
+    hessian = np.zeros((img.ndim, img.ndim) + scaled_img.shape, dtype=scaled_img.dtype)
+    for index, (ax0, ax1) in enumerate(combinations_with_replacement(range(img.ndim), 2)):
         element = hessian_elements[index]
         hessian[ax0, ax1, ...] = element
         if ax0 != ax1:
@@ -264,30 +261,6 @@ def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma):
     return vesselness
 
 
-def convert_frangi_scales(scales_um, px_size):
-    """
-    Compute the Frangi filter scales in pixel.
-
-    Parameters
-    ----------
-    scales_um: list (dtype=float)
-        Frangi filter scales [μm]
-
-    px_size: int
-        isotropic pixel size [μm]
-
-    Returns
-    -------
-    scales_px: numpy.ndarray (dtype=int)
-        Frangi filter scales [px]
-    """
-
-    scales_um = np.asarray(scales_um)
-    scales_px = scales_um / px_size
-
-    return scales_px
-
-
 def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True):
     """
     Apply 3D Frangi filter to microscopy volume image.
@@ -329,7 +302,7 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True
     n_scales = len(scales_px)
     if n_scales == 1:
         enhanced_img, fiber_vec, eigenval \
-            = compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma,  dark=dark)
+            = compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark)
 
     # parallel scaled vesselness analysis
     else:
@@ -381,17 +354,13 @@ def reject_vesselness_background(vesselness, eigen2, eigen3, dark):
 
     Returns
     -------
-    vesselness: numpy.ndarray (axis order=(Z,Y,X))
+    vesselness: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
         masked Frangi's vesselness likelihood image
     """
 
-    if dark:
-        vesselness[eigen2 < 0] = 0
-        vesselness[eigen3 < 0] = 0
-    else:
-        vesselness[eigen2 > 0] = 0
-        vesselness[eigen3 > 0] = 0
-    vesselness[np.isnan(vesselness)] = 0
+    bg_msk = np.logical_or(eigen2 < 0, eigen3 < 0) if dark else np.logical_or(eigen2 > 0, eigen3 > 0)
+    bg_msk = np.logical_or(bg_msk, np.isnan(vesselness))
+    vesselness[bg_msk] = 0
 
     return vesselness