Functional imaging based on tomographic X-ray imaging relies on the reconstruction of a temporal sequence of images which accurately reproduces the time attenuation curves of the tissue. The main constraints of these techniques are temporal resolution and dose. Using current techniques the data acquisition has to be performed fast so that the dynamic attenuation values can be regarded as static during the scan. Due to the relatively high number of repeated scans the dose per single scan has to be low yielding a poor signal-to-noise ratio (SNR) in the reconstructed images. In a previous publication a temporal interpolation scheme in the projection data space was relaxing the temporal resolution constraint. The aim of this contribution is the improvement of the SNR. A temporal smoothing term is introduced in the temporal interpolation scheme such that only the physiologic relevant bandwidth is considered. A significant increase of the SNR is achieved. The obtained level of noise only depends on the total dose applied and is independent of the number of scans and the SNR of a single reconstructed image. The approach might be the first step towards using slowly rotating CT systems for perfusion imaging like C-arm or small animal CT scanners.

1. Introduction

The application of tomographic X-ray techniques to carry out perfusion studies relies on the reconstruction of a sequence of images which accurately reproduces the dynamics of the flow of an injected contrast agent. The temporal evolution of the attenuation coefficient corresponding to a voxel is denoted as time attenuation curve (TAC) or perfusion signal and is proportional to the concentration of contrast agent in the region occupied by said voxel. Hence, the temporal evolution of the gray value of a voxel is proportional to the temporal evolution of the average concentration of contrast agent within the voxel. TACs are used as input for algorithms based on a tracer kinetic model which compute functional parameter maps. Several approaches have been presented in the literature for this purpose. We refer the reader to [1] for an overview of these methods. The obtained maps provide the physician with physiological information related to blood supply to the tissues in the ROI which is crucial for the diagnosis.

One of the main constraints posed by these techniques is the temporal resolution; that is, the data acquisition-reconstruction process has to be capable of reproducing the perfusion signals. For this reason, up to now, data acquisition is performed fast enough to assume that the attenuation values are constant during data acquisition. Another important issue are dose considerations. Due to the relatively high number of repeated scans in a perfusion study the dose per single scan has to be low yielding a poor signal to noise ratio in the reconstructed images. The noise level can even be higher than the level of contrast enhancement caused by contrast agent flow. Under such conditions, the detection of perfusion signals becomes a challenging issue. This aspect will become even more critical with the introduction of large area detectors where a larger area of the patient is exposed.

Currently, only CT scanners allow a sampling rate high enough to assume constant attenuation during acquisition and therefore the only clinical application of tomographic X-ray techniques to perform perfusion studies is perfusion Computed Tomography (perfusion CT). This technique has already found its way into clinical routine where it represents, together with MRI, one of the primary imaging techniques for the diagnosis of patients with symptoms of stroke. It presents several advantages over other techniques to perform perfusion studies as MRI which include the widespread availability of CT scanners, an easier access to the patient and lower cost [2]; moreover, several studies indicate that the potential of this technique is not yet fully exploited [35]. In this paper we propose an acquisition-reconstruction scheme which tackles both the temporal resolution and the noise problem making it suitable for perfusion studies with other tomographic X-ray techniques.

In [6] we introduced an algorithm with increased temporal resolution under the assumption of no motion, which weakened the requirement of fast data acquisition. This algorithm (TIA-TFDK: Temporal Interpolation Approach Tent Feldkamp Davis Kress) is based on a convolution and backprojection approach for circular trajectories and reconstructs an intermediate volume image denoted as partial block backprojection (PBB) from only that part of the projection data that are acquired in a time period of a given length. The final volume image at a certain time is computed by accumulating a complete set of PBBs of the same time, thus reducing data inconsistencies. For this purpose a temporal spline interpolation approach is used to estimate the value of a PBB at any desired time. Image acquisition has to be fast enough such that the PBBs satisfy Nyquist's sampling condition. However, the TIA-TFDK algorithm provides a fixed temporal resolution for a given rotation time of the scanner and it does not take noise into account. In order to tackle this two aspects, we introduce in this a simple modification of the TIA-TFDK algorithm which exploits the fact that perfusion signals are essentially low pass. An estimation of the maximum frequency of the fastest perfusion signal in the region of interest is used as a priori information. By substituting the interpolation step for a more flexible estimation scheme which allows to choose the temporal bandwidth of the output sequence, higher frequency components containing mainly noise are suppressed and thus the SNR is improved. In our modification, which we call Temporal Smoothing Approach-TFDK (TSA-TFDK), we propose to use smoothing splines as a natural evolution from the interpolation splines scheme used in [6]. With the proposed algorithm the temporal bandwidth is an additional parameter which is independent of the rotation time. We also justify theoretically that under certain assumption (no motion and only quantum noise) the image quality can be shown to be the same independently of the rotation time.

This algorithm opens the possibility to use other tomographic X-ray techniques for perfusion studies as, for example, small animal CT scanners and represents a first step toward the application of C-arm systems for perfusion studies. However, the assumption of no motion limits the possible clinical applications to fixed patients.

We start with the setting of our dynamic reconstruction problem. We then describe the low-noise estimation of time dependent projections in the presence of noise when the only a priori information available is the signal's bandwidth. We describe a practical implementation of this low-noise estimation with polynomial spline filters. Subsequently, we apply these ideas to the dynamic reconstruction problem and describe the TSA-TFDK approach for dynamic reconstruction. Finally, we present a numerical simulation and provide an illustrative example of the effect of the noise reduction on clinical data.

2. Problem Setting

