Ensemble Kalman Filters for reliability estimation in perfusion inference
IInternational Journal for Uncertainty Quantification , x(x): 1–19 (2017October 23, 2018)
ENSEMBLE KALMAN FILTERS FOR RELIABILITYESTIMATION IN PERFUSION INFERENCE
Peter Zaspel ˚ *Address all correspondence to: Peter Zaspel, University of Basel, 4051 Basel, Switzerland, E-mail: [email protected] Original Manuscript Submitted: mm/dd/yyyy; Final Draft Received: mm/dd/yyyyWe consider the solution of inverse problems in dynamic contrast–enhanced imaging by means of Ensemble KalmanFilters. Our quantity of interest is blood perfusion, i.e. blood flow rates in tissue. While existing approaches to com-pute blood perfusion parameters for given time series of radiological measurements mainly rely on deterministic,deconvolution–based methods, we aim at recovering probabilistic solution information for given noisy measurements.To this end, we model radiological image capturing as sequential data assimilation process and solve it by an EnsembleKalman Filter. Thereby, we recover deterministic results as ensemble–based mean and are able to compute reliability in-formation such as probabilities for the perfusion to be in a given range. Our target application is the inference of bloodperfusion parameters in the human brain. A numerical study shows promising results for artificial measurementsgenerated by a Digital Perfusion Phantom.
KEY WORDS:
Medical Imaging, Stochastic Modeling, Inverse Problems, Ensemble Kalman Filter, Dy-namic Contrast-Enhanced Imaging, Perfusion, Inference
1. INTRODUCTION
Medical imaging by x-rays, magnetic resonance imaging (MRI) and computed tomography (CT) has considerablychanged medical diagnosis throughout the last decades. Often, contrast agents , i.e. specific liquid chemicals, areinjected into the patients blood circulation during the imaging process. This leads to contrast–enhanced images ofhigher contrast in some regions of the human body. In this work, we study inverse problems for a specific MRI or CTimaging task. That is, we aim at recovering a quantity of interest in medical imaging, which is derived by dynamic contrast–enhanced (DCE) imaging. In dynamic contrast–enhanced imaging, a time-dependent series of radiologicalimages of a part of the patients body (e.g. the brain) is taken immediately after injecting a contrast agent into thepatient’s blood circulation. By observing the time-dependent concentration evolution of the contrast agent inside thepatient’s tissue, it is possible to recover information about the blood flow rates, i.e. the perfusion .The outcome of the image acquisition process is a time-discrete series of tree-dimensional (space-discrete) con-centration images c of a part of the patient body. The actual perfusion evaluation is a post-processing step, beingpreceded by image de-noising and motion compensation. Currently, blood perfusion is computed independently perdiscrete tissue volume element, i.e. voxel , thus spatial information is mostly neglected. Variants of the indicator-dilution theory [1–3] describe the concentration of a contrast agent in tissue at a given point in time as the result of aconvolution in time of the (known) time-dependent arterial or blood circulation inflow concentration c art with an un-known tissue-dependent kernel function k . Blood perfusion is computed as a weighted maximum or point evaluationof the unknown kernel function.Current state of the art methods aim at recovering the unknown time-dependent kernel function for given dis-crete measurements of the contrast agent in tissue. The kernel function is either modeled as a parametrized analytic function [2,4,5] or discretized as a fully unknown function [1,2,6]. Then, most approaches rely on a deterministic reconstruction of the kernel function, involving the solution of a deconvolution problem with regularization. a r X i v : . [ m a t h . NA ] O c t P. Zaspel
One drawback of the use of motion compensation, de-noising and deterministic deconvolution lies in the loss ofinformation on the quality of a computed solution. That is, probabilistic information about the measurement accuracyand errors in space and time with their influence on the exactness on the computed quantity of interest are neglectedor even lost.In this work, we propose an approach to infer perfusion in the discussed application case while keeping the prob-abilistic information on the solution. Thereby, we overcome the discussed drawback of knowledge loss. To achievethis, we model the inference problem as a sequential data assimilation problem: First, the unknown kernel function k is described as an unknown system state, for which a predictive time-discrete stochastic system state model isintroduced. In this model, the kernel function is represented as a random variable. Then, a time-discrete stochasticobservation model describes the relationship of the current approximation of k and the noisy measurements deliv-ered by medical imaging. Finally, the well-known Ensemble Kalman Filter (EnKF) [7–10] is used to compute anensemble–based approximation of the posterior probability density function (PDF) of k given the system state modeland the (noisy) measurements. Based on this PDF, means, cumulative distribution functions, etc. can be computed.The whole sequential data assimilation methodology is applied to (noisy) artificial measurement data generated by a Digital Perfusion Phantom [11–13], i.e. a forward model describing the mapping of perfusion information to med- ical images. Note that we stick to the use of a Digital Perfusion Phantom, instead of using exemplary patient datasince, first, it allows to artificially create arbitrary amounts of radiological images and, second, a scanner- and patient-independent way to analyse perfusion estimation methods is highly desired in radiology. Certainly, applying theproposed methodology (together with radiologists) to real patient data is future work.The data assimilation problem that we will model in Section 3 will stick to Gaussian random fields and a linearforward model. This is a strong simplification, allowing to analytically solve the data assimilation problem by theoriginal
Kalman filter [14], However, we use the EnKF, which is a generalization of this Kalman filter, usually beingapplied in a non-linear, non-Gaussian setting. For more details on the connection between EnKF and the Kalman filterand a convergence analysis in the linear Gaussian case see [15,16], while EnKFs for inverse problems are discussedin [9,17,18] and many extensions and alternatives for the EnKF are, e.g., developed in [19–24]. We decided to usethe EnKF here, since we consider this work as a starting point for much more involved approaches for the predictionof perfusion. In fact, the rather simple linear forward model from the indicator dilution theory should be replaced bymore complex or even PDE-based models, which will certainly no longer be linear. Moreover, we expect that a moreinvolved evolution model, cf. Section 3.3.1, with non-Gaussian noise might become valid in real application cases.Therefore, we here already introduce the more involved EnKF framework, while considering other forward modelsand non-Gaussian noise as future work.To the best of our knowledge, we consider the discussed work to be a new contribution to the field. Nevertheless,there has been previous work on the use of Ensemble Kalman Filters in the application scenario. In [25], the authorsconcentrate on the introduction of a tissue model that includes space-dependent information. To achieve this, a bloodflow model is combined with an EnKF. Preliminary results for this approach are given. In contrast, we focus heredirectly on the mathematical setting based on the indicator-dilution theory that is well-known and, therefore, wellaccepted by radiologist. Hence, our methodology is considered as an extension to the existing standard methodologyintroducing the opportunity to derive statistical information on the computed solution. In addition to the differentobjective compared to [25], we also perform a large number of parameter studies and convergence tests, which arecrucial to understand the properties of the method.This article is organized as follows. In Section 2, we give a mathematical model for the radiological imagingand perfusion extraction mechanism. Section 3 outlines our numerical approach based on sequential data assimilationusing EnKF. Numerical results are given in Section 4 while Section 5 summarizes the discussed work.
2. MODELING RADIOLOGICAL IMAGING AND PERFUSION EXTRACTION
In the following, we start by giving an abstract model for the transport of contrast agent. Then, measurements bye.g. MRI are abstractly modeled. A concrete model for the contrast agent distribution is given by the indicator- dilution theory. Finally, our quantity of interest, i.e. blood perfusion, is introduced and the deterministic inference problem is summarized.
International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference Shall D tiss Ă R be the tissue domain in the human body for which we want to derive information by dynamiccontrast–enhanced imaging. We study contrast agent transport / concentrations in a time interval r , T s Ă R with T being the final time. The inflow concentration of the contrast agent (at some arterial inlet) is a function c art : r , T s Ñ R ě . The time-continuous contrast agent concentration in tissue can be modeled as a function c : D tiss ˆ r , T s Ñ R ě . Both are related to each other by an (unknown) operator B , with c p¨ , t q “ B r c art sp t q (1)that models the function of the human body with respect to contrast agent transport. Appropriate measurement devices (CT, MRI, . . . ) usually have a cuboidal measurement domain. Therefore, we startby limiting D tiss to D meas “ Ś d “ r , a d s with a “ p a , a , a q J P R describing the size of the measurement domain. For simplicity, we assume the measurement domain and the area of interest to match exactly, i.e. D meas “ D tiss , excluding cases in which some part of the measurement domain does not contain valid tissue. Moreover, D meas is simplified as being stationary in time, i.e. the measurement device (or the patient) does not move or movements areconsidered as measurement error.The finite spatial resolution N D P N of the measurement device leads to a decomposition of D meas into N voxel “ ś d “ N p d q D volume elements or voxels of volume V voxel “ ś d “ a d { N p d q D for which we obtain aver-aged (constant) measurements. We introduce a measurement operator Ψ that gives for a given exact contrast agentconcentration c and a chosen point in time t P r , T s a measurement vector c p t q P R N voxel ě as c p t q “ Ψ r c sp t q : “ Θ r c sp t q ` E r c sp t q . Here, Θ is a noise-free measurement-operator and is usually a volumetric average over each voxel being equivalentto a piece-wise constant approximation in space. E abstractly models a (potentially non-linear) additive error (noise,movements, technical problems, . . . ).To reflect time-discrete measurements, we introduce N obs ordered, pair-wise different discrete observation times t obsi P r , T s , i P t , . . . , N obs u , at which measurements or observations are done, giving the observation matrix C : “ p c ji q j “ ,...,N voxel ,i “ ,...,N obs composed of observation vectors c i as C “ p c | . . . | c N obs q with c i “ Ψ r c sp t obsi q : “ Θ r c sp t obsi q ` E r c sp t obsi q . (2) The indicator-dilution theory (IDT) [1] provides a model for the time evolution of the contrast agent concentration ina reference voxel D voxel Ă D tiss with volume V voxel , given the arterial inflow c art . While, in this standard model,the contrast agent’s concentrations are assumed to be constant in each voxel, we first want to formulate the IDT asa space-continuous model and then move over to a discrete description as consequence of a measurement process.Our continuous version of the indicator-dilution theory–based transport model replaces B in eq. (1) with the modeloperator B IDT given via c p x , t q “ B IDT r c art , k sp t q : “ ż T c art p τ q k p x , t ´ τ q d τ , p x , t q P D tiss ˆ r , T s . (3) Kernel k : D tiss ˆ r , T s Ñ R fully characterizes the properties of the tissue at point x . In order to have an well- defined integrand, we assume k p¨ , t q “ t ă
0. Note that the model operator B IDT is actually independent of the spatial position.
Volume x, Issue x, 2017October 23, 2018
P. Zaspel
We now apply the measurement operator Ψ to eq. (3) obtaining Ψ r c sp t q “ Θ r B IDT r c art , k ss p t q ` E r B IDT r c art , k ss p t q“ ż T c art p τ q k p t ´ τ q d τ ` E r B IDT r c art , k ss p t q , where k “ p k , . . . , k N voxel q J is a vector of univariate kernel functions k j : r , T s Ñ R . Since we are interested intime-discrete observations, we limit our discussion to observation times t obsi yielding c i “ Ψ r c sp t obsi q “ ż T c art p τ q k p t obsi ´ τ q d τ ` e i p c art , k q , with the abbreviation e i p c art , k q : “ E r B IDT r c art , k ss p t obsi q . For a single voxel j P t , . . . , N voxel u , we obtain c obsj,i “ ż T c art p τ q k j p t obsi ´ τ q d τ ` e j,i p c art , k q . In case of e j,i p c art , k q “
0, this boils down to the classical indicator-dilution-theory model given on a reference voxel j . Obviously, this model is independent of the spatial position of the voxel j . The classical theory further introducesa mean density ρ j P R ě in a voxel j , which becomes of interest in the following subsection. The inference task discussed in this article is to compute a time-stationary perfusion (blood flow) information p P R N voxel given the (assumed to be exactly known) inflow concentration c art and the observation matrix C . Formally,the blood perfusion in a given voxel j can be evaluated as quantity of interest of the computed response function k j as p j : “ p p k j q : “ ρ j k j p q . From a mathematical point of view, this quantity has nice properties, since it is just a point evaluation of the responsefunction. In practice [1], perfusion is however often evaluated as ˜ p j : “ ˜ p p k j q : “ ρ j max t Pr ,T s k j p t q . For simplicity and since we use just artificial input data, we stick to the first version of this quantity of interest.
To summarize this section, we formulate the deterministic problem that we aim to solve: For given measurement time T P R , arterial inflow c art : r , T s Ñ R ě , measurement/observation times t obs ă t obs ă . . . ă t obsN obs and obser-vation matrix C or vectors c i P R N voxel , i P t , . . . , N obs u , we aim at computing a vector k “ p k , . . . , k N voxel q T of kernel functions k j : r , T s Ñ R and the derived quantity of interest p “ p p , . . . , p N voxel q J with p j “ k j p q{ ρ j such that c j,i « ż T c art p τ q k j p t obsi ´ τ q d τ ` e j,i p c art , k q , j P t , . . . , N voxel u , i P t , . . . , N obs u . (4)Clearly, this problem is underdetermined with the given requirements. Furthermore, we have not specified the natureof the error term, yet. This is why we used the notion “ « ”. A much clearer idea of the concept of a solution to this problem is given in the next section, where we reformulate the problem as Bayesian sequential data assimilation problem. International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference
3. NUMERICAL APPROACH BY SEQUENTIAL DATA ASSIMILATION
In this section, we first introduce a discretization for the model discussed in the last section. This is necessary, sincewe will use its discretized version in context of sequential data assimilation, afterwards. An approximation to thesolution of the assimilation problem is derived by the Ensemble Kalman Filter that is briefly introduced as final partof this section.
We start by discretizing eq. (4) for fixed i P t , . . . , N obs u and fixed j P t , . . . , N voxel u . Numerical quadrature usinga rectangular rule gives c j,i « ∆ τ N q ´ ÿ q “ c art p τ q q k j p t obsi ´ τ q q ` e j,i p c art , k q , with N q equidistant abscissas τ q : “ q ¨ ∆ τ and ∆ τ : “ TN q . In the original problem setting, the observation times t obsi can be chosen arbitrarily. However, we here introduce a simplification, in which we assume the observation times tobe given for a fixed time step size ∆ t obs . Moreover, this time step size shall be a multiple of ∆ τ , i.e. t obsi : “ i ∆ t obs , ∆ t obs : “ s ¨ ∆ τ , s P N . Thereby, we obtain c j,i « ∆ τ N q ´ ÿ q “ c art p τ q q k j p t obsi ´ τ q q ` e j,i p c art , k q“ ∆ τ N q ´ ÿ q “ c art p q ∆ τ q k j pp i s ´ q q ∆ τ q ` e j,i p c art , k q . Since we now only need c art and k j being evaluated at multiples of ∆ τ , we can replace them by vector c art “ ` c art, , . . . , c art,N q ´ ˘ J such that c art,q : “ c art p q ∆ τ q and matrix K P R N q ˆ N voxel with K : “ p k q,j q q,j such that k q,j : “ k j p q ∆ τ q , yielding c j,i « ∆ τ N q ´ ÿ q “ c art,q k p i s ´ q q ,j ` e j,i p c art , k q . With the extension of k j p t q “ t ă j , i ) a(degenerated) matrix H j,i P R ˆ N q such that c j,i « H j,i k j ` e j,i p c art , k q , (5)where the k j P R N q are the column vectors of matrix K , i.e. K “ p k | . . . | k N q ´ q .Following the nomenclature of [14], we next reformulate the deterministic inference problem from Section 2.5as a sequential data assimilation problem. To this end, we first translate the involved quantities into random variablesas in a Bayesian inference problem. Thereafter, we introduce the basic concepts of sequential data assimilation.Since the problem decouples for all voxels j P t , . . . N voxel u , we keep j fixed for the rest of this section. Let be p Ω , F , P q a probability space. In Bayesian inference we want to gain information on a system state variable for given observation(s). In our context, the state variable is the time-continuous kernel function k j . However, forsimplicity and since we deal with discrete data anyway, we infer the discrete k j P R N q from (5), instead. Therefore, we introduce a new random variable k j : Ω Ñ R N q , (6) Volume x, Issue x, 2017October 23, 2018
P. Zaspel replacing the time-discrete deterministic solution vector k j . ˚ Moreover, we introduce a random variable e j,i p c art q :Ω Ñ R , replacing the error term used before. Note that we assume k j and e j,i to be independent random variables.This is a rather strong simplification, since we initially modeled e j,i p c art , k q to be a potentially non-linear error inthe observation data, which itself is given in the indicator-dilution-theory by the arterial inflow c art and the tissueproperties modeled by kernel k . That is, we – at this point – decouple the error in the observation from the specificpatient tissue. This decoupling is reflected by the new notation e j,i p c art q . Finally, we also consider each observation c j,i as random variable c j,i : Ω Ñ R , which is usually called observed variable . Using eq. (5), c j,i is defined as c j,i p ω q : “ H j,i k j p ω q ` e j,i p c art qp ω q , @ ω P Ω . (7)We will call matrix H j,i (linear) forward map . The aim of inference is to find a reference trajectory k refj , being arealization of k j such that the (measured) observations fit to the observed variable. Sequential data assimilation relies on an evolution model and a forward model to obtain k refj . The models run ondifferent time scales. The evolution model is a stochastic difference equation implying a certain predicted evolution ofthe system state variable over many small time steps. The forward model defines a relationship between the referencetrajectory (which is to be found) and the observed data at the observation times. In context of sequential data assimilation for dynamic processes, it is usually assumed that the coupling between themeasured observations and the system state variable is time-local. That is, a new observation at time t obs only affectsthe system state variable for times t ě t obs . In our application, this is different, since the forward model, i.e. theindicator-dilution theory, is a non-local operator in time. Therefore, we need an evolution model that allows to doglobal updates to the system state variable k j . The probably most simplistic approach to model this type of globalupdates is given by the evolution model k p l ` q j “ k p l q j ` ? ∆ τ n p l q , l P t , . . . , N q ´ u , (8)Here ´ k p q j , . . . , k p N q ´ q j ¯ , is a sequence of random variables of the type given in eq. (6) for time steps l ∆ τ . Weassume k p q j „ N p , σ Σ n q , corresponding to a zero initial guess for the kernel function with Gaussian noise with acovariance matrix Σ n . σ P R is a scaling coefficient. The n p l q s are a sequence of independent identically distributedrandom variables with n p n q : Ω Ñ R N q drawn as n p l q „ N p , Σ n q . Σ n P R N q ˆ N q will be chosen using a Gaussiancovariance kernel such that Σ n : “ p σ l,l q N q ´ l,l “ with σ l,l : “ α e ´ } τ l ´ τ l (cid:96) and a parametrization in the scale α P R and the correlation length (cid:96) P R .Analyzing this evolution model, we can state that it can be understood as Euler-Maruyama-based discretizationof the system of stochastic ordinary differential equations d k j “ d W t , (9)where W t is a vector of correlated univariate Wiener processes. Moreover, we observe that this evolution model, incontrast to the standard setting of dynamical processes, now only takes the role of coupling the time-discrete values in k . This coupling is imposed by the covariance of the noise term. In fact, as we will see in Section 4.5, the correlationlength in the Gaussian covariance kernel will have a regularizing influence on the inferred solution. Note that thechoice of Gaussian noise might lead to a locally negative kernel K , while this kernel is supposed to be positive. The choice of a better noise distribution is future work. ˚ We use sans serif letters to indicate that a given quantity is a random variable.
International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference A rather natural question in context of the proposed application is, whether it would preferable to apply anEnKF-based approach for direct inversion, cf. [9,17,18] to the measurement matrix C , avoiding sequential dataassimilation. In fact, we prefer sequential data assimilation, since it allows to treat the given inference problem as atime-dependent problem. In the real application case of DCE imaging, one objective of researchers is to find meansto effectively control the image capturing process, in terms of a feedback loop. In that context, it is important to beable to continuously monitor the achieved approximation of the perfusion information during the image capturingprocess. Based on that monitoring, one might be able to select the next observation time or the required quality forthe next observation. This, however, cannot be done in direct inversion of the problem. We choose eq. (7) as our forward model, i.e. we get the forward model with respect to the reference trajectory k refj c j,i “ H j,i k refj,i s ` e j,i , i P t , . . . , N obs u . (10)The e j,i are sequences of i.i.d. random variables for growing observation time index i following e j,i „ N p , σ e q , for all i P t , . . . , N obs u , with σ e P R the observation error variance. Note that this choice of the distribution of e j,i is afurther simplification over the simplification that has been made in Section 3.2. There, we decoupled the observationerror from the kernel function k describing the tissue properties. Here, we further decouple the observation errorfrom the arterial inflow and make it a purely data-independent error that is furthermore only modeled as normallydistributed. It is very clear, that this is a very strong simplification. Finding a much better, maybe imaging devicedependent, error is future work.Due to ∆ t obs “ s ∆ τ , k refj,i s is the unknown reference trajectory evaluated at observation time t obsi . To be concise, we here only briefly summarize the general idea of the actual assimilation task with notation from[14]. Further details can be found e.g. in [14].Let π k p i s q j p k j q be the probability density function of the random variable k p i s q j at time t obsi for i P t , . . . , N obs u .Then, sequential data assimilation computes posterior PDFs π k p i s q j p k j | c j, : i q , i P t , . . . , N obs u , i.e. probability density functions of the random variables k p i s q j with an instance k j conditioned to the observa-tions c j, , . . . , c j,i that are instances of c j, , . . . , c j,i . This is done using an iterative approach. It is started with π k p q j p k j | c j, : q being the PDF of k p q j . Then, for a given PDF π k pp i ´ q s q j p k j | c j, : i ´ q , it iteratively1. computes the density π k p i s q j p k j | c j, : i ´ q and thereby solves a prediction problem for the given evolution modeleq. (8),2. applies Bayes theorem π k p i s q j p k j | c j, : i q “ π c j,i p c j,i | k j q π k p i s q j p k j | c j, : i ´ q ş R Nq π c j,i p c j,i | k j q π k p i s q j p k j | c j, : i ´ q d k j in an update step to compute π k p i s q j p k j | c j, : i q .In other words, the idea is to start from knowledge (encoded in π k pp i ´ q s q j p k j | c j, : i ´ q ) at an observation time step t obsi ´ . Then, knowledge for a new observation time step is forecasted / predicted using only the evolution model eq. (8). This forecast is finally corrected using the information given by observation c j,i . The unknown reference trajectory is ultimately given as mean of the marginal PDF π k p i s q j p k j | c j, : N obs q . Volume x, Issue x, 2017October 23, 2018
P. Zaspel
The EnKF is a Monte-Carlo–type implementation of the above discussed iterative data assimilation task. Insteadof explicitly computing the posterior PDFs π k p i s q j p k j | c j, : i ´ q and π k p i s q j p k j | c j, : i q , the EnKF constructs an ensem-ble of realizations of random variables representing these PDFs in an empirical sense. In that context, forecast and analysis ensembles are distinguished. As we will see, the computation of the forecast ensemble corresponds to ap-proximating π k p i s q j p k j | c j, : i ´ ,j q , while the computation of the analysis ensemble corresponds to the approximationof π k p i s q j p k j | c j, : i,j q .Shall N e be the size of the ensembles. Then, the EnKF algorithm starts by drawing N e samples k p q , j , . . . , k p q ,N e j of the (initial) system state according to the PDF of k p q j . The algorithm consists of two main steps which are iterativelydone for i P t , . . . , N obs u . In the forecast step, the ensemble is propagated over s steps of the evolution model in eq. (8) to reach the nextobservation time step t obsi “ i s ∆ τ . To achieve this, realizations n p l q ,m P R N q for m P t , . . . , N e u are drawni.i.d. from n p n q in each of the s steps. Then the propagation equation reads for n “ , . . . , s as k p i p s ´ q` l q ,mj “ k p i p s ´ q`p l ´ qq ,mj ` ? ∆ τ n p l q ,m , m P t , . . . , N e u . The newly generated ensemble is the forecast ensemble ´ k f,mj ¯ N e m “ with k f,mj : “ k i s,m j . We further computethe empirical forecast mean k fj : “ N e N e ÿ m “ k f,mj P R N q (11)and the empirical forecast covariance (matrix) Σ f k j : “ N e ´ N e ÿ m “ ´ k f,mj ´ k fj ¯ ´ k f,mj ´ k fj ¯ J P R N q ˆ N q . (12) In the analysis step, the Kalman filter [14,26] is applied to the forecast ensemble to compute an analysis ensemble ´ k a,mj ¯ N e m “ representing the PDF π k p i s q j p k j | c j, : i q , which is conditioned to the new observation c j,i . As part of theKalman filter, the forward model eq. (10) with k refj,i s being replaced by k p i s q j is evaluated. Here, we use a linearforward map H j,i . Moreover all involved random variables are Gaussian. Therefore, it can be shown that the analysisensemble follows a Gaussian distribution, too and thus it can be fully characterized by the empirical analysis mean k aj and the empirical analysis covariance Σ a k j .Based on this observation, the core idea of the Kalman filter is to compute the empirical analysis mean asminimization problem k aj “ arg min k j P R Nq ˜››› k j ´ k fj ››› ´ Σ f k j ¯ ´ ` } H j,i k j ´ c j,i } σ ´ e ¸ . Given the linearity of H j,i , the minimum can be exactly computed as k aj “ k fj ´ U j,i p H j,i k j ´ c j,i q , International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference where U i,j is the Kalman (update) matrix U j,i “ Σ k fj H J j,i p H j,i Σ f k j H J j,i ` σ e q ´ . Instead of explicitly computing the empirical analysis mean and covariance (the latter by an analogous update idea),the analysis part of the Ensemble Kalman Filter (with perturbed observations) [14, Chapter 7] directly updates theforecast ensemble by k a,mj “ k f,mj ´ U j,i p H j,i k f,mj ` e j,i,m ´ c j,i q , m P t , . . . , N e u , where t e j,i,m u N e m “ are realizations of e j,i . If required, empirical versions of the analysis mean and analysis covariancecan be computed analogously to eq. (11) and eq. (12). Finally, the next forecast step is initialized with k p i s q ,mj “ k a,mj ,that is, the analysis ensemble replaces the system state for t obsi . For i “ N obs the algorithm terminates with an analysis ensemble, representing the posterior PDF π k p Nobs ¨ s q j p k j | c j, : N obs q . The reference trajectory is extracted as empirical mean k j : “ N e ř N e m “ k p N obs s q ,mj . The(mean) perfusion p j can be derived as p j “ ρ j k j | t “ . Empirical covariances are extracted as discussed before.Moreover, in case cumulative distribution functions or other probabilistic quantities shall be extracted, a kernel-density estimator (such as ksdensity in Matlab) is applied to the generated ensemble.
4. NUMERICAL RESULTS
In this section, we demonstrate the beforehand introduced numerical method for artificial test data. To this end, wefirst introduce the source of this test data, which is a
Digital Perfusion Phantom . Then, we study the numericalproperties of our method in terms of convergence, parameter dependence and input dependence in a single-voxelscenario. Finally we solve the perfusion inference problem for a slice of a full (artificial) DCE imaging brain data set.
Digital Perfusion Phantoms (DPP) [11–13] allow to artificially generate DCE image data for perfusion analysis.Thereby new algorithms can be tested on such data without the additional constraints of true patient data. PerfusionPhantoms basically solve the forward problem, which involves to transform perfusion information into contrast agentconcentrations. In our work, we use the
Digital Brain Perfusion Phantom package [13], which is a Matlab imple-mentation of the model introduced in [11]. The software provides a radiological image of a brain. A user interface,see Figure 1, allows to mark regions of reduced and strongly reduced perfusion. It is possible to control the obser-vation snapshot time step size (i.e. ∆ t obs ) of the artificial radiological imaging process. The measurement time is T “
49. The resolution of the artificially generated data is N D “ p , , q . The arterial input function isprovided as discrete evaluations c art p t inputi q with t inputi “ i . The Perfusion Phantom package uses a piecewisecubic spline interpolant through this data as exact c art , see Figure 2(a). During the artificial imaging process, eachsnapshot (i.e. c i ) is written in a separate file. A baseline for the radiological images is written, too. It contains themeasurement data without contrast agent concentrations. In our examples, we always subtract this baseline data fromthe artificial measurements to obtain just the necessary concentration information.We perform a major part of our numerical tests on a single reference voxel which has been chosen arbitrarily as p , , q . The observation data for that single voxel is stored with a time step size of ∆ t obs “ .
25. This datais interpolated by a piecewise cubic spline to obtain measurement data at arbitrary points in time for our initial tests,cf. Figure 2(b). Towards the end of this section, results for a full slice p¨ , ¨ , q of the full data set are discussed.We start by showing a series of numerical results obtained for given artificial input without noise. These results will give an insight into the choice of the different parameters of the method and into the convergence properties of the method. Noisy data is discussed afterwards. Volume x, Issue x, 2017October 23, 20180
P. Zaspel
Figure 1:
The source of our artificial measurements is the Digital Brain Perfusion Phantom package [13]. A Matlab implementationof this work is available. It allows to mark brain regions with reduced and severely reduced perfusion, here shown with the colorsyellow and red. Given this data, artificial DCE imaging data is are created.
Let us first have a look at the evolution of the analysis ensemble during the sequential data assimilation process. Wehave chosen an observation time step size of ∆ t obs “ .
25, a quadrature step size of ∆ τ “ . s “
4) andan ensemble size of N e “ k j has a magnitude of about 10 ´ .Therefore, the scaling α of the covariance matrix Σ n should be relative to a standard deviation of 10 ´ . With thisin mind, we set α “ p ´ q . relative variance of 0 . k j and the scaling α . This is considered future work.The correlation length is set to (cid:96) “
2. For the covariance of the initial state k p q j , we impose an additional scaling of σ “ . σ e “ . k aj , i.e. the prediction for the unknown kernel function,for different observation times during the operation of the EnKF. Note that a scaled evolution of k aj at t “ p . Therefore, discussing numerical results for k aj is equivalent to discussingresults for p . The major information gain for the predicted result is in time interval r , s . This is the time intervalin which the concentration at the arterial inlet grows. Afterwards, the data assimilation process only gains very littlemore information and converges towards the final result. Next, we discuss the convergence of the empirical mean of k j with respect to the ensemble size N e . In the following,we will always concentrate on the last analysis ensemble obtained after assimilating the observation for t obs “
49. Toshorten notation, we skip additional indices, indicating this and call the empirical mean of this analysis ensemble k j .Our convergence study with respect to the ensemble size uses the same parameters as in the previous paragraph.However, this time, we change the size of the ensemble. In Figure 4(a), we show the empirical mean k j for ensemble sizes N e P t , , , , u . The convergence in the error of the empirical mean is shown in Figure 4(b). Here, we define the solution for N e “ (cid:96) error International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference t c a r t (a) arterial input function t c t i ss (b) concentration measurement in single voxel Figure 2:
We use idealized concentration functions for one tissue voxel in order to test the implemented numerical method. ¨ ´ time t m ea n e s ti m a t e k a j t obs “ . t obs “ . t obs “ . t obs “ . t obs “ . t obs “ . t obs “ . t obs “ . Figure 3:
During the sequential data assimilation process, the analysis ensemble and thereby the empirical mean of the kernelfunction gets continuously updated, here shown for different update time steps. } k j ´ k overkillj } (cid:96) } k overkillj } (cid:96) towards this solution. The results indicate a convergence order of approximately . This is the ex-pected order of convergence, since we use a Monte Carlo-type estimator. Note that an ensemble size of about 20,which is often used for Ensemble Kalman Filters, seems not to be enough in this application. In that case, we observea highly oscillatory result with a strong overshooting for the initial peak of the mean estimate (which will be theperfusion estimate). ∆ τ In the following, we have a look at convergence with respect to the quadrature and evolution model step size ∆ τ .Here, we do not use an overkill solution. To achieve this, we (discretely) fold the empirical mean k j against the(discretized) arterial input function c art , i.e. we transfer the prediction for k j into observation space. In observationspace, we compare against the analytically given artificial measurement result c . Our numerical study uses a variationof the sub-step number s , i.e. we change ∆ τ while keeping all other parameters as in Section 4.2.In Figure 5(a), we visually compare the results obtained for an increasing number of sub-steps s (i.e. decreasing ∆ τ ). The convergence plot in Figure 5(b) further shows the error reduction in the relative (cid:96) norm for decreasing ∆ τ if we compare the convolved mean estimate k j with the real observation data. The results indicate approximately first Volume x, Issue x, 2017October 23, 20182
P. Zaspel ´ . . ¨ ´ time t m ea n e s ti m a t e k j N e “ N e “ N e “ N e “ N e “ (a) ensemble estimates for k j ´ ensemble size M ‘ e r o r o f k j measured errorconvergence order (b) convergence wrt. overkill solution Figure 4:
With growing ensemble size, the empirical estimate for the mean of the response / kernel function gets more accurate(left) and converges with roughly order (right). order convergence. In fact, parameter ∆ τ influences the Euler-Maruyama approximation of the continuous stochas-tic differential equation eq. (9) and the quadrature of the convolution integral. While the Euler-Maruyama method isknown to have halve order convergence, the rectangular rule is convergent of second order for sufficiently smooth inte-grands. The observed convergence behavior strongly depends on the dominance of one of the errors (time-integration,quadrature). The observed first order seems to indicate that the quadrature error for the convolution integral is dom-inant. Nevertheless, full second order convergence is not achieved. This observation is clearly a pre-asymptotic andstrongly problem-dependent result. Our next study shall give an insight into the influence of the system state model, more specifically the influence of thecorrelation length (cid:96) of the random variable n on the inferred solution. To study the influence of the correlation length,we keep the parameters as in Section 4.2 and apply different correlation lengths (cid:96) P t . , . , u . The results of thisnumerical study are given in Figure 6. Here, the inferred kernel function k j is shown for different correlation lengths.For growing correlation length the result gets less noisy. Hence, a larger correlation length has a regularizing effecton the solution. Since, in general, we seek for smooth solutions, we always choose (cid:96) “ Our final test with noise-free model data on a single voxel highlights the influence of a change of the observation timestep size ∆ t obs , i.e. a change in the number of observations that are made during the imaging process. To test this, wetake the same parameters as in Section 4.2, but change the observation time step size as ∆ t obs P t . , . , . , . u while keeping ∆ τ constant. The quantity that we study is the computed probability density function for k | t “ , hencea scaled version of p . We use the kernel density estimator ksdensity in Matlab to reconstruct a continuous PDFfor the ensemble data.The results of this study can be seen in Figure 7. Here, we make two observations. First, the mean of the PDFstill changes for growing number of measurements, converging towards a true solution. Second, and more important,we observe a variance reduction if we increase the number of measurements. This type of information would not be available in classical inverse approaches for compute perfusion estimation. That is, we can now obtain confidence information for our solution. International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference t c a r t ˚ k j s “ s “ s “ s “ s “ c tiss (a) k j in observation space ´ . ´ ´ . ´ ´ time step size ∆ τ e rr o r i n c a r t ˚ k j w r t . c measured errorconvergence order 1convergence order (b) convergence Figure 5:
A smaller time step size for the quadrature / system state model leads to convergence of k j in observation space towardsthe measurement concentration c . ´ . . . . . ¨ ´ time t m ea n e s ti m a t e k j ‘ “ . ‘ “ . ‘ “ Figure 6:
Longer correlation lengths (cid:96) impose a higher smoothness on the ensemble estimate.
Until now, we considered noise-free input data. Instead, we now discuss the same one-voxel input as before, but addartificial noise as c noisyj,i “ c j,i ` w j,i , i P t , . . . , N obs u , where the w j,i are realizations of i.i.d. random variables w j,i : Ω Ñ R , w j,i „ N p , σ w q with σ w P R the varianceof the noise.We use a series of test cases with σ w “ α rel and α rel P t ´ , ´ , ´ , ´ , ´ u . Hence, α rel correspondsto the relative variance with respect to the magnitude of the measurements. We keep a major part of the parametersfrom Section 4.2. However, we change the fixed observation error variance σ e , to a problem-adapted one, namely,i.e. σ e “ σ w . Note that in practice, one would empirically estimate the noise in the measurement data and wouldset σ e accordingly. Another change concerns the number of samples N e in the EnKF. As our experiments showed,the size of the ensemble has to be increased for higher variances. This is well covered by classical Monte-Carlo theory. Therefore, we set N e “ α rel P t ´ , ´ , ´ u while we use N e “ N e “ α rel “ ´ and α rel “ ´ , respectively. Volume x, Issue x, 2017October 23, 20184
P. Zaspel . . . . . ¨ ´ k | t “ e s ti m a t e d P D F f o r k | t “ ∆ t obs “ . ∆ t obs “ . ∆ t obs “ . ∆ t obs “ . Figure 7:
The more observation samples are taken, the more reliable the estimate of the solution. Hence, the estimated PDF for k | t “ shows a smaller variance for smaller observation time steps ∆ t obs . In Figure 8, we give two examples of noisy inputs for α rel “ . α rel “ . c art ˚ k j arealso given in Figure 8. In fact, the reconstructed solution is almost not influenced for α rel “ . k | t “ for growing noise in the data.Since we appropriately account for the noise in the input, the mean is almost identical up to α rel “ . e j,i , which acts here as a regularization for the noisy input. Nonetheless,as long as the observation error variance is set in the range of the input noise variance, the variance in the solutioncorrectly represents the variance coming from the noise in the input. We finally apply the beforehand studied method to a full application problem given by the Digital Brain PerfusionPhantom introduced in Section 4.1. We use the slice p¨ , ¨ , q with the choice of regions with reduced and severelyreduced perfusion as in Figure 1. We discuss a result for a large observation time-step size ∆ t obs “ .
0, a highlyresolved quadrature with ∆ τ “ . σ w “ . σ e “ σ w . All other parameters are kept as in Section 4.2. Storing and computing the ensembles for the discussed test cases is a rather challenging task. Just considering theanalysis ensemble for a single slice, we need to store for each of the 256 ˆ
256 voxels 5000 realizations of discretekernel functions k j given via N q “
785 double precision values leading to a total storage requirement of256 ˆ ˆ ˆ ˆ « . All our calculations are done in Matlab. We always compute 8 rows of the final 256 ˆ
256 slice at the same timeand reuse the random input for each voxel in order to reduce the runtime. Note that especially sampling from n p n q is very computationally demanding. In order to do the calculations, we need constant access to way more than 64GBytes of RAM. Due to storage und memory requirements, we use nodes of the cluster Rhea at Oak Ridge National
Lab to compute the full problem. Each node has 128 GBytes of RAM and a dual Intel® Xeon® E5-2650 CPU with
16 cores. To compute 8 lines, i.e. results for 8 ˆ “ International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference t c on t r a s t a g e n t c on ce n t r a ti on c noisyj for α rel « . c art ˚ k j (a) t c on t r a s t a g e n t c on ce n t r a ti on c noisyj for α rel “ . c art ˚ k j (b) Figure 8:
Even for stronger noise on the input data the inference of k j is acceptable, as long as the measurement variance σ e is chosen appropriately. Here, we compare the noisy input c noisy j k j in observation space for α rel “ . α rel “ . that Matlab uses approximately 14 cores of the full machine. The total computing time (with respect to one node ofRhea) is thereby roughly 104 hours or about 4 . In our application examples, we consider the approximation of probabilistic quantities of interest in connection withthe perfusion p j : “ p p k j q “ ρ j k p q . Besides of the mean p j we are especially interested in probabilities for thecorresponding random variable p j to be in a given range. To be more specific, we compute the probabilities P p p j ă q , P p ď p j ă q , P p p j ě q , noting that the underlying perfusion p j lies in the interval r , s in the case of the Digital Brain Perfusion Phantomdata that we consider. These quantities give probabilities for low, medium and high perfusion in some region of thebrain. Given the final analysis ensemble for k j , it is easy to compute the above quantities by using the kernel densityestimator ksdensity . The latter one can compute a cumulative distribution function (CDF) for each voxel, whichis finally evaluated appropriately. An important advantage of the use of a Digital Perfusion Phantom is the existence of a reference solution to com-pare with. The DPP software that we use stores the reference solution together with the other generated data. InFigure 10(b), we show the reference solution for our full application test case. The approximated result of our appli-cation example study, i.e. p , is shown in Figure 10(a). As expected from our single-voxel study, the inferred perfusionmatches the exact perfusion result well. Note that this is the case even though we add a considerable amount of noiseon the measurements. As discussed before, we can use the ensemble-based estimate of the posterior probability density function to extract a wide range of probabilistic information on the inferred solution. This is the main result of this work. To
Volume x, Issue x, 2017October 23, 20186
P. Zaspel . . ¨ ´ ¯ k j e s ti m a t e d P D F f o r k | t “ α rel « . α rel « . α rel « . α rel “ . α rel “ . Figure 9:
The proposed method is pretty robust with respect to noise. This can be seen, if we study the estimated probability densityfunctions for k | t “ . With growing noise variance, the empirical PDF estimate still recovers the mean appropriately. Extreme noisevariances degenerate the result, as expected.
50 100 150 200 25050100150200250 x y (a) Approximated perfusion p
50 100 150 200 25050100150200250 x y (b) Noise-free reference perfusion given by the DPP Figure 10:
Our approximation method recovers the reference perfusion result (right) as mean of the ensemble in the EnsembleKalman Filter. Both results match well, even though we introduced a considerable amount of artificial noise. exemplify this, we compute space-dependent probabilities for low (Figure 11), medium (Figure 12(a)) and high(Figure 12(b)) perfusion ranges, cf. Section 4.8.2. In case of Figure 11, we e.g. can now easily identify ranges of lowperfusion and even can give a probability for this result.In general, we claim that this probability information or derived probabilistic quantities (variance, percentiles,etc.) can give domain-experts in radiology a much clearer information on the reliability of the inferred estimates.
5. SUMMARY
In this work, we have discussed the use of Ensemble Kalman Filters for sequential data assimilation in order to inferprobabilistic information on (blood) perfusion in tissue for given measurements from dynamic contrast–enhancedimaging. The deterministic inference of perfusion is well-known in the field of radiological imaging. However, to thebest of the author’s knowledge, the new contribution is the approximation of PDFs for the perfusion given (noisy) mea-surements. EnKF are well-known in inference for dynamical systems and partial differential equations with stochasticcoefficients. Hence, modeling the dynamic contrast–enhanced imaging process as sequential data assimilation in a
Bayesian context was the main contribution of the work. Given the ensemble-based approximation of the PDF, we could compute probabilistic quantities such as probabilities for perfusion parameter ranges.
International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference
50 100 150 200 25050100150200250 x y . . . . Figure 11:
The advantage of the proposed method is that we are now also able to compute probabilistic information for thesolution, here shown by plotting the probability P p p j ă q . Hence, the depicted results give the space-dependent probability forlow ( ă
10) perfusion.
50 100 150 200 25050100150200250 x y . . . . (a) P p ď p j ă q
50 100 150 20050100150200250 x y . . . . (b) P p p j ě q Figure 12:
Based on the results of the EnKF, it is easily possible to identify regions of high probability to have medium (left) andhigh (right) perfusion.
The new approach was first investigated for a single-voxel example with respect to convergence and parameterinfluence. Afterwards, it was applied to artificial application data generated by a Digital Perfusion Phantom, i.e. amodel for deriving DCE image data for given perfusion data. Overall, the effectiveness of the method could bedemonstrated, showing empirical convergence results and appropriate approximations of probabilistic information.The use of realistic patient data, refined problem-adapted covariance kernels, advanced filtering techniques and anefficient parallel implementation are future work.
ACKNOWLEDGEMENTS
This work is funded by the Swiss National Science Foundation (SNF) under project number 407540 167186. Fur-thermore, this research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National
Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-
AC05-00OR22725.
Volume x, Issue x, 2017October 23, 20188
P. Zaspel
The author also likes to thank Wolfram Stiller and Christian Weis of the department of
Diagnostic and Inter-ventional Radiology of the
University Medical Center Heidelberg and Holger Fr¨oning of the
Institute of ComputerEngineering at University of Heidelberg for fruitful initial discussions on the application background.
References
1. Brix, G., Griebel, J., Kiessling, F., and Wenz, F., Tracer kinetic modelling of tumour angiogenesis based on dynamic contrast-enhanced CT and MRI measurements,
European journal of nuclear medicine and molecular imaging , 37(1):30–51, 2010.2. Fieselmann, A., Kowarschik, M., Ganguly, A., Hornegger, J., and Fahrig, R., Deconvolution-based CT and MR brain perfusionmeasurement: theoretical model revisited and practical implementation details,
Journal of Biomedical Imaging , 2011:14,2011.3. Sourbron, S., A tracer-kinetic field theory for medical imaging,
IEEE transactions on medical imaging , 33(4):935–946, 2014.4. Tofts, P.S., Modeling tracer kinetics in dynamic Gd-DTPA MR imaging,
Journal of Magnetic Resonance Imaging , 7(1):91–101, 1997.5. Parker, G.J., Roberts, C., Macdonald, A., Buonaccorsi, G.A., Cheung, S., Buckley, D.L., Jackson, A., Watson, Y., Davies, K.,and Jayson, G.C., Experimentally-derived functional form for a population-averaged high-temporal-resolution arterial inputfunction for dynamic contrast-enhanced mri,
Magnetic resonance in medicine , 56(5):993–1000, 2006.6. Østergaard, L., Weisskoff, R.M., Chesler, D.A., Gyldensted, C., and Rosen, B.R., High resolution measurement of cerebralblood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis,
Magnetic reso-nance in medicine , 36(5):715–725, 1996.7. Evensen, G.,
Data assimilation: the ensemble Kalman filter , Springer Science & Business Media, 2009.8. Stuart, A.M., Inverse problems: a Bayesian perspective,
Acta Numerica , 19:451–559, 2010.9. Iglesias, M.A., Law, K.J., and Stuart, A.M., Ensemble Kalman methods for inverse problems,
Inverse Problems , 29(4):045001,2013.10. Ernst, O.G., Sprungk, B., and Starkloff, H.J., Analysis of the ensemble and polynomial chaos Kalman filters in Bayesianinverse problems,
SIAM/ASA Journal on Uncertainty Quantification , 3(1):823–851, 2015.11. Riordan, A.J., Prokop, M., Viergever, M.A., Dankbaar, J.W., Smit, E.J., and de Jong, H.W., Validation of CT brain perfusionmethods using a realistic dynamic head phantom,
Medical physics , 38(6):3212–3221, 2011.12. Pianykh, O.S., Digital perfusion phantoms for visual perfusion validation,
American Journal of Roentgenology , 199(3):627–634, 2012.13. Manhart, M.
Digital Brain Perfusion Phantom Documentation . provided on the data web page of the pattern recognition labat FAU Erlangen-N¨urnberg, Germany (last check: Feb. 23, 2017).14. Reich, S. and Cotter, C.,
Probabilistic forecasting and Bayesian data assimilation , Cambridge University Press, 2015.15. Majda, A.J. and Tong, X.T., Performance of ensemble Kalman filters in large dimensions,
Communications on Pure andApplied Mathematics , 71(5):892–937, 2018.16. Tong, X.T., Performance Analysis of Local Ensemble Kalman Filter,
Journal of Nonlinear Science , 28(4):1397–1442, 2018.17. Schillings, C. and Stuart, A.M., Analysis of the ensemble Kalman filter for inverse problems,
SIAM Journal on NumericalAnalysis , 55(3):1264–1290, 2017.18. Schillings, C. and Stuart, A.M., Convergence analysis of ensemble Kalman inversion: the linear, noisy case,
Applicable Anal-ysis. An International Journal , 97(1):107–123, 2018.19. Anderson, J.L., An ensemble adjustment Kalman filter for data assimilation,
Monthly weather review , 129(12):2884–2903,2001.20. Anderson, J.L., A local least squares framework for ensemble filtering,
Monthly Weather Review , 131(4):634–642, 2003.21. Tippett, M.K., Anderson, J.L., Bishop, C.H., Hamill, T.M., and Whitaker, J.S., Ensemble square root filters,
Monthly WeatherReview , 131(7):1485–1490, 2003.22. Zhang, F., Zhang, M., and Hansen, J.A., Coupling ensemble Kalman filter with four-dimensional variational data assimilation,
Advances in Atmospheric Sciences , 26(1):1–8, 2009.
International Journal for Uncertainty Quantification nKF for reliability estimation in perfusion inference
23. Reich, S., A nonparametric ensemble transform method for Bayesian inference,
SIAM Journal on Scientific Computing ,35(4):A2013–A2024, 2013.24. Reich, S. and Cotter, C.J., Ensemble filter techniques for intermittent data assimilation,
Large Scale Inverse Problems. Com-putational Methods and Applications in the Earth Sciences , 13:91–134, 2013.25. Nævdal, G., Sævareid, O., and Lorentzen, R.J., Data assimilation using MRI data, In
Proceedings of the ECCOMAS Congress2016 , 2016.26. Kalman, R.E. , A new approach to linear filtering and prediction problems,
Journal of basic Engineering , 82(1):35–45, 1960., 82(1):35–45, 1960.