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

@@ -223,7 +223,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 +292,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):