State Estimation with Structural Priors in fMRI

Journal of Mathematical Imaging and Vision, Jul 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10851-017-0749-x.pdf

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 .


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10851-017-0749-x.pdf

Ville-Veikko Wettenhovi, Ville Kolehmainen, Joanna Huttunen, Mikko Kettunen, Olli Gröhn, Marko Vauhkonen. State Estimation with Structural Priors in fMRI, Journal of Mathematical Imaging and Vision, 2017, 1-15, DOI: 10.1007/s10851-017-0749-x