# Data-Driven Latent Space Representation for Robust Bipedal Locomotion
Learning

Guillermo A. Castillo1, Bowen Weng1, Wei Zhang2, and Ayonga Hereid3 This work
was supported in part by the National Science Foundation under grant
FRR-21441568.1Electrical and Computer Engineering, Ohio State University,
Columbus, OH, USA<EMAIL_ADDRESS>Institute of
Robotics, Southern University of Science and Technology (SUSTech), China;
<EMAIL_ADDRESS>and Aerospace Engineering, Ohio State
University, Columbus, OH, USA<EMAIL_ADDRESS>
###### Abstract

This paper presents a novel framework for learning robust bipedal walking by
combining a data-driven state representation with a Reinforcement Learning
(RL) based locomotion policy. The framework utilizes an autoencoder to learn a
low-dimensional latent space that captures the complex dynamics of bipedal
locomotion from existing locomotion data. This reduced dimensional state
representation is then used as states for training a robust RL-based gait
policy, eliminating the need for heuristic state selections or the use of
template models for gait planning. The results demonstrate that the learned
latent variables are disentangled and directly correspond to different gaits
or speeds, such as moving forward, backward, or walking in place. Compared to
traditional template model-based approaches, our framework exhibits superior
performance and robustness in simulation. The trained policy effectively
tracks a wide range of walking speeds and demonstrates good generalization
capabilities to unseen scenarios.

##  I Introduction

Bipedal robots have long been a subject of fascination and research within the
field of robotics due to their potential for versatile and agile locomotion in
complex environments, closely mimicking the morphology and capabilities of
humans. However, despite significant advancements in control techniques and
hardware, developing a robust and efficient bipedal locomotion control system
remains a challenge, largely due to the high dimensionality, underactuation,
and highly nonlinear and hybrid dynamics of bipedal locomotion.

Conventional methods for bipedal walking often involve solving optimization
problems using the robot’s full-order [1] or reduced-order model [2, 3] to
find feasible trajectories that enable stable walking gaits. The full-order
model captures all the complexities and details of the robot’s dynamics, but
it can be computationally demanding and not suitable for real-time control
[4]. On the other hand, reduced-order template models (such as linear inverted
pendulum (LIP) and its variants [2, 5, 6]) simplify the dynamics of the
system, making it easier to plan trajectories for the robot’s center of mass
and end-effector. However, these reduced-order models often require strict
constraints to account for the mismatch between the reduced and full-order
states of the robot.

Figure 1: An overview of the overall structure and flow of the proposed
learning-based framework. In the pre-training phase, the autoencoder learns a
latent space that captures the dynamics of the full-order system. During the
RL training, the policy maps the latent representation to a set of task space
actions that are translated into task space trajectories. Finally, a whole-
body task space feedback controller computes the motor torque to track the
desired task-space trajectories.

The evolution of modern control theory has seen an infusion of machine
learning and RL techniques, particularly with the growing abundance and
accessibility of data. These data-driven approaches offer novel ways to
address challenges in control design, often sidestepping traditional, more
rigid, and computationally expensive methods. There is a growing interest in
using reinforcement learning based approaches that allow exploiting data from
simulation to train controllers in a model-free fashion [7, 8, 9, 10].
However, similar to model-based techniques, they are highly dependent on the
quality of the state representation provided to the learning algorithm. As an
alternative, more complex frameworks have been proposed to combine learning
algorithms with model-based controllers. In [11], an HZD-based approach is
used to learn a policy that satisfies Control Barrier Functions (CBF) defined
on the reduced-order dynamics. In our previous work [12, 13], a cascade
structure is implemented to compensate the learned trajectories with feedback
regulators to increase the robustness of the walking gait. A hierarchical
structure that combines a template-based RL policy with a model-based low-
level controller is proposed in [14].

An effective state representation that can accurately capture the complex
dynamics of the whole system can significantly enhance the learning process,
enabling more efficient learning and better transferability of control
strategies. Unsupervised learning, dimensionality reduction, and
representation learning methods can be employed to extract relevant features
from high-dimensional sensory data, enabling the development of more efficient
and interpretable state representations [15, 16]. Dai et al. [17] learn the
step-to-step residual dynamics via an adaptive control approach, which is then
used to design a foot-stepping controller for a bipedal robot in simulation.
More complex learned residual dynamic models are also combined with Model
Predictive Control (MPC) for agile systems [18]. These techniques require
previous knowledge of the dynamics of the system’s model. Several works also
have exploited latent representations through reinforcement learning for
locomotion. Peng et al. [19] combine techniques from adversarial imitation
learning and unsupervised reinforcement learning to develop skill embeddings
that produce locomotion behaviors. Starke et al. [20] extract a multi-
dimensional phase space from the full-order state motion data, which
effectively clusters animations and produces a manifold with better temporal
and spatial alignment. However, these are end-to-end frameworks, which makes
it difficult to establish a relationship between the latent space and the
control actions of the policy. Moreover, these approaches have been mostly
used to control animated characters in simulation, and there are no studies
focused on implementation for actual bipedal robots.

In this paper, we propose a novel data-driven framework for bipedal walking
that combines a learned low dimensional state representation of bipedal
locomotion with a robust gait planner using RL, as depicted in Fig. 1. Our
framework uses an autoencoder to extract an effective reduced-dimensional
state representation of the full-order system dynamics. We then integrate this
reduced-order latent space with reinforcement learning and a task space
feedback controller to train robust locomotion policies. This paper makes two
key contributions. First, we demonstrate that the complex dynamics of bipedal
robots can be effectively captured using a low-dimensional learned latent
space. This allows for a more compact representation of the system’s behavior.
Second, we show that the learned latent variables can be leveraged to design a
robust locomotion policy using RL. By bridging the gap between state
representation learning and learning-based control policies, this work enables
leveraging existing locomotion data to develop more effective and adaptable
frameworks for versatile and robust bipedal locomotion.

##  II Preliminaries and Problem Formulation

###  II-A Hybrid System Model of Bipedal Locomotion

The bipedal locomotion problem can be characterized as a hybrid system
determined by a collection of phases of continuous dynamics with discrete
events between the transitions of the continuous phases. Formally, the hybrid
system model for biped locomotion can be defined as [21]:

$\displaystyle\Sigma:\left\\{\begin{array}[]{lcl}\dot{x}=f(x)+g(x)u+\mathbf{\omega}(x,u)&&x\in\mathcal{X}\setminus\mathcal{H}\\\
x^{+}=\Delta(x^{-})&&x^{-}\in\mathcal{H},\end{array}\right.$ (3)

where $f(\cdot)$ and $g(\cdot)$ are vector fields,
$x=(q,\dot{q})\in\mathcal{X}\subseteq\mathbb{R}^{n}$ represents the robot
states with $q$ being the vector of generalized coordinates,
$u\in\mathcal{U}\subseteq\mathbb{R}^{m}$ is a vector of actuator inputs, and
$\mathbf{\omega}\in\Omega\subseteq\mathbb{R}^{w}$ captures external
disturbances and model uncertainties. The reset map
$\Delta:\mathcal{H}\to\mathcal{X}$ denotes the mapping between the post-impact
states $x^{+}$ immediately after impacts and the pre-impact states $x^{-}$
right before impacts, and $\mathcal{H}$ is the switching surface corresponding
to swing foot impacts.

For a typical humanoid robot with high degrees of freedom (DoF), the dimension
of the robot states $x$ is too large to be effectively used for feedback
motion planning. Low-dimensional models have become a powerful tool for motion
planning of bipedal locomotion, given their potential to characterize the
dynamics of bipedal walking into simple linear or nonlinear models. However,
existing reduced-order template models, such as the Linear Inverted Pendulum
(LIP), make assumptions that limit the full capabilities of walking robots.
These assumptions include a constant center of mass (CoM) height and zero
angular momentum about the CoM during the walking gait. While previous
research has explored alternative template models for bipedal locomotion, this
study aims to investigate whether available locomotion data can be utilized to
identify an effective reduced-dimensional state representation of bipedal
locomotion for motion planning purposes.

###  II-B Data-driven Low Dimensional Latent Space

Autoencoders are a great tool to harness high dimensional gait data to extract
a reduced dimensional latent representation of the system that captures the
essence of the full-order robot dynamics. An autoencoder works by compressing
the input data into a low-dimensional latent space through a feature-
extracting function in a specific parameterized closed form, such as neural
networks [16]. This function, called encoder is determined by

$\displaystyle z=h(x,\theta_{e}),$ (4)

where $z$ is the latent variable or representation encoded from the input $x$,
and $\theta_{e}$ is the vector of parameters for the encoder neural network.
Another closed-form parameterized function called the decoder, maps from the
latent space back to the full-order state. This function is defined by

$\displaystyle\hat{x}=d(z,\theta_{d}),$ (5)

where $z$ is the encoded latent variable, $\hat{x}$ is the reconstruction of
the original input data $x$, and $\theta_{d}$ is the vector of parameters for
the decoder neural network.

The strength of autoencoders lies in their ability to preserve most of the
crucial information from the original data in the latent representation, even
though it is of much lower dimensionality. This is ensured by training the
autoencoder to minimize the reconstruction loss $\mathcal{L}$, which measures
the error between the original data and its reconstruction from the latent
space. In summary, autoencoder training consists of finding a value of the
parameter vectors $\theta_{e},\theta_{d}$ that minimize the reconstruction
error:

$\displaystyle\mathcal{J}_{\mathrm{AE}}(\theta_{e},\theta_{d})=\sum_{x^{(i)}\in\mathbf{X}}\mathcal{L}\left(x^{(i)},d\left(h\left(x^{(i)},\theta_{e}\right),\theta_{d}\right)\right),$
(6)

where $x^{(i)}$ is a training sample containing the vectors of position and
velocity of the generalized coordinates, and $\mathbf{X}$ is the training set
containing all the data samples. $h$ and $d$ are the encoding and decoding
functions introduced in (4) and (5). Once the autoencoder has been trained,
the resulting latent representation can then be utilized as part of the state
for the high-level planner policy, offering a compact yet expressive state
space for learning and control.

##  III Method

The proposed data-driven framework is shown in Fig. 1. The gait policy will be
learned via RL by utilizing an effective latent representation of the full-
order robot’s dynamics to a set of task space commands that generate desired
trajectories for the robot. Then, a whole body task space controller (TSC)
from [14] is employed to accurately track the desired trajectories to realize
stable locomotion.

The proposed framework is tested with the robot Digit, which is a 3D fully
actuated bipedal robot with $30$ DoF and $20$ actuated joints built by the
company Agility Robotics. Each leg has six actuated joints corresponding to
the motors located on the robot’s hip, knee, and ankle and three passive
joints corresponding to the robot’s tarsus, shin spring and heel spring
joints. In addition, it has four actuated joints per arm corresponding to the
shoulder and elbow joints. Since the spring joints are very stiff, we
considered them as fixed joints in this work. Therefore, the vector of
generalized coordinates for Digit is defined by

$\displaystyle q=(p,q_{\phi},q_{j}),$ (7)

where $p=(p_{x},p_{y},p_{z})$ is the position of the robot’s base,
$q_{\phi}=(q_{x},q_{y},q_{z},q_{w})$ is the quaternion representation of the
robot’s orientation, and $q_{j}=(q_{1},\dots,q_{n_{j}})$ is the vector of the
robot’s joints with $n_{j}=24$. Therefore, $q\in
SE(3)\times\mathbb{R}^{24}\subset\mathbb{R}^{31}$,
$\dot{q}\in\mathbb{R}^{30}$, and $x\in\mathbb{R}^{61}$.

###  III-A Reduced Dimensional Latent Space Representation

To collect the locomotion dataset to train the autoencoder, we use the
hierarchical controller proposed in [14]. The dataset is collected by
performing walking gaits at various velocities. Specifically, the velocities
$v_{x}$ and $v_{y}$ are varied within the ranges of $[-0.5,1.0]$ m/s and
$[-0.2,0.2]$ m/s, respectively, with a step size of $0.1$ m/s, resulting in a
total of $16$ different walking gaits. Each walking gait has a duration of
$10$ seconds, with the data being collected at a frequency of $50$ Hz.
Therefore, the complete locomotion dataset $\mathbf{X}$ consists of $40,000$
samples of the robot’s full-order states, i.e.,
$\mathbf{X}=\\{x^{(i)}|i\in[1,40000]\\}$.

Remark:  Note that any locomotion controller could be used to collect the gait
data. For instance, many commercial robots are equipped with proprietary
controllers that could be used to collect this data even though they are a
black box from the user’s perspective. Moreover, we do not make any
assumptions about the distribution of the gait data.

In this work, we use an encoder parameterized by a fully connected neural
network with three hidden layers of $128$, $64$, and $32$ units, respectively,
and ReLU activation functions. The input of the encoder is the robot’s full
order states. However, since the inputs of the encoder need to be bounded, the
absolute base position of the robot cannot be directly used. This is because
the absolute position can grow unbounded as the robot moves in the sagittal or
frontal plane. To address this, the base position and orientation of the robot
are transformed from the world frame to the stance foot frame, which ensures
that the inputs to the encoder remain bounded and allows for effective
learning of the latent representation of the robot’s dynamics.

The autoencoder is trained using Adam optimizer [22] with a learning rate of
$0.001$ and a batch size $B$ of $128$. The reconstruction loss is computed
with the mean squared error (MSE) between the original values of the gait
dataset $x^{(i)}\in\mathbf{X}$ and their reconstructed values
$\hat{x}^{(i)}\in\mathbf{\hat{X}}$, given as

$\displaystyle\mathcal{L}=\frac{1}{B}\left(x^{(i)}-\hat{x}^{(i)}\right)^{2}.$
(8)

The autoencoder is trained for 400 episodes in a 12-core CPU machine with an
NVIDIA RTX 2080 GPU. The training takes about 10 minutes using the described
locomotion dataset $\mathbf{X}$ and PyTorch.

The selection of the dimension of the latent variable $z$ involves a trade-
off. On one hand, a smaller dimension is desirable as it reduces the number of
inputs for the RL policy. However, the dimension should also be large enough
for the autoencoder to accurately reconstruct the full-order state of the
system. In our study, we investigated this trade-off and discovered that even
for a complex system like a humanoid robot, a latent variable dimension of
$N=2$ is sufficient to capture the dynamics of the full-order systems
effectively. This finding is supported by the results shown in Fig. 2, where
increasing the dimension of the latent variable does not lead to significant
improvements in reconstruction quality. Even in cases where performance is
slightly degraded, the reconstructed variable still follows the same pattern
as the original state. Interestingly, our findings align with those reported
in [20], where the authors demonstrated that periodic movements in bipedal
locomotion can be represented using fewer than five phase variables. This
similarity may be attributed to the symmetric and periodic nature of bipedal
walking.

Figure 2: Reconstruction of the robot’s state with different dimensions $N$
for the latent variable $z$. The plots show the position (left column) and
velocity (right column) for the robot’s base $x$ coordinate (top), robot’s
base $y$ coordinate (middle), and left knee joint (bottom), respectively.

###  III-B RL-based Gait Policy using Learned Latent Variables

Given the reduced dimensional latent space, we train an RL policy for robust
locomotion based on our previous work in [14]. The states of the policy are
defined as

$\displaystyle s=(z,e_{\bar{v}},v^{d},a_{k-1}),$ (9)

where $z=(z_{1},\ldots,z_{N})\in\mathbb{R}^{N}$ is the encoded latent state,
$e_{\bar{v}}=(e_{\bar{v}_{x}},e_{\bar{v}_{y}})$ is the error between the
average velocity, $\bar{v}=(\bar{v}_{x},\bar{v}_{y})$, and the desired
velocity, $v^{d}=(v_{x}^{d},v_{x}^{d})$, of the robot, and $a_{k-1}$ is the
last action of the planner policy.

The selection of the action space in this work is designed to exploit the
natural nonlinear dynamics of the biped robot and enhance the robustness of
the policy under various challenging scenarios. Since Digit is a fully
actuated system during the single support phase, we include the instantaneous
base velocity as part of the action $a$. We choose to control instantaneous
velocity over the base position because controlling the position during the
stance phase through the TSC is more challenging than velocities in practical
hardware implementation. This is caused by the noisy and inaccurate estimation
of the position and the limited amount of torque in the robot’s ankles.
Moreover, sudden motions of the foot position due to terrain irregularities
could produce big errors in the base position tracking that could result in
aggressive control maneuvers. Therefore, controlling the instantaneous
velocity is a less aggressive strategy for the TSC and provides some degree of
damping to the ankle motion that helps with the stability of the walking gait
under irregular terrains. Thus the action $a\in\mathcal{A}$ of the policy for
Digit is chosen to be:

$\displaystyle
a=(p^{x}_{\text{sw},T},p^{y}_{\text{sw},T},v_{x}^{d},v_{y}^{d}),$ (10)

which corresponds to the landing position of the swing foot in $x$ and $y$
coordinates with respect to the robot’s base and an offset to the
instantaneous velocity of the robot’s base in $x$ and $y$ coordinates, as
illustrated in Fig. 3. The trajectory generation module transforms the policy
action $a$ into smooth task-space trajectories for the robot’s base and end-
effectors. Specifically, the trajectory for the relative swing foot positions,
$p^{x}_{\text{sw}}$ and $p^{y}_{\text{sw}}$, are generated using a minimum
jerk trajectory connecting initial foot positions with target foot positions
from the policy action. In particular, the initial foot positions will be
computed at every touch-down event and kept constant throughout the current
step.

Figure 3: Visualization of the policy actions and the trajectories generated
for the task space controller. Additionally, we keep the torso straight up and
the base at a constant height.

The neural network (NN) chosen to parameterize the gait policy is a feed-
forward network with two hidden layers, each layer with $128$ units. The
hidden layers use the ReLU activation function, and the output layer is
bounded by the $\mathrm{Tanh}$ activation function and a scaling factor to
constrain the maximum value of the policy commands within the feasible
physical limits of the robot hardware. Finally, we implemented a model-based
TSC controller following the structure to track the task space trajectories,
as described in our previous work [14].

Figure 4: Two-dimensional principal component analysis (PCA) of the learned
latent manifolds at different walking speeds.

###  III-C Learning Procedure of the RL Gait Policy

To train the RL policy, we use the Proximal Policy Optimization algorithm [23]
with input normalization, fixed covariance, and parallel experience
collection. We follow the algorithm implementation described in [8].

For each episode during the RL training, the initial state of the robot is set
randomly from a normal distribution about an initial pose corresponding to the
robot standing in the double support phase. An episode will be terminated
early if the torso pitch and roll angles exceed $1$ rad or the base height
falls below $0.8$ m. The reward function adopted in this work is designed to
(1) keep track of desired velocities, (2) minimize the angular momentum around
the center of mass, denoted as $L_{\textrm{CoM}}$, and (3) reduce the
variation of the policy actions between each iteration. More specifically,

$\displaystyle\mathbf{r}=\mathbf{w}^{T}[r_{v},r_{L_{\textrm{CoM}}},r_{a}]^{T},$
(11)

with

$\displaystyle r_{v_{x}}$
$\displaystyle=\exp{({-\left\lVert\bar{v}-v^{d}\right\rVert}^{2})},$ (12)
$\displaystyle r_{L_{\textrm{CoM}}}$ $\displaystyle=\exp{({-\left\lVert
L_{\textrm{CoM}}\right\rVert}^{2})},$ (13) $\displaystyle r_{a}$
$\displaystyle=\exp{({-\left\lVert a_{k}-a_{k-1}\right\rVert}^{2})},$ (14)

and the weights are chosen as $\mathbf{w}^{T}=[0.6,0.3,0.1]$.

One iteration step of the policy corresponds to the interaction of the
learning agent with the environment. The RL policy takes the reduced order
state $s$ and computes an action $a$ that is converted in desired task-space
trajectories $y^{d}$ at the time $t_{k}$. The reference trajectories are then
sent to the task-space controller, which sends torque commands to the robot.
This workflow is depicted in Fig. 1. The feedback control loop runs at a
frequency of $1$ kHz, while the high-level planner policy runs at $50$ Hz. The
maximum length of each episode is $600$ iteration steps, which corresponds to
$12$ seconds of simulated time.

##  IV Simulation results

In this section, we show the performance of the learned planner policy under
different testing scenarios with the bipedal robot Digit. Moreover, we analyze
the latent representation generated during the testing of the RL policy.

###  IV-A Latent State Representation

To visualize the learned latent manifolds, we applied principal component
analysis (PCA) to reduce the dimensionality of the latent space to 2D. Fig. 4
shows the 2D representation of the latent space for dimensions $N=2$, $N=4$,
and $N=8$. The data used for visualization correspond to the robot Digit
walking with the learned RL policy within the range of
$v_{x}^{d}\in[-0.75,1.4]$ m/s. In all three cases, the latent space exhibits a
well-distributed and disentangled representation of the data, with each
walking speed corresponding to a specific area in the 2D plane. For
comparison, Fig. 4 also includes the 2D PCA visualization of the full-order
state of the robot. It is evident that data points for different speeds
overlap with no clear structure. This highlights that the latent space, even
with a very low dimensionality $N$, effectively captures the distribution of
the data in the full-order system. This demonstrates the potential of
exploiting the latent representation for control purposes.

Figure 5: Comparison of velocity tracking performance of the learned policy
with different dimensions of the latent space $N$ with a model-based ALIP
planner described in [5].

###  IV-B Tracking Performance of Different Velocity Profiles

We evaluated the performance of the learned gait policies with latent space
dimensions with $N=2$, $N=4$, and $N=8$ in tracking a velocity profile in
different directions using the Digit robot. As shown in Fig. 5, the policy
successfully tracks walking speeds in the range $v_{x}^{d}\in[-0.75,1.4]$ m/s,
even with aggressive changes in the velocity profile. Importantly, the data
collected to train the autoencoder was within the range of
$v_{x}^{d}\in[-0.5,1.0]$. These results show the effectivity of the learned
latent states to control the walking velocity even with a low dimension $N$.
Furthermore, the results demonstrate that the latent representation captures
the dynamic nature of the walking gait and can be generalized to scenarios
outside the training distribution. Additional tests about this aspect are
presented in Sec. IV-C.

For comparison, Fig. 5 also shows the speed-tracking performance of a model-
based ALIP planner based on the design in [5]. Our policy outperforms the
ALIP-based planner in terms of velocity tracking and offers a wider range of
admissible walking speeds. The ALIP planner fails to maintain a stable walking
gait for speeds higher than 1.2 m/s. In addition, Fig. 6 illustrates the
correspondence between the walking speeds and control actions of the RL
policies with different $N$s and the ALIP planner. Specifically, we focus on
the swing foot landing position with respect to the robot’s base
$p^{sw}_{x,T}$ when the robot is walking at $0.7$ and $1.0$ m/s.
Interestingly, the actions of the RL policy exhibit similar patterns across
different latent space dimensions $N$. This finding aligns with the results
presented in Sec. IV-A, where we showed that even a small $N$ is sufficient to
characterize the dynamics of the full-order system fully. When provided with a
larger latent space dimension, the policy learns to disregard less important
states in the latent space, resulting in similar policy actions for different
$N$ values.

Figure 6: Comparison of the desired swing foot landing positions (the policy
action) from different latent space dimensions $N$ and the ALIP planner.

Furthermore, the actions of the latent space-based policies exhibit similar
patterns to those of the ALIP planner. This is intriguing because the latent
space used by the RL policy does not have a direct physical interpretation.
Nevertheless, the RL policy learns to behave in a manner similar to the
template model-based planner. These behaviors are not enforced during
training, suggesting that the latent space naturally captures the dynamics of
both the template model and the full-order system.

###  IV-C Policy Generalization to Out-of-Distribution Scenarios

In addition to evaluating the policy’s performance on data it was trained on,
we also assess its ability to handle out-of-distribution scenarios. To do
this, we conduct tests where the policy is instructed to maintain a constant
speed while the height of the base is varied. It is important to note that the
training of the autoencoder and the RL policy did not include any data with
varying base heights, as the locomotion data used for training latent spaces
was collected with a fixed base height of 1 m. However, as depicted in Fig. 7,
the policy demonstrates successful tracking of the desired walking speed
regardless of the different base heights commanded. An interesting future work
direction may be the exploration of the robust generalization of the latent
space to generate transitions between different locomotion tasks such as
walking, jogging, sitting, and jumping.

Figure 7: Generalization of the latent space and trained policy to out-of-
distribution data.

Furthermore, we conducted tests to evaluate the policy’s robustness against
external disturbances in the forward and backward directions. The disturbances
ranged from $-100$ N to $60$ N, with durations ranging from $0.1$ s to $1.5$
s. As illustrated in Fig. 8, the policy demonstrated effective reactions to
the different disturbances, successfully maintaining stability and tracking
the desired walking speeds without falling.

Figure 8: Robustness test to external disturbances.

##  V Conclusion

In this work, we present a novel data-driven learning framework to realize
robust bipedal locomotion. The design of the high-level RL gait policy takes
data-driven reduced dimensional latent variables as input states and generates
a set of task space commands, including the robot’s step length with respect
to the base and instantaneous velocity offset of the robot’s base. The latent
representation of the full-order state is obtained using an autoencoder
trained with supervised learning from locomotion data collected with existing
locomotion controllers. Our work shows that the learned latent representation
manifold has a disentangled structure that is directly correlated with the
speed of the walking robot. The insightful choice of the RL state and action
spaces results in a compact policy that learns effective strategies for robust
and dynamic locomotion in simulation. Future work will focus on implementing
and validating the proposed framework on Digit, further demonstrating its
effectiveness and adaptability for real-world bipedal locomotion.

## References

  * [1] A. D. Ames, “Human-inspired control of bipedal robots via control lyapunov functions and quadratic programs,” in _Proceedings of the 16 th international conference on Hybrid systems: computation and control_, C. Belta and F. Ivancic, Eds., ACM. ACM, 2013, pp. 31–32.
  * [2] S. Kajita, F. Kanehiro, K. Kaneko, K. Fujiwara, K. Harada, K. Yokoi, and H. Hirukawa, “Biped walking pattern generation by using preview control of zero-moment point,” in _IEEE International Conference on Robotics and Automation_ , vol. 2, Sept. 2003, pp. 1620–1626.
  * [3] P. M. Wensing and D. E. Orin, “High-speed humanoid running through control with a 3D-SLIP model,” _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pp. 5134–5140, Nov. 2013.
  * [4] A. Hereid, C. M. Hubicki, E. A. Cousineau, and A. D. Ames, “Dynamic humanoid locomotion: a scalable formulation for HZD gait optimization,” _IEEE Transactions on Robotics_ , vol. 34, no. 2, pp. 370–387, Apr. 2018.
  * [5] Y. Gong and J. W. Grizzle, “Zero Dynamics, Pendulum Models, and Angular Momentum in Feedback Control of Bipedal Locomotion,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 144, no. 12, p. 121006, 10 2022.
  * [6] X. Xiong, J. Reher, and A. D. Ames, “Global position control on underactuated bipedal robots: Step-to-step dynamics approximation for step planning,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 2825–2831.
  * [7] Z. Xie, G. Berseth, P. Clary, J. Hurst, and M. van de Panne, “Feedback control for cassie with deep reinforcement learning,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 1241–1246.
  * [8] J. Siekmann, S. Valluri, J. Dao, L. Bermillo, H. Duan, A. Fern, and J. Hurst, “Learning memory-based control for human-scale bipedal locomotion,” in _Robotics science and systems_ , 2020.
  * [9] Z. Li, X. Cheng, X. B. Peng, P. Abbeel, S. Levine, G. Berseth, and K. Sreenath, “Reinforcement learning for robust parameterized locomotion control of bipedal robots,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.
  * [10] K. Green, Y. Godse, J. Dao, R. L. Hatton, A. Fern, and J. Hurst, “Learning spring mass locomotion: Guiding policies with a reduced-order model,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3926–3932, 2021.
  * [11] I. D. J. Rodriguez, N. Csomay-Shanklin, Y. Yue, and A. D. Ames, “Neural gaits: Learning bipedal locomotion via control barrier functions and zero dynamics policies,” in _Learning for Dynamics and Control Conference_. PMLR, 2022, pp. 1060–1072.
  * [12] G. A. Castillo, B. Weng, W. Zhang, and A. Hereid, “Robust feedback motion policy design using reinforcement learning on a 3d digit bipedal robot,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2021, pp. 5136–5143.
  * [13] G. Castillo, B. Weng, W. Zhang, and A. Hereid, “Reinforcement learning-based cascade motion policy design for robust 3d bipedal locomotion,” _IEEE Access_ , vol. 10, pp. 20 135–20 148, 2022.
  * [14] G. A. Castillo, B. Weng, S. Yang, W. Zhang, and A. Hereid, “Template model inspired task space learning for robust bipedal locomotion,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2023.
  * [15] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” in _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2014.
  * [16] Y. Bengio, A. Courville, and P. Vincent, “Representation learning: A review and new perspectives,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 35, no. 8, pp. 1798–1828, 2013.
  * [17] M. Dai, X. Xiong, and A. D. Ames, “Data-driven adaptation for robust bipedal locomotion with step-to-step dynamics,” _arXiv preprint arXiv:2209.08458_ , 2023.
  * [18] T. Salzmann, E. Kaufmann, J. Arrizabalaga, M. Pavone, D. Scaramuzza, and M. Ryll, “Real-time neural mpc: Deep learning model predictive control for quadrotors and agile robotic platforms,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 4, pp. 2397–2404, 2023.
  * [19] X. B. Peng, Y. Guo, L. Halper, S. Levine, and S. Fidler, “ASE: large-scale reusable adversarial skill embeddings for physically simulated characters,” _ACM Trans. Graph._ , vol. 41, no. 4, pp. 94:1–94:17, 2022.
  * [20] S. Starke, I. Mason, and T. Komura, “Deepphase: Periodic autoencoders for learning motion phase manifolds,” _ACM Trans. Graph._ , vol. 41, no. 4, jul 2022.
  * [21] J. W. Grizzle, C. Chevallereau, R. W. Sinnet, and A. D. Ames, “Models, feedback control, and open problems of 3D bipedal robotic walking,” _Automatica_ , vol. 50, no. 8, pp. 1955–1988, 2014.
  * [22] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2015.
  * [23] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017.


# Learning the Domain Specific Inverse NUFFT for Accelerated Spiral MRI using
Diffusion Models

###### Abstract

Deep learning methods for accelerated MRI achieve state-of-the-art results but
largely ignore additional speedups possible with noncartesian sampling
trajectories. To address this gap, we created a generative diffusion model-
based reconstruction algorithm for multi-coil highly undersampled spiral MRI.
This model uses conditioning during training as well as frequency-based
guidance to ensure consistency between images and measurements. Evaluated on
retrospective data, we show high quality (structural similarity >0.87) in
reconstructed images with ultrafast scan times (0.02 seconds for a 2D image).
We use this algorithm to identify a set of optimal variable-density spiral
trajectories and show large improvements in image quality compared to
conventional reconstruction using the non-uniform fast Fourier transform. By
combining efficient spiral sampling trajectories, multicoil imaging, and deep
learning reconstruction, these methods could enable the extremely high
acceleration factors needed for real-time 3D imaging.

Index Terms—  Accelerated MRI, Spiral MRI, Deep Learning, Image Reconstruction

##  1 Introduction

Despite the numerous advantages of MRI, inherent physical and hardware
constraints cap acquisition speed and lead to long scan times. This creates
myriad downstream obstacles, including low patient compliance, inefficient
resource allocation, and image motion artifacts, among others. For this
reason, methods to accelerate MRI acquisition have been and continue to be an
active area of research. Acceleration can be achieved through a combination of
using more efficient scanning sequences and reducing the number of
measurements made in the frequency space of the image, k-space.

Considering the former, successful techniques for faster scanning include
radial and spiral imaging methods, which exploit the unequal distribution of
information across k-space by sampling lower frequencies more densely [1, 2].

Considering the latter, acquiring fewer measurments in k-space rewards a
decrease in scan time proportional the the ratio of undersampling, but comes
at the cost of image quality. Sampling below the Nyquist limit introduces
ambiguities during reconstruction which manifest as artifacts in the final
image. Reconstruction algorithms must therefore leverage additional
information, including multicoil data in the case of parallel imaging, and
sparse priors, in the case of compressed sensing, in order to resolve these
ambiguities [3, 4].

A third more recent approach to undersampled MRI reconstruction lies in deep
learning methods, which essentially learn a set of image priors and use these
to regularize solutions to the ill-posed reconstruction problem [5]. Within
this category, diffusion models stand out for producing state-of-the art
results on image reconstruction tasks for faster scanning, motion correction,
noise reduction, and others [6, 7, 8, 9]. Despite this, the vast majority of
deep learning approaches, and to our knowledge, all diffusion-based approaches
to image reconstruction, focus on Cartesian-sampled MRI, missing out on
potential acceleration gains attained by using more efficient non-Cartesian
sampling trajectories. This work addresses this gap by introducing a diffusion
model-based method for trajectory-agnostic undersampled MRI reconstruction of
multicoil spiral MRI. Contributions:

  * •

Creation of a novel multi-conditioning strategy for solving the inverse non-
uniform fast Fourier transform (nufft) using a learned conditional score
function and weak frequency guidance during sampling

  * •

Efficient hyperparameter search of the joint trajectory-reconstruction space
and identification of optimal sampling trajectories

  * •

Retrospective acquisition and reconstruction of a 2D, 256x256 pixel, 22x22 cm2
image with a readout duration of 0.02 seconds.

Fig. 1: Example trajectories (A) and the corresponding readout gradients in
$k_{x}$ and $k_{y}$ (B). All trajectories shown cover the frequency space of a
256x256 image and have a readout duration of 10.0 ms.

##  2 Background

Canonical score-based models consider the mapping between a known distribution
of independently and identically distributed samples of gaussian noise and an
observed, but unknown distribution of data $p(x)$. These distributions bookend
a Langevin diffusion process described by the following stochastic
differential equation representing the trajectory of a sample from our data
distribution into a sample from our noise distribution:

$dx=f(x,t)dt+g(t)d\mathrm{w}.$ (1)

Here, functions $f(\cdot,t)$ and $g(\cdot)$ are the drift and diffusion
coefficiets of $x(t)$ respectively, and w is a standard Wiener process, or
Brownian motion.

In order to generate a novel sample from our data distribution, we can
generate a random noise vector and attempt to solve the reverse-time SDE, but
this is generally intractable. However, we can approximate this process by
estimating the noise-conditioned score function, $\nabla_{x}\log p_{t}(x)$,
which computes the likelihood of a sample $x$ existing between the noise and
image distributions. With this, the reverse-time SDE becomes:

$dx=\left[f(x,t)-g(t)^{2}\nabla_{x}\log
p_{t}(x)\right]dt+g(t)d\bar{\mathrm{w}}$ (2)

The score function can be trained using a score matching with Langevin
dynamics algorithm [10, 11].

##  3 Methods

This research study was conducted retrospectively using human subject data
made available in open access by [12]. Ethical approval was not required.

Fig. 2: Given measurements $y_{0}$, reconstruction follows a modified
diffusion sampling process. At each timestep, a noisy latent $x_{t}$ is
concatenated with a prior $p_{0}$ and passed to the denoising model to obtain
$\tilde{x}_{t-1}$. To enforce consistency with $y_{0}$, we compute a frequency
gradient $\nabla y_{t-1}$ and solve for the image gradient using a modified
iterative inverse nufft (section 4). A weighted sum of $x_{t-1}$ and $\nabla
x_{t-1}$ yields the corrected image $x_{t-1}$. This is repeated until $t=0$.

###  3.1 Data

We use the NYU FastMRI dataset [12] consisting of 6970 fully sampled 2D brain
scans on hardware ranging from 4 to 24 coils. For training and testing, we
consider axial T2 weighted turbo spin echo sequences characterized by the
following sequence parameters: scan time=140 s, TR=6 s, TE=113 ms, slices=30,
slice thickness=5 mm, field of view=22 cm, matrix size=320x320. Effective scan
time for a 2D slice at 2562 resolution is $140\mathrm{s}/320*256/30\approx
3.7\mathrm{s}$. As this data is initially acquired using Cartesian sequences,
we simulated spiral acquisition by retrospectively interpolating in k-space to
attain complex-valued measurements along generated spiral trajectories.

###  3.2 Generating spiral trajectories

Following Kim et al. [13], we consider spiral trajectories of the form

$k(\tau)=\int_{0}^{\tau}\frac{1}{\rho(\phi)}d\phi
e^{j\omega\tau}\;\;\approx\;\;\lambda\tau^{\alpha}e^{j\omega\tau}.$ (3)

Here, $\rho$ denotes sampling density, $\tau$ is a function of time, $\phi$ is
angular position, $\omega=2\pi n$ is frequency, with $n$ the number of turns
in k-space, $\lambda$ a scaling factor equal to $matrix\;size/(2*FOV)$, and
$\alpha$ is a bias term for oversampling the center of k-space relative to the
edges. Solving this parametric equation under the constraints of capped
gradient slew rate and capped gradient amplitude yields gradients ($g_{x}(t)$
and $g_{y}(t)$) as well as a spiral trajectory in the $k_{x}$,$k_{y}$-plane
(figure 1). In doing so, we can tune sampling parameters to control for
factors such as read out duration and dwell time, while varying the number of
interleaves and ratio of low-to-high frequency oversampling.

Fig. 3: Representative reconstruction results for a single 2D 16 coil image.
Retrospective k-space data was sampled with an optimized 23 interleave
sequence with a total readout duration of 0.02 s. Rows 1 and 2 show the RSS-
reconstructed images and log-scaled k-space magnitudes for the ground truth,
inverse nufft, and proposed model reconstructions. Below are the individual
coil magnitude and phase images for the fully sampled image, the inverse nufft
reconstructions, and the model predictions.

Fig. 4: We performed a grid hyperparameter search over a 2D trajectory space.
We fixed readout duration at 0.02 seconds and varied the number of interleaves
from 1 to 125 and alpha from 1 to 4. Based on structural similarity of the
model-reconstructed images, we found multiple trajectories that yield improved
image quality. In comparison, the naive Archimedean spiral, corresponding to 1
interleave and $\alpha=1$, performs very poorly.

Fig. 5: (A) Effect of sampling trajectory optimization, model reconstruction
without frequency guidance, and model reconstruction with frequency guidance.
For the non-optimized trajectory, we used a single interleave Archimedean
spiral with a readout duration of 0.02 s. The optimized trajectory uses a 23
interleave, $\alpha=1.23$ sequence with an identical readout duration. (B)
Snapshots of the image latent $x_{t}$ and the gradient signal $\nabla x_{t}$
taken during a diffusion sampling process.

###  3.3 Image reconstruction is inverse problem solving

MRI undersampled acquisition amounts to measuring an unknown signal $x$
through some imperfect sampling function $A$: $y=Ax+\epsilon$. Here, $y$ is
the measured multicoil k-space data, and $A$ is the non-uniform fourier
transform. $\epsilon$ is measurement noise and exists in the same domain as
the $y$; in MRI, noise is gaussian-distributed across the real and imaginary
components of $y$ for each coil.

Reconstruction is an ill-posed inverse problem of recovering an image signal
$x$ from a set of incomplete k-space measurements $y$. As $x$ and $y$ exist in
different domains, $x$ is hidden behind a sampling operator $A$. Solving this
problem necessitates prior knowledge. In our case, we learn an underlying
conditional distribution of images and seek to reconstruct samples from this
distribution consistent with the measurements. Information is supplied in two
forms: first, we learn a conditional score function $\nabla_{x}\log
p_{t}(x_{t}|\tilde{A}^{-1}y_{0})$, where $y_{0}$ is the measurement in
frequency space and $\tilde{A}^{-1}$ is an approximate inverse of $A$, in our
case the inverse nufft solved iteratively using conjugate gradients. We find
that adding this supervision during training helps to constrain the model when
faced with a large number of input image channels and the periodic ambiguity
inherent when operating on complex numbers.

Second, we use frequency space gradients to weakly guide the sampling process.
At each time step during sampling, we compute the forward nufft of an
uncorrected $\tilde{x}_{t-1}$ and take a difference between that and the
measurements $y_{0}$. To minimize this difference, the approximate gradient in
image space is calculated by solving a modified approximate inverse nufft
$\tilde{A_{t}}^{-1}$, which corrects for low frequency biases and applies time
step-dependent noising determined by the noise schedule $\sigma(t)$, which is
necessary for sampling with langevin diffusion.

$\tilde{A_{t}}^{-1}x_{t}=\frac{\tilde{A}^{-1}x_{t}}{c_{1}e^{-c_{2}r^{2}}}+\mathcal{N}(0,\sigma(t)^{2})$
(4)

Following [11], we choose a linear noise schedule and observe that the
underlying ordinary differential equation describing transit from latent to
image is locally linear, so summation of $\tilde{x}_{t-1}$ and $\nabla
x_{t-1}$ to obtain a frequency-corrected image $x_{t-1}$ is akin to gradient
descent (figure 2).

In practice, due to the non-invertibility of the nufft, imperfections in the
approximate inverse nufft bleed into the final image reconstruction,
introducing artifacts and reducing quality. To avoid this, we anneal the
guidance signal following an empirically chosen linear schedule
$\gamma(t)=\beta(1-t)$, ensuring strong guidance at the outset of sampling and
minimal artifacts at the end of sampling. A consequence of this choice is that
we do not strongly enforce that $Ax_{0}=y_{0}$ at time 0.

##  4 Results

Model reconstruction performance was evaluated on a held-out test dataset.
Test trajectories have a fixed readout duration of 0.02 seconds, in which time
the measurements needed to reconstruct a 256x256 pixel, 22x22 cm2 2D image are
acquired. Reconstructed image quality was scored using structural similarity
(SSIM) (figure 3).

To investigate the effect choosing different scanning trajectories has on the
quality of reconstructed images, we also perform a grid hyperparameter search
of spiral trajectories with a fixed readout duration of 0.02 seconds (figure
4) and varying $\alpha$ and interleaves. Surprisingly, the common ’naive’
trajectory, a single interleave Archimedean spiral, corresponding to
$\alpha=1$, performs very poorly when sampled below the Nyquist limit.
Trajectories which perform better tended to lie along two logarithmic curves
roughly characterized by $\alpha=1.33\log(0.39\,interleaves)$ and
$\alpha=0.87\log(0.54\,interleaves)$.

Finally, we conduct an ablation study to disentangle the effects of optimal
sampling trajectory without model reconstruction, model reconstruction without
frequency guidance, and model reconstruction with frequency guidance (figure
5). We find that all three contribute to noticeable increases in image
quality, both visual, and quantitative based on SSIM. The combination of
choosing an optimal trajectory, performing model reconstruction with
conditioning, and using annealed frequency guidance results in large
improvements in image quality, up to and exceeding a 0.15 boost in SSIM.

##  5 Discussion

While initial results are promising, the main limitation of this project is
the reliance on retrospective, Cartesian-sampled data. Implementing the
sequences outlined in this work will likely require customizing spiral
sequences so as to match the contrast and signal of the original Cartesian
sequences, which will constrain the space of realizable trajectories. Until a
dataset of raw non-Cartesian MRI data becomes available, this will continue to
be an obstacle. For a similar reason, it is difficult to make head-to-head
comparisons between the original sequence and the proposed sequences without
prospective validation. For this reason, the immediate task is to acquire
prospective data with custom sequences and use it to validate image
reconstruction. Currently, a concrete direct comparison is between the
proposed sequences and their Nyquist-sampled counterparts, which run roughly
3x longer.

Apart from the short-term task of matching sequence parameters between spiral
and Cartesian sequences, the choice of spiral sequence leaves a considerable
amount of flexibility even within the space of optimized interleave and
density pairs identified above. Variation in number of interleaves, and by
extension the duration of a single interleave, allows for tailoring of
sequence contrast, signal, and speed to task-specific requirements. An area in
which these sequences could provide additional benefit, even outside of sheer
acceleration, would be in imaging tissues with a very short T2*, as
acceleration within an interleave allows for proportionally more data
acquisition to occur before signal has decayed.

##  6 Conclusion

Here we introduce a new method and show preliminary results for reconstructing
spiral MRI using a diffusion model. Combining multicoil imaging, spiral
scanning, and undersampling enables dramatically faster imaging speeds.
Applications of this work are widespread; in addition to the numerous typical
benefits associated with faster scanning, including better patient compliance
and fewer motion artifacts, these methods have the potential to reach the
extremely high acceleration factors necessary to achieve high resolution real-
time 3D imaging.

##  7 Acknowledgments

No funding was received for conducting this study. The authors have no
relevant financial or non-financial interests to disclose.

## References

  * [1] Mark J Blum, Michael Braun, and Dov Rosenfeld,  “Fast magnetic resonance imaging using spiral trajectories,”  in Medical Imaging. SPIE, 1987, vol. 767, pp. 40–46.
  * [2] Stefanie Winkelmann, Tobias Schaeffter, Thomas Koehler, Holger Eggers, and Olaf Doessel,  “An optimal radial profile order based on the golden ratio for time-resolved mri,”  IEEE transactions on medical imaging, vol. 26, no. 1, pp. 68–76, 2006.
  * [3] Daniel K Sodickson and Warren J Manning,  “Simultaneous acquisition of spatial harmonics (smash): fast imaging with radiofrequency coil arrays,”  Magnetic resonance in medicine, vol. 38, no. 4, pp. 591–603, 1997.
  * [4] Michael Lustig, David Donoho, and John M Pauly,  “Sparse mri: The application of compressed sensing for rapid mr imaging,”  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007.
  * [5] Julio A Oscanoa, Matthew J Middione, Cagan Alkan, Mahmut Yurt, Michael Loecher, Shreyas S Vasanawala, and Daniel B Ennis,  “Deep learning-based reconstruction for cardiac mri: A review,”  Bioengineering, vol. 10, no. 3, pp. 334, 2023.
  * [6] Yang Song, Liyue Shen, Lei Xing, and Stefano Ermon,  “Solving inverse problems in medical imaging with score-based generative models,”  arXiv preprint arXiv:2111.08005, 2021.
  * [7] Patricia M Johnson, Michael P Recht, and Florian Knoll,  “Improving the speed of mri with artificial intelligence,”  in Seminars in musculoskeletal radiology. Thieme Medical Publishers, 2020, vol. 24, pp. 012–020.
  * [8] Zhuo-Xu Cui, Chentao Cao, Jing Cheng, Sen Jia, Hairong Zheng, Dong Liang, and Yanjie Zhu,  “Spirit-diffusion: Self-consistency driven diffusion model for accelerated mri,”  arXiv preprint arXiv:2304.05060, 2023.
  * [9] Asad Aali, Marius Arvinte, Sidharth Kumar, and Jonathan I Tamir,  “Solving inverse problems with score-based generative priors learned from noisy data,”  arXiv preprint arXiv:2305.01166, 2023.
  * [10] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole,  “Score-based generative modeling through stochastic differential equations,”  arXiv preprint arXiv:2011.13456, 2020.
  * [11] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine,  “Elucidating the design space of diffusion-based generative models,”  Advances in Neural Information Processing Systems, vol. 35, pp. 26565–26577, 2022.
  * [12] Florian Knoll, Jure Zbontar, Anuroop Sriram, Matthew J Muckley, Mary Bruno, Aaron Defazio, Marc Parente, Krzysztof J Geras, Joe Katsnelson, Hersh Chandarana, et al.,  “fastmri: A publicly available raw k-space and dicom dataset of knee images for accelerated mr image reconstruction using machine learning,”  Radiology: Artificial Intelligence, vol. 2, no. 1, pp. e190007, 2020.
  * [13] Dong-hyun Kim, Elfar Adalsteinsson, and Daniel M Spielman,  “Simple analytic variable density spiral design,”  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, vol. 50, no. 1, pp. 214–219, 2003.