From 989464f137e8b5edec7da47f821d210c83097592 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 22 Feb 2024 13:37:03 +0100 Subject: [PATCH 1/5] Feat: Improve Cerebellum --- .../layer_based_direction_vectors.py | 2 +- atlas_direction_vectors/cerebellum.py | 33 +++++++++++-------- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py b/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py index 0decece..9e795a1 100644 --- a/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py +++ b/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py @@ -15,6 +15,7 @@ This script is used to compute the mouse isocortex and the mouse thalamus directions vectors. """ + from __future__ import annotations import enum @@ -331,7 +332,6 @@ def compute_layered_region_direction_vectors( layer_to_weight, border_region_mask, shading_complement, - expansion_width, ) # remove the grown regions into the void from the dilation and the external_id values diff --git a/atlas_direction_vectors/cerebellum.py b/atlas_direction_vectors/cerebellum.py index a979a66..67d346e 100644 --- a/atlas_direction_vectors/cerebellum.py +++ b/atlas_direction_vectors/cerebellum.py @@ -4,9 +4,12 @@ and high values where fiber tracts are outgoing. The direction vectors are given by the gradient of this scalar field. """ + import logging from typing import TYPE_CHECKING, Optional +from functools import partial +from joblib import Parallel, delayed import numpy as np from atlas_commons.typing import FloatArray @@ -49,10 +52,16 @@ def compute_direction_vectors(region_map: "RegionMap", annotation: "VoxelData") and parent_id not in subregion_ids ): subregion_ids.append(parent_id) - L.info("Computing direction vectors for region %s.", region_map.get(parent_id, "name")) - subregion_direction_vectors = cereb_subregion_direction_vectors( - parent_id, region_map, annotation - ) + L.info( + "Computing direction vectors for regions %s.", + [region_map.get(parent_id, "name") for parent_id in subregion_ids], + ) + + f = partial(cereb_subregion_direction_vectors, region_map=region_map, annotation=annotation) + with Parallel(n_jobs=-1, verbose=10) as parallel: + for subregion_direction_vectors in parallel( + delayed(f)(parent_id) for parent_id in subregion_ids + ): # Assembles subregion direction vectors. subregion_mask = np.logical_not(np.isnan(subregion_direction_vectors)) direction_vectors[subregion_mask] = subregion_direction_vectors[subregion_mask] @@ -91,23 +100,19 @@ def cereb_subregion_direction_vectors( "layers": { "names": [ name + ", molecular layer", - name + ", Purkinje layer", name + ", granular layer", "cerebellum related fiber tracts", ], - "queries": [acronym + "mo", acronym + "pu", acronym + "gr", "cbf"], + "queries": [acronym + "mo", acronym + "gr", "cbf"], "attribute": "acronym", "with_descendants": True, }, } region_to_weight = { - acronym + "mo": 1.0 if "mo" not in weights else weights["mo"], - acronym + "pu": 0.0 if "pu" not in weights else weights["pu"], - acronym + "gr": -1.0 if "gr" not in weights else weights["gr"], - "cbf": -5.0 if "cbf" not in weights else weights["cbf"], - "outside_of_brain": 3.0 - if "outside_of_brain" not in weights - else weights["outside_of_brain"], + acronym + "mo": weights.get("mo", 1.0), + acronym + "gr": weights.get("gr", -1.0), + "cbf": weights.get("cbf", -5.0), + "outside_of_brain": weights.get("outside_of_brain", 3), } return compute_layered_region_direction_vectors( @@ -117,5 +122,5 @@ def cereb_subregion_direction_vectors( region_to_weight=region_to_weight, shading_width=4, expansion_width=8, - has_hemispheres=False, + has_hemispheres=True, ) From 592a88fdc975cbc260cc01c63094fef0f24a8b4f Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 29 Feb 2024 11:06:27 +0100 Subject: [PATCH 2/5] option to enfore x/z directions --- atlas_direction_vectors/algorithms/utils.py | 41 +++++++++++++++++-- .../app/orientation_field.py | 11 ++++- atlas_direction_vectors/cerebellum.py | 4 +- 3 files changed, 49 insertions(+), 7 deletions(-) diff --git a/atlas_direction_vectors/algorithms/utils.py b/atlas_direction_vectors/algorithms/utils.py index 3362212..421ba37 100644 --- a/atlas_direction_vectors/algorithms/utils.py +++ b/atlas_direction_vectors/algorithms/utils.py @@ -1,11 +1,16 @@ """Low-level tools for the computation of direction vectors""" +import logging +from typing import Union + import numpy as np # type: ignore from atlas_commons.utils import FloatArray, NumericArray, normalize from scipy.ndimage import gaussian_filter # type: ignore from scipy.ndimage import generate_binary_structure # type: ignore from scipy.signal import correlate # type: ignore +L = logging.getLogger(__name__) + def compute_blur_gradient( scalar_field: FloatArray, gaussian_stddev=3.0, gaussian_radius=None @@ -45,26 +50,53 @@ def compute_blur_gradient( def _quaternion_from_vectors( # pylint: disable=invalid-name - s: NumericArray, t: NumericArray + s: NumericArray, t: NumericArray, align_to: Union[NumericArray, None] = None ) -> NumericArray: """ Returns the quaternion (s cross t) + (s dot t + |s||t|). This quaternion q maps s to t, i.e., qsq^{-1} = t. + If align_to is specified (either `x`, or `z`), we will rotate it around its axis ([0, 1, 0]) + so that it is maximally aligned with `align_to`. + Args: s: numeric array of shape (3,) or (N, 3). t: numeric array of shape (N, 3) if s has two dimensions and its first dimension is N. + align_to: if not None, the returned quaternion is aligned with the given axis. Returns: Numeric array of shape (N, 4) where N is the first dimension of t. This data is interpreted as a 1D array of quaternions with size N. A quaternion is a 4D vector [w, x, y, z] where [x, y, z] is the imaginary part. """ w = np.matmul(s, np.array(t).T) + np.linalg.norm(s, axis=-1) * np.linalg.norm(t, axis=-1) - return np.hstack([w[:, np.newaxis], np.cross(s, t)]) + quaternions = np.hstack([w[:, np.newaxis], np.cross(s, t)]) + + if align_to is not None: + target_dir = {"x": np.array([1.0, 0.0, 0.0]), "z": np.array([0.0, 0.0, 1.0])}[align_to] + + import pyquaternion as pyq + from joblib import Parallel, delayed + from scipy.optimize import minimize + + def rotate_q(_q): + q = pyq.Quaternion(_q) + + def cost(a): + q_rot = pyq.Quaternion(axis=[0, 1, 0], angle=a) + return -(q * q_rot).rotate(target_dir).dot(target_dir) + + angle = minimize(cost, 0.0).x + return (q * pyq.Quaternion(axis=[0, 1, 0], angle=angle)).q.tolist() + + L.info("We are aligning %s quaternions to %s", len(quaternions), align_to) + with Parallel(n_jobs=-1, batch_size=10000, verbose=10) as parallel: + quaternions = np.array(parallel(delayed(rotate_q)(q) for q in quaternions)) + + return quaternions -def vector_to_quaternion(vector_field: FloatArray) -> FloatArray: +def vector_to_quaternion(vector_field: FloatArray, align_to: Union[str, None] = None) -> FloatArray: """ Find quaternions which rotate [0.0, 1.0, 0.0] to each vector in `vector_field`. @@ -75,6 +107,7 @@ def vector_to_quaternion(vector_field: FloatArray) -> FloatArray: Arguments: vector_field: field of floating point 3D unit vectors, i.e., a float array of shape (..., 3). + align_to: if not None, the returned quaternion is aligned with the given axis. Returns: numpy.ndarray of shape (..., 4) and of the same type as the input @@ -88,7 +121,7 @@ def vector_to_quaternion(vector_field: FloatArray) -> FloatArray: quaternions = np.full(vector_field.shape[:-1] + (4,), np.nan, dtype=vector_field.dtype) non_nan_mask = ~np.isnan(np.linalg.norm(vector_field, axis=-1)) quaternions[non_nan_mask] = _quaternion_from_vectors( - [0.0, 1.0, 0.0], vector_field[non_nan_mask] + [0.0, 1.0, 0.0], vector_field[non_nan_mask], align_to=align_to ) return quaternions diff --git a/atlas_direction_vectors/app/orientation_field.py b/atlas_direction_vectors/app/orientation_field.py index fd015fd..ab2067e 100644 --- a/atlas_direction_vectors/app/orientation_field.py +++ b/atlas_direction_vectors/app/orientation_field.py @@ -17,7 +17,9 @@ convention used by the placement algorithm, see https://bbpteam.epfl.ch/documentation/projects/placement-algorithm/latest/index.html. """ + import logging +from typing import Union import click import voxcell # type: ignore @@ -45,15 +47,22 @@ "NaN vectors indicate out-of-domain voxels but also voxels for which an orientation could" " not be derived.", ) +@click.option( + "--align-to", + type=str, + required=False, + help=("Direction to align the direction vectors to."), +) @log_args(L) def cmd( direction_vectors_path: str, output_path: str, + align_to: Union[str, None] = None, ) -> None: """Turn direction vectors into quaternions interpreted as 3D orientations.""" direction_vectors = voxcell.VoxelData.load_nrrd(direction_vectors_path) - quaternions = vector_to_quaternion(direction_vectors.raw) + quaternions = vector_to_quaternion(direction_vectors.raw, align_to=align_to) voxcell.OrientationField( quaternions, direction_vectors.voxel_dimensions, direction_vectors.offset ).save_nrrd(output_path) diff --git a/atlas_direction_vectors/cerebellum.py b/atlas_direction_vectors/cerebellum.py index 67d346e..ed6322d 100644 --- a/atlas_direction_vectors/cerebellum.py +++ b/atlas_direction_vectors/cerebellum.py @@ -6,12 +6,12 @@ """ import logging -from typing import TYPE_CHECKING, Optional from functools import partial +from typing import TYPE_CHECKING, Optional -from joblib import Parallel, delayed import numpy as np from atlas_commons.typing import FloatArray +from joblib import Parallel, delayed from atlas_direction_vectors.algorithms.layer_based_direction_vectors import ( compute_layered_region_direction_vectors, From 5cd82ff194b03f745a4c9492d3745f5086a379eb Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 29 Feb 2024 11:27:04 +0100 Subject: [PATCH 3/5] fix --- .../algorithms/layer_based_direction_vectors.py | 1 + setup.py | 1 + 2 files changed, 2 insertions(+) diff --git a/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py b/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py index 9e795a1..79cee74 100644 --- a/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py +++ b/atlas_direction_vectors/algorithms/layer_based_direction_vectors.py @@ -332,6 +332,7 @@ def compute_layered_region_direction_vectors( layer_to_weight, border_region_mask, shading_complement, + expansion_width, ) # remove the grown regions into the void from the dilation and the external_id values diff --git a/setup.py b/setup.py index 471c5b3..742b69d 100644 --- a/setup.py +++ b/setup.py @@ -25,6 +25,7 @@ "numpy>=1.15.0", "scipy>=1.10.0", "voxcell>=3.0.0", + "joblib>=1.3" ], extras_require={ "tests": [ From cbf2d348789a3d53c4671652e96bc242203c7d2d Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 29 Feb 2024 11:34:39 +0100 Subject: [PATCH 4/5] fix --- atlas_direction_vectors/algorithms/utils.py | 7 +++---- setup.py | 3 ++- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/atlas_direction_vectors/algorithms/utils.py b/atlas_direction_vectors/algorithms/utils.py index 421ba37..f83ff46 100644 --- a/atlas_direction_vectors/algorithms/utils.py +++ b/atlas_direction_vectors/algorithms/utils.py @@ -4,9 +4,12 @@ from typing import Union import numpy as np # type: ignore +import pyquaternion as pyq from atlas_commons.utils import FloatArray, NumericArray, normalize +from joblib import Parallel, delayed from scipy.ndimage import gaussian_filter # type: ignore from scipy.ndimage import generate_binary_structure # type: ignore +from scipy.optimize import minimize from scipy.signal import correlate # type: ignore L = logging.getLogger(__name__) @@ -75,10 +78,6 @@ def _quaternion_from_vectors( # pylint: disable=invalid-name if align_to is not None: target_dir = {"x": np.array([1.0, 0.0, 0.0]), "z": np.array([0.0, 0.0, 1.0])}[align_to] - import pyquaternion as pyq - from joblib import Parallel, delayed - from scipy.optimize import minimize - def rotate_q(_q): q = pyq.Quaternion(_q) diff --git a/setup.py b/setup.py index 742b69d..65215cd 100644 --- a/setup.py +++ b/setup.py @@ -25,7 +25,8 @@ "numpy>=1.15.0", "scipy>=1.10.0", "voxcell>=3.0.0", - "joblib>=1.3" + "joblib>=1.3", + "pyquaternion>=0.9.9", ], extras_require={ "tests": [ From 0c3e742559d29490dce619aab1a702359174da3a Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 29 Feb 2024 18:27:52 +0100 Subject: [PATCH 5/5] fix test --- ...ebellum_shading_gradient_orientations.nrrd | Bin 1232 -> 1023 bytes tests/test_cerebellum.py | 9 ++++----- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/tests/data/cerebellum_shading_gradient_orientations.nrrd b/tests/data/cerebellum_shading_gradient_orientations.nrrd index caf1faf064d06283730d79bf62aaddb66715eb61..bef222f4d59ff3dc8c6641a802f94c90323c55e7 100644 GIT binary patch delta 647 zcmV;20(kw<3I7L>IuJB1FfuJNIUq4NIx;spGB>eGP62;CXjDNIh1WTb5)Ez{ak&Y%=`%$>Qw^Uv@;Nr~*=c>tn_=qs z${>^$JHt~;OucE_0QdW6Q~#>>Kv8|Fhr#IY?a_bj+!gig?4RJ>iIJB7>Ry$pCzqz; z`h*ugO}IA(dgGf5KT^3p`-{q}jQiz|9npNkZF_t=*xO^wC5THjxMjrUCfuga6!}aw zeCEh!ZsN1k`c~x&QzO}v&|cUXzFKAq8o%OgyQ-;w_29#hYX8|4bUyMUx_VLNE;ND~ zec6BVUyUZ*zOAq0gVVdiu7tbq4vo(Yd`snSnix`tzdrP{+>e$ME_r)281s#}1aXN5 zw~V;lgxmC)BAZK zUultFf$|F)`IRWY(*3McN`DrpKZ8bpmZ)n#ON;(2P=5xE{wz^{mKOb4p#BUR{aK>^ zEG_!8K>Zms`m;p++2z+MZNCDvUxCJcCD48)E%qxw`xR*HR|4%<(qg{?v|oY7ekIU; hB`x+VK>HPF>{kNqS1!L!PEKoH{{U{Xx|#qU004G2LmU7A delta 858 zcmey*et~m>m4KnHp^>hEnS!B-m4T_1fx$*EKgRm#jPr>S$Nxtw8LEkUZ)!0#&N+4A zw_U#8lM2OlnPpFR^vImP_2xtU^9h-ocP1)d5%aawQQdvx_rCq}c$Dzwko@TSCZ zu|%5q?jQGd+<0H!djDwg-ml&Ped`s-KEbsX{)0WO^ z*O{+ue?4sH?dSET)8G3~d@r7|Gy1#fc9EW`=Fe4kO0PG)p3XbH@TqOs_RPFnb{^N3 z+HYZg+P9`-)2g+?{nzcvdu#8A&%K)DpT4!a$g%kN6E4Gw?{43pHa=UlS~V%=kI46j zz0Z#H$ZYJZznT4ZmVVu7k(0*qKc7lHfJ6H`|K{jE$*EUN=eCSKantGCHUD!D|K2WN zxp~6(wX@XQub%m4chqHPN>-PR)xYVpuZTYV87`Fi>Xe)K&2^QD-h1Ax(=V>p%U-uj zTw6V5UCHuynd|KSh0gKzDrA0|R~>y~WB9qQW!HG#f9X$(lZaMS>+P&aeH-}hZj}p< zZ0_C5IcI}seeY`dQZM|jdhWA5UthOe-+k-px5CSDy26*Ug8mlq7QkSyd|v%Iv7&7g zHn)OYteE~-{rQJ~Z(pA{_x$xz;nY0GdFADSg5~QyCM;g(_uVRK@tsl~zj+?-esAm# z)Qc5Uo^2MjJ#q50!z*)>Jd>72*SHlQTQ~ht?%{2W|3e?vf2y&Ijnlg1Ci8`_Q)%Bj zE|+uXCGY*6mt?o0a`_H3^GBuox<7rNVVV>&?fIR;P7r{wr56Zbtkbu z{edM=p<(K9?ieUc&HrhA2yfia@=xf8_XB-a5S!~y=mUL_P*6QXz37kC2jW>#WK};* zKM>DUuleWZhkqC5i+)LO`sZ_DKIboE$NFUqeu@P zWcrfM`gi66dts{z;7Aw445RCG59Ty0M<)n>Od~ ssXPb;K((GX)~7Up6hP#g{(2p%Q$0MN7nq#?ef!NjskAeZL5_g|0NfJ0%m4rY diff --git a/tests/test_cerebellum.py b/tests/test_cerebellum.py index 6f6feb8..5637543 100644 --- a/tests/test_cerebellum.py +++ b/tests/test_cerebellum.py @@ -183,17 +183,16 @@ def _check_vectors_direction_dominance( def test_compute_cerebellum_direction_vectors(region_map, annotation): res = tested.compute_direction_vectors(region_map, annotation) _check_vectors_defined_in_regions( - res, region_map, annotation, ["FLgr", "FLpu", "FLmo"] + ["LINGgr", "LINGpu", "LINGmo"] + res, region_map, annotation, ["FLgr", "FLmo"] + ["LINGgr", "LINGmo"] ) _check_vectors_direction_dominance( - res, region_map, annotation, ["FLgr", "FLpu", "FLmo"], [-1.0, 0.0, 0.0] + res, region_map, annotation, ["FLgr", "FLmo"], [-1.0, 0.0, 0.0] ) _check_vectors_direction_dominance( - res, region_map, annotation, ["LINGgr", "LINGpu", "LINGmo"], [-1.0, 0.0, 0.0] + res, region_map, annotation, ["LINGgr", "LINGmo"], [-1.0, 0.0, 0.0] ) - expected_direction_vectors = VoxelData.load_nrrd( DATA_DIR / "cerebellum_shading_gradient_orientations.nrrd" ) @@ -205,4 +204,4 @@ def test_cereb_subreg_direction_vectors(region_map, annotation): res = tested.cereb_subregion_direction_vectors( region_map.find("FL", "acronym").pop(), region_map, annotation ) - _check_vectors_defined_in_regions(res, region_map, annotation, ["FLgr", "FLpu", "FLmo"]) + _check_vectors_defined_in_regions(res, region_map, annotation, ["FLgr", "FLmo"])