A source-detector arrangement with a cylindrical detector as shown in Figure 1 rotates around the 𝑧 axis with a constant angular speed 𝜔 on a circular trajectory. The angular speed 𝜔 determines the rotation time 𝑇2𝜋, 𝑇2𝜋=2𝜋/𝜔. The plane containing the source trajectory is denoted as 𝑥𝑦 plane. We assume that the source is situated at the angular position 𝛼=0 at 𝑡=0; hence, the cone-beam projection 𝑃𝛼(𝛾,𝜑,𝑡𝛼) is acquired at 𝑡𝛼=𝛼/𝜔. 𝛾 and 𝜑 denote the fan and cone-angles, respectively, and 𝛾(𝛼,𝐱), and 𝜑(𝛼,𝐱), the fan and cone-angles which determine the ray passing through 𝐱, where 𝐱3 denotes the spatial coordinate, for a given angular position 𝛼. The object is represented by a time dependent spatial distribution of an attenuation coefficient 𝜇(𝐱,𝑡), where 𝑡 denotes time, and is located within a cylinder of radius 𝐱max, that is, 𝜇(𝐱,𝑡)=0for𝐱𝐱max.(1) We concentrate in changes of the attenuation coefficient caused by the change in the composition of a substance as in the case of contrast agent flow. We assume that neither motion nor deformation occurs.

The source rotates continuously during a total acquisition time of 𝑇tot and during this time a total dose 𝐷tot is applied. We assume that the total dose is distributed uniformly among all projections acquired during the acquisition time. We use as an estimation of the dose applied the mAs product, which is a reasonable measure to compare reconstruction algorithms in terms of dose efficiency as long as the input data are obtained with the same scanner [7]. During its continuous rotation the source may be switched off during regular periods of time.

We also assume that the only source of noise is quantum noise. This is a reasonable assumption for objects with small sections, for example, in neurological applications.

Our goal is the estimation of 𝜇(𝐱,𝑡) at temporal increments of 𝑇fr during a total time of 𝑇tot.

3. Low-Noise Estimation in Dynamic X-Ray Tomographic Imaging

3.1. Low-Noise Estimation of Time Dependent Projections

Let us now consider the problem of the estimation of the continuous projection value 𝑝(𝑡) from the noisy measurements 𝑦[𝑛] taken every 𝑇𝑠.

The continuous projection value 𝑝(𝑡) is an accumulation of time attenuation curves along the direction of the x-ray beam. In order to find an appropriate method to estimate 𝑝(𝑡) it is convenient first to characterize TACs in order to be able to extract the maximum information from the measurements. The model, however, should be as general as possible in order for it to be valid both in pathologic and nonpathologic cases. The morphology of perfusion signals depends on the patient, on the tissue, and on the injection parameters (amount of contrast agent, injection rate, etc.) so that it is very variable. The task to find a model that represents their variety is a very challenging one. The model we propose is motivated by the linear systems approach for the modeling of tracer kinetics [8]. According to this approach, the time-attenuation curve in a tissue is obtained by convolution of the time-attenuation curve of the input artery with the impulse response of the tissue which is essentially a low-pass filter. The flow through several tissues is modeled by successive convolutions. Additionally, prior to reaching tissue, venously injected contrast agent goes through the cardio-pulmonary system which has a very strong low-pass filtering effect. As a consequence of this, TACs are very smooth curves. In Fourier domain, this translates into a Fourier transform that can be neglected over a certain threshold which we denote as 𝜈max. Moreover, the injection profile is usually a rectangular pulse; thus, the power spectrum of TACs is concentrated around 0 Hz. As an example, typical values of 𝜈max measured in the brain at the arteria cerebri anterior with an injection rate of 20 mL/s and an injected volume of 50 mL are 0.15 Hz [9].

According to the model presented in Appendix A, the power spectral density of quantum noise is constant over all frequencies. The natural conclusion is then that suppressing frequency components over |𝜈max|, the signal is preserved and only noise is eliminated, increasing thus the signal to noise ratio. The maximum enhancement of the SNR is obtained by using a continuous ideal low-pass filter with cut-off frequency 𝜈max for the estimation. This principle works only if the sampling condition is fulfilled

𝜈max<12𝑇𝑠,(2) which ensures that the repetitions of the frequency spectrum (light shaded spectra in Figure 2) do not overlap.

