Merge branch 'mic-devel' of https://github.com/lens-biophotonics/fibe…

…r-orientation into mic-devel

foa3d/frangi.py

@@ -127,7 +127,7 @@ def compute_scaled_hessian(img, sigma=1, trunc=4):
     Parameters
     ----------
     img: numpy.ndarray (axis order=(Z,Y,X))
-        microscopy volume image
+        3D microscopy image
 
     sigma: int
         spatial scale [px]
@@ -179,7 +179,7 @@ def compute_scaled_orientation(scale_px, img, alpha=0.001, beta=1, gamma=None, d
         spatial scale [px]
 
     img: numpy.ndarray (axis order=(Z,Y,X))
-        microscopy volume image
+        3D microscopy image
 
     alpha: float
         plate-like score sensitivity
@@ -263,12 +263,12 @@ def compute_scaled_vesselness(eigen1, eigen2, eigen3, alpha, beta, gamma):
 
 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.
+    Apply 3D Frangi filter to 3D microscopy image.
 
     Parameters
     ----------
     img: numpy.ndarray (axis order=(Z,Y,X))
-        microscopy volume image
+        3D microscopy image
 
     scales_px: int or numpy.ndarray (dtype=int)
         analyzed spatial scales [px]
@@ -288,10 +288,10 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True
 
     Returns
     -------
-    enhanced_img: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
+    frangi_img: numpy.ndarray (axis order=(Z,Y,X), dtype=float)
         Frangi's vesselness likelihood image
 
-    fiber_vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
+    fbr_vec: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
         3D fiber orientation map
 
     eigenval: numpy.ndarray (axis order=(Z,Y,X,C), dtype=float)
@@ -301,12 +301,12 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True
     # single-scale vesselness analysis
     n_scales = len(scales_px)
     if n_scales == 1:
-        enhanced_img, fiber_vec, eigenval \
+        frangi_img, fbr_vec, eigenval \
             = compute_scaled_orientation(scales_px[0], img, alpha=alpha, beta=beta, gamma=gamma, dark=dark)
 
     # parallel scaled vesselness analysis
     else:
-        with Parallel(n_jobs=n_scales, backend='threading', max_nbytes=None) as parallel:
+        with Parallel(n_jobs=n_scales, prefer='threads', require='sharedmem') as parallel:
             par_res = \
                 parallel(
                     delayed(compute_scaled_orientation)(
@@ -317,19 +317,19 @@ def frangi_filter(img, scales_px=1, alpha=0.001, beta=1.0, gamma=None, dark=True
             enhanced_img_tpl, eigenvectors_tpl, eigenvalues_tpl = zip(*par_res)
             eigenval = np.stack(eigenvalues_tpl, axis=0)
             eigenvec = np.stack(eigenvectors_tpl, axis=0)
-            enhanced_img = np.stack(enhanced_img_tpl, axis=0)
+            frangi_img = np.stack(enhanced_img_tpl, axis=0)
 
             # get max scale-wise vesselness
-            best_idx = np.argmax(enhanced_img, axis=0)
+            best_idx = np.argmax(frangi_img, axis=0)
             best_idx = np.expand_dims(best_idx, axis=0)
-            enhanced_img = np.take_along_axis(enhanced_img, best_idx, axis=0).squeeze(axis=0)
+            frangi_img = np.take_along_axis(frangi_img, 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)
-            fiber_vec = np.take_along_axis(eigenvec, best_idx, axis=0).squeeze(axis=0)
+            fbr_vec = np.take_along_axis(eigenvec, best_idx, axis=0).squeeze(axis=0)
 
-    return enhanced_img, fiber_vec, eigenval
+    return frangi_img, fbr_vec, eigenval
 
 
 def reject_fiber_background(vesselness, eigen2, eigen3, dark):

foa3d/pipeline.py

@@ -657,7 +657,7 @@ def odf_analysis(fbr_vec_img, iso_fbr_img, px_sz_iso, save_dir, tmp_dir, img_nam
 
 
 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):
+                              ram=None, jobs=4, backend='threads', is_tiled=False, verbose=10):
     """
     Apply 3D Frangi filtering to basic TPFM image slices using parallel threads.
 
@@ -737,7 +737,7 @@ def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name, ts_msk
 
     # parallel Frangi filter-based fiber orientation analysis of microscopy image slices
     start_time = perf_counter()
-    with Parallel(n_jobs=batch_sz, backend=backend, verbose=verbose, max_nbytes=None) as parallel:
+    with Parallel(n_jobs=batch_sz, prefer=backend, verbose=verbose, require='sharedmem') as parallel:
         parallel(
             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,

foa3d/slicing.py

@@ -201,7 +201,7 @@ def compute_overlap(smooth_sigma, frangi_sigma_um, px_rsz_ratio=None, trunc=2):
         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) / px_rsz_ratio
     if smooth_sigma is not None:
         sigma = np.concatenate((smooth_sigma, sigma))
 
@@ -247,9 +247,11 @@ def compute_slice_shape(img_shp, item_sz, px_sz=None, slc_sz=1e6, ovlp=0):
     if len(img_shp) == 4:
         item_sz *= 3
 
-    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)
+    tot_ovlp = 2 * ovlp
+    side_xyz = np.floor((slc_sz / item_sz)**(1/3)).astype(int)
+    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
@@ -365,7 +367,7 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi
     ns = len(frangi_sigma_um)
 
     # initialize slice batch size
-    batch_sz = np.min([jobs // ns, num_cpu]).astype(int) + 1
+    batch_sz = np.min([jobs, num_cpu]).astype(int) + 1
 
     # get pixel resize ratio
     px_rsz_ratio = np.divide(px_sz, px_sz_iso)
@@ -387,7 +389,7 @@ def config_frangi_batch(px_sz, px_sz_iso, img_shp, item_sz, smooth_sigma, frangi
     return batch_sz, in_slc_shp, in_slc_shp_um, px_rsz_ratio, ovlp, ovlp_rsz
 
 
-def crop(img, rng, ovlp, flip=()):
+def crop(img, rng, ovlp=None, flip=()):
     """
     Shrink image slice at total volume boundaries, for overall shape consistency.