·Home ·Table of Contents ·Computer Processing and Simulation

Dynamic Image Reconstruction: General estimation Principles for Dynamic tomography V. M. Artemiev, A. O. Naumov Institute of Applied Physics, Akademicheskaya str.16, 220072, Minsk, Belarus. G.-R.Tillack Federal Institute for Materials Research and Testing (BAM), Unter den Eichen 87,12205 Berlin, Germany. Contact

Key words: dynamic image, stochastic estimation, dynamic reconstruction, Kalman filter

1. Introduction

In recent years the computerized tomography (CT), e.g. the X-ray CT, is rapidly involved into the problem of process tomography [1]. Here the problem is faced to provide information about the internals of a process, in particular its dynamical behavior. These properties can be interpreted as dynamic field or image. However the state of the art methods applied for process tomography (compare [1]) are based on the classical reconstruction principles developed for standard CT applications [2]. These approaches assume that the image properties are constant during the data acquisition. This requirement restricts the applicability of these approaches to slow varying processes or high speed acquisition systems. The intend of the paper is to discuss a new approach to real-time optimal linear reconstruction of dynamic images as a new framework for process tomography. This technique is called dynamic image reconstruction.
The following assumptions are made: The image is supposed to be - in a statistical sense - a stationary Gaussian random field discrete in space and time. Additionally the acquisition task is assumed to be linear and includes additive noise and blurring. Accordingly the process image and the acquisition task can be represented as a Markovian random sequence. The general principle for the reconstruction of the dynamic image properties is given by Bayesian approach to the Markovian stochastic estimation theory[3,4]. This approach leads to Kalman filter type reconstruction algorithm in the time domain.

2. Gaussian dynamic image model

It is supposed that the dynamic image can be represented as a Gaussian random field discrete in space and in time. The domain of the field is assumed to be a hyper parallelepiped described by an equidistant grid of image elements. The total number of image elements is given by S defining the image dimension. The discrete time coordinate is noted by k = 0,1,2,... and the image value for the element i=1,2,...,S is noted as x_i(k) . The assemblage of these values can be represented as S -dimensional image vector x(k) = [x₁(k),...,x_i(k),...x_s(k)]^T

A dynamic Gaussian random field with the above properties can be represented by the following S-dimensional linear stochastic difference equation (SDE)

x(k+1)=Ax(k)+u+w(k)

(2.1)

with the S ´ S dimensional constant matrix A and the S dimensional constant vector u . The S dimensional random vector w(k) represents Gaussian distributed white noise with zero mean and constant covariance matrix

given by the probability density distribution (PDD) function

(2.2)

The SDE image model (2.1) define the image field as a multidimensional Markovian random sequence. Assuming Gaussian initial conditions the process x(k) remains Gaussian for all time. Additionally it is supposed that the process
x(k) tends to a stationary regime, i.e. the field PDD function tends to a stationary Gaussian distribution for k(r)¥

The Markovian random sequence (2.1) is a priori completely defined by two functions [3,4]:

•  the image prior PDD function p_pr(x,k) and
• the transition probability (TP) p_pr(x, k+1 | x',k). Here the prior PDD function is chosen in the stationary regime

(2.3)

The S dimensional vector m_o gives the image field prior mean value and the S ´ S matrix p_o describes the prior covariance matrix in the space domain. The corresponding TP function is given by

(2.4)

In the next step the matrices A,Q as well as the vector u of the SDE image model (2.1) have to be determined for the given prior m_o and P_o.
From eq. (2.1) the following expression for the image mean value m_o(k)=áx(k)ñ is found

(2.5)

For the stationary regime k(r)µ follows

(2.6)

Hence the stationary field mean value m_o is determined by the choice of the constant vector u if the matrix A is known.
The SDE for the centralized image vector follows from eq.(2.1) and eq.(2.5)

(2.7)

The space-time S´S covariance matrix P₀(k+1)|k)can be introduced by the expression

(2.8)