The variance of the noise after the estimation 𝜎2Est is obtained by integrating the power spectral density of the noise in the frequency band ]𝜈max,𝜈max[. Since the power spectral density of noise is constant equal to 𝜀2𝑇𝑠 and the frequency response of the ideal low-pass filter is constant in ]𝜈max,𝜈max[, we get

𝜎2Est=2𝜈max𝜀2𝑇𝑠.(3) Note that in order to come to (3) we have only used the hypotheses of quantum white noise (see Appendix A), bandlimited signal, and constant total dose. In the following we consider the influence of dose settings in the noise variance.

According to the definition of dose, the mAs product dose is linearly accumulated. Hence, the total dose applied 𝐷tot is

𝐷tot=𝐷𝑝𝑁rot𝑁𝛼,(4) where 𝐷𝑝 is the dose applied to measure each projection, 𝑁rot is the number of rotations with activated X-ray source, and 𝑁𝛼 is the number of projections acquired per rotation. Every projection is measured once during a rotation, hence 𝑁rot, 𝑇𝑠, and the total acquisition time are related by

𝑁rot=𝑇tot𝑇𝑠.(5) Note that 𝑇𝑠 need not to be 𝑇2𝜋 but might be a multiple of it. Since the total dose is distributed uniformly among all projections acquired, the dose applied to measure a projection is

𝐷𝑝=𝐷tot𝑁rot𝑁𝛼.(6) It is a well-known fact that for quantum noise the mean square value is inversely proportional to the dose applied 𝜀21/𝐷𝑝. Using (5) and (6), we get

𝜀2𝑇tot𝑁𝛼𝑇𝑠𝐷tot.(7) Finally, using (3) we get

𝜎2Est𝜈max𝐷tot𝑇tot𝑁𝛼.(8) This expression can be interpreted as follows. The estimation with a continuous low-pass filter exploits the redundancy in the signal to reduce noise, and 1/𝜈max is a good indicator for signal redundancy. Hence, for a given total dose 𝐷tot distributed during a total acquisition time 𝑇tot and acquiring 𝑁𝛼 projections per rotation, the higher the redundancy is (i.e., the lower 𝜈max), the stronger the achieved noise reduction will be. Note that 𝜎2Est does not depend on 𝑇𝑠. Thus the same level of noise will be obtained by acquiring a few scans with a high dose per view or a higher amount of scans with lower dose per view, as long as condition (2) is fulfilled.

3.2. Low-Noise Estimation with Polynomial Spline Filters

As stated in the introduction, we use the TIA-TFDK algorithm from [6] as a starting point for our new approach. Hence, the most natural choice to approximate the ideal low-pass filters described in the previous section is to use polynomial spline filters, since this will imply a minor modification in the algorithm.

An ideal low-pass filter has an infinite support in temporal domain and decays very slowly (1/|𝑡|) so that samples that lie far from the value to be estimated will still contribute significantly to it. This is very inefficient for practical purposes. Polynomial spline filters are appropiate to approximate ideal low-pass filters for our purposes for two main reasons. First, from [6] as a starting point for our new approach, hence this is a natural choice since it implies a minor modification in the algorithm. Additionally, polynomial splines provide a viable compromise between accuracy and efficiency. The filter kernels are quite local and converge rapidly to the ideal low-pass kernel when the order 𝑛 of the polynomial splines increases. Given a set of samples, by properly choosing the coefficients of the polynomial spline basis functions (B-splines), estimations of any degree of smoothing and, in the extreme, of no smoothing at all (interpolation) can be performed [10, 11]. In this section we describe how temporal smoothing with polynomial splines can be interpreted as a low-pass filtering in time.

We denote with 𝑦[𝑘] the temporal noisy samples of a given signal 𝑝(𝑡) taken every 𝑇𝑠. Throughout this section 𝑡 denotes time normalized by the sampling interval 𝑇𝑠, that is, 𝑡=𝑡/𝑇𝑠; the corresponding dimensionless frequency is denoted by 𝜈=𝜈𝑇𝑠. We approximate the continuous signal 𝑝(𝑡) by a polynomial spline function of order 𝑛: 𝑠𝑛𝜆(𝑡), where 𝜆 indicates the degree of smoothing. According to Schönberg's theorem (see [10]), 𝑠𝑛𝜆(𝑡) can be expressed as 𝑠𝑛𝜆𝑡=𝑘𝑐𝜆[𝑘]𝛽𝑛𝑡𝑘,(9) where 𝛽𝑛(𝑡) are the polynomial B-splines and 𝑐𝜆[𝑘] the B-spline coefficients. The polynomial B-splines are uniquely determined by its order 𝑛 so that, for a given order 𝑛, we only need to find the coefficients 𝑐𝜆[𝑘]. For this purpose we minimize the following functional:

𝐹𝐜𝝀,𝐲=𝑘𝑦[𝑘]𝑠𝑛𝜆[𝑘]2+𝜆𝐷𝐿𝑠𝑛𝜆(𝑡)2𝐿2,(10) where 𝑛=2𝐿1, 𝐷𝐿 denotes the 𝐿th derivative and 𝐿2 the 𝐿2 norm. The first term in (10) forces the estimated function to be close to the sampled values at the sampling points. The second term is a regularity constraint which favors a smooth estimation of the signal. It is controlled by the smoothing parameter 𝜆. For 𝜆=0, 𝐹{𝐜𝝀,𝐲}=0 if 𝑠𝑛𝜆[𝑘]=𝑦[𝑘], which is the interpolation condition. For larger values of the smoothing parameter, the smoothness constraint might determine a curve not passing through the sampled values. If we denote 𝑏𝑛[𝑘]=𝛽𝑛(𝑘), it can be shown (see [11]) that the coefficients of the 𝑛th-order polynomial spline that minimize (10) can be computed as

𝐶𝜆𝑒𝑖2𝜋𝜈=̂𝑌𝑒𝑖2𝜋𝜈𝐵𝑛𝑒𝑖2𝜋𝜈+𝜆𝑒𝑖2𝜋𝜈+2𝑒𝑖2𝜋𝜈𝐿,(11) where 𝐶𝜆(𝑒𝑖2𝜋𝜈),̂𝑌(𝑒𝑖2𝜋𝜈)- and 𝐵𝑛(𝑒𝑖2𝜋𝜈) denote the Fourier Transform of the discrete sequences 𝑐𝜆[𝑘], 𝑦[𝑘], and 𝑏𝑛[𝑘], respectively. Substituting the spline coefficients 𝐶𝜆(𝑒𝑖2𝜋𝜈) of (11) in the Fourier Transform of (9) we obtain

̂𝑠𝑛𝜆𝜈=̂𝑌𝑒𝑖2𝜋𝜈sinc𝜈𝑛+1𝐵𝑛𝑒𝑖2𝜋𝜈+𝜆𝑒𝑖2𝜋𝜈+2𝑒𝑖2𝜋𝜈𝐿̂𝜂𝑛𝜆(𝜈),(12) where we have used that ̂𝛽𝑛(𝜈)=(sinc(𝜈))𝑛+1.

Equation (12) shows that the smoothing operation can be interpreted as the discrete convolution of the samples 𝑦[𝑘] with a continuous polynomial spline low-pass filter 𝜂𝑛𝜆(𝑡). When 𝜆=0, 𝜂𝑛0(𝑡) is the polynomial spline interpolator of order 𝑛 and it verifies lim𝑛𝜂𝑛0(𝑡)=sinc(𝑡) [12].

We denote as cut-off frequency the frequency 𝜈𝑐 at which the frequency response of the filter falls to half of the maximum, that is, ̂𝜂𝜆(𝜈𝑐)=0.5. Figure 3 shows the frequency responses of different polynomial spline filters 𝜂𝑛𝜆(𝑡) for different values of 𝑛 and 𝜆. By properly choosing these parameters we can obtain good approximations to the ideal low-pass filters used in Appendix A. In practice, 𝑛 is responsible for the sharpness of the edges of the frequency response, that is, how close the filter is to the ideal low-pass filter (Figure 3(a)) and 𝜆 for the position of the cut-off frequency (Figure 3(b)).

Even if 𝜆 is responsible for the position of the cut-off frequency, this dependency is different for different orders 𝑛. This is illustrated in Figure 4. For 𝜆=0, 𝜈𝑐=0.5 for all 𝑛 values. The cut-off frequency decreases slower for lower orders until 𝜆1. After this value, lower orders decrease much faster. For high values of 𝜆 the decay is very slow and therefore large increases in the smoothing parameter provide only very small reduction of the cut-off frequency. The value of 𝜆, the normalized cut-off frequency 𝜈𝑐 and the spline of order 𝑛 are related by [13]


Contrary to the frequency response of the ideal low-pass filter, the frequency response of polynomial spline filters starts decaying already for frequencies under 𝜈𝑐. The proportion of the frequency interval [0,𝜈𝑐] which can be reproduced accurately depends mainly on the order of the splines 𝑛. For the purpose of this paper we will use 𝑛=9 as an approximation for the ideal low-pass filter. Figure 3 shows that for 𝑛=9 the frequency response is almost constant until approximately 0.8𝜈𝑐. As a consequence, we will use hereafter the following expression to calculate the cut-off frequency to ensure that the frequency range from 0 to 𝜈max remains unaltered after filtering


In Section 3.1, we provided an expression for the estimation of the variance of the filtered sequence. We will derive now an expression for the case where polynomial spline filtering is used instead of ideal filtering. The power spectrum of a filtered stationary process is obtained by integrating the product of the power spectrum of the process with the square of the absolute value of the frequency response of the filter [14]. Since we used an ideal low-pass filter, the integral from 𝜈max to 𝜈max could simply be substituted by 2𝜈max times the value of the power spectral density. If we use a spline low-pass filter this is no longer true. We now have for(3)

𝜎2Spl=𝜀2𝑇𝑠+||̂𝜂𝑛𝜆𝜈||2𝑑𝜈.(15) The value of the integral depends on 𝑛 and 𝜆. In Figure 3, we show the frequency response of 𝜂9𝜆(𝑡) for different values of 𝜆. Qualitatively it is clear that all these filters are very close to an ideal low-pass filter with the corresponding cut-off frequency. Hence, the value of the integral must be close to 2𝜈𝑐. In Figure 5 we show the ratio

+||̂𝜂9𝜆(𝜈)||2𝑑𝜈2𝜈𝑐(16) for values of 𝜈𝑐[0.02,0.5]. The values are all between 0.914 and 0.95. We would like to have a simple rule of thumb to estimate the value of the variance of a sequence filtered with a low-pass spline filter. For this purpose, we propose to use the average value of the curve in Figure 5, 0.92 as a representative value and then add the factor 0.92/0.8 to (3) as a correction factor for filtering with splines. This yields


4. Temporal Smoothing Approach-TFDK

In this section we discuss the modification of the TIA-TFDK algorithm to exploit the ideas presented in the previous sections. For the sake of simplicity we initially assume that the source is always switched on during acquisition time, that is, it operates in continuous acquisition mode.

The TIA-TFDK algorithm can be seen as a dynamic extension of the T-FDK algorithm for cone-beam reconstruction on circular trajectories described in [15]. We first introduce this algorithm to which we have added the temporal dependence of the object. Unclear relationships between time variables will be explained further on. The T-FDK algorithm consists in first rebinning the cone-beam projections to a fan-parallel beam (see Figure 6)

𝑃𝛼𝛾,𝜙,𝑡𝛼𝑃𝑏𝛽𝑢,𝑣,𝑡𝛽,(18) and then performing filtering and backprojection on the rebinned projections:

𝑄𝛽𝑢,𝑣,𝑡𝛽=𝑢max𝑢max𝑃𝑏𝛽𝑠,𝑣,𝑡𝛽𝑤(𝑠,𝑣)𝑔(𝑢𝑠)𝑑𝑠,𝜇𝑟𝐱,𝑡𝜋=122𝜋0𝑄𝛽𝑢(𝐱,𝛽),𝑣(𝐱,𝛽),𝑡𝛽𝑑𝛽,(19) where 𝑤(𝑢,𝑣) is a weighting function which only depends on the detector coordinates and (𝑢(𝐱,𝛽),𝑣(𝐱,𝛽)) are the coordinates of the detector pixel where the ray passing through 𝐱 intersects the detector for the projection at projection angle 𝛽 (see Figure 6).

Both the rebinning and the backprojection steps include approximations since in both, projections acquired at different times are used to compute a value which is associated to a unique time. During the rebinning step, in order to calculate the rebinned value at projection angle 𝛽, cone-beam projections acquired during the time interval 𝑡𝛼[𝑡𝛽𝛾max,𝑡𝛽+𝛾max] are used, where 𝛾max is the maximum fan-angle. Hence, the time uncertainty introduced in the rebinning step depends on the maximum fan-angle, that is, on the diameter of the object, whereas the time uncertainty due to the backprojection step depends on the length of the backprojection interval. We can therefore reduce the uncertainty due to backprojection by dividing the angular interval [0,2𝜋] in 𝑁 subintervals of length 2𝜋/𝑁. We denote the backprojection over these subintervals as partial block backprojections (PBBs) as in [16]. Since every PBB occupies the same amount of memory as the reconstructed volume, taking a lot of PBB intervals has a direct impact on the amount of memory required for the reconstruction. On the other hand, as shown in [6], the uncertainty due to rebinning does not depend on the number of PBBs per rotation 𝑁; there is a certain threshold 𝑁max, which only depends on the maximum fan-angle, over which the error due to rebinning predominates and makes it not worth increasing 𝑁, that is, to reduce the length of the backprojection interval. For perfusion signals, [6] shows that with the typical fields of view of most current CT scanners and 𝑁𝑁max, the errors due to the uncertainties introduced by rebinning and backprojection can be neglected when compared to the errors due to resolution or noise. In the TSA-TFDK, the reconstruction of a PBB over several rotations can be interpreted as the sampling of time dependent PBB with a sampling interval 𝑇𝑠=𝑇2𝜋. Furthermore, in the rebinned geometry, for small cone-angles, equivalent rays are acquired every 𝑇𝑠=𝑇𝜋. That is to say, reconstructing PBBs over several rotations provides a series of noisy samples acquired every half-rotation.

In Appendix A we assumed that the variances of the projection values do not depend on time and derived a model for the temporal behavior of noise in projections based on this assumption. It should be noted, however, that if this assumption holds for projection values, it automatically holds for filtered projection values and also for linear combinations thereof, for example, to PBBs. Similarly, the statistical independence of the fluctuations at different times is guaranteed for filtered projections and linear combinations of projections as long as the projection sets are not overlapping. As a consequence, the proposed model for the temporal behavior of noise in projections, as a stationary process in time with statistically independent time samples can be also applied to filtered projections and partial block backprojections. Therefore, a low-noise estimation is obtained from the series of noisy samples if these are filtered along the time axis with a continuous polynomial spline filter with cut-off frequency 𝜈𝑐=𝜈max/0.8. According to the sampling rate of 1/𝑇𝜋 obtained when combining the time-series of the 𝑗th PBB with the time series of the 𝑗+𝑁/2th PBB (computed with equivalent rays), accurate reconstruction is possible as long as 𝜈𝑐=𝜈max0.812𝑇𝜋𝑇2𝜋0.8𝜈max.(20) By filtering the PBBs this way, we can estimate their values at any time and thus reconstruct an image frame at any desired reconstruction time by accumulating 𝑁/2 PBBs. Contrary to the TIA-TFDK algorithm, the temporal estimation is carried out using polynomial spline filters with a normalized cut-off frequency 𝜈𝑐0.5, that is, the PBB samples are smoothed along the time axis.

Note that (20) is based on the assumption that we use polynomial splines of order 𝑛=9. If other values of 𝑛 are used, (20) should be adapted accordingly (see Section 3.2).

A detailed description of the steps of the algorithm is provided in Appendix B.

Since typically 𝜈max<0.15 Hz (see [9]), the condition (14) translates into the rotation time having to fulfill 𝑇2𝜋<5.33 s, which is achieved by a wide variety of X-ray tomographic imaging systems. For the choice of 𝑁 the most challenging case according to (20) is when 𝑇2𝜋𝜈max=0.8, for which we provided the following empirical formula in [6]:

𝑁max𝜋610.02+𝑟2/4Γ2,(21) where 𝑟 is the radius of the object and Γ is the radius of the circular trajectory of the source. This 𝑁max gives the value over which the uncertainty introduced by rebinning predominates over the uncertainty produced by backprojection.

A especially interesting case results from the application of the algorithm to scanners with high rotational speeds (e.g., 2 rotations per second) because then using the lowest rotation time of the scanner 𝑇min2𝜋, 𝑁, can be reduced to 1 without significant loss of accuracy. Thus, the temporal estimation is carried out as a postprocessing of the reconstructed sequence.

5. Numerical Example

In this section we present an example of low-noise acquisition and reconstruction with the TSA-TFDK algorithm and compare it to a perfusion protocol which is representative of current clinical routine perfusion CT. During this section, we will refer to this protocol as the reference protocol. For the acquisition, 800 projections per rotation were simulated with a cylindrical scanner of 256 channels. A spherical phantom containing 6 spherical inserts with time-dependent attenuation was used. In this phantom, the inserts follow the temporal law

𝜇𝑖(𝑡)=𝐶𝑖𝑡𝑝1𝑝2𝑒(𝑡𝑝1)/𝑝3,(22) with 𝑝1=5 s, 𝑝2=2.3, and 𝑝3=3 s. This curve fulfills 𝜈max0.15 Hz and is representative of the temporal evolution of the concentration of contrast agent in vessels and tissue after the injection of a bolus of contrast agent. All inserts follow the same law, except for the amplitude. 𝐶𝑖 (𝑠2) determines the amplitude and was chosen in such a way that the maximum values of the enhancement were 10,18,26,34,42, and 50 HU.

We first simulated acquisition and reconstruction with the reference protocol, consisting in full-scan T-FDK reconstructions [15]. 𝑁sc=40 scans were simulated with a rotation time of 𝑇2𝜋=0.5 s during 40 s. The source was switched off every second rotation. Quantum noise was added to the projections. The corresponding values for the rotation time and smoothing parameter were calculated using (B.1) and (B.2) and are given in Table 1. The lowest sampling rate implies a slow rotation scanning, in which the signal is not oversampled and therefore 𝜆=0. In this case, TSA-TFDK reduced to the TIA-TFDK algorithm as described in [6]. The highest sampling rate implies a fast rotation scanning such that the signal is strongly oversampled and therefore a high value of 𝜆 is needed to adjust the cut-off frequency to 𝜈max. The noise parameter for each simulation was adjusted according to the number of scans in such a way that the accumulated total dose was kept constant.

Figures 7(a) and 7(b) show a frame of the reconstructed sequence with the reference protocol (left) and with TSA-TFDK with fast scanning (right). The noise reduction effect can be clearly observed. In Figure 8, we show the temporal evolution of the mean of the reconstructed attenuation values within the insert with maximum enhancement 18 HU. The curves of the TSA-TFDK sequences are clearly smoother than the curves obtained with the reference protocol, but the temporal resolution is preserved. The results with fast and slow scanning are qualitatively similar, except for some oscillations towards the end of the sequence in the slow scanned sequence, probably as a consequence of a light aliasing.

In Table 2, we show the value of the standard deviation measured within the inserts. As expected, the standard deviation in the TSA-TFDK reconstructions is lower than in the reference protocol sequence. The TSA-TFDK reconstructions of slow and fast scans exhibit a similar noise level. Indeed, according to (A.5) and (17) the ratio of the variances should be 𝑟𝜎212.3𝜈max=2.898.(23) With the values from Table 2 we can calculate the measured reduction of the variance as 4.48/2.62=1.71 and (1.71)2=2.924, which shows that TSA-TFDK behaves as expected.

6. Results on Clinical Data

In this section, we show the effect of the low-noise reconstruction on a clinical data set. For this purpose, we present an example where the temporal smoothing is performed on the reconstructed images of a neurological clinical dataset. While this could be seen as a particular case of the TSA-TFDK algorithm for 𝑁=1, we would like to remark that the TSA-TFDK has been developed for slow rotating scanners. This example is presented to illustrate the effect of the low-noise reconstruction and compare to current methods.

6.1. Data and Method

In the perfusion protocol data were acquired for 40 s, during which the source rotated at 0.5 s/rot and was switched on every two rotations. Hence, 40 images were generated which corresponds to a rate of 1 image/s. The scan was performed with a tube voltage of 120 kVp and at a tube current of 220 mA. Since the rotation time was 𝑇2𝜋=0.5 s, the dose per reconstructed image in mAs was 𝐷2𝜋=110 mAs. Every image frame was reconstructed from projection data of a full-scan with the reconstruction algorithm provided by the commercial scanner. At this point, we would like to point out that this study does not depend on the reconstruction algorithm, since data are processed after reconstruction. The X-ray beam was collimated to obtain a slice width of 10 mm.

Since we had already a reconstructed sequence at our disposal, we used it to estimate a value of 𝜈max instead of using an a priori general estimation for it. Note that this was possible because the sampling interval was 𝑇𝑠=1 s which satisfies the sampling condition since 𝜈max0.5 Hz [9]. For the purpose of estimating 𝜈max, we used an ROI within the arteria cerebri anterior. The temporal evolution of the concentration of contrast agent in this artery can be assumed to be the fastest in the dataset. Additionally, the arteria cerebri anterior is approximately orthogonal to the slice plane so that partial volume effects are avoided. With this TAC we estimated 𝜈max as the value for which [𝜈max,𝜈max] contains 99% of the energy of the signal. We obtained a value of 𝜈max0.0966 Hz.

According to this estimation, the signal is clearly oversampled and there is leeway for low-pass filtering to increase the SNR as discussed in Section 3.2. Particularly, using (14) we get

𝜈𝑐=0.0966Hz×𝑇𝑠0.80.12(24) for which 𝜆=15.82 (see (13)).

As we saw in Section 3.1, this low-noise level can also be achieved by concentrating the same total dose to fewer rotations. A simple way to simulate sequences with different sampling rates is to downsample the original sequence by a factor of 𝐾, that is, to keep one every 𝐾 samples of the original sequence. The maximum value of 𝐾 can be obtained from (14) taking into account that 𝜈𝑐0.5,

𝜈𝑐=𝜈max𝐾𝑇𝑠0.80.5𝐾4.14.(25) Hence, according to the analysis of Section 3.1, we should be able to obtain complete sequences with the same low level of noise from sequences sampled every 1, 2, 3, or 4 seconds. Using (13) and (25) we get 𝜈𝑐 and 𝜆 for every sequence, these values are shown in Table 3.

Note that since we do not increase the dose by the same downsampling factor, the total dose applied 𝐷tot is reduced by the downsampling factor 𝐾. The only way to avoid this would be to perform the acquisition of every sequence independently so that the dose per rotation is increased by the corresponding value for every sequence. Since this would be highly invasive for the patient, we use here a more simple way to verify the validity of our approach. We simply exploit the fact that the variance of the noise is inversely proportional to the dose applied as shown in (8). Thus according to our model, if we take as a reference the sequence with 𝐾=1, the variance of the other sequences will be 𝜎2𝐾=𝐾𝜎2Spl2.3𝜀2𝐾𝑇𝑠𝜈max.(26) The variance of the noise of the original data is 𝜎2𝑂𝑟=𝜀2. Hence, we can express the ratio between the variance of the noise in the original sequence and that of the downsampled and filtered sequences as 𝜎2𝑂𝑟𝜎2𝐾12.3𝐾𝑇𝑠𝜈max.(27)

6.2. Results

In Figure 9, we show some examples of time-attenuation curves of the original (gray) and the low-noise sequence 1 with 𝐾=1 (black). The curves show the temporal evolution of the average value of the pixels within a ROI in different tissues. As can be seen, the enhancement curve for the arteria cerebri anterior in the low-noise sequence is smooth but can still follow the changes in the original curve. The same can be observed for the sinus sagittalis. In such large vessels, the enhancement due to contrast agent flow is higher than in tissue or in small vessels, and the original SNR is high enough to compute physiological parameters without smoothing. In other regions, enhancement is larger than noise but of the same order of magnitude. An example for such regions is shown in curve Figure 9(c)) where the TAC of an ROI within a tumor is shown. The TAC of the low-noise estimation clearly eliminates noise and delivers a curve which is physiologically more plausible than the curve obtained from the original dataset. Finally, Figure 9(d)) shows the TAC of an ROI within gray matter.

