Effect of errors in the system matrix on maximum a posteriori image reconstruction

Statistically based iterative image reconstruction methods have been developed for emission tomography. One important component in iterative image reconstruction is the system matrix, which defines the mapping from the image space to the data space. Several groups have demonstrated that an accurate system matrix can improve image quality in both single photon emission computed tomography (SPECT) and positron emission tomography (PET). While iterative methods are amenable to arbitrary and complicated system models, the true system response is never known exactly. In practice, one also has to sacrifice the accuracy of the system model because of limited computing and imaging resources. This paper analyses the effect of errors in the system matrix on iterative image reconstruction methods that are based on the maximum a posteriori principle. We derived an analytical expression for calculating artefacts in a reconstructed image that are caused by errors in the system matrix using the first-order Taylor series approximation. The theoretical expression is used to determine the required minimum accuracy of the system matrix in emission tomography. Computer simulations show that the theoretical results work reasonably well in low-noise situations.     

Tomographic motion detection and correction directly in sinogram space

Patient motion, especially respiratory motion, results in various artefacts such as blurring and streaks in tomographic images. The interplay of the movement of the beam aperture and variations of organ anatomy during delivery can create 'hot' and 'cold' spots throughout the field in intensity-modulated radiation therapy (IMRT). Detection and correction of patient motion is extremely important in tomographic imaging and IMRT. Tomographic projection data (sinogram) encode not only the patient anatomy information, but also the intra-scanning motion information. In this paper, we developed an algorithm to detect and correct the in-plane respiratory motion directly in sinogram space. The respiratory motion is modelled as time-varying scaling along the x and y directions. Its effects on the sinogram are discussed. Based on the traces of some nodal points in the sinogram, the intra-scanning motion is determined. The motion correction is also implemented in sinogram space. The motion-corrected sinogram is used for reconstruction by the filtered back-projection (FBP) method. Computer simulations validate the motion detection and correction algorithm. The reconstructed images from the motion-corrected sinogram eliminate the majority of the artefacts. The method could be applied to projection data used in CT and ECT, as well as in tomotherapy delivery modification and dose reconstruction.     

基于投影矩阵广义逆的CT重建迭代算法

目的在CT检查时,有限角度投影和稀疏矩阵投影能够减少X射线的剂量,然而这会导致投影数据不足,给图像重建带来一定的困难。为了克服这一难题得到较好的重建图像,本文提出一种基于计算投影矩阵广义逆的CT迭代重建算法。方法该算法在计算过程中,将重建图像表示为投影矩阵以及其广义逆的乘积。首先使用一阶迭代计算广义逆矩阵,但是由于投影矩阵和其广义逆矩阵都比较大,在迭代过程中以投影和滤波反投影代替。然后通过不同的算法分别对平行束投影、有限角度投影、稀疏矩阵投影的数据进行重建,并对重建结果的均方差、通用图像质量指标以及图像互信息进行比较。结果本文提出的方法重建出图像的均方差、通用图像质量指标和图像互信息更优,而且重建时间较短。结论该方法能够在没有未知图像先验结构信息和伪影猜想的情况下有效地提高重建图像的质量,而且该算法不需要计算投影过程,重建过程简单易行。     

迭代成像中一种快速求解投影矩阵的方法

迭代成像中投影矩阵的获取速度直接影响算法的效率.本文提出了一种快速求解投影矩阵的方法,它尤其适用于扇形扫描投影数据的重建.利用扇形射束与重建图像像素的几何对称性一次获取多个投影矩阵的向量以及实时获取投影矩阵向量方式,减少了获得投影矩阵花费时间,节省了存储矩阵所占用的空间.通过模拟投影数据进行图像重建,其结果表明,该法较全部求出投影矩阵法极大地减少了内存空间的占用,较数据转换法节省了转换时间,从而提高了扇形投影数据的图像重建效率.     

Three-dimensional image reconstruction for PET by multi-slice rebinning and axial image filtering

A fast method is described for reconstructing volume images from three-dimensional (3D) coincidence data in positron emission tomography (PET). The reconstruction method makes use of all coincidence data acquired by high-sensitivity PET systems that do not have inter-slice absorbers (septa) to restrict the axial acceptance angle. The reconstruction method requires only a small amount of storage and computation, making it well suited for dynamic and whole-body studies. The method consists of three steps: (i) rebinning of coincidence data into a stack of 2D sinograms; (ii) slice-by-slice reconstruction of the sinogram associated with each slice to produce a preliminary 3D image having strong blurring in the axial (z) direction, but with different blurring at different z positions; and (iii) spatially variant filtering of the 3D image in the axial direction (i.e. 1D filtering in z for each x-y column) to produce the final image. The first step involves a new form of the rebinning operation in which multiple sinograms are incremented for each oblique coincidence line (multi-slice rebinning). The axial filtering step is formulated and implemented using the singular value decomposition (SVD). The method has been applied successfully to simulated data and to measured data for different kinds of phantom (multiple point sources, multiple discs, a cylinder with cold spheres, and a 3D brain phantom).     

