#### State Estimation with Structural Priors in fMRI

State Estimation with Structural Priors in fMRI
Ville-Veikko Wettenhovi 0 1
Ville Kolehmainen 0 1
Joanna Huttunen 0 1
Mikko Kettunen 0 1
Olli Gröhn 0 1
Marko Vauhkonen 0 1
0 A.I. Virtanen Institute for Molecular Sciences, University of Eastern Finland , 70211 Kuopio , Finland
1 Department of Applied Physics, University of Eastern Finland , 70211 Kuopio , Finland
We propose a state estimation approach for functional magnetic resonance imaging (fMRI). In state estimation, time-dependent image reconstruction problem is modeled by separate state evolution and observation models, and the objective is to estimate the time series of system states, given the models and the time-dependent measurement data. Our method computes the state estimates by using the Kalman filter (KF) and Kalman smoother (KS) algorithms. We propose to complement the state estimation formulation with a structural prior which can be derived from the anatomical MRI, acquired as a part of the fMRI protocol. Two different constructions of the structural prior are considered. The first one is a structured smoothness prior where the state observation matrix is augmented with a spatially weighted regularization matrix which promotes structural similarity of the gradient of the unknown image with the gradient of the anatomical image. The second approach is based on applying structured total variation denoising to the KF estimate at each time step of the Kalman recursions. The proposed approaches are evaluated using simulated and experimental, radially sampled, small animal fMRI data from a rat brain. In our method, the state estimates are updated after each new spoke of radial data becomes available, leading to faster frame rate compared to the conventional image reconstruction approaches. The results are compared to a sliding window method and a conventional reconstruction which produces new image only after a full circle of k-space spokes becomes available. The results suggest that the state estimation approach with the structural prior can improve both spatial and temporal resolution of fMRI.
Functional magnetic resonance imaging; Kalman filter; Structural prior; Total variation; Radial MRI; GPU
1 Introduction
Functional magnetic resonance imaging (fMRI) is an
imaging technique that can be used to noninvasively monitor brain
activity. The working principle of fMRI is based on the fact
that neuronal activity is coupled with cerebral blood flow
and volume. When a certain part of the brain is activated,
the blood flow and amount of oxygen to that area increase
more than the metabolic need for oxygen, enabling the related
brain activity to be indirectly evaluated by the blood oxygen
level-dependent (BOLD) contrast [
2
].
In a typical fMRI scan, the MRI scanner collects a time
series of k-space data, which is used to reconstruct a time
series of images. The image reconstruction is traditionally
carried out using static image reconstruction methods, based
on assumption of time invariance of the unknowns during
the acquisition of the data for the selected k-space
trajectory. Typically, when the k-space trajectory is a fully
sampled Cartesian trajectory, the reconstruction is carried
out using inverse FFT. When employing non-Cartesian
sampling schemes, the usual reconstruction is regridding where
the non-Cartesian data are interpolated onto a Cartesian grid
before the inverse FFT [
16
]. To obtain sufficient spatial
resolution with the classical methods, enough k-space samples
per image frame are necessary in order to fulfill the Nyquist
criterion. This comes at the expense of temporal resolution
resulting in a loss of accuracy when studying a physiological
process with changes faster than the sampling time for each
image frame.
One approach to partially alleviate the loss of temporal
resolution is to use a sliding window method [
6
], where the last
available measurement data are combined with a time
window of previous data in order to obtain a sufficiently sampled
frame of k-space data. This approach, however, is
suboptimal as the reconstruction problem becomes inconsistent in
the sense that changes might have occurred in the target
during the time window of the augmented data frame.
A natural approach to improve the temporal resolution is
to use undersampled k-space data. When using undersampled
data, it is preferable to use non-Cartesian acquisition, such
as spiral or radial sampling, as it allows for denser sampling
of the center of the k-space [
1
]. However, if undersampled
data are used, the image reconstruction problem becomes
ill-posed. This can lead to severe aliasing artifacts when
conventional reconstruction methods, such as inverse FFT
or regridding, are used. A well-known static reconstruction
approach that allows for the use of drastically undersampled
data is compressed sensing (CS) [
4
]. CS image
reconstruction has been extensively used in static MRI, in cardiac MRI
[
9,18,25,26,29
] and recently also in fMRI [
14,17,28,51
].
In this paper, we propose a state estimation approach to
fMRI. In the state estimation paradigm, the image
reconstruction problem is considered explicitly as a time-dependent
problem where the image can be updated after a new
observation, such as a single spoke of radial data, becomes available.
A state evolution model is used to model the unknowns
(time series of fMRI images) as a time-dependent process,
while the relation between the unknown fMRI image and
k-space measurements at each time instant is modeled by a
separate observation equation. The objective is to estimate
the sequence of states (fMRI images) given the models and
the time series of k-space data. Since the state estimation
approach allows an update of the state estimate after each new
measurement becomes available, it avoids the data
inconsistency which is present in the conventional static approaches
due to longer time window of acquiring data for a single
image in the fMRI time series.
One of the most commonly used methods for
computing the state estimates is the Kalman filter (KF) [
21
], which
can additionally be smoothed by using the Kalman smoother
(KS). In dynamic cardiac MRI, KF has been utilized in
[
8,27,32,39,40,42
] and KS in [32]. In fMRI, KF has been
mainly used to estimate certain parameters from the
reconstructed data [
10,24,36
]; KF, enhanced with CS, has been
used in fMRI image reconstruction [
49
].
Further, we propose to complement the state estimation
approach with a structural prior, which is designed to promote
structural similarity of the gradient of the unknown image
with the gradient of the anatomical MRI. This anatomical
image is usually taken as part of the fMRI measurement
protocol, but used only for visualization purposes. Two different
reconstructions are considered. The first one is a structured
smoothness prior, obtained by augmenting the state
observation matrix with a spatially weighted difference matrix,
where the weighting is constructed based on the
magnitude of the gradient of the anatomical MRI image [
19
]. The
approach corresponds to augmenting the observation
updating density of Bayes filtering with a Gaussian prior model
for the spatial features of the unknown state [
20
]. The
second approach is an approximation, where the spatial prior
information is utilized by applying a structurally guided total
variation (TV) denoising to the Kalman filter estimates at
each time step. The rationale behind this heuristic
approximation is to obtain faster computational times compared to
the augmented model. However, the approximation does not
have known convergence results, while the augmented model
has. In recent years, there have been several contributions to
include an anatomical prior into a static image
reconstruction problem, see, e.g., [
5,7,44,50
]. In this paper, the idea
of structured regularization is extended to the dynamic state
estimation framework.
The state estimates for both approaches are computed
using the Kalman filter and the Rauch–Tung–Striebel (RTS)
fixed-interval Kalman smoother. In this study, the state
estimation approach is evaluated using simulated and
experimental, sparsely collected radial fMRI data from a rat brain.
In the computations, the state estimates are updated after each
new spoke of radial data becomes available. The approaches
are compared with the sliding window method using a full
circle of n spokes as well as with a conventional
frame-byframe method using the k-space data in frames of n-spokes.
The sliding window and conventional estimates are
reconstructed with linear least squares (LS).
2 Theory
2.1 Forward Model and Conventional Reconstruction
Procedure
Let z j ∈ CM denote the k-space data for a single spoke in the
radial sampling of the k-space, and let H j denote the forward
matrix modeling the mapping from the image space to the
j th spoke of data in the data space. Assuming that the radial
sampling pattern consists of n equally spaced spokes around
the full circle, the data and forward matrix for one traditional
fMRI image frame can be written as
⎛ z1⎞
z = ⎜⎜ z.2⎟⎟
⎜ .. ⎟
⎝ ⎠
zn
⎛ H1⎞
⎜ H2
H = ⎜⎜ ... ⎟⎟⎟ ,
⎝ ⎠
Hn
(1)
respectively, leading to an observation equation
z = Hf + v
for a single frame of the sampling scheme. In (2), v models the
measurement noise and f ∈ CNpix is the vector representation
of the complex-valued Npix = N × N image.
Assume that the fMRI experiment consists of Nframes
frames of data of the form (2), and let z( ) denote the th frame
of data in the time series {z( ) ∈ CnM , = 1, . . . , Nframes}.
In the conventional setup to fMRI reconstruction, the
respective image frames of the time series {f( ) ∈ CNpix , =
1, . . . , Nframes} are reconstructed from the data vector z( )
using conventional static methods, such as the regridding
algorithm or LS. When calculating the LS solution, we need
to solve the following minimization problem
f ( )
ˆ
= mf(il)n z( ) − Hf (l) 2 ,
which can be solved iteratively by, for example, the LSQR
algorithm [
31
].
We remark that when working with radially sampled data,
one can carry out the reconstruction either by using the raw
k-space data directly or by using data which have been
transformed to complex-valued Radon transform formalism by
1D Fourier transform of each measured spoke of k-space
data. In the former case, H corresponds to a complex-valued
(non-uniform) FFT matrix and in the latter to a real-valued
Radon transform operator.
2.2 State Estimation Representation
In the state estimation approach, the fMRI problem is
modeled by the pair of equations
f t = Ft−1f t−1 + wt−1
zt = Ht f t + vt ,
where (4a) is the state evolution model and (4b) is the
equation for the observation zt at time index t . In (4a), f t ∈ CNpix
is system state at time t , Ft−1 is the state transition matrix,
and wt−1 is the process noise. Ht in (4b) is the observation
matrix. In this study, zt ∈ CM corresponds to the k-space
data for a single radial spoke. Therefore, the state estimation
approach provides n images during the cycle of one
conventional frame of (2) and t = 1, . . . , T , where T = n Nframes,
images from the whole fMRI experiment with Nframes
conventional image frames.
The Kalman filter (KF) is a linear recursive state
estimation algorithm that can be used to estimate the states of a
dynamic system of the form (4). Under the standard
assumptions that vt and wt in (4a)–(4b) are Gaussian zero-mean
(2)
(3)
(4a)
(4b)
noise processes with known covariances and mutually
independent in the sense that vk ⊥ v j , wk ⊥ w j for k = j
and wk ⊥ vk, the Kalman filter produces the minimum
mean-square-error estimate of ft based on the
measurements z0, . . . , zT . The derivation of the Kalman filter can be
found from [
38
]. In the statistical setup, the KF produces the
means and covariances of the time evolution and observation
updating probability densities for a Gaussian Bayes filtering
problem, see [20, Chapter 4]. The standard KF equations are
given by
fˆ0+ = E (f 0)
P0+ = E [(f 0 − fˆ0+)(f0 − fˆ0 ) ]
+ T
fˆt− = Ft−1fˆt+−1
Pt− = Ft−1Pt+−1FtT−1 + Qt−1
Kt = Pt−HtT(Ht Pt−HtT + Rt )−1
fˆt+ = fˆt− + Kt (zt − Ht f ˆt−)
Pt+ = (I − Kt Ht )Pt−,
(5a)
(5b)
(5c)
(5d)
(5e)
where fˆt− is the (complex-valued) a priori estimate, Pt− is the
a priori error covariance of the estimate, Qt is the covariance
of the process noise, fˆt+ is the a posteriori estimate, Pt+ is
the a posteriori error covariance of the estimate, Kt is the
Kalman gain, and Rt is the covariance of the observation
noise. The error covariances Pt represent the uncertainty in
the a priori and a posteriori estimates given the model and the
data; large values mean high uncertainty, while smaller values
mean more reliable estimate [
38
]. The covariance Qt reflects
the uncertainty in the state evolution model, and it
determines the weighting that is given for measurements; large
covariance values give more weight to the current
measurement, while smaller covariance values emphasize the past
estimates. Small covariance values lead to smoother
realizations of the time series, but can also cause the estimate to lag
behind, while larger values can lead to noisy images.
Since the observation noise vt in MRI is thermal noise, it
can be modeled well as zero-mean white Gaussian process
with covariance Rt [
39
]. Furthermore, we assume that there
is no correlation in the measurement noise between the time
steps and that the variance of the observation noise is the same
with all k-space samples, leading to a covariance matrix of
the form Rt = σvI.
The state evolution model can be used to incorporate
temporal prior information and models about the unknowns into
the image reconstruction problem. In this study, we employ
a simple random walk formulation, where Ft−1 = I and
the process noise wt is modeled as zero-mean Gaussian
process with diagonal covariance Qt = σw2 I. In other words,
we use random walk formulation to model the unknown
process where the variance σw2 is used to control how much the
process is allowed to change between consecutive time steps.
2.3 Structural Priors
In this section, we describe how the anatomical image, which
is acquired as part of the fMRI experiment, is used to
incorporate spatial prior information about the target into the image
reconstruction problem. The rationale in the structural prior is
to promote structural similarity of the unknown fMRI images
with the gradient of the anatomical image. In the
following, we denote the magnitude of the anatomical image by
f˜ref = |f ref |.
2.3.1 Structured Smoothness Prior
Structured smoothness prior can be obtained by
augmenting the observation equation in the state estimation model
(4b) by a linear regularization term, leading to the following
augmented observation equation
zt
ακ Lf˜ref
=
Ht
ακ L
f t + vt ,
where α > 0 is a regularization parameter, that controls the
strength of spatial regularization, L is an Npix × Npix first
order difference matrix, f˜ref is the anatomical MRI, which has
been (if needed) resized to the same resolution with the fMRI
problem, and κ is a spatially varying weighting functional
κ = κ x + κ y
where κ x is
κx = exp ⎝ −
and i denotes the pixel numbers. κ y is the same, but for the
y-direction. The coefficient κ controls the weighting of the
difference matrix in a given pixel, and C is a positive
constant that can be considered as an edge threshold parameter; if
the edge magnitude in the anatomical image is clearly larger
than C , the weighting κ for the smoothness regularization
in that particular location will be significantly smaller than
one, implying locally less penalty for an edge in the unknown
image at that location, and edges that have magnitude smaller
than C in the anatomical image are considered as noise. For
details on constructing a similar type of structured
regularization for a static inverse problem, see [
19,22,43
].
To make Eq. (6) compatible with the Kalman filter
recursions (5a)–(5e), we write
z˜t = H˜ t f t + v˜ t ,
is a weighting matrix, which implements a structured TV
using the parallel level sets approach [
7,23
], and Ω the whole
(6)
(7)
(8)
where
z˜t =
H˜ t =
zt
ακ Lf˜ref
αHκtL .
R˜ t =
R0t 0I .
Using these notations, z˜t replaces zt and H˜ t replaces Ht
in (5c)–(5e). Rt in (5c) becomes R˜ t which is
(10)
(11)
(12)
The use of the augmented model (6) corresponds to
augmenting the observation updating probability density of a
Gaussian Bayes filtering problem with a Gaussian prior
model for the spatial features of the state f t . The augmented
Kalman filter algorithm produces the means and covariances
for the evolution and modified observation updating
densities. For a detailed derivation and interpretation of the
augmented model, see [20, Chapter 4.4] and [
37
].
2.3.2 Structured Total Variation Denoising
As the second approach for incorporating structural prior
information into the state estimation, we consider a heuristic
approximation, where a structurally guided TV denoising is
applied to the Kalman filter estimate (5d) at each time step
t of the Kalman recursion [utilizing the original observation
equation in (4)]. The rationale in the approach is to utilize
the spatial prior information in an approximate but
computationally efficient way. However, a downside of the approach
is that it does not have theoretical results for convergence.
TV was first introduced for image denoising by Rudin et
al. in 1992 [
35
], and it has since then been used extensively
in image denoising and regularization of image
reconstruction problems. TV denoising has also been previously
implemented in fMRI [
28
], but did not include anatomical
information.
The functional we minimize in the denoising step is a
locally and directionally weighted TV functional
Ψ (uts ) =
Ω
(∇uts )TD(∇uts ) + β dx ,
where ut := fˆt+ is the KF estimate from (5d) at time t and
D = I − λννT
(13)
(14)
(15)
(16)
When calculating the KS estimate, the KF recursion (5a)–
(5e) is calculated first through t = 1, . . . , T (forward pass
through the data) and the smoother is calculated as a
backward pass. This requires saving the smoother gains (17a)
during the forward pass.
3 Materials and Methods
3.1 Estimates
The following estimates are computed and compared in the
results section:
(AKF)
(AKS)
(TV-KF)
(TV-KS)
(LS)
(SW)
Augmented Kalman filter with structured
regularization. Equations (5a)–(5e), using the observation
Eq. (9).
Kalman smoother applied to the results of AKF.
Kalman filter (5a)–(5e) with application of
structured TV denoising to the KF (5d) estimate at each
time step. Observation equation given by (4)
Kalman smoother applied to the results of TV-KF.
Conventional static frame-by-frame estimation
based on successive data frames of form (2).
Reconstructions are computed using the LS
method (3).
Sliding window using a full frame of n-spokes,
such that the last measured spoke of data (say with
index n − k over the rotation of n spokes) is
augmented atop of the previous n −1 spokes. The data
and observation matrix are of the form
image domain. Choosing D = I would convert the functional
to conventional (isotropic) TV. λ in (12) is defined as
(18)
and ν is defined as
∇f˜ref
where G(ut ) is a function that defines the domain of ut ,
s = 0, . . . , S is the iteration index of the interior point
algorithm with a total of S iterations and γ is a regularization
parameter that controls the strength of the (structured) TV
denoising. In this study, the function G(ut ) is selected as
G(ut ) = (1/2) ut 2, leading to gradient descent iteration
uts+1 = uts − γ W∇Ψ (uts )
of Ψ (ut ), where W a weighting matrix. The weighting matrix
W in (11) is selected as the a posterior error covariance
Pt+ in (5e) and has been included to balance the weighting
of the TV denoising based on the magnitude of the errors
of the state estimate; the denoising smooths more the areas
of high uncertainty, while having less of an effect on areas
of higher certainty. The inclusion of the weighting matrix
improves the contrast of the KF estimates and was validated
experimentally.
With the selected G, we do not impose constraints to the
domain of the estimate. By using different choices of G(ut ),
one could select G(ut ) such that it imposes, for example,
a positivity constraint to the solution uts+1. In this case we
could choose G(ut ) = utT ln ut −utTe, where e is a unit vector
[
34
].
2.4 Kalman Smoother
The Kalman smoother (KS) implemented in this work is
the Rauch-Tung-Striebel (RTS) fixed-interval smoother. The
RTS smoother is defined with the following equations
K f s,t = Pt+FtT(Pt−+1)−1
fˆ f s,t = fˆt+ + K f s,t (fˆ f s,t+1 − fˆt++1),
where K f s,t is the smoother gain and t = T − 1, . . . , 0.
(17a)
(17b)
⎛
zn−k ⎞
z = ⎜⎜⎜⎜⎜⎜ zn−z...nk+1⎟⎟⎟⎟⎟⎟ , H = ⎜⎜⎜⎜⎜⎜
⎜⎜⎜⎜ z...1 ⎟⎟⎟⎟ ⎜⎜⎜⎜
⎝ ⎠ ⎝
zn−k−1
⎛
Hn−k ⎞
Hn−k+1⎟
... ⎟⎟⎟
Hn ⎟⎟ .
H1 ⎟⎟
... ⎠⎟⎟
Hn−k−1
The sliding window reconstructions are computed
by the least squares approach using data of the
form (18).
For an fMRI experiment of Nframes conventional frames
with n spokes in each frame, the Kalman estimates produce
a time series of n Nframes images, the sliding window a time
series of n(Nframes − 1) + 1 images and the LS a time series
of Nframes images from the same time interval. First sliding
window estimate is made after n measurements have been
collected, causing the sliding window estimate to be n − 1
shorter than the Kalman estimates.
We remark that due to the high amount of time points in
the fMRI time series, the RTS smoother is not calculated for
the entire duration nor after every time point in the estimates
AKS and TV-KS. Instead, smoother gain matrices (17a) are
calculated for all time points, but only h gain matrices are
stored in memory at every time point; previous gain matrices
are removed when the next ones are calculated. Furthermore,
for increased computational speed, after each smoother
calculation of (17b) there are b time points that skip the smoother
calculations. Selection of the parameters is a compromise
between reconstruction properties and computational cost;
increasing both h and b would give a temporally smoother
estimate at the cost of increased memory usage and
computational cost.
3.2 Experimental fMRI Data from a Wistar Rat
Experimental fMRI data were acquired in a 9.4 T
horizontal magnet (Agilent, Palo Alto, USA) using a volume coil
transmit/quadrature surface coil receive pair (Rapid Biomed,
Rimpar, Germany). For animal experiments, a gradient
echobased radial pulse sequence (repetition time 20 ms, echo time
10 ms, flip angle 30◦, field-of-view 32 × 32 mm, slice
thickness 1.5 mm, number of points in each spoke 64, 52 spokes
collected in sequential order, time per image 1.04 s) was used
to repeatedly collect data before and during nicotine injection
for 18 min.
There were a total of 1000 frames with 52 spokes in the
data. In these computations, we used Nframes = 650 frames
of data by cutting a part of the initial baseline and end of the
experiment away. Thus, while the LS approach produced a
time series of 650 images, the state estimation produced a
time series of T = n Nframes = 52 · 650 =33,800 images and
the sliding window a time series of T = n(Nframes − 1) +
1 =33,749 images from the same time interval.
An anatomical image was collected from the same slice
using a gradient echo pulse sequence with otherwise identical
acquisition parameters but using a full Cartesian sampling to
obtain a 128 × 128 anatomical image as the reference image.
For the reconstructions, the anatomical image was resized
to the same pixel resolution (Npix = 642) with the fMRI
and scaled to magnitude values between [
0, 1
]. The resized
magnitude image is shown in Fig. 1.
All animal procedures were approved by the Animal
Ethics Committee of the Provincial Government of
Southern Finland and conducted in accordance with the
guidelines set by the European Community Council Directives
2010/63/EU. Adult male Wistar rat weighing 385 g was
anesthetized with isoflurane for surgery. The femoral artery was
cannulated to allow the monitoring of the blood gases and
pH. The femoral vein was cannulated for the
administration of nicotine. After surgery, the anesthesia was changed
to urethane (1.25 g/kg, i.p.). This ensured that the anesthesia
level during functional imaging was not at the surgical level
to hinder the detection of activation [
15
]. Rat was allowed
to breathe normally a mixture of 70% N2–30% O2, and the
breathing rate was monitored. Rat was fixed in nonmagnetic
stereotaxic frame with a bite bar and earplugs for the MRI
experiments. Temperature was maintained at approximately
37◦C with water circulation.
During the fMRI experiment, a dose of 0.25 mg/kg
nicotine tartrate salt (0.081 mg/kg free base) was injected
intravenously. Nicotine was chosen because it is widely used in
pharmacological studies and is known to produce large
cortical responses [
11
] under these experimental conditions [
30
].
For the analysis of the results, a region of interest (ROI)
from the cortex was selected. The ROI is highlighted with a
red line in Fig. 1.
3.3 Simulated Data
The 128 × 128 anatomical image of the Wistar rat brain was
used as the baseline target of the fMRI simulation. A region
of interest (ROI), highlighted in red in the left image of Fig. 2,
was defined for the simulation.
Radial k-space data were simulated by using the
MATLAB codes by Guerquin-Kern [
12
] and Guerquin-Kern et al.
[
13
]. Fifty-one equiv-spaced spokes were used in one frame
for a total of Nframes = 50 conventional measurement frames,
leading to times series of 2550 true states. The states were
obtained by adding a simulated, time-varying response with
peak intensity of 0.1 to the pixel values inside the ROI. The
right image in Fig. 2 shows the mean of the time series of
true images in the ROI.
Gaussian noise with standard deviation (STD) of 0.005
was added to simulate physiological noise in the phantom
images. Random noise was added to simulated raw k-space
data using the code package’s noise simulator, where the
noise power is computed relatively to the highest k-space
frequencies sampled. The average STD of the noise was 0.001,
leading to average SNR of 4.08 in the simulated fMRI data.
Resized version of the baseline phantom was used as the
anatomical prior f˜ref and normalized to have values between
[
0, 1
].
3.4 Kalman Filter Computations
Computations were done in MATLAB (2015a, The
MathWorks Inc., Natick, MA) with Nvidia Geforce GTX Titan X
and Nvidia Tesla K40c using single precision.
With both simulated data and experimental data, the
TVKF and TV-KS estimates were computed with Algorithm
1 while the AKF and AKS estimates were computed with
Algorithm 2.
The structured TV prior in (11) was discretized by
Ψ (ut ) =
Npix
k=1
(∇uts )kTDk (∇uts )k + β,
where the gradient was calculated with the forward difference
approximation
(∇us )k = (usk − uis , usk − usj )T,
where i corresponds to the neighboring pixel of pixel k in
x -direction and j to the neighboring pixel in y-direction.
In both the simulated and experimental case, the number
of iteration steps in the TV denoising (16) was S = 10.
3.4.1 Simulated Test Case
The reconstructions were computed using resolution Npix =
642. The simulated raw k-space data were transformed to the
Radon transform formalism by 1D Fourier transform of each
simulated spoke. Therefore, the matrix H in (2) corresponds
to a matrix implementing the Radon transform (or parallel
ray x-ray tomography transform).
The value of the variance of the observation noise, σv2 =
1.33 × 10−6, was set equal to the variance of the simulated
(19)
(20)
noise that was added to the measurements. Process noise
variance, σw2 = 1 × 10−5, was chosen by test runs, such that
the reconstruction errors were smallest, and kept constant
with respect to time. A complex-valued image using the first
51 spokes reconstructed with the LS approach was used as
the initial value for the KF. The initial value of the error
covariance P0 was diagonal with the values being 0.01% of
the variance of the initial estimate. This was selected so that
the diagonal values would be close to the final values, or
steady-state values, that can be obtained either by iterating
Eqs. (5b), (5c) and (5e) which can be done without the data,
or by running a test run of the whole filter with the data to be
reconstructed. In the Kalman smoother, the number of stored
smoother gains were h = 153 (three full k-space cycles)
and the number of time points that were skipped after each
smoother calculation was b = 3. With these values, a good
compromise was achieved with computational efficiency and
temporal smoothness.
The parameters of the structural priors were selected
manually (by test runs) for visually optimal image quality. The
parameter values for the structured smoothness
regularization (AKF) in Eqs. (6–8) were α = 0.02 and C = 0.01
The parameter γ in the structured TV denoising (16) was
also selected manually for (visually) optimal image quality,
leading to choice γ = 0.25.
We remark that the structured TV denoising was applied
separately for the real and imaginary parts of the Kalman
filter estimate (5d), using the real and imaginary parts of the
anatomical image as the reference images, respectively. Same
value of γ was used for the denoising of real and imaginary
parts of the KF estimate.
3.4.2 Experimental Data
The radial fMRI data were transformed to the Radon
transform formalism by 1D Fourier transform of each measured
k-space spoke.
The variance of the measurement noise was not available
from a calibration or phantom measurement, but was selected
manually, leading to choice σ 2
v = 5 × 10−2. The variance
of the state noise σw2 was set to the same value as in the
simulated test case. The initial values for the estimate and
error covariance were also chosen similarly to the simulated
case.
The regularization parameter of the structural prior was
selected manually for (visually) optimal estimates. The
regularization parameter in (6) was α = 1, while the edge
weighting parameter C was the same as in the simulated case.
In the TV-KF, different values of γ were employed for the
denoising of the real and imaginary parts of the Kalman filter
estimate; for the real part, the value was γr = 4.7 and for
the imaginary part γi = 3.2. Different values were selected
due to different levels of noise in the real and imaginary parts
Original
LS
TV-KF
TV-KS
SW
AKF
AKS
LS
TV-KF
TV-KS
SW
AKF
AKS
of the images. All other parameters were the same as in the
simulated case.
3.5 Image Fidelity Measures
In the simulated test case, the fidelity of the reconstructed
images with respect to the true target is assessed using the
relative L 2 norm error
Err(f ) =
f − f true ,
f true
peak signal-to-noise ratio (PSNR) [
45
] and mean
structural similarity index (SSIM) [
46
] over the whole image
domain. The reconstruction quality of the simulated
activation response inside the ROI is assessed using the L 2 norm
error and contrast-to-noise ratio
C N R = | A − Abase| ,
σv
(21)
where A is the peak (or minimum if the signal decreases)
signal amplitude after the stimulus, Abase is the mean
amplitude of the baseline signal (signal before stimulus), and σv
is the STD of the initial baseline. Other means of calculating
the CNR exist for fMRI, but (21) is commonly used [
47
].
1000
1500
Time point
1000
1500
Time point
2000
2500
2000
2500
1000
1500
Time point
2000
2500
KF, TV-KS, AKF and AKS. Bottom relative error for the entire image
Ω for LS, SW, TV-KF, TV-KS, AKF and AKS
Table 1 Integrated values of the relative L2 norm error in the ROI and
the whole domain Ω from Fig. 3, mean SSIM and PSNR values for Ω,
and mean CNR values for the ROI in the simulated case
LS
SW
TV-KF
TV-KS
AKF
AKS
ROI
Ω
L2
CNR
L2
SSIM
PSNR
Fig. 4 Reconstructions from simulated data. Rows from top to bottom
1 Left side shows the true target image with the temporal variations of
the indicated column of pixels shown on the right. The vertical line
on the temporal variation image shows the time of the image. 2 Same
as above, but with SW reconstruction. 3 LS reconstruction. 4 TV-KF
reconstruction. 5 TV-KS reconstruction. 6 AKF reconstruction. 7 AKS
reconstruction. All images have the color scale of the original image
of one column of image pixels for the entire time series.
The sliding window estimate is 50 time points shorter in
length than the KF and KS estimates, since the SW time
series can be started from the first frame of n = 51
spokes.
Figure 3 and Table 1 reveal that the state estimation
estimates have smaller reconstruction errors, both in the whole
domain and in the ROI, compared to the LS approach or the
sliding window method. Also, the other fidelity measures
for the state estimation approaches are better than those for
Algorithm 1 Kalman filter and smoother algorithms
initialize fˆ0+ = E (f0), P0+ = E [(f0 − fˆ0)(f0 − fˆ0)T], ut0 = fˆt+, γ > 0,
0 < h ≤ T , δ = 0, b ≥ 0
1: for t ← 0, . . . , T do
2: Compute the KF estimate by using (5a), (5b), (5c), (5d), and (5e)
by using (4). Compute the smoother gain with (17a) and store it
3: for s ← 0, . . . , S do
4: Set ut0 = fˆt+.
5: Compute the TV prior with (19) and then calculate the prior
enhanced estimate with (16).
6: end for
7: if t ≥ h ∧ δ ≥ b || t = T then
8: Compute the KS estimate with (17b)
9: δ = 0
10: else
11: δ = δ + 1
12: end if
13: end for
the LS or SW. Based on the fidelity measures, the Kalman
smoothers give the best estimates with the TV-KS yielding
slightly better results than AKS.
The rapid and periodic changes shown in Fig. 3 with
the TV-KF and AKF are caused by the periodic collection
of the radial data, with the largest drop occurring when
moving from the last spoke of a circle to the first spoke
of the next circle. However, as can be seen, the backward
run by the smoother removes this periodic oscillation from
the error. Furthermore, the smoother estimates also produce
Algorithm
algorithms
initialize fˆ+
0 = E (f0), P0+ = E [(f0 − fˆ0)(f0 − fˆ0)T], α > 0, 0 < h ≤
T , δ = 0, b ≥ 0
1: for t ← 0, . . . , T do
2: Compute the AKF estimate by using (5a), (5b), (5c), (5d), and
(5e) by using (9). Compute the smoother gain with (17a) and store
it
3: if t ≥ h ∧ δ ≥ b || t = T then
4: Compute the KS estimate with (17b)
5: δ = 0
6: else
7: δ = δ + 1
8: end if
9: end for
temporally smoother estimates when compared to the other
methods.
The computational times of one time point with the current
implementations were on average 0.0354 s for the TV-KF,
0.4312 s for the TV-KS, 0.4632 s for AKF and 0.8552 s for
AKS on GTX Titan X and 0.0605 s for LS/SW on Intel Xeon
E5507.
4.2 Experimental Data
Figure 5a shows the mean of the ROI, indicated by a line in
Fig. 1, as a function of time. Figure 5b shows a close-up of
Fig. 6 Reconstructions from experimental data. Rows from top to
bottom 1 Left side shows the LS reconstruction with the temporal variations
of the indicated row of pixels shown on the right. The vertical line on
the temporal variation image shows the time of the image. 2 Same as
above, but with SW reconstruction. 3 TV-KF reconstruction. 4 TV-KS
reconstruction. 5 AKF reconstruction. 6 AKS reconstruction
the activation region of Fig. 5a, the time interval shown with
the black vertical lines in the top figure. Figure 6 shows a
snapshot of the reconstructions together with time series of
one row of pixels.
The results in Figs. 5a, b and 6 conform to the findings of
the simulated test case; the state estimation approach yields
better image quality than the LS or SW. The spatial details
in the state estimates are more clear, and the images are less
noisy compared to the conventional approaches. The KF
estimates also produce temporally much smoother estimates.
All the methods detect well the nicotine induced cortical
response in the ROI. However, the response is detected
Table 2 Mean CNR values calculated with (21) for the ROI in
experimental case
ROI
LS
5 Discussion and Conclusions
In this paper, we proposed a state estimation approach for
fMRI. The proposed method utilizes a structural prior from
the anatomical image which is acquired as part of the fMRI
measurement protocol. Two different constructions of the
structural prior were considered, the first one being a
structured smoothness regularization obtained by augmenting the
observation equation with a structurally weighted
regularization matrix and the second a heuristic approach where the
Kalman filter estimates at each time step were denoised by a
structurally guided TV denoising.
The state estimation approach was evaluated using
simulated and experimental small animal imaging fMRI data from
a rat brain and is shown to improve reconstruction quality
compared to a conventional frame-by-frame reconstruction
or the sliding window method. Overall, the Kalman smoother
using the structurally guided TV denoising provided the
best estimates of all the methods. While the approach using
the structurally guided TV denoising is heuristic and does
not have theoretical convergence results, the finding that
it produced the best estimates in a computationally
efficient form for all test cases considered suggests that it
could be a feasible choice in practical applications. The
augmented KF, on the other hand, has a sound theoretical
justification in the statistical Bayes filtering framework with
the algorithm converging to the means and covariances of
evolution updating density and spatially regularized
observation updating density in a Gaussian Bayes filtering problem
[
20, 37
].
We also experimented the KF without any spatial
regularization and with non-structured regularization. When no
regularization was employed, the reconstructed images were
highly noisy and clearly inferior compared to augmented KF
or the TV-KF, implying that the use of spatial regularization
in state estimation approach for fMRI is highly beneficial.
We also studied structured regularization against their
nonstructured versions, finding out that the structurally guided
regularization improves the details of the images compared
to non-structurally guided regularization. To compare the KF
estimates with a sparsity promoting regularized least squares
estimation approach, we computed regularized LS estimates
using the structured TV functional as the regularization
functional and the same data that were used in the LS and SW
approaches. While the regularization improved the
reconstruction fidelity over the non-regularized LS or SW, the
Kalman filter estimates were more accurate than the TV
regularized reconstructions.
The construction of the structural prior is based on the
assumption that the reference image, and the fMRI images are
co-registered. In small animal imaging this seems reasonable,
since the animal is anesthetized during the experiment and
kept in fixed position with a support inside the MRI scanner.
If needed, the co-registration of the images can be verified by
checking the co-registration of the anatomical image against
a standard reconstruction of the fMRI. Also, recent results
in [
23
] indicate that the structurally guided regularization
can be an improvement over conventional regularization even
in cases where there is uncertainty about locations of the
boundaries.
The structurally guided weighting for the structurally
guided TV denoising in the TV-KF algorithm can be
constructed either by using the magnitude of the anatomical
image, or real and imaginary parts of the anatomical image
separately for the real and imaginary parts of the unknown,
respectively. Whenever the fMRI and anatomical reference
are acquired using similar measurement and k-space
sampling protocol, either choice produces roughly the same
results. However, when different samplings are employed,
such as the experimental case in this study with
Cartesian sampling for the reference and radial sampling for
fMRI, the magnitude image has to be used for the
weighting due to different spatial features of the real and imaginary
parts of the images obtained with different sampling
protocols.
In this work, the state evolution was evaluated as a
simple random walk process. A topic of future research is the
improvement of the methods by using more advanced state
evolution models for encoding temporal prior information.
These could be, for example, based on simple models of the
physiological signals [
33
], or if slow changes are expected,
based on kinematic models [
41
] where rate of change of
the signal derivatives are controlled by using a higher-order
time series model. These approaches, however, have the
complication that the computational demand increases due to
increased number of unknown state variables.
In this study, the state estimation approach was evaluated
using time steps of single k-space spoke. The time
resolution and sampling can in principle be selected quite freely;
for example, one could use more than one spoke per time
instant or use the data in the sliding window manner. For
radial sampling, an interesting choice could be to sample the
spokes using the golden angle rule [
48
] instead of sampling
the spokes in small incremental angles as this would make the
sampling of the k-space more uniform and possibly prevent
the periodic component of the reconstruction errors which
can now be seen with TV-KF and AKF.
Acknowledgements This work was supported by the Academy of
Finland (Project 250215 Finnish Center of Excellence in Inverse Problems
Research) and by Jane and Aatos Erkko foundation (Project 64741).
The authors wish to thank Nvidia for the Tesla K40c used in this study
through the Hardware Grant Program.
Open Access This article is distributed under the terms of the Creative
Commons Attribution 4.0 International License (http://creativecomm
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit
to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made.
Ville-Veikko Wettenhovi is a
Ph.D. student at the University of
Eastern Finland. He received his
M.Sc. in medical physics from
the University of Eastern Finland
in 2015. His current research
interests are state estimation
methods, medical imaging,
especially MRI and PET imaging,
and time-varying reconstruction.
Ville Kolehmainen is a Profes
sor of Computational Physics at
the University of Eastern
Finland. He received Ph.D. degree
in applied physics from the
University of Kuopio in 2001.
His research interests includes
computational inverse problems
with applications to tomography
imaging problems.
Joanna Huttunen is a
postdoctoral researcher in A.I.
Virtanen Institute for Molecular
Sciences at the University of
Eastern Finland (UEF). She received
her Ph.D. degree in biomedical
imaging in 2013 from the UEF.
Her research interests include
functional MRI and
electrophysiological measurements in
animal models.
Mikko Kettunen is a Research
Director at the A.I. Virtanen
Institute for Molecular Sciences,
University of Eastern Finland.
He received his PhD degree in
Biomedical NMR from the
University of Kuopio, Finland in
2002. From 2003 to 2013, he
worked as post-doctoral/senior
scientist at the Cancer Research
UK Cambridge Institute,
University of Cambridge, UK. His
research interests include
preclinical brain and cancer MRI.
Olli Gröhn is a Professor of
Biomedical NMR and Director
of Biomedical Imaging Unit that
serves as national infrastructure
for experimental MRI in Finland.
He received PhD in biomedical
NMR from University of Kuopio
in 2000. He has 20-year
experience in development of MRI
technique. His research group is
focusing on technical
development of novel MRI techniques,
which can be applied in several
different neurological diseases.
The special interest areas
currently include fMRI and resting state fMRI combined with brain
stimulation.
1. Bernstein , M.A. , King , K.F. , Zhou , X.J.: Handbook of MRI Pulse Sequences . Elsevier, Burlington, MA ( 2004 )
2. Buxton , R.B.: Introduction to Functional Magnetic Resonance Imaging: Principles and Techniques, 2nd edn . Cambridge University Press, Cambridge ( 2009 )
3. Byrne , C. : Block-iterative interior point optimization methods for image reconstruction from limited data . Inverse Probl . 16 ( 5 ), 1405 - 1419 ( 2000 )
4. Candes , E. , Wakin , M.: An introduction to compressive sampling . IEEE Signal Process. Mag . 25 ( 2 ), 21 - 30 ( 2008 )
5. Chan , C. , Fulton , R. , Barnett , R. , Feng , D.D. , Meikle , S. : Postreconstruction nonlocal means filtering of whole-body PET with an anatomical prior . IEEE Trans. Med. Imaging 33 ( 3 ), 636 - 650 ( 2014 )
6. d'Arcy , J.A. , Collins, D.J. , Rowland , I.J. , Padhani , A.R. , Leach , M.O. : Applications of sliding window reconstruction with Cartesian sampling for dynamic contrast enhanced MRI . NMRI Biomed . 15 ( 2 ), 174 - 183 ( 2002 )
7. Ehrhardt , M.J. , Markiewicz , P. , Liljeroth , M. , Barnes , A. , Kolehmainen , V. , Duncan , J.S. , Pizarro , L. , Atkinson , D. , Hutton , B.F. , Ourselin , S. , Thielemans , K. , Arridge , S.R.: PET reconstruction with an anatomical MRI prior using parallel level sets . IEEE Trans. Med. Imaging 35 ( 9 ), 2189 - 2199 ( 2016 )
8. Feng , X. , Salerno , M. , Kramer , C.M., Meyer, C.H.: Kalman filter techniques for accelerated Cartesian dynamic cardiac imaging . Magn. Reson. Med . 69 ( 5 ), 1346 - 1356 ( 2013 )
9. Gamper , U. , Boesiger , P. , Kozerke , S. : Compressed sensing in dynamic MRI . Magn. Reson. Med . 59 ( 2 ), 365 - 373 ( 2008 )
10. Gössl , C. , Auer , D.P. , Fahrmeir , L. : Dynamic models in fMRI . Magn. Reson. Med . 43 ( 1 ), 72 - 81 ( 2000 )
11. Gozzi , A. , Schwarz , A. , Reese , T. , Bertani , S. , Crestan , V. , Bifone , A. : Region-specific effects of nicotine on brain activity: a pharmacological MRI study in the drug-naive rat . Neuropsychopharmacology 31 ( 8 ), 1690 - 1703 ( 2006 )
12. Guerquin-Kern , M. : Matlab Framework for MRI Simulation and Reconstruction ( 2014 ). http://bigwww.epfl.ch/algorithms/ mri-reconstruction/
13. Guerquin-Kern , M. , Lejeune , L. , Pruessmann , K. , Unser , M. : Realistic analytical phantoms for parallel magnetic resonance imaging . IEEE Trans. Med. Imaging 31 ( 3 ), 626 - 636 ( 2012 )
14. Holland, D.J., Liu, C. , Song , X. , Mazerolle , E.L. , Stevens , M.T. , Sederman , A.J. , Gladden , L.F. , D'Arcy , R.C.N. , Bowen , C.V. , Beyea , S.D.: Compressed sensing reconstruction improves sensitivity of variable density spiral fMRI . Magn. Reson. Med . 70 ( 6 ), 1634 - 1643 ( 2013 )
15. Huttunen , J. , Gröhn , O. , Penttonen , M. : Coupling between simultaneously recorded bold response and neuronal activity in the rat somatosensory cortex . NeuroImage 39 ( 2 ), 775 - 785 ( 2008 )
16. Jackson , J.I. , Meyer, C.H. , Nishimura , D.G. , Macovski , A. : Selection of a convolution function for Fourier inversion using gridding . IEEE Trans. Med. Imaging 10 ( 3 ), 473 - 478 ( 1991 )
17. Jeromin , O. , Pattichis , M.S. , Calhoun , V.D.: Optimal compressed sensing reconstructions of fMRI using 2D deterministic and stochastic sampling geometries . Biomed. Eng. Online 11 ( 25 ), 1 - 36 ( 2012 )
18. Jung , H. , Sung , K. , Nayak , K.S. , Kim , E.Y. , Ye , J.C. : k-t FOCUSS: a general compressed sensing framework for high resolution dynamic MRI . Magn. Reson. Med . 61 ( 1 ), 103 - 116 ( 2009 )
19. Kaipio , J.P. , Karjalainen , P.A. , Somersalo , E. , Vauhkonen , M. : State estimation in time-varying electrical impedance tomography . Ann. N. Y. Acad. Sci . 873 ( 1 ), 430 - 439 ( 1999 )
20. Kaipio , J.P. : Statistical and Computational Inverse Problems. Springer, New York ( 2005 )
21. Kalman , R.E.: A new approach to linear filtering and prediction problems . J. Basic Eng . 82 ( 1 ), 35 - 45 ( 1960 )
22. Kervinen , M. , Vauhkonen , M. , Kaipio , J.P. , Karjalainen , P.A. : Time-varying reconstruction in single photon emission computed tomography . Int. J. Imaging Syst. Techol . 14 ( 4 ), 186 - 197 ( 2004 )
23. Kolehmainen , V. , Ehrhardt , M.J. , Arridge , S. : Incorporating structural prior information and sparsity into EIT using parallel level sets . Inverse Probl ( 2016 , under review)
24. Li , L. , Yan , B. , Tong , L. , Wang , L. , Li , J. : Incremental activation detection for real-time fMRI series using robust Kalman filter . Comput. Math. Methods Med . 2014 , 759805 ( 2014 ). doi: 10 .1155/ 2014 /759805
25. Lustig , M. , Donoho , D. , Pauly , J.M. : Sparse MRI : the application of compressed sensing for rapid MR imaging . Magn. Reson. Med . 58 ( 6 ), 1182 - 1195 ( 2008 )
26. Lustig , M. , Donoho , D. , Santos , J. , Pauly , J.: Compressed sensing MRI . IEEE Signal Process. Mag . 25 ( 2 ), 72 - 82 ( 2008 )
27. Majumdar , A. , Ward , R.K. , Aboulnasr , T. : Compressed sensing based real-time dynamic MRI reconstruction . IEEE Trans. Med. Imaging 31 ( 12 ), 2253 - 2266 ( 2012 )
28. Michel , V. , Gramfort , A. , Varoquaux , G. , Eger , E. , Thirion , B. : Total variation regularization for fMRI-based prediction of behavior . IEEE Trans. Med. Imaging 30 ( 7 ), 1328 - 1340 ( 2011 )
29. Otazo , R. , Kim , D. , Axel , L. , Sodickson , D.K. : Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI . Magn. Reson. Med . 64 ( 3 ), 767 - 776 ( 2010 )
30. Paasonen , J. , Salo , R. , Shatillo , A. , Forsberg , M. , Närvänen , J. , Juttunen , J.K. , Gröhn , O. : Comparison of seven different anesthesia protocols for nicotine pharmacologic magnetic resonance imaging in rat . Eur. Neuropsychopharmacol . 26 ( 3 ), 518 - 531 ( 2016 )
31. Paige , C.C. , Saunders , M.A. : LSQR: an algorithm for sparse linear equations and sparse least squares . ACM Trans. Math. Softw . 8 ( 1 ), 43 - 71 ( 1982 )
32. Park , S. , Park , J.: Accelerated dynamic cardiac MRI exploiting sparse-Kalman-smoother self-calibration and reconstruction (k-t SPARKS) . Phys. Med. Biol . 60 , 3655 - 3671 ( 2015 )
33. Prince , S. , Kolehmainen , V. , Kaipio , J.P. , Franceschini , M.A. , Boas , D. , Arridge , S.R. : Time series estimation of biological factors in optical diffusion tomography . Phys. Med. Biol . 48 , 1491 - 1504 ( 2003 )
34. Qranfal , J. , Tanoh , G.: Regularized Kalman filtering for dynamic SPECT . J. Phys.: Conf. Ser . 124 , 012042 ( 2008 )
35. Rudin , L. , Osher , S. , Fatemi , E.: Nonlinear total variation based noise removal algorithms . Physica D 60 ( 1-4 ), 259 - 268 ( 1992 )
36. Särkkä , S. , Solin , A. , Nummenmaa , A. , Vehtari , A. , Auranen , T. , Vanni , S. , Lin , F.H. : Dynamic retrospective filtering of physiological noise in BOLD fMRI: DRIFTER . NeuroImage 60 ( 2 ), 1517 - 1527 ( 2012 )
37. Sbarbaro , D. , Vauhkonen , M. , Johansen , T.A. : State estimation and inverse problems in electrical impedance tomography: observability, convergence and regularization . Inverse Probl . 31 ( 4 ), 1 - 27 ( 2015 )
38. Simon , D. : Optimal State Estimation: Kalman, H∞ and Nonlinear Approaches . Wiley, New York ( 2006 )
39. Sümbül , U. , Santos , J.M. , Pauly , J.M. : Improved time series reconstruction for dynamic magnetic resonance imaging . IEEE Trans. Med. Imaging 28 ( 7 ), 1093 - 1104 ( 2009 )
40. Sümbül , U. , Santos , J.M. , Pauly , J.M.: A practical acceleration algorithm for real-time imaging . IEEE Trans. Med. Imaging 28 ( 12 ), 2042 - 2051 ( 2009 )
41. Tossavainen , O.P. , Vauhkonen , M. , Kolehmainen , V. , Kim , K.Y.: Tracking of moving interfaces in sedimentation processes using electrical impedance tomography . Chem. Eng. Sci . 61 , 7717 - 7729 ( 2006 )
42. Vaswani , N.: LS-CS-residual (LS-CS): compressive sensing on least squares residual . IEEE Trans. Signal Process . 58 ( 8 ), 4108 - 4120 ( 2010 )
43. Vauhkonen , M. , Karjalainen , P.A. , Kaipio , J.P.: A Kalman filter approach to track fast impedance changes in electrical impedance tomography . IEEE Trans. Biomed. Eng . 45 ( 4 ), 486 - 493 ( 1997 )
44. Vunckx , K. , Atre , A. , Baete , K. , Reilhac , A. , Deroose , C.M. , Laere , K.V. , Nuyts , J. : Evaluation of three MRI-based anatomical priors for quantitative PET brain imaging . IEEE Trans. Med. Imaging 31 ( 3 ), 599 - 612 ( 2012 )
45. Wang , Z. , Bovik , A.C. : Mean squared error: love it or leave it? IEEE Signal Process . Mag. 26 ( 9 ), 98 - 117 ( 2009 )
46. Wang , Z. , Bovik , A.C. , Sheikh , H.R. , Simoncelli , E.P.: An optimal radial profile order based on the golden ratio for time-resolved MRI . IEEE Trans. Image Process . 13 ( 4 ), 600 - 612 ( 2004 )
47. Welvaert , M. , Rosseel , Y. : On the definition of signal-to-noise ratio and contrast-to-noise ratio for fMRI data . PLoS One 8 ( 11 ), e77089 ( 2013 )
48. Winkelmann , S. , Schaeffter , T. , Koehler , T. , Eggers , H. , Doessel , O. : An optimal radial profile order based on the golden ratio for time-resolved MRI . IEEE Trans. Med. Imaging 26 ( 1 ), 68 - 76 ( 2007 )
49. Yan , S. , Nie , L. , Wu , C. , Guo , Y. : Linear dynamic sparse modelling for functional MR imaging . Brain Inform . 1 ( 1 ), 11 - 18 ( 2014 )
50. Yendiki , A. , Reuter , M. , Wilkens , P. , Rosas , H.D. , Fischl , B. : Joint reconstruction of white-matter pathways from longitudinal diffusion MRI data with anatomical priors . NeuroImage 127 , 277 - 286 ( 2016 )
51. Zong , X. , Lee , J. , Poplawsky , A.J. , Kim , S.G. , Ye , J.C. : Compressed sensing fMRI using gradient-recalled echo and EPI sequences . NeuroImage 92 , 312 - 321 ( 2014 ) Marko Vauhkonen received his PhD in physics in 1997 at the University of Kuopio, Finland. After graduation, he worked as a researcher and research director at the same university , until moving to Germany in 2006 . There he worked as a MarieCurie Research Fellow for two years at the Philips Research GmbH, Aachen . During 2008 - 2009 Vauhkonen worked as a CTO in a spin-off company Numcore Ltd. until starting as a professor at the University of
Kuopio (currently University of Eastern Finland), Department of Applied Physics, in 2009 . His research interests include inverse problems, time-varying reconstruction, process tomography and medical imaging such as PET, SPECT and MRI. He has published around 100 scientific journal articles .