While in the original curve the enhancement cannot really be perceived, the TAC of the low-noise sequence clearly shows an increase with time in the concentration of contrast agent. On the other hand, the TAC of the low-noise sequence contains low frequency oscillations which are not of physiological nature. These are due to low-frequency noise in the frequency band [𝜈max,𝜈max].

We would like now to verify if the reduction of the noise caused by the temporal smoothing is in accordance with the predictions of our model. Noise is measured as the variance of pixel values in a certain region of a homogeneous object. The variance might be due not only to quantum noise but also to pixels in the region of interest that correspond to different tissues and have therefore a different temporal behavior. Since our model concerns the temporal behavior of the statistical fluctuations caused by quantum noise only, we must segment a region within a tissue where all points have the same temporal behavior. This is not practicable in clinical images since we do not know a priori the exact temporal behavior in any region, except for in those regions where the attenuation is constant because contrast agent does not reach them. For this reason we used a segmented ROI within the ventricular system to estimate the variance. This delivered a set of 2300 pixels with the same temporal behavior. We then used the whole 2300×40 pixels to estimate the variance in each reconstructed sequence. The ratio of the variance of the original sequence to the variance of the postprocessed sequences is shown in the second column of Table 4. The third column shows the values estimated using (27). The estimated values match approximately the measured values. According to the proposed model for the temporal behavior of noise, the reduction of the variance should be inversely proportional to the sampling rate since the total dose is decreased by downsampling. This is approximately in accordance with the values in Table 4. The effect on image quality is shown in Figure 10. The frame shown corresponds to 𝑡=19 s. Sequence 4 has samples at 𝑡=16 s and 𝑡=20 s so that the frame at 𝑡=19 s has been computed at the filtering step. The image quality of the frame estimated from sequence 4 appears to be equivalent to the original frame although four times fewer input data were used for the computation of the sequence. Sequence 1, with strong temporal smoothing, presents substantially reduced noise while visually preserving spatial resolution. The reconstructed sequences were used as input for the Perfusion CT software (Siemens AG, Healthcare Sector, Forchheim, Germany) that computes the functional parameter maps. This software first segments vessels and bone, performs a strong spatial smoothing, and subsequently computes the functional parameters. The segmented regions are excluded from the functional maps and represented in black. Figure 11 shows cerebral blood flow maps computed from the original sequence (a), sequence 4 (b), and sequence 1 (c). Map (a) presents many small isolated segmented regions compared to map (c). These correspond to areas where the noise level led to a wrong classification as vessels. Most of them disappear in map (c). Finally, map (c) is smoother in space which is physiologically more plausible. The quality of maps (a) and (b) is equivalent although (b) was computed from a sequence reconstructed with four times fewer data, equivalent to four times less dose.