Scatter correction for positron emission mammography

In this paper we present a scatter correction method for a regularized list mode maximum likelihood reconstruction algorithm for the positron emission mammograph (PEM) that is being developed at our laboratory. The scatter events inside the object are modelled as additive Poisson random variables in the forward model of the reconstruction algorithm. The mean scatter sinogram is estimated using a Monte Carlo simulation program. With the assumption that the background activity is nearly uniform, the Monte Carlo scatter simulation only needs to run once for each PEM configuration. This saves computation time. The crystal scatters are modelled as a shift-invariant blurring in image domain because they are more localized. Thus, the useful information in the crystal scatters can be deconvolved in high-resolution reconstructions. The propagation of the noise from the estimated scatter sinogram into the reconstruction is analysed theoretically. The results provide an easy way to calculate the required number of events in the Monte Carlo scatter simulation for a given noise level in the image. The analysis is also applicable to other scatter estimation methods, provided that the covariance of the estimated scatter sinogram is available.     

Sparse angular CT reconstruction using non-local means based iterative-correction POCS

In divergent-beam computed tomography (CT), sparse angular sampling frequently leads to conspicuous streak artifacts. In this paper, we propose a novel non-local means (NL-means) based iterative-correction projection onto convex sets (POCS) algorithm, named as NLMIC-POCS, for effective and robust sparse angular CT reconstruction. The motivation for using NLMIC-POCS is that NL-means filtered image can produce an acceptable priori solution for sequential POCS iterative reconstruction. The NLMIC-POCS algorithm has been tested on simulated and real phantom data. The experimental results show that the presented NLMIC-POCS algorithm can significantly improve the image quality of the sparse angular CT reconstruction in suppressing streak artifacts and preserving the edges of the image.     

Single scan parameterization of space-variant point spread functions in image space via a printed array: the impact for two PET/CT scanners

Incorporation of a resolution model during statistical image reconstruction often produces images of improved resolution and signal-to-noise ratio. A novel and practical methodology to rapidly and accurately determine the overall emission and detection blurring component of the system matrix using a printed point source array within a custom-made Perspex phantom is presented. The array was scanned at different positions and orientations within the field of view (FOV) to examine the feasibility of extrapolating the measured point source blurring to other locations in the FOV and the robustness of measurements from a single point source array scan. We measured the spatially-variant image-based blurring on two PET/CT scanners, the B-Hi-Rez and the TruePoint TrueV. These measured spatially-variant kernels and the spatially-invariant kernel at the FOV centre were then incorporated within an ordinary Poisson ordered subset expectation maximization (OP-OSEM) algorithm and compared to the manufacturer's implementation using projection space resolution modelling (RM). Comparisons were based on a point source array, the NEMA IEC image quality phantom, the Cologne resolution phantom and two clinical studies (carbon-11 labelled anti-sense oligonucleotide [(11)C]-ASO and fluorine-18 labelled fluoro-l-thymidine [(18)F]-FLT). Robust and accurate measurements of spatially-variant image blurring were successfully obtained from a single scan. Spatially-variant resolution modelling resulted in notable resolution improvements away from the centre of the FOV. Comparison between spatially-variant image-space methods and the projection-space approach (the first such report, using a range of studies) demonstrated very similar performance with our image-based implementation producing slightly better contrast recovery (CR) for the same level of image roughness (IR). These results demonstrate that image-based resolution modelling within reconstruction is a valid alternative to projection-based modelling, and that, when using the proposed practical methodology, the necessary resolution measurements can be obtained from a single scan. This approach avoids the relatively time-consuming and involved procedures previously proposed in the literature.     

Differential attenuation method for simultaneous estimation of activity and attenuation in multiemission single photon emission computed tomography

