Compute fiber density image

foa3d/odf.py

@@ -7,7 +7,8 @@
 from foa3d.utils import create_memory_map, normalize_image, transform_axes
 
 
-def compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, vxl_side, odf_norm, odf_deg=6, vxl_thr=0.5, vec_thr=0.000001):
+def compute_odf_map(fbr_vec, px_sz, odf, odi, fbr_dnst, vec_tensor_eigen, vxl_side, odf_norm,
+                    odf_deg=6, vxl_thr=0.5, vec_thr=0.000001):
     """
     Compute the spherical harmonics coefficients iterating over super-voxels
     of fiber orientation vectors.
@@ -17,6 +18,9 @@ def compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, vxl_side, odf_norm, odf
     fbr_vec: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float)
         fiber orientation vectors
 
+    px_sz: numpy.ndarray (shape=(3,), dtype=float)
+        adjusted isotropic pixel size [μm]
+
     odf: NumPy memory-map object (axis order=(X,Y,Z,C), dtype=float32)
                 initialized array of ODF spherical harmonics coefficients
 
@@ -35,6 +39,8 @@ def compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, vxl_side, odf_norm, odf
             odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
                 orientation dispersion index anisotropy
 
+    fbr_dnst
+
     vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
         initialized array of orientation tensor eigenvalues
 
@@ -74,14 +80,17 @@ def compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, vxl_side, odf_norm, odf
                 zerovec = np.count_nonzero(np.all(vec_vxl == 0, axis=-1))
                 sli_vxl_size = np.prod(vec_vxl.shape[:-1])
 
+                # local fiber density estimate
+                zv, yv, xv = z // vxl_side, y // vxl_side, x // vxl_side
+                fbr_dnst[zv, yv, xv] = (sli_vxl_size - zerovec) / (sli_vxl_size * np.prod(px_sz))
+
                 # skip boundary voxels and voxels without enough non-zero orientation vectors
                 if sli_vxl_size / ref_vxl_size > vxl_thr and 1 - zerovec / sli_vxl_size > vec_thr:
 
                     # flatten orientation supervoxel
                     vec_arr = vec_vxl.ravel()
 
                     # compute ODF and orientation tensor eigenvalues
-                    zv, yv, xv = z // vxl_side, y // vxl_side, x // vxl_side
                     odf[zv, yv, xv, :] = fiber_vectors_to_sph_harm(vec_arr, odf_deg, odf_norm)
                     vec_tensor_eigen[zv, yv, xv, :] = compute_vec_tensor_eigen(vec_arr)
 
@@ -94,7 +103,7 @@ def compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, vxl_side, odf_norm, odf
     # transform axes
     odf = transform_axes(odf, swapped=(0, 2), flipped=(1, 2))
 
-    return odf
+    return odf, fbr_dnst
 
 
 @njit(cache=True)
@@ -262,8 +271,11 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False):
             odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
                 initialized array of orientation dispersion anisotropy parameters
 
+    fbr_dnst: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
+        initialized fiber density map
+
     bg_mrtrix: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8)
-                initialized background for ODF visualization in MRtrix3
+        initialized background for ODF visualization in MRtrix3
 
     vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
         initialized fiber orientation tensor eigenvalues
@@ -285,6 +297,9 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False):
     vec_tensor_eigen = create_memory_map(vec_tensor_shape, dtype='float32',
                                          name=f'tensor_tmp{odf_scale}', tmp_dir=tmp_dir)
 
+    # fiber density memory map
+    fbr_dnst = create_memory_map(odi_shp, dtype='float32', name=f'fbr_dnst_tmp{odf_scale}', tmp_dir=tmp_dir)
+
     # create ODI memory maps
     odi_tot = create_memory_map(odi_shp, dtype='float32', name=f'odi_tot_tmp{odf_scale}', tmp_dir=tmp_dir)
     if exp_all:
@@ -301,7 +316,7 @@ def init_odf_arrays(vec_img_shp, tmp_dir, odf_scale, odf_deg=6, exp_all=False):
     odi['odi_tot'] = odi_tot
     odi['odi_anis'] = odi_anis
 
-    return odf, odi, bg_mrtrix, vec_tensor_eigen
+    return odf, odi, fbr_dnst, bg_mrtrix, vec_tensor_eigen
 
 
 def mask_orientation_dispersion(vec_tensor_eigen, odi):

foa3d/output.py

@@ -196,7 +196,7 @@ def save_frangi_arrays(save_dir, img_name, out_img, px_sz, ram=None):
     print_flsh(f"\nFrangi filter arrays saved to: {save_dir}\n")
 
 
-def save_odf_arrays(save_dir, img_name, odf_scale_um, px_sz, odf, bg, odi_pri, odi_sec, odi_tot, odi_anis):
+def save_odf_arrays(save_dir, img_name, odf_scale_um, px_sz, odf, bg, fbr_dnst, odi_pri, odi_sec, odi_tot, odi_anis):
     """
     Save the output arrays of the ODF analysis stage to TIF and Nifti files.
     Arrays tagged with 'mrtrixview' are preliminarily transformed
