.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/07-tours/06-navigator-motion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_07-tours_06-navigator-motion.py: ================================= Rigid head motion from navigators ================================= In this tour you measure how a head moves during a scan, from navigators. If the head moves, the anatomy is in a different place for the readouts that encode it, and the image is blurred and ghosted. Prospective correction moves the imaging field of view with the head before each readout, and for that it needs the rigid pose of the head, six degrees of freedom, measured again and again during the scan. A navigator is a short, low-resolution acquisition interleaved with the imaging readouts; registering it against the first navigator of the scan tells you how the head has moved since [#ehman]_. A 2D navigator plane, registered against its reference, measures three degrees of freedom: two shifts in the plane and one rotation about its normal. Three orthogonal planes give nine measurements of the six degrees of freedom of the pose, which a least-squares fit combines. Over a scan, a Kalman filter smooths the poses, as in PROMO [#promo]_, which tracks the head with three orthogonal spiral navigators. Prerequisites: lesson 7, :doc:`../04-non-cartesian/01-trajectories-and-transforms`. **Learning objectives** - Relate each navigator plane to the three degrees of freedom it measures, and the rigid pose to the nine measurements of three orthogonal planes. - Reconstruct radial navigator planes with :func:`~bartorch.tools.reconstruct_navigator` and measure the pose with :class:`~bartorch.tools.NavigatorMotionTracker`. - Check a measured pose by comparing the navigator of the moved head with the navigator of the head moved by that pose. - Set the filter's process and measurement noise against the precision of the navigator and the dynamics of the motion. .. GENERATED FROM PYTHON SOURCE LINES 40-85 .. GENERATED FROM PYTHON SOURCE LINES 86-107 The head -------- Your head is the BrainWeb T1-weighted volume [#brainweb]_, downsampled to 2 mm isotropic and padded to :math:`112^3`. Array axis 0 runs inferior to superior (S/I), axis 1 posterior to anterior (A/P) and axis 2 left to right (L/R). You state poses in that frame: a rotation vector in radians about the centre of the volume, printed in degrees, and a translation in millimetres. ``move(rotvec, translation_mm)``, hidden on this page, returns the volume after such a motion. The navigator ------------- Each navigator is three central planes, axial, coronal and sagittal, each acquired along 96 golden-angle radial spokes of 56 samples. The trajectory is in grid units of a :math:`56^2` matrix, so the spokes reach half the imaging resolution and you reconstruct each plane at 4 mm over the full field of view. :func:`~bartorch.tools.reconstruct_navigator` grids a plane by the density-compensated adjoint NUFFT, with weights from :func:`bartorch.estimate_density`. .. GENERATED FROM PYTHON SOURCE LINES 107-123 .. code-block:: Python import numpy as np import torch import bartorch import bartorch.tools as bt from scipy.spatial.transform import Rotation NAV, SPOKES = 56, 96 golden = np.pi * (3 - np.sqrt(5)) angle = np.arange(SPOKES)[:, None] * golden radius = np.arange(NAV) - NAV / 2 spokes = np.stack([np.cos(angle) * radius, np.sin(angle) * radius], axis=-1) trajectory = torch.tensor(spokes.reshape(-1, 2), dtype=torch.float32) density = bartorch.estimate_density(trajectory, (NAV, NAV)) .. GENERATED FROM PYTHON SOURCE LINES 124-127 ``AXES`` tells the tracker, for each plane, which direction in the volume its image's rows and columns run along. That is what relates each plane's in-plane measurement to the 3D pose. .. GENERATED FROM PYTHON SOURCE LINES 127-171 .. code-block:: Python AXES = [ ((0, 1, 0), (0, 0, 1)), # axial: rows A/P, columns L/R ((1, 0, 0), (0, 0, 1)), # coronal: rows S/I, columns L/R ((1, 0, 0), (0, 1, 0)), # sagittal: rows S/I, columns A/P ] SPACING_MM = VOXEL_MM * SIZE / NAV def navigator(head, noise=0.0, generator=None): """The three navigator planes of ``head``, with k-space noise.""" middle = SIZE // 2 traj = torch.nn.functional.pad(trajectory, (0, 1))[None] planes = [] for plane in (head[middle], head[:, middle], head[:, :, middle]): plane = torch.tensor(plane, dtype=torch.complex64) samples = bartorch.nufft(plane, traj).reshape(1, -1) if noise: samples = samples + noise * torch.randn( samples.shape, dtype=torch.complex64, generator=generator, ) nav = bt.reconstruct_navigator( samples, trajectory[None], (NAV, NAV), density=density, ) planes.append(nav[0]) return planes reference = navigator(volume) .. image-sg:: /auto_examples/07-tours/images/sphx_glr_06-navigator-motion_001.png :alt: axial, coronal, sagittal :srcset: /auto_examples/07-tours/images/sphx_glr_06-navigator-motion_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 172-186 These are the three reference planes at 4 mm resolution. Rigid registration relies on the outline of the head and the large structures, which this resolution retains. Pose estimation --------------- :class:`~bartorch.tools.NavigatorMotionTracker` takes the first navigator you give it as the reference and returns the pose of every later one relative to it. It registers each plane in 2D with :class:`~bartorch.tools.RigidRegistration`, and solves for the pose that fits the three planes in the least-squares sense. You move the head by a few degrees and millimetres along every axis and see what the tracker measures. .. GENERATED FROM PYTHON SOURCE LINES 186-210 .. code-block:: Python rotvec = np.array([0.02, -0.03, 0.05]) # rad translation = np.array([3.0, -2.0, 4.0]) # mm tracker = bt.NavigatorMotionTracker(measurement_noise=1e-3) tracker.track(reference, AXES, spacing=SPACING_MM) moved = navigator(move(rotvec, translation)) pose = tracker.track(moved, AXES, spacing=SPACING_MM) rotation = Rotation.from_matrix(np.asarray(pose.matrix)[:3, :3]) measured_rotvec = rotation.as_rotvec() measured_translation = np.asarray(pose.translation) .. rst-class:: sphx-glr-script-out .. code-block:: none truth measured rotation S/I 1.15° 1.46° rotation A/P -1.72° -1.45° rotation L/R 2.86° 2.51° shift S/I 3.0 mm 3.37 mm shift A/P -2.0 mm -2.66 mm shift L/R 4.0 mm 3.65 mm .. GENERATED FROM PYTHON SOURCE LINES 211-216 In a real scan the true pose is unknown, so you check a measured pose on the navigator itself: move the head by the measured pose, acquire its navigator, and compare that with the navigator after the real motion. Without any correction, the difference is that of the motion; with the measured pose, what is left is the error of the pose. .. GENERATED FROM PYTHON SOURCE LINES 216-243 .. code-block:: Python realigned = navigator(move(measured_rotvec, measured_translation)) .. image-sg:: /auto_examples/07-tours/images/sphx_glr_06-navigator-motion_002.png :alt: axial, coronal, sagittal :srcset: /auto_examples/07-tours/images/sphx_glr_06-navigator-motion_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none axial: NRMSE against the moved navigator, reference 0.335, measured pose 0.051 coronal: NRMSE against the moved navigator, reference 0.332, measured pose 0.049 sagittal: NRMSE against the moved navigator, reference 0.421, measured pose 0.067 .. GENERATED FROM PYTHON SOURCE LINES 244-260 The top row is the motion as the navigator sees it: strong edges all around the head. In the bottom row, after the measured pose, only a faint difference at the edges of the head is left. The pose is off by a fraction of a degree and of a millimetre, a small fraction of the 4 mm navigator voxel. That error has two sources: the coarse resolution of the navigator, and the through-plane motion, which each plane sees as a change of the anatomy in it rather than as a rigid motion within it. Motion over a scan ------------------ Over a scan, the head nods, a rotation about the left-right axis, by up to 2.3° with a period of 8 s, and drifts by 3 mm along the superior-inferior axis over 12 s. You acquire a navigator every 0.5 s, and each one carries complex Gaussian k-space noise, so each measured pose carries registration error. .. GENERATED FROM PYTHON SOURCE LINES 260-273 .. code-block:: Python DT, COUNT = 0.5, 24 seconds = DT * np.arange(1, COUNT + 1) nod = 0.04 * np.sin(2 * np.pi * seconds / 8) # rad, about L/R drift = 3.0 * seconds / seconds[-1] # mm, along S/I truth = np.zeros((COUNT, 6)) truth[:, 2], truth[:, 3] = nod, drift generator = torch.Generator().manual_seed(0) scan = [] for a, d in zip(nod, drift): head = move((0.0, 0.0, a), (d, 0.0, 0.0)) scan.append(navigator(head, noise=0.5, generator=generator)) .. GENERATED FROM PYTHON SOURCE LINES 274-280 The tracker filters the poses with a constant-velocity extended Kalman filter, and two numbers set it. ``process_noise`` is the variance of the acceleration the filter allows between navigators, per pose coordinate, in rad²/s⁴ and mm²/s⁴. ``measurement_noise`` is the variance it assigns to each measured coordinate. You filter the same scan with three values of ``process_noise``; the largest one leaves the filter at the measurements. .. GENERATED FROM PYTHON SOURCE LINES 280-325 .. code-block:: Python def track(process_noise): tracker = bt.NavigatorMotionTracker( process_noise=process_noise, measurement_noise=1e-4, ) tracker.track(reference, AXES, spacing=SPACING_MM) poses = [] for planes in scan: poses.append(tracker.track(planes, AXES, dt=DT, spacing=SPACING_MM)) return poses traces = {noise: track(noise) for noise in (1e2, 1e-3, 1e-4)} .. image-sg:: /auto_examples/07-tours/images/sphx_glr_06-navigator-motion_003.png :alt: nod, about L/R, drift, along S/I :srcset: /auto_examples/07-tours/images/sphx_glr_06-navigator-motion_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none process_noise 100 rms error: rotation S/I 0.47°, A/P 0.36°, L/R 0.40°; shift S/I 0.19 mm, A/P 0.42 mm, L/R 0.19 mm process_noise 0.001 rms error: rotation S/I 0.44°, A/P 0.29°, L/R 0.36°; shift S/I 0.16 mm, A/P 0.41 mm, L/R 0.15 mm process_noise 0.0001 rms error: rotation S/I 0.43°, A/P 0.24°, L/R 0.60°; shift S/I 0.15 mm, A/P 0.42 mm, L/R 0.13 mm .. GENERATED FROM PYTHON SOURCE LINES 326-354 A large ``process_noise`` leaves the filter at the measurements (the dots). A small one makes it trust its constant-velocity prediction, which smooths the registration error but lags a change of direction: at ``1e-4`` the error of the nod grows while that of the other coordinates falls. An error common to every navigator, such as the offset of the A/P shift, is a bias of the measurement, and no value of ``process_noise`` reduces it. The value that gives the smallest error depends on the motion and on the navigator's precision, so you set it against a motion trace of your own application. As a spec --------- What this tour built, stated the way you would ask an agent for it: .. code-block:: text Track rigid head motion from navigators with bartorch. Each navigator is three orthogonal central planes (axial, coronal, sagittal), each 96 golden-angle radial spokes of 56 samples, reconstructed on a 56^2 grid with bt.reconstruct_navigator and bartorch.estimate_density weights. Pass the plane orientations as row and column directions in the volume (S/I, A/P, L/R) and the pixel spacing in mm to bt.NavigatorMotionTracker; the first navigator is the reference. Report each pose as a rotation vector and a translation in mm, check it by re-navigating the head moved by the measured pose, and filter a 24-navigator scan (one every 0.5 s) with process_noise 1e2, 1e-3 and 1e-4, reporting the RMS error per coordinate against the true motion. .. GENERATED FROM PYTHON SOURCE LINES 356-372 References ---------- .. [#ehman] Ehman RL, Felmlee JP. Adaptive technique for high-definition MR imaging of moving structures. *Radiology* 173(1):255-263 (1989). https://doi.org/10.1148/radiology.173.1.2781017 .. [#promo] White N, Roddey C, Shankaranarayanan A, Han E, Rettmann D, Santos J, Kuperman J, Dale A. PROMO: Real-time prospective motion correction in MRI using image-based tracking. *Magn Reson Med* 63(1):91-105 (2010). https://doi.org/10.1002/mrm.22176 .. [#brainweb] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, Evans AC. Design and construction of a realistic digital brain phantom. *IEEE Trans Med Imaging* 17(3):463-468 (1998). https://doi.org/10.1109/42.712135 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.928 seconds) .. _sphx_glr_download_auto_examples_07-tours_06-navigator-motion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 06-navigator-motion.ipynb <06-navigator-motion.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 06-navigator-motion.py <06-navigator-motion.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 06-navigator-motion.zip <06-navigator-motion.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_