A penalized weighted least squares reconstruction algorithm is described that simultaneously estimates activity and attenuation distributions from emission sinogram data alone. This estimation technique is based on differential attenuation information and is applicable to any single photon emission computed tomography imaging isotope with emissions at two or more distinct energies, after accurate compensation for Compton scatter. A rotation-based forward projector is used to efficiently model photon attenuation at multiple emission energies, as well as distance-dependent spatial resolution. The algorithm was tested using simulated scatter-free 201T1 projection data from a single-slice numerical cardiac phantom with and without cold myocardial defects. Poisson noise was added to the projection data to mimic clinically realistic count densities. The activity estimates resulting from the proposed method had fewer artifacts and were substantially more accurate than images reconstructed with filtered backprojection without compensation for attenuation. Several techniques were employed to reduce the time required for the iterative routine to converge and to reduce the sensitivity of the solution to noise in the projection data. These included: (1) a preconditioning image variable transformation; (2) a coarse-to-fine grid initialization schedule; and (3) a convex hull image mask determined directly from the data. The combined effect of these techniques substantially reduced the compute time required for the reconstruction.     

GPU-based iterative cone-beam CT reconstruction using tight frame regularization

The x-ray imaging dose from serial cone-beam computed tomography (CBCT) scans raises a clinical concern in most image-guided radiation therapy procedures. It is the goal of this paper to develop a fast graphic processing unit (GPU)-based algorithm to reconstruct high-quality CBCT images from undersampled and noisy projection data so as to lower the imaging dose. For this purpose, we have developed an iterative tight-frame (TF)-based CBCT reconstruction algorithm. A condition that a real CBCT image has a sparse representation under a TF basis is imposed in the iteration process as regularization to the solution. To speed up the computation, a multi-grid method is employed. Our GPU implementation has achieved high computational efficiency and a CBCT image of resolution 512 × 512 × 70 can be reconstructed in ～5 min. We have tested our algorithm on a digital NCAT phantom and a physical Catphan phantom. It is found that our TF-based algorithm is able to reconstruct CBCT in the context of undersampling and low mAs levels. We have also quantitatively analyzed the reconstructed CBCT image quality in terms of the modulation-transfer function and contrast-to-noise ratio under various scanning conditions. The results confirm the high CBCT image quality obtained from our TF algorithm. Moreover, our algorithm has also been validated in a real clinical context using a head-and-neck patient case. Comparisons of the developed TF algorithm and the current state-of-the-art TV algorithm have also been made in various cases studied in terms of reconstructed image quality and computation efficiency.     

GPU-accelerated 3D Bayesian image reconstruction from Compton scattered data

This paper describes the development of fast Bayesian reconstruction methods for Compton cameras using commodity graphics hardware. For fast iterative reconstruction, not only is it important to increase the convergence rate, but also it is equally important to accelerate the computation of time-consuming and repeated operations, such as projection and backprojection. Since the size of the system matrix for a typical Compton camera is intractably large, it is impractical to use a conventional caching scheme that stores the pre-calculated elements of a system matrix and uses them for the calculation of projection and backprojection. In this paper we propose GPU (graphics processing unit)-accelerated methods that can rapidly perform conical projection and backprojection on the fly. Since the conventional ray-based backprojection method is inefficient for parallel computing on GPUs, we develop voxel-based conical backprojection methods using two different approximation schemes. In the first scheme, we approximate the intersecting chord length of the ray passing through a voxel by the perpendicular distance from the center to the ray. In the second scheme, each voxel is regarded as a dimensionless point rather than a cube so that the backprojection can be performed without the need for calculating intersecting chord lengths or their approximations. Our simulation studies show that the GPU-based method dramatically improves the computational speed with only minor loss of accuracy in reconstruction. With the development of high-resolution detectors, the difference in the reconstruction accuracy between the GPU-based method and the CPU-based method will eventually be negligible.     

Pragmatic fully 3D image reconstruction for the MiCES mouse imaging PET scanner

We present a pragmatic approach to image reconstruction for data from the micro crystal elements system (MiCES) fully 3D mouse imaging positron emission tomography (PET) scanner under construction at the University of Washington. Our approach is modelled on fully 3D image reconstruction used in clinical PET scanners, which is based on Fourier rebinning (FORE) followed by 2D iterative image reconstruction using ordered-subsets expectation-maximization (OSEM). The use of iterative methods allows modelling of physical effects (e.g., statistical noise, detector blurring, attenuation, etc), while FORE accelerates the reconstruction process by reducing the fully 3D data to a stacked set of independent 2D sinograms. Previous investigations have indicated that non-stationary detector point-spread response effects, which are typically ignored for clinical imaging, significantly impact image quality for the MiCES scanner geometry. To model the effect of non-stationary detector blurring (DB) in the FORE+OSEM(DB) algorithm, we have added a factorized system matrix to the ASPIRE reconstruction library. Initial results indicate that the proposed approach produces an improvement in resolution without an undue increase in noise and without a significant increase in the computational burden. The impact on task performance, however, remains to be evaluated.     

