Compute local angular disarray image

foa3d/odf.py

@@ -6,6 +6,28 @@
 from foa3d.utils import normalize_image, transform_axes
 
 
+def compute_disarray(fiber_vec):
+    """
+    Compute local angular disarray, i.e. the misalignment degree of nearby orientation vectors
+    with respect to the mean direction (Giardini et al., Front. Physiol. 2021).
+
+    Parameters
+    ----------
+    fiber_vec: numpy.ndarray (shape=(N,3), dtype=float)
+        array of fiber orientation vectors
+        (reshaped super-voxel of shape=(Nz,Ny,Nx), i.e. N=Nz*Ny*Nx)
+
+    Returns
+    -------
+    disarray: float32
+        local angular disarray
+    """
+    fiber_vec = np.delete(fiber_vec, np.all(fiber_vec == 0, axis=-1), axis=0)
+    disarray = (1 - np.linalg.norm(np.mean(fiber_vec, axis=0))).astype(np.float32)
+
+    return disarray
+
+
 @njit(cache=True)
 def compute_fiber_angles(fiber_vec, norm):
     """
@@ -38,7 +60,7 @@ def compute_fiber_angles(fiber_vec, norm):
     return phi, theta
 
 
-def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, vec_tensor_eigen, vxl_side, odf_norm,
+def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, vec_tensor_eigen, vxl_side, odf_norm,
                     odf_deg=6, vxl_thr=0.5, vec_thr=-1):
     """
     Compute the spherical harmonics coefficients iterating over super-voxels
@@ -64,6 +86,9 @@ def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, vec_ten
     odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
         orientation dispersion index anisotropy
 
+    disarray: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
+        local angular disarray
+
     vec_tensor_eigen: NumPy memory-map object (axis order=(Z,Y,X,C), dtype=float32)
         initialized array of orientation tensor eigenvalues
 
@@ -114,6 +139,7 @@ def compute_odf_map(fiber_vec, odf, odi_pri, odi_sec, odi_tot, odi_anis, vec_ten
                     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)
+                    disarray[zv, yv, xv] = compute_disarray(vec_arr)
 
     # compute slice-wise dispersion and anisotropy parameters
     compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot, odi_anis)

foa3d/pipeline.py

@@ -190,6 +190,9 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, ram=None)
     odi_anis: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
         initialized array of orientation dispersion anisotropy parameters
 
+    disarray: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32)
+        local angular disarray
+
     vec_tensor_eigen: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X,C), dtype=float32)
         initialized array of fiber orientation tensor eigenvalues
     """
@@ -224,7 +227,11 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, ram=None)
     odi_anis, _ = init_volume(odi_shape, dtype='uint8',
                               name='odi_anis_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram)
 
-    return odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, vec_tensor_eigen
+    # create disarray map
+    disarray, _ = init_volume(odi_shape, dtype='float32',
+                              name='disarray{0}'.format(odf_scale), tmp=tmp_dir, ram=ram)
+
+    return odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray, vec_tensor_eigen
 
 
 def init_volume(shape, dtype, chunks=True, name='tmp', tmp=None, mmap_mode='r+', ram=None):
@@ -623,19 +630,20 @@ def odf_analysis(fiber_vec_img, iso_fiber_img, px_size_iso, save_dir, tmp_dir, i
     odf_scale = int(np.ceil(odf_scale_um / px_size_iso[0]))
 
     # initialize ODF analysis output volumes
-    odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, tensor \
+    odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray, tensor \
         = init_odf_volumes(fiber_vec_img.shape[:-1], tmp_dir, odf_scale=odf_scale, odf_degrees=odf_deg, ram=ram)
 
     # generate downsampled background for MRtrix3 mrview
     bg_img = fiber_vec_img if iso_fiber_img is None else iso_fiber_img
     generate_odf_background(bg_img, bg_mrtrix, vxl_side=odf_scale)
 
     # compute ODF coefficients
-    odf = compute_odf_map(fiber_vec_img, odf, odi_pri, odi_sec, odi_tot, odi_anis, tensor, odf_scale, odf_norm,
-                          odf_deg=odf_deg)
+    odf = compute_odf_map(fiber_vec_img, odf, odi_pri, odi_sec, odi_tot, odi_anis, disarray, tensor,
+                          odf_scale, odf_norm, odf_deg=odf_deg)
 
     # save memory maps to file
-    save_odf_arrays(odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, px_size_iso, save_dir, img_name, odf_scale_um)
+    save_odf_arrays(odf, bg_mrtrix, odi_pri, odi_sec, odi_tot, odi_anis, disarray,
+                    px_size_iso, save_dir, img_name, odf_scale_um)
 
 
 def parallel_frangi_on_slices(img, cli_args, save_dir, tmp_dir, img_name,
@@ -889,7 +897,7 @@ def save_frangi_arrays(fiber_dset_path, fiber_vec_img, fiber_vec_clr, frac_anis_
         save_array('soma_msk_{0}'.format(img_name), save_dir, soma_msk, px_size)
 
 
-def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_dir, img_name, odf_scale_um):
+def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, disarray, px_size, save_dir, img_name, odf_scale_um):
     """
     Save the output arrays of the ODF analysis stage to TIF and Nifti files.
     Arrays tagged with 'mrtrixview' are preliminarily transformed
@@ -916,6 +924,8 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_
     odi_anis: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
         orientation dispersion anisotropy parameter
 
+    disarray: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32)
+
     px_size: numpy.ndarray (shape=(3,), dtype=float)
         pixel size (Z,Y,X) [μm]
 
@@ -939,3 +949,4 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_
     save_array('odi_sec_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_sec, px_size, odi=True)
     save_array('odi_tot_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_tot, px_size, odi=True)
     save_array('odi_anis_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, odi_anis, px_size, odi=True)
+    save_array('disarray_sv{0}_{1}'.format(odf_scale_um, img_name), save_dir, disarray, px_size, odi=True)