@@ -223,6 +223,9 @@ def save_odf_arrays(save_dir, img_name, odf_scale_um, px_sz, odf, bg, odi_pri, o
     bg: NumPy memory-map object (axis order=(X,Y,Z), dtype=uint8)
         background for ODF visualization in MRtrix3
 
+    fbr_dnst: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
+        fiber orientation density [1/μm³]
+
     odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
         primary orientation dispersion parameter
 
@@ -244,6 +247,9 @@ def save_odf_arrays(save_dir, img_name, odf_scale_um, px_sz, odf, bg, odi_pri, o
     save_array(f'bg_mrtrixview_sv{sbfx}', save_dir, bg, fmt='nii')
     save_array(f'odf_mrtrixview_sv{sbfx}', save_dir, odf, fmt='nii')
 
+    # save fiber density
+    save_array(f'fbr_dnst_sv{sbfx}', save_dir, fbr_dnst, px_sz)
+
     # save total orientation dispersion
     save_array(f'odi_tot_sv{sbfx}', save_dir, odi_tot, px_sz)
 

foa3d/pipeline.py

@@ -633,14 +633,14 @@ def odf_analysis(fbr_vec, iso_fbr, px_sz, save_dirs, img_name, odf_scale_um, odf
     odf_scale = int(np.ceil(odf_scale_um / px_sz[0]))
 
     # initialize ODF analysis output arrays
-    odf, odi, bg_mrtrix, vec_tensor_eigen = \
+    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)
 
     # generate downsampled background for MRtrix3 mrview
     generate_odf_background(bg_mrtrix, fbr_vec, vxl_side=odf_scale, iso_fbr=iso_fbr)
 
     # compute ODF coefficients
-    odf = compute_odf_map(fbr_vec, odf, odi, vec_tensor_eigen, odf_scale, odf_norm, odf_deg=odf_deg)
+    odf, dnst = compute_odf_map(fbr_vec, px_sz, odf, odi, dnst, vec_tensor_eigen, odf_scale, odf_norm, odf_deg=odf_deg)
 
     # save memory maps to TIFF or NIfTI files
-    save_odf_arrays(save_dirs['odf'], img_name, odf_scale_um, px_sz, odf, bg_mrtrix, **odi)
+    save_odf_arrays(save_dirs['odf'], img_name, odf_scale_um, px_sz, odf, bg_mrtrix, dnst, **odi)

foa3d/printing.py

@@ -35,7 +35,7 @@ def color_text(r, g, b, text):
     return clr_text
 
 
-def print_flsh(string_to_print=""):
+def print_flsh(string_to_print="", end='\n'):
     """
     Print string and flush output data buffer.
 
@@ -44,11 +44,14 @@ def print_flsh(string_to_print=""):
     string_to_print: str
         string to be printed
 
+    end: str
+        string appended after the last value, default a newline
+
     Returns
     -------
     None
     """
-    print(string_to_print, flush=True)
+    print(string_to_print, flush=True, end=end)
 
 
 def print_blur(sigma_um, psf_fwhm):
@@ -184,7 +187,7 @@ def print_frangi_info(in_img, frangi_cfg, in_slc_shp_um, tot_slc_num):
     if gamma is None:
         gamma = 'auto'
 
-    print_flsh(color_text(0, 191, 255, "\n3D Frangi Filter\n") + "\nSensitivity\n" +
+    print_flsh(color_text(0, 191, 255, "\n\n\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")
 
@@ -223,7 +226,7 @@ def print_frangi_progress(info, is_valid, verbose=1):
         prog_prc = 100 * slc_cnt / tot_slc
         _, hrs, mins, secs = elapsed_time(start_time)
         print_flsh(f"[Parallel(n_jobs={batch_sz})]:\t{slc_cnt}/{tot_slc} done\t|\t" +
-                   f"elapsed: {hrs} hrs {mins} min {int(secs)} s\t{prog_prc}%")
+                   f"elapsed: {hrs} hrs {mins} min {int(secs)} s\t{prog_prc:.1f}%")
 
 
 def print_image_info(in_img):
@@ -292,8 +295,8 @@ def print_odf_info(odf_scales_um, odf_degrees):
     None
     """
     print_flsh(
-        color_text(0, 191, 255,
-                   f"\n3D ODF Analysis\n\nResolution   [μm]: {odf_scales_um}\nExpansion degrees: {odf_degrees}\n"))
+        color_text(0, 191, 255, "\n3D ODF Analysis") +
+        f"\n\nResolution   [μm]: {odf_scales_um}\nExpansion degrees: {odf_degrees}\n")
 
 
 def print_output_res(px_sz_iso):

foa3d/slicing.py

@@ -4,6 +4,7 @@
 import psutil
 
 from numba import njit
+from foa3d.printing import print_flsh
 from foa3d.utils import get_available_cores
 
 
@@ -523,42 +524,41 @@ def generate_slice_lists(in_slc_shp, img_shp, batch_sz, px_rsz_ratio, ovlp=0, ms
     slc_per_dim = np.floor(np.divide(img_shp, in_slc_shp)).astype(int)
 
     # initialize empty range lists
-    in_pad_lst = list()
-    in_rng_lst = list()
-    bc_rng_lst = list()
-    out_rng_lst = list()
-    for z, y, x in product(range(slc_per_dim[0]), range(slc_per_dim[1]), range(slc_per_dim[2])):
+    # 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]))):
 
         # index ranges of analyzed image slices (with padding)
-        in_rng, pad = compute_slice_range((z, y, x), in_slc_shp, img_shp, slc_per_dim, ovlp=ovlp)
-        in_rng_lst.append(in_rng)
-        in_pad_lst.append(pad)
+        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)
 
         # output index ranges
-        out_rng, _ = compute_slice_range((z, y, x), out_slc_shp, out_img_shp, slc_per_dim)
-        out_rng_lst.append(out_rng)
+        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((z, y, x), in_slc_shp, img_shp, slc_per_dim)
-            bc_rng_lst.append(bc_rng)
+            bc_rng, _ = compute_slice_range(zyx, in_slc_shp, img_shp, slc_per_dim)
+            slc_rng['bc_rng'].append(bc_rng)
         else:
-            bc_rng_lst.append(None)
+            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(in_rng_lst)
+    tot_slc = len(slc_rng['in_rng'])
 
     # adjust slice batch size
     if batch_sz > tot_slc:
         batch_sz = tot_slc
 
-    # fill slice range dictionary
-    slc_rng = dict()
-    slc_rng['in_rng'] = in_rng_lst
-    slc_rng['in_pad'] = in_pad_lst
-    slc_rng['out_rng'] = out_rng_lst
-    slc_rng['bc_rng'] = bc_rng_lst
-
     return slc_rng, out_slc_shp, tot_slc, batch_sz