7. Conclusion

In this paper we have presented the dynamic reconstruction algorithm TSA-TFDK, which is a further development of the previously published TIA-TFDK algorithm [6]. A temporal sequence of image frames is generated. In the presented modification a smoothing term is incorporated into the temporal estimation step reducing the noise level. The continuous low-pass filter covers only the relevant bandwidth of the dynamic process and higher frequencies are strongly attenuated. No relevant information is lost. Image properties are independent of the sampling rate and the rotation time of the scanner as long as Nyquist's sampling condition is fulfilled. The obtained level of noise only depends on the total dose applied. As an example the TSA-TFDK allows to reconstruct a dynamic perfusion process with a physiologic relevant maximum frequency of 𝜈max=0.15 Hz from data acquired by scanning with any rotation time up to 5.33 s. In a simulation study similar temporal resolution and noise level are achieved by acquiring with slow and fast rotation time. We have illustrated the effect of the low-pass filtering using a clinical data set, and the results obtained were in accordance with the predictions obtained from a model for the temporal behavior of noise in dynamic X-ray imaging presented. The proposed approach opens the possibility to use other tomographic X-ray techniques such as small animal imaging CT scanners for perfusion studies. It might even be a starting point for the development of dynamic reconstruction algorithms for C-arm systems by adapting the estimation step to the irregular sampling scheme caused by the back and forth movement. However, the assumption of no motion limits the possible clinical applications to fixed patients.


