From 4d2b2845459721f16a95457e2caf272a7ba2148a Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 12:24:44 -0400 Subject: [PATCH 1/7] asym odf clean up --- scilpy/reconst/aodf.py | 44 ++++++++++++++++++++++++ scripts/scil_compute_asym_odf_metrics.py | 35 ------------------- 2 files changed, 44 insertions(+), 35 deletions(-) create mode 100644 scilpy/reconst/aodf.py diff --git a/scilpy/reconst/aodf.py b/scilpy/reconst/aodf.py new file mode 100644 index 000000000..6d4eb6d24 --- /dev/null +++ b/scilpy/reconst/aodf.py @@ -0,0 +1,44 @@ +# -*- coding: utf-8 -*- + + +# +# aodf := Assymetric Orientation Distribution Function (aODF) +# + +import numpy as np + +from dipy.reconst.shm import sph_harm_ind_list + +def compute_asymmetry_index(sh_coeffs, order, mask): + _, l_list = sph_harm_ind_list(order, full_basis=True) + + sign = np.power(-1.0, l_list) + sign = np.reshape(sign, (1, 1, 1, len(l_list))) + sh_squared = sh_coeffs**2 + mask = np.logical_and(sh_squared.sum(axis=-1) > 0., mask) + + asi_map = np.zeros(sh_coeffs.shape[:-1]) + asi_map[mask] = np.sum(sh_squared * sign, axis=-1)[mask] / \ + np.sum(sh_squared, axis=-1)[mask] + + # Negatives should not happen (amplitudes always positive) + asi_map = np.clip(asi_map, 0.0, 1.0) + asi_map = np.sqrt(1 - asi_map**2) * mask + + return asi_map + + +def compute_odd_power_map(sh_coeffs, order, mask): + _, l_list = sph_harm_ind_list(order, full_basis=True) + odd_l_list = (l_list % 2 == 1).reshape((1, 1, 1, -1)) + + odd_order_norm = np.linalg.norm(sh_coeffs * odd_l_list, + ord=2, axis=-1) + + full_order_norm = np.linalg.norm(sh_coeffs, ord=2, axis=-1) + + asym_map = np.zeros(sh_coeffs.shape[:-1]) + mask = np.logical_and(full_order_norm > 0, mask) + asym_map[mask] = odd_order_norm[mask] / full_order_norm[mask] + + return asym_map diff --git a/scripts/scil_compute_asym_odf_metrics.py b/scripts/scil_compute_asym_odf_metrics.py index 1e5de5015..9d68708ec 100755 --- a/scripts/scil_compute_asym_odf_metrics.py +++ b/scripts/scil_compute_asym_odf_metrics.py @@ -109,41 +109,6 @@ def _build_arg_parser(): return p -def compute_asymmetry_index(sh_coeffs, order, mask): - _, l_list = sph_harm_ind_list(order, full_basis=True) - - sign = np.power(-1.0, l_list) - sign = np.reshape(sign, (1, 1, 1, len(l_list))) - sh_squared = sh_coeffs**2 - mask = np.logical_and(sh_squared.sum(axis=-1) > 0., mask) - - asi_map = np.zeros(sh_coeffs.shape[:-1]) - asi_map[mask] = np.sum(sh_squared * sign, axis=-1)[mask] / \ - np.sum(sh_squared, axis=-1)[mask] - - # Negatives should not happen (amplitudes always positive) - asi_map = np.clip(asi_map, 0.0, 1.0) - asi_map = np.sqrt(1 - asi_map**2) * mask - - return asi_map - - -def compute_odd_power_map(sh_coeffs, order, mask): - _, l_list = sph_harm_ind_list(order, full_basis=True) - odd_l_list = (l_list % 2 == 1).reshape((1, 1, 1, -1)) - - odd_order_norm = np.linalg.norm(sh_coeffs * odd_l_list, - ord=2, axis=-1) - - full_order_norm = np.linalg.norm(sh_coeffs, ord=2, axis=-1) - - asym_map = np.zeros(sh_coeffs.shape[:-1]) - mask = np.logical_and(full_order_norm > 0, mask) - asym_map[mask] = odd_order_norm[mask] / full_order_norm[mask] - - return asym_map - - def main(): parser = _build_arg_parser() args = parser.parse_args() From dbd302203795bb73c09fe9c955d0dae9fc9fd426 Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 13:22:01 -0400 Subject: [PATCH 2/7] wip --- scripts/scil_compute_asym_odf_metrics.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/scripts/scil_compute_asym_odf_metrics.py b/scripts/scil_compute_asym_odf_metrics.py index 9d68708ec..dac39da66 100755 --- a/scripts/scil_compute_asym_odf_metrics.py +++ b/scripts/scil_compute_asym_odf_metrics.py @@ -35,6 +35,8 @@ from scilpy.reconst.multi_processes import peaks_from_sh from scilpy.reconst.utils import get_sh_order_and_fullness +from scilpy.reconst.aodf import (compute_asymmetry_index, + compute_odd_power_map) from scilpy.io.utils import (add_processes_arg, add_sh_basis_args, assert_inputs_exist, From 17763d8d2499f1a5408cadf0738842f48b9d7792 Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 13:30:13 -0400 Subject: [PATCH 3/7] pep8 de marde --- scilpy/reconst/aodf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scilpy/reconst/aodf.py b/scilpy/reconst/aodf.py index 6d4eb6d24..f9aa6d2f6 100644 --- a/scilpy/reconst/aodf.py +++ b/scilpy/reconst/aodf.py @@ -9,6 +9,7 @@ from dipy.reconst.shm import sph_harm_ind_list + def compute_asymmetry_index(sh_coeffs, order, mask): _, l_list = sph_harm_ind_list(order, full_basis=True) From cbfd672c50a75728eca6835d2801dd8274e1f661 Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 14:34:35 -0400 Subject: [PATCH 4/7] addressing Charles comments --- scilpy/reconst/aodf.py | 41 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/scilpy/reconst/aodf.py b/scilpy/reconst/aodf.py index f9aa6d2f6..25760debe 100644 --- a/scilpy/reconst/aodf.py +++ b/scilpy/reconst/aodf.py @@ -1,8 +1,6 @@ # -*- coding: utf-8 -*- - - # -# aodf := Assymetric Orientation Distribution Function (aODF) +# aodf := Asymmetric Orientation Distribution Function (aODF) # import numpy as np @@ -11,6 +9,25 @@ def compute_asymmetry_index(sh_coeffs, order, mask): + """ + Compute asymmetry index (ASI) (Cetin-Karayumak et al, 2018) from + asymmetric ODF volume expressed in full SH basis. + + Parameters + ----------------- + sh_coeffs: ndarray (x, y, z, ncoeffs) + Input spherical harmonics coefficients. + order: int > 0 + Max SH order. + mask: ndarray (x, y, z), bool + Mask inside which ASI should be computed. + + Returns + ----------- + asi_map: ndarray (x, y, z) + Asymmetry index map. + """ + _, l_list = sph_harm_ind_list(order, full_basis=True) sign = np.power(-1.0, l_list) @@ -30,6 +47,24 @@ def compute_asymmetry_index(sh_coeffs, order, mask): def compute_odd_power_map(sh_coeffs, order, mask): + """ + Compute odd-power map (Poirier et al, 2020) from + asymmetric ODF volume expressed in full SH basis. + + Parameters + ----------------- + sh_coeffs: ndarray (x, y, z, ncoeffs) + Input spherical harmonics coefficients. + order: int > 0 + Max SH order. + mask: ndarray (x, y, z), bool + Mask inside which odd-power map should be computed. + + Returns + ----------- + odd_power_map: ndarray (x, y, z) + Odd-power map. + """ _, l_list = sph_harm_ind_list(order, full_basis=True) odd_l_list = (l_list % 2 == 1).reshape((1, 1, 1, -1)) From 54e5d3d12812ecc4345d221e1d97e6701046a19c Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 14:35:06 -0400 Subject: [PATCH 5/7] pep8 --- scilpy/reconst/aodf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scilpy/reconst/aodf.py b/scilpy/reconst/aodf.py index 25760debe..bd60c49ad 100644 --- a/scilpy/reconst/aodf.py +++ b/scilpy/reconst/aodf.py @@ -64,7 +64,7 @@ def compute_odd_power_map(sh_coeffs, order, mask): ----------- odd_power_map: ndarray (x, y, z) Odd-power map. - """ + """ _, l_list = sph_harm_ind_list(order, full_basis=True) odd_l_list = (l_list % 2 == 1).reshape((1, 1, 1, -1)) From 74612fc8841a328d9f0519bad14a0e7b65b31ec6 Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 14:52:33 -0400 Subject: [PATCH 6/7] correcting style --- scilpy/reconst/aodf.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scilpy/reconst/aodf.py b/scilpy/reconst/aodf.py index bd60c49ad..4d75adbf2 100644 --- a/scilpy/reconst/aodf.py +++ b/scilpy/reconst/aodf.py @@ -14,7 +14,7 @@ def compute_asymmetry_index(sh_coeffs, order, mask): asymmetric ODF volume expressed in full SH basis. Parameters - ----------------- + ---------- sh_coeffs: ndarray (x, y, z, ncoeffs) Input spherical harmonics coefficients. order: int > 0 @@ -23,7 +23,7 @@ def compute_asymmetry_index(sh_coeffs, order, mask): Mask inside which ASI should be computed. Returns - ----------- + ------- asi_map: ndarray (x, y, z) Asymmetry index map. """ @@ -52,7 +52,7 @@ def compute_odd_power_map(sh_coeffs, order, mask): asymmetric ODF volume expressed in full SH basis. Parameters - ----------------- + ---------- sh_coeffs: ndarray (x, y, z, ncoeffs) Input spherical harmonics coefficients. order: int > 0 @@ -61,7 +61,7 @@ def compute_odd_power_map(sh_coeffs, order, mask): Mask inside which odd-power map should be computed. Returns - ----------- + ------- odd_power_map: ndarray (x, y, z) Odd-power map. """ From f8292a26a8e74160ad6d3e506be6573e8ac72fa4 Mon Sep 17 00:00:00 2001 From: mdesco Date: Mon, 23 Oct 2023 14:53:53 -0400 Subject: [PATCH 7/7] last clean up --- scripts/scil_compute_asym_odf_metrics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/scil_compute_asym_odf_metrics.py b/scripts/scil_compute_asym_odf_metrics.py index dac39da66..40e5661e6 100755 --- a/scripts/scil_compute_asym_odf_metrics.py +++ b/scripts/scil_compute_asym_odf_metrics.py @@ -31,7 +31,6 @@ from dipy.data import get_sphere, SPHERE_FILES from dipy.direction.peaks import reshape_peaks_for_visualization -from dipy.reconst.shm import sph_harm_ind_list from scilpy.reconst.multi_processes import peaks_from_sh from scilpy.reconst.utils import get_sh_order_and_fullness