Progressive cross-section display of 3D medical images

The paper presents a hierarchical coding algorithm for 3D medical images based upon hierarchical interpolation with radial basis function networks. By using the properties of the Kronecker product, the computation of the network parameters and the 3D image reconstruction are efficiently done in (L⁴) computation time and O(L³) storage space, when applied to 3D images of size (L×L×L). A further reduction in processing time is accomplished by using sparse matrix techniques. The salient features of the proposed coding method are that arbitrary cross-section images can be progressively displayed without reconstruction of the whole 3D image; the first image reconstruction starts as soon as the first data transmission has been completed; no expanding procedure is required in 3D image reconstruction, and the blocking effects are not apparent even in the lowest-resolution image. Experimental results using two 3D MRI images, of size (128×18×64) and with 8-bit grey levels, show that the coding performance is better than that of the 3D DCT coding by about 0.25 bits pixel⁻¹ at higher bit rates, and that the new cross-section display method synthesises the coarsest (finest) section image about six (three) times faster than the standard method that requires the whole 3D image reconstruction.     

A practical method for three-dimensional reconstruction of joints using a C-arm system and shift-and-add algorithm

Currently, radiography with C-arm systems is playing a major role in the assessment of arthritis. However, the radiographic two-dimensional projection images of joints often interfere with physicians' efforts to better understand and measure the structure changes of joints due to the overlap of bone structures at different depths. An accurate, low-cost, and practical three-dimensional (3D) reconstruction approach of joints will be beneficial in diagnosing arthritis. Toward this end, a novel method is developed in this paper based on a C-arm system. The idea is to apply the shift-and-add algorithm (commonly used in digital tomosynthesis) on the segmented projection images at multiple angles, which results in accurate reconstruction of the 3D structures of joints. The method provides a new solution to precisely distinguish objects from blurring background. The proposed method has been tested and evaluated on simulated cylinders, a chicken bone phantom with known structure, and an in vivo human index finger. The results are demonstrated and discussed.     

CT image registration in sinogram space

Object displacement in a CT scan is generally reflected in CT projection data or sinogram. In this work, the direct relationship between object motion and the change of CT projection data (sinogram) is investigated and this knowledge is applied to create a novel algorithm for sinogram registration. Calculated and experimental results demonstrate that the registration technique works well for registering rigid 2D or 3D motion in parallel and fan beam samplings. Problem and solution for 3D sinogram-based registration of metallic fiducials are also addressed. Since the motion is registered before image reconstruction, the presented algorithm is particularly useful when registering images with metal or truncation artifacts. In addition, this algorithm is valuable for dealing with situations where only limited projection data are available, making it appealing for various applications in image guided radiation therapy.     

Transaxial system models for jPET-D4 image reconstruction

A high-performance brain PET scanner, jPET-D4, which provides four-layer depth-of-interaction (DOI) information, is being developed to achieve not only high spatial resolution, but also high scanner sensitivity. One technical issue to be dealt with is the data dimensions which increase in proportion to the square of the number of DOI layers. It is, therefore, difficult to apply algebraic or statistical image reconstruction methods directly to DOI-PET, though they improve image quality through accurate system modelling. The process that requires the most computational time and storage space is the calculation of the huge number of system matrix elements. The DOI compression (DOIC) method, which we have previously proposed, reduces data dimensions by a factor of 1/5. In this paper, we propose a transaxial imaging system model optimized for jPET-D4 with the DOIC method. The proposed model assumes that detector response functions (DRFs) are uniform along line-of-responses (LORs). Then each element of the system matrix is calculated as the summed intersection lengths between a pixel and sub-LORs weighted by a value from the DRF look-up-table. 2D numerical simulation results showed that the proposed model cut the calculation time by a factor of several hundred while keeping image quality, compared with the accurate system model. A 3D image reconstruction with the on-the-fly calculation of the system matrix is within the practical limitations by incorporating the proposed model and the DOIC method with one-pass accelerated iterative methods.     

Iterative image reconstruction in helical cone-beam x-ray CT using a stored system matrix approach

