diff --git a/foa3d/frangi.py b/foa3d/frangi.py index 0efd4c2..1349695 100644 --- a/foa3d/frangi.py +++ b/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): diff --git a/foa3d/pipeline.py b/foa3d/pipeline.py index e5cdb92..5176c29 100644 --- a/foa3d/pipeline.py +++ b/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, diff --git a/foa3d/slicing.py b/foa3d/slicing.py index cb55070..36e9b65 100644 --- a/foa3d/slicing.py +++ b/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.