A. Temporal Sampling of Noise in Dynamic X-Ray Tomographic Imaging

We consider the temporal evolution of a single noisy projection value which we denote by 𝑦(𝑡). This value is determined by the angular position of the source, the fan-angle and the cone-angle:

𝑦(𝑡)=𝑃𝛼(𝛾,𝜑,𝑡).(A.1) The temporal evolution 𝑦(𝑡) is modeled as a continuous deterministic signal 𝑝(𝑡) contaminated by additive noise 𝜀(𝑡)

𝑦(𝑡)=𝑝(𝑡)+𝜀(𝑡).(A.2) The noise is assumed to be a continuous random process with zero mean and mean square value 𝜀2.

For the following analysis we consider quantum as the main source of noise in X-ray projections and neglect the contribution of electronic noise. This is a reasonable assumption for applications where the maximum width traversed by the X-ray beam is small as, for example, neurological applications. According to the Poisson law the mean number of quanta detected by a detector pixel fluctuates around a mean value 𝑁 with a variance equal to the mean. As a consequence of the process of contrast agent flow, the mean number of detected quanta depends on time and therefore so does the noise variance. However, the increase in attenuation due to contrast agent flow is small compared to the overall attenuation of the passing X-ray beam. Hence, it is reasonable to assume that the noise variance does not depend on time (is constant) [9] and thus 𝜀(𝑡) is a wide-sense stationary process in time [17]. We use hereafter the term stationary process to refer to a wide-sense stationary process. Moreover, for a given projection value, the fluctuations around the mean at two time instants 𝑡1 and 𝑡2 are statistically independent.