with the time delay parameter l =0,1,2,.... For l=0 the covariance matrix (2.8) transfers to the space covariance matrix P_o(k)=P_o( k | k ).
As for eq. (2.7) the equation for the image vector are given by

(2.9)

By substituting the equations for the image vectors into eq. (2.9) the following recursive equation is obtained

(2.10)

with A¹ denoting the l-th power of the matrix A. By substituting Eqs. (2.10) and (2.7) into expression (2.8) the equation for the space-time covariance matrix is obtained.

In the steady state and .The S´S matrix F_o(l) is the stationary space-time covariance matrix that satisfies

(2.11)

For l=0 the matrix F_o(0) transforms to P_o with

(2.12)

After multiplying matrix A¹ by eq.(2.12) and substituting the result into eq. (2.11) the final equation for the stationary space-time covariance matrix is obtained

(2.13)

Eqs.(2.6), (2.12) and (2.13) define the mean value m_o, space covariance matrix P_o and stationary space-time covariance matrix F_o(l) for given A,Q and u .
The following inverse equations are valid to obtain the model parameters A , Q and u for given P_o,F_o(l) and m_o :

(2.14)

It is important to note that the solution of eq. (2.14) is unique for arbitrary mean image vector m_o and space covariance matrix P_o but only for one chosen time delay parameter l .

The stochastic boundary conditions for the field model (2.1) can be specified by choosing the mean image vector m_oand the space covariance matrix P_o on the boundary of the hyper parallelepiped under consideration.

3. The observability of the reconstruction as an estimation problem

Observations in terms of projection data are used for the image reconstruction. The observation z_i(k) is given by the sum of S_k < S components of the vector x(k) and is commonly called ray sum[2]

(3.1)

with the 1´S row matrix H₁(k) having S-S_k elements equal to zero. For the following time k + 1 the observation is given by another set of S_k+1 < S components of the vector x(k+1) with z₁(k +1)=H₁(k +1)x(k + 1). The observations z₁(k) for k=1,2,... are scalars and the corresponding observation matrices have the properties

(3.2)

Here the coefficient c_k give the square sum of the matrix elements of matrix H₁(k). If these elements hold only the values zero or one c_k gives the number of non-vanishing elements of H₁(k). The coefficients c_kl describe the degree of the matrices' intersection. If the coefficient c_kl is equal to zero the intersection of the matrices is the empty set.
The m-dimensional vector z_m(k) holds the simultaneous observation of m ray sums(3.1) defined by

zm(k)=Hm(k)x(k)

(3.3)

with the m ´ S observation matrix H_m(k)which is built from the row matrices H_li(k) ,i=1,2,...,m , i.e.

Hm(k)=[H11(k),...,H1i(k),...,H1m(k)]T

(3.4)

If for each k number of row matrices m=M and if

(3.5)

with the unit vector 1, then the observation z_m(k) is called projection.
The first step of the reconstruction like an estimation procedure is given by the evaluation of the observability properties. The reconstruction process is fully observable at the time interval [0,k] if the image x(k) can be exactly estimated with observations z₁^k =[z^_m(1),z^_m(2),...,z^_m(k)]^T . The observability criteria consists in calculating the rank of the
(k ´m) ´ S observation matrix [3,4]

(3.6)

The dynamic image reconstruction process is fully observable if Rank(H₁^k)=S.

It is supposed that the number of projections is equal toN and that each projection consists of M ray sums. Then the rank of the observation matrix (3.6) can not be greater than NM. In practical applications it usually turns out that Rank(H₁^N)<S and the dynamic image reconstruction is an ill-posed problem [3,4]. This is the main difficulty one has to face compared to standard estimation problems.
Usually the observation is corrupted by additive noise v(k)and linear image blurring given by y(k)=D(k)x(k) with the blurring matrix D(k). As a result the observation takes the form

(3.7)

The noise v(k) is supposed to be m-dimensional Gaussian distributed random vector with zero-mean and the m ´ m-dimensional covariance matrix

(3.8)

represented by the PDD function

(3.9)

For equation (3.7) and noise PDD function (3.9) the following likelihood function is valid

(3.10)