We present a stored system matrix (SM) approach for iterative x-ray CT image reconstruction with helical cone-beam geometry. Because of the symmetry of a helical source trajectory, it is sufficient to calculate and store the SM entries for one transaxial slice only and for all source positions illuminating the slice. This is made possible by (1) selecting the reconstruction slice thickness to be an integer multiple of the source translation per projection view, and (2) discretizing the 3D reconstruction volume on a rotated stack of slices. Using the proposed method, the memory requirement for reconstructing a full field-of-view of clinical scanners is manageable on current computing platforms. The same storage principle can be generalized and applied to volume-of-interest (VOI) image reconstruction for helical cone-beam CT. In this case, the stored SM entries correspond to a partial- or full-ring region on one transaxial slice, and for all source positions illuminating the ring. The size and location of the ring depend on the size and the location of the VOI and the scan geometry. We demonstrate by both computer simulations and clinical patient data the speed and efficacy of iterative image reconstruction using the stored SM approach.     

Reduction of motion blurring artifacts using respiratory gated CT in sinogram space: a quantitative evaluation

Techniques have been developed for reducing motion blurring artifacts by using respiratory gated computed tomography (CT) in sinogram space and quantitatively evaluating the artifact reduction. A synthetic sinogram was built from multiple scans intercepting a respiratory gating window. A gated CT image was then reconstructed using the filtered back-projection algorithm. Wedge phantoms, developed for quantifying the motion artifact reduction, were scanned while being moved using a computer-controlled linear stage. The resulting artifacts appeared between the high and low density regions as an apparent feature with a Hounsfield value that was the average of the two regions. A CT profile through these regions was fit using two error functions, each modeling the partial-volume averaging characteristics for the unmoving phantom. The motion artifact was quantified by determining the apparent distance between the two functions. The blurring artifact had a linear relationship with both the speed and the tangent of the wedge angles. When gating was employed, the blurring artifact was reduced systematically at the air-phantom interface. The gated image of phantoms moving at 20 mm/s showed similar blurring artifacts as the nongated image of phantoms moving at 10 mm/s. Nine patients were also scanned using the synchronized respiratory motion technique. Image artifacts were evaluated in the diaphragm, where high contrast interfaces intercepted the imaging plane. For patients, this respiratory gating technique reduced the blurring artifacts by 9%-41% at the lung-diaphragm interface.     

An optimization-based method for geometrical calibration in cone-beam CT without dedicated phantoms

In this paper we present a new method for the determination of geometrical misalignments in cone-beam CT scanners, from the analysis of the projection data of a generic object. No a priori knowledge of the object shape and positioning is required. We show that a cost function, which depends on the misalignment parameters, can be defined using the projection data and that such a cost function has a local minimum in correspondence to the actual parameters of the system. Hence, the calibration of the scanner can be carried out by minimizing the cost function using standard optimization techniques. The method is developed for a particular class of 3D object functions, for which the redundancy of the fan beam sinogram in the transaxial midplane can be extended to cone-beam projection data, even at wide cone angles. The method has an approximated validity for objects which do not belong to that class; in that case, a suitable subset of the projection data can be selected in order to compute the cost function. We show by numerical simulations that our method is capable to determine with high accuracy the most critical misalignment parameters of the scanner, i.e., the transversal shift and the skew of the detector. Additionally, the detector slant can be determined. Other parameters such as the detector tilt, the longitudinal shift and the error in the source-detector distance cannot be determined with our method, as the proposed cost function has a very weak dependence on them. However, due to the negligible influence of these latter parameters in the reconstructed image quality, they can be kept fixed at estimated values in both calibration and reconstruction processes without compromising the final result. A trade-off between computational cost and calibration accuracy must be considered when choosing the data subset used for the computation of the cost function. Results on real data of a mouse femur as obtained with a small animal micro-CT are shown as well, proving the capability of the proposed calibration method. In principle, the method can be adapted to other cone-beam imaging modalities (e.g., single photon emission computed tomography).     

Statistical inversion for medical x-ray tomography with few radiographs: II. Application to dental radiology

Diagnostic and operational tasks in dental radiology often require three-dimensional information that is difficult or impossible to see in a projection image. A CT-scan provides the dentist with comprehensive three-dimensional data. However, often CT-scan is impractical and, instead, only a few projection radiographs with sparsely distributed projection directions are available. Statistical (Bayesian) inversion is well-suited approach for reconstruction from such incomplete data. In statistical inversion, a priori information is used to compensate for the incomplete information of the data. The inverse problem is recast in the form of statistical inference from the posterior probability distribution that is based on statistical models of the projection data and the a priori information of the tissue. In this paper, a statistical model for three-dimensional imaging of dentomaxillofacial structures is proposed. Optimization and MCMC algorithms are implemented for the computation of posterior statistics. Results are given with in vitro projection data that were taken with a commercial intraoral x-ray sensor. Examples include limited-angle tomography and full-angle tomography with sparse projection data. Reconstructions with traditional tomographic reconstruction methods are given as reference for the assessment of the estimates that are based on the statistical model.