In dynamic acquisition the projection value 𝑦(𝑡) is measured once during every rotation. This can be seen as a temporal sampling of the signal 𝑝(𝑡) contaminated by noise 𝜀(𝑡). Let us now analyze the sampling of stationary process 𝜀(𝑡). We denote by 𝜀1[𝑛]=𝜀(𝑛𝑇𝑠) the discrete sequence of noise samples obtained by sampling 𝜀(𝑡) with a sampling rate of 1/𝑇𝑠. As mentioned above, the fluctuations at two different times are statistically independent and therefore also uncorrelated. This can be expressed by the discrete autocorrelation function 𝑅𝜀1 as 𝑅𝜀1[𝑛]=𝐸𝜀1[𝑘]𝜀1[𝑘+𝑛]=𝜀2𝛿[𝑛],(A.3) where 𝐸 denotes the expectation value and 𝛿[𝑛] is the Kronecker symbol (𝛿[0]=1 and 𝛿[𝑛]=0 if 𝑛0). The discrete power spectral density 𝑆𝜀1 is defined as (see [18])

𝑆𝜀1𝑒𝑖2𝜋𝜈=𝑛𝑅𝜀1[𝑛]𝑒𝑖2𝜋𝜈𝑛=𝜀2=const(A.4) where 𝜈=𝜈𝑇𝑠 is the normalized frequency with 𝜈 as the physical frequency. The notation (𝑒𝑖2𝜋𝜈) indicates that it is the Fourier Transform of a discrete sequence and is therefore periodic (we follow the notational conventions of [18]). 𝑆𝜀1(𝑒𝑖2𝜋𝜈) owes its name to the fact that integrating it over one period of the normalized frequency, ]1/2,1/2[, yields the mean square value or variance of the noise. Downsampling 𝜀1[𝑛] by a factor of 𝐾 (or sampling 𝜀(𝑡) with sampling interval 𝐾𝑇𝑠) would deliver a sequence with the same discrete autocorrelation as in (A.3) and therefore the same discrete power spectral density as in (A.4).

Definition (A.4), however, does not provide any insight into the physical frequency of the underlying continuous process. For this reason we define the discrete physical power spectral density 𝑆𝜀1 of the samples of a continuous process sampled every 𝑇𝑠 by introducing the change of variables 𝜈=𝜈𝑇𝑠 in (A.4). Thus, we get

𝑆𝜀1𝑒𝑖2𝜋𝜈𝑇𝑠=𝑇𝑠𝑛𝑅𝜀1[𝑛]𝑒𝑖2𝜋𝜈𝑛𝑇𝑠=𝜀2𝑇𝑠.(A.5) Integrating (A.5) over one period on the physical frequency axis we obtain the mean square value or variance of the discrete process

𝜀21=1/2𝑇𝑠1/2𝑇𝑠𝑆𝜀1𝑒𝑖2𝜋𝜈𝑇𝑠𝑑𝜈=𝜀2.(A.6) Therefore 𝑆𝜀1(𝑒𝑖2𝜋𝜈𝑇𝑠) describes the distribution of power density over one period of the physical frequency axis. In order to simplify nomenclature, we will use hereafter the wording power spectral density to refer to the physical power spectral density.

As a consequence of (A.5) and (A.6), the discrete power spectral density of the discrete sequence 𝜀𝐾[𝑛], obtained by downsampling 𝜀1[𝑛] by a factor of 𝐾 (or sampling 𝜀(𝑡) with sampling interval 𝐾𝑇𝑠), is increased by a factor of 𝐾 with respect to the discrete power spectral density of 𝜀1[𝑛] (see Figure 12)



The TSA-TFDK is summarized in the following steps.

Algorithm TSA-TFDK
(1) Choose rotation time such that 𝑇2𝜋0.8𝜈max.(B.1)
(2) Calculate 𝜆 for 𝑛=9𝜆=𝜋𝜈max𝑇2𝜋0.810𝜋10.(B.2)
(3) Rebin cone-beam projections to fan-parallel beam 𝑃𝛼𝛾,𝜙,𝑡𝛼𝑃𝑏𝛽𝑢,𝑣,𝑡𝛽.(B.3)
(4) Divide every rotation in 𝑁 sets of projections and reconstruct 𝑁 partial block backprojections with the T-FDK algorithm 𝑄𝛽𝑢,𝑣,𝑡𝛽+𝑘𝑇2𝜋=𝑢max𝑢max𝑃𝑏𝛽𝑠,𝑣,𝑡𝛽+𝑘𝑇2𝜋𝑤(𝑠,𝑣)𝑔(𝑢𝑠)𝑑𝑠,𝜁𝑗𝐱,𝑡(𝜋/𝑁)(2𝑗+1)+𝑘𝑇2𝜋=12(2𝜋/𝑁)(𝑗+1)(2𝜋/𝑁)𝑗𝑄𝛽𝑢(𝐱,𝛽),𝑣(𝐱,𝛽),𝑡𝛽+𝑘𝑇2𝜋𝑑𝛽,(B.4) where for all 𝑗=0,,𝑁1, for all 𝑘.
(5) Combine the 𝑗th and the (𝑗+𝑁/2)th sequences of partial block backprojections to a unique sequence 𝜁𝑗𝐱,𝑡(𝜋/𝑁)(2𝑗+1)+𝑘𝑇𝜋=𝜁𝑗𝐱,𝑡(𝜋/𝑁)(2𝑗+1)+𝑘𝑇2𝜋,𝑘even,𝜁𝑗+(𝑁/2)𝐱,𝑡(𝜋/𝑁)(2𝑗+1)+𝑘𝑇2𝜋𝑇𝜋,𝑘odd,(B.5) where for all 𝑗=0,,𝑁/21, for all 𝑘.
(6) Estimate the values of the partial reconstructions at the output frame times with polynomial spline filtering 𝜁𝑗(𝐱,𝑡)=𝑘𝜁𝑗𝐱,𝑡(𝜋/𝑁)(2𝑗+1)+𝑘𝑇𝜋×𝜂𝑛𝜆𝑡𝑡(𝜋/𝑁)(2𝑗+1)𝑘𝑇𝜋𝑇𝜋,(B.6) where for all 𝑗=0,,𝑁/21.
(7) Accumulate the estimated 𝑁/2 partial block backprojections to obtain the time dependent attenuation map 𝜇𝑟(𝐱,𝑡)=2𝑁/21𝑗=0𝜁𝑗(𝐱,𝑡).(B.7)


P. Montes gratefully acknowledges the financial support by Siemens AG Healthcare Sector, Forchheim, Germany. The authors would like thank Dr. Peter Schramm (Department of Neuroradiology of the University of Heidelberg Medical School) for providing the clinical dataset used in Section 6.