Research Article | Open Access
Pau Montes, Günter Lauritsch, "Low-Noise Dynamic Reconstruction for X-Ray Tomographic Perfusion Studies Using Low Sampling Rates", International Journal of Biomedical Imaging, vol. 2009, Article ID 108028, 13 pages, 2009. https://doi.org/10.1155/2009/108028
Low-Noise Dynamic Reconstruction for X-Ray Tomographic Perfusion Studies Using Low Sampling Rates
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.
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  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 ; moreover, several studies indicate that the potential of this technique is not yet fully exploited [3–5]. 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  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 . 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 , . The plane containing the source trajectory is denoted as plane. We assume that the source is situated at the angular position at ; 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 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 , that is, 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 and during this time a total dose 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 . 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 during a total time of .
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 . 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 . Moreover, the injection profile is usually a rectangular pulse; thus, the power spectrum of TACs is concentrated around Hz. As an example, typical values of measured in the brain at the arteria cerebri anterior with an injection rate of mL/s and an injected volume of mL are 0.15 Hz .
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 , 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 for the estimation. This principle works only if the sampling condition is fulfilled
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 is obtained by integrating the power spectral density of the noise in the frequency band . Since the power spectral density of noise is constant equal to and the frequency response of the ideal low-pass filter is constant in , we get
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 is
where is the dose applied to measure each projection, 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 , , and the total acquisition time are related by
Note that need not to be 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
Finally, using (3) we get
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 is a good indicator for signal redundancy. Hence, for a given total dose distributed during a total acquisition time and acquiring projections per rotation, the higher the redundancy is (i.e., the lower ), the stronger the achieved noise reduction will be. Note that 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  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 () 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  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 ), can be expressed as 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:
where , denotes the th derivative and the 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 , 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 ) that the coefficients of the th-order polynomial spline that minimize (10) can be computed as
where we have used that .
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 , is the polynomial spline interpolator of order and it verifies .
We denote as cut-off frequency the frequency at which the frequency response of the filter falls to half of the maximum, that is, . 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 , for all values. The cut-off frequency decreases slower for lower orders until . 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 
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 which can be reproduced accurately depends mainly on the order of the splines . For the purpose of this paper we will use as an approximation for the ideal low-pass filter. Figure 3 shows that for the frequency response is almost constant until approximately . As a consequence, we will use hereafter the following expression to calculate the cut-off frequency to ensure that the frequency range from to 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 . Since we used an ideal low-pass filter, the integral from to could simply be substituted by 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)
The value of the integral depends on and . In Figure 3, we show the frequency response of 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 . In Figure 5 we show the ratio
for values of . The values are all between and . 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, as a representative value and then add the factor 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 . 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)
and then performing filtering and backprojection on the rebinned projections:
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 are used, where 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 in subintervals of length . We denote the backprojection over these subintervals as partial block backprojections (PBBs) as in . 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 , the uncertainty due to rebinning does not depend on the number of PBBs per rotation ; there is a certain threshold , 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,  shows that with the typical fields of view of most current CT scanners and , 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 . 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 . According to the sampling rate of obtained when combining the time-series of the th PBB with the time series of the th PBB (computed with equivalent rays), accurate reconstruction is possible as long as 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 PBBs. Contrary to the TIA-TFDK algorithm, the temporal estimation is carried out using polynomial spline filters with a normalized cut-off frequency , that is, the PBB samples are smoothed along the time axis.
A detailed description of the steps of the algorithm is provided in Appendix B.
Since typically Hz (see ), the condition (14) translates into the rotation time having to fulfill 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 , for which we provided the following empirical formula in :
where is the radius of the object and is the radius of the circular trajectory of the source. This 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 , can be reduced to 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, projections per rotation were simulated with a cylindrical scanner of channels. A spherical phantom containing spherical inserts with time-dependent attenuation was used. In this phantom, the inserts follow the temporal law
with s, and s. This curve fulfills 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. () determines the amplitude and was chosen in such a way that the maximum values of the enhancement were and HU.
We first simulated acquisition and reconstruction with the reference protocol, consisting in full-scan T-FDK reconstructions . scans were simulated with a rotation time of s during 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 . In this case, TSA-TFDK reduced to the TIA-TFDK algorithm as described in . 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 . 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 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.
(a) Reference protocol
(a) Max. enh.: 18 HU, fast scanning
(b) Max. enh.: 18 HU, slow scanning
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 With the values from Table 2 we can calculate the measured reduction of the variance as and , 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 , 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 s, during which the source rotated at s/rot and was switched on every two rotations. Hence, images were generated which corresponds to a rate of image/s. The scan was performed with a tube voltage of kVp and at a tube current of mA. Since the rotation time was s, the dose per reconstructed image in mAs was 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 mm.
Since we had already a reconstructed sequence at our disposal, we used it to estimate a value of instead of using an a priori general estimation for it. Note that this was possible because the sampling interval was s which satisfies the sampling condition since Hz . For the purpose of estimating , 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 as the value for which contains % of the energy of the signal. We obtained a value of Hz.
for which (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
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 , , or 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 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 , the variance of the other sequences will be The variance of the noise of the original data is . 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
In Figure 9, we show some examples of time-attenuation curves of the original (gray) and the low-noise sequence 1 with (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.
(a) Arteria cerebri anterior
(b) Sinus sagittalis
(d) 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 .
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 pixels with the same temporal behavior. We then used the whole 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 s. Sequence 4 has samples at s and s so that the frame at 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.
In this paper we have presented the dynamic reconstruction algorithm TSA-TFDK, which is a further development of the previously published TIA-TFDK algorithm . 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 Hz from data acquired by scanning with any rotation time up to 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:
The temporal evolution is modeled as a continuous deterministic signal contaminated by additive noise
The noise is assumed to be a continuous random process with zero mean and mean square value .
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)  and thus is a wide-sense stationary process in time . 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 and 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 the discrete sequence of noise samples obtained by sampling with a sampling rate of . 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 as where denotes the expectation value and is the Kronecker symbol ( and if ). The discrete power spectral density is defined as (see )
where is the normalized frequency with as the physical frequency. The notation indicates that it is the Fourier Transform of a discrete sequence and is therefore periodic (we follow the notational conventions of ). owes its name to the fact that integrating it over one period of the normalized frequency, , yields the mean square value or variance of the noise. Downsampling 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 of the samples of a continuous process sampled every by introducing the change of variables in (A.4). Thus, we get
Integrating (A.5) over one period on the physical frequency axis we obtain the mean square value or variance of the discrete process
Therefore 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 by a factor of (or sampling with sampling interval ), is increased by a factor of with respect to the discrete power spectral density of (see Figure 12)
The TSA-TFDK is summarized in the following steps.
(1) Choose rotation time such that
(2) Calculate for
(3) Rebin cone-beam projections to fan-parallel beam
(4) Divide every rotation in sets of projections and reconstruct partial block backprojections with the T-FDK algorithm where for all for all .
(5) Combine the th and the th sequences of partial block backprojections to a unique sequence where for all , for all .
(6) Estimate the values of the partial reconstructions at the output frame times with polynomial spline filtering where for all .
(7) Accumulate the estimated partial block backprojections to obtain the time dependent attenuation map
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.
- K. Miles, P. Dawson, and M. Blomley, Functional Computed Tomography, Isis Medical Media, Oxford, UK, 1997.
- K. Sartor, Ed., Neuroradiologie, Thieme, Stuttgart, Germany, 2001.
- M. Koenig, M. Kraus, C. Theek, E. Klotz, W. Gehlen, and L. Heuser, “Quantitative assessment of the ischemic brain by means of perfusion-related parameters derived from perfusion CT,” Stroke, vol. 32, no. 2, pp. 431–437, 2001.
- M. Wintermark, M. Reichhart, O. Cuisenaire et al., “Comparison of admission perfusion computed tomography and qualitative diffusion- and perfusion-weighted magnetic resonance imaging in acute stroke patients,” Stroke, vol. 33, no. 8, pp. 2025–2031, 2002.
- P. Schramm, P. D. Schellinger, E. Klotz et al., “Comparison of perfusion computed tomography and computed tomography angiography source images with perfusion-weighted imaging and diffusion-weighted imaging in patients with acute stoke of less than 6 hours' duration,” Stroke, vol. 35, no. 7, pp. 1652–1658, 2004.
- P. Montes and G. Lauritsch, “A temporal interpolation approach for dynamic reconstruction in perfusion CT,” Medical Physics, vol. 34, no. 7, pp. 3077–3092, 2007.
- W. A. Kalender, Computed Tomography: Fundamentals, System Technology, Image Quality, Applications, Publicis MCD, Munich, Germany, 2000.
- K. L. Zierler, “Equations for measuring blood flow by external monitoring of radioisotopes,” Circulation Research, vol. 16, pp. 309–321, 1965.
- P. Montes, Dynamic cone-beam reconstruction for perfusion computed tomography, Ph.D. thesis, Interdisciplinary Center for Scientific Computing, Ruprecht-Karls-Universität Heidelberg, Heidelberg, Germany, 2006, http://www.ub.uni-heidelberg.de/archiv/7020/.
- M. Unser, A. Aldroubi, and M. Eden, “B-spline signal processing: part I—theory,” IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821–833, 1993.
- M. Unser, A. Aldroubi, and M. Eden, “B-spline signal processing: part II—efficient design and applications,” IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834–848, 1993.
- A. Aldroubi, M. Unser, and M. Eden, “Cardinal spline filters: stability and convergence to the ideal sinc interpolator,” Signal Processing, vol. 28, no. 2, pp. 127–138, 1992.
- M. Unser and T. Blu, “Self-similarity: part I—splines and operators,” IEEE Transactions on Signal Processing, vol. 55, no. 4, pp. 1352–1363, 2007.
- H. H. Barret and K. J. Myers, Foundations of Image Science, John Wiley & Sons, New York, NY, USA, 2004.
- M. Grass, Th. Köhler, and R. Proksa, “3D cone-beam CT reconstruction for circular trajectories,” Physics in Medicine and Biology, vol. 45, no. 2, pp. 329–347, 2000.
- P. Grangeat, A. Koenig, T. Rodet, and S. Bonnet, “Theoretical framework for a dynamic cone-beam reconstruction algorithm based on a dynamic particle model,” Physics in Medicine and Biology, vol. 47, no. 15, pp. 2611–2625, 2002.
- A. Papoulis, Probability, Random Variables and Stochastic Processes, Series in Systems Science, McGraw-Hill, New York, NY, USA, 1965.
- A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, Upper Saddle River, NJ, USA, 1998.
Copyright © 2009 Pau Montes and Günter Lauritsch. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.