Skip to content

Commit

Permalink
Compute local angular disarray image and improve orientation dispersi…
Browse files Browse the repository at this point in the history 
…on indices
  • Loading branch information
@msorelli
msorelli committed Jan 29, 2024
1 parent 9211a7e commit 50903e5
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 29 deletions.
52 changes: 37 additions & 15 deletions foa3d/odf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand Down Expand Up @@ -38,8 +60,8 @@ 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,
odf_deg=6, vxl_thr=0.5, vec_thr=-1):
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=0.000001):
"""
Compute the spherical harmonics coefficients iterating over super-voxels
of fiber orientation vectors.
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -152,16 +178,16 @@ def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot,
Returns
-------
odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
odi_pri: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
primary orientation dispersion index
odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
odi_sec: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
secondary orientation dispersion index
odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
odi_tot: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
total orientation dispersion index
odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=uint8)
odi_anis: NumPy memory-map object (axis order=(Z,Y,X), dtype=float32)
orientation dispersion index anisotropy
"""

Expand All @@ -170,21 +196,17 @@ def compute_orientation_dispersion(vec_tensor_eigen, odi_pri, odi_sec, odi_tot,
k3 = np.abs(vec_tensor_eigen[..., 0])
diff = np.abs(vec_tensor_eigen[..., 1] - vec_tensor_eigen[..., 0])

# primary dispersion (1.2732395447351628 = 4/π)
odi_pri[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, k2))) \
.astype(np.uint8)
# primary dispersion (0.3183098861837907 = 1/π)
odi_pri[:] = (0.5 - 0.3183098861837907 * np.arctan2(k1, k2)).astype(np.float32)

# secondary dispersion
odi_sec[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, k3))) \
.astype(np.uint8)
odi_sec[:] = (0.5 - 0.3183098861837907 * np.arctan2(k1, k3)).astype(np.float32)

# total dispersion
odi_tot[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, np.sqrt(np.abs(np.multiply(k2, k3)))))) \
.astype(np.uint8)
odi_tot[:] = (0.5 - 0.3183098861837907 * np.arctan2(k1, np.sqrt(np.abs(np.multiply(k2, k3))))).astype(np.float32)

# dispersion anisotropy
odi_anis[:] = (255 * (2 - 1.2732395447351628 * np.arctan2(k1, diff))) \
.astype(np.uint8)
odi_anis[:] = (0.5 - 0.3183098861837907 * np.arctan2(k1, diff)).astype(np.float32)


@njit(cache=True)
Expand Down
39 changes: 25 additions & 14 deletions foa3d/pipeline.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""
Expand All @@ -215,16 +218,20 @@ def init_odf_volumes(vec_img_shape, tmp_dir, odf_scale, odf_degrees=6, ram=None)
name='tensor_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram)

# create ODI memory maps
odi_pri, _ = init_volume(odi_shape, dtype='uint8',
odi_pri, _ = init_volume(odi_shape, dtype='float32',
name='odi_pri_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram)
odi_sec, _ = init_volume(odi_shape, dtype='uint8',
odi_sec, _ = init_volume(odi_shape, dtype='float32',
name='odi_sec_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram)
odi_tot, _ = init_volume(odi_shape, dtype='uint8',
odi_tot, _ = init_volume(odi_shape, dtype='float32',
name='odi_tot_tmp{0}'.format(odf_scale), tmp=tmp_dir, ram=ram)
odi_anis, _ = init_volume(odi_shape, dtype='uint8',
odi_anis, _ = init_volume(odi_shape, dtype='float32',
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):
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -904,18 +912,20 @@ def save_odf_arrays(odf, bg, odi_pri, odi_sec, odi_tot, odi_anis, px_size, save_
bg: NumPy memory-map object or HDF5 dataset (axis order=(X,Y,Z), dtype=uint8)
background for ODF visualization in MRtrix3
odi_pri: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
odi_pri: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32)
primary orientation dispersion parameter
odi_sec: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
odi_sec: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32)
secondary orientation dispersion parameter
odi_tot: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
odi_tot: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32)
total orientation dispersion parameter
odi_anis: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=uint8)
odi_anis: NumPy memory-map object or HDF5 dataset (axis order=(Z,Y,X), dtype=float32)
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]
Expand All @@ -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)

0 comments on commit 50903e5

Please sign in to comment.