with the normalization coefficient q₁(k).
For the dynamic image model (2.1) and the observation (3.7) the image reconstruction is reduced to the well known optimal stochastic estimation problem as formulated in [3,4]. This estimation problem is ill-posed and a standard approach to solve it is the Bayesian approach here applied to dynamic image reconstruction.

4. The Bayesian approach for the dynamic image reconstruction.

The dynamic image vector x (k) defined by eq.(2.1) together with the observation z_m(k) given by eq. (3.7) describe a joint Markovian random sequence. The posterior PDD function p _ps(x, k | z ₁ ^k) for the current image vector satisfies the following recursive equation [3,4]

(4.1)

with the normalization coefficient q_ps(k) . The likelihood function p₁(z , k | x , k) is defined by eq. (3.10) and the transition probability p_pr(x,k|x',k-1) follows from eq. (2.4). The initial condition for eq. (4.1) consists in the image prior PDD function given by eq. (2.3). Assuming a stationary Gaussian image model and a linear observation the posterior PDD function (4.1) remains Gaussian for any k and has the form:

(4.2)

with S-dimensional posterior mean value vector

(4.3)

and S ´S posterior covariance matrix

(4.4)

For the Gaussian assumption the optimal image vector estimate is given by the posterior mean value vector (4.3) together with the posterior error covariance matrix (4.4). The optimal estimation satisfies the recursive Kalman filter (KF) algorithm which consists of the following three equations:
• the equation for the optimal estimation

  (4.5)

• the Riccaty equation for the error covariance matrix

  (4.6)

• and the equation for the optimal gain matrix

  (4.7)

The notation defines a one-step prediction of the estimation

(4.8)

and the P(k + 1 | k) is an one-step prediction for the error covariance matrix

P(k+1|k)=AP(k)AT+Q

(4.9)

For Markovian models one-step predictors (4.8) and (4.9) are optimal. For other models higher order predictors may become optimal.
The initial condition for eq.(4.5) is the prior image mean value m_o given by eq.(2.6) while for the Riccaty equation (4.6) it is given by the prior covariance matrix P_o from eq.(2.12). The image model parameters A,Q and u are introduced into KF algorithm by the expressions(4.8) and(4.9).
The equations for the error covariance matrix(4.6) and the optimal gain matrix (4.7) are independent on the z_m(k) observations . Therefore both matrices can be calculated in advance for later use. Accordingly it is possible to investigate the reconstruction performance

(4.10)

in advance and evaluate the influence of the observation matrix H_m(k).
There are basically three different cases to be discussed: The first one is the scalar case (m = 1 ) for sequent ray-sum observations z₁(k)with statistical independent noise components v(k) having variances r(k) . As a result in eqs.(4.6) and (4.7) becomes a scalar with trivial inverse operation. For this case the calculation expenses are sufficiently decreased such that it is sometimes efficient to represent the m-dimensional observation vector as an artificial sequence of scalar observations.
In the second case the observation vector z_m(k) is equal to the total projection with m =M like a number of parallel working detectors. Here the inversion of the high-dimensional matrix in eqs. (4.6) and (4.7) has to be performed leading to large numerical expenses.
The third case supposes that at the current time k all N projections with dimension M are observed simultaneously giving the observation vector z_MN(k) . The image reconstruction becomes a one-step procedure with the following quasi-static KF algorithm

(4.11)

The static random image is a partial case of the above dynamic image(2.1) with A =I,Q = 0 ,u=0 given by the SDE
x(k +1)=x(k) The Bayesian approach for the reconstruction of random static images was studied for a long period of time after the publication[2] where the one-step optimal linear reconstruction algorithm was found which is similar to proposed KF algorithm(4.11).

5. The dynamic image reconstruction for simulated data.

The performance of the derived dynamic image reconstruction procedure is investigated for a non-homogeneous random solid object moving with constant velocity. Fig.1 shows the geometry of the object together with the X-ray data acquisition system. The object has a square cross section in the(r₁,r₂) -plane with random distribution of material properties. The object moves with given velocity along the axis. As a result the image values in the (r₁,r₂) -plane are given by a random process assumed to be Gaussian and stationary. The radiation source takes N_o discrete positions arranged on a circle with its center on the r₃axis. The image projections are measured using detector elements placed on the circle with its center at the source position. The projection data acquisition procedure is implemented at discrete time k=1,2,... within the acquisition interval Dt. It is supposed that during this interval the image values are constant, i.e. the dynamic properties of the random object are within this period of time negligible, and the n₀£N₀ projections are valid. The remaining time at the interval between k and k+1 is applied for the image reconstruction.

Fig 1: Geometrical setup.

The image vector is supposed to haveS=250 x 250 image elements, the number of detector elements is chosen to M = 400, and the total number of source positions is N_o=50. The image is assumed to be a Markovian random sequence with zero-mean value, axial symmetric covariance matrix with covariance radius of 30 image elements. The relaxation time depends on the motion velocity and is chosen to 200 time intervals.
The image reconstruction is carried out sequentially applying the KF algorithm (4.5)-(4.9) for the case m = M. Two reconstruction variants are investigated. First it is supposed that the image reconstruction procedure is applied to the quasi-static regime using 50 projections acquired during the time interval between the reconstruction steps kand k+1. The reconstruction performance, given by the criteria(4.10), is shown in Fig.2 with representing the number of projections actually used for the reconstruction. Each reconstruction step is performed without taking into account any dynamic properties of the image, i.e. the time covariance matrix (2.11) is neglected. This procedure yields the reconstruction error of at the end of every reconstruction step.

Fig 2: Reconstruction error for quasi-static regime

Fig 3: Reconstruction error for dynamic regime

In the second variant the dynamic reconstruction technology is employed. Here only 10 projections are included in each reconstruction step within the same time interval as used in the quasi-static regime. To account for the dynamic properties of the image in this case the time covariance matrix (2.11) is employed. Fig.3 shows the performance of the dynamic reconstruction. At the end of intervals the same accuracy is obtained as in the previous case. Thus the dynamic reconstruction requires a smaller number of projections per reconstruction step to reach the same accuracy because of considering the time covariance matrix as a kind of prior knowledge about the dynamic properties of the image to complete the previously introduced spatial prior.
The reconstruction results are shown in fig.4. Here fig. 4a gives a part of the original image represented over a total time of 300 intervals:k=1,...,300. . Fig4b represents the result of the dynamic reconstruction for the same period of time considering a total number of 3000 projections.

Fig 4: Reconstruction result for dynamic regime (b) in comparison to original image (a)

6. Conclusions

The proposed dynamic reconstruction for process tomography solves the task of reconstruction of time dependent image values from projection data. The dynamic image model is supposed to be represented as a Gaussian random field discrete in space and time. The projection model is assumed to be linear taking into account additive noise and image blurring. Therefore the resulting reconstruction problem is ill-posed. Accordingly the stochastic estimation theory employing the Bayesian approach is apply to regularize the above problem. As a result a discrete-time Kalman filter algorithm is found for the optimal linear dynamic reconstruction. The presented examples show the principle capability of the derived algorithm to reconstruct the dynamic properties of the image under investigation from restricted number of projections. Additional investigations have to be carried out to evaluate the capability of the dynamic reconstruction applied to industrial problems in the sense of process tomography. Furthermore the generalization of the proposed approach to nonlinear algorithms is under research.

7. References

1. Williams R A and Beck M S (Eds.) 1995 Process Tomography - Principles, Techniques and Applications (Oxford: Butterworth-Heinemann)
2. Herman G T 1980 Image Reconstruction from Projections. The Fundamentals of Computerized Tomography (New York: Academic Press)
3. Sage A P and Melse I L 1971 Estimation Theory with Applications to Communication and Control (New York: McGraw-Hill)
4. Jarlikov M S and Mironov M A 1993 Markov Estimation Theory for Random Processes (Moscow: Radio i Sviaz, in Russian)

© AIPnD , created by NDT.net

|Home|    |Top|