Review of Parallel Computing Techniques for Computed Tomography Abstract Image Reconstructi

Review of Parallel Computing Techniques for Computed TomographyImage ReconstructionJun Ni1, 3, Xiang Li2, Tao He3, Ge Wang1,21Medical Imaging High Performance Computing Lab, Department of Radiology,2Department of Biomedical Engineering, 3Department of Computer ScienceThe University of Iowa, Iowa City, IA 52242{jun-ni, xiang-li, tao-he, ge-wang}AbstractAfter we briefly review representative analytic and iterative reconstruction algorithms for X-ray computed tomography (CT), we address the need for faster reconstruction by parallel computing techniques. For a decent, a cone-beam reconstruction usually takes hours on a regular PC, since most of algorithms take more than 60 iterations even longer. In order to speedup the performance, people introduce various acceleration methodologies including algorithm improvements, chip utilization, and parallel computing technique. This paper focuses on the speedup the computation using parallel computing. The first generation of parallel computing systems was based on a centralized parallel configuration. The second generation of such systems employed a cluster of general-purpose computers that are connected by a fast local area network (LAN). Hereby, we highlight distributed parallel computing techniques: from a locally distributed client-server topology to a peer-to-peer (P2P) enhanced network model. With the P2P technology, the client would be directly connected to all other computing peers seamlessly, forming a virtual parallel computer. There are multiple Internet connections between the client and other computing peers. This way, a single failure of node wouldn’t cause the entire failure of computation. Finally, we state that by integrating the large-scale geographically distributed systems such as Grid computing technology the future of the CT reconstruction will be highly parallel, efficient, scalable over the Internet, so will be other biomedical imaging tasks.1. IntroductionX-ray computed tomography (CT) is one of the most important non-invasive medical imaging techniques [1]. X-ray CT reconstructs a cross-sectional image by computing the X-ray absorption coefficient distribution of an object from projection data, which records the relative number of photons passing through the object. Hence, X-ray CT is regarded as transmission CT. Another imaging modality is emission CT such aspositron emission tomography (PET) [2-4] and single photon emission computed tomography (SPECT) [5-7] where the distribution of injected radioactive chemicals is estimated. No matter it is emission or transmission CT, the principles for image reconstruction remains the same.The development of X-ray CT technology is closely related to the evolution of the detector design and the scanning mode. The first generation of X-ray CT scanners used the parallel-beam geometry. (Fig. 1) The next generation systems were in the fan-beam geometry, which may be further divided into sub-categories [8, 9]. For 3D image reconstruction, an image volume was traditionally reconstructed by stacking 2D cross-sectional images. This method resulted in poor resolution in the axial direction. The modern scanning mode is to let the gantry rotate continuously while a patient table is simultaneously translated [10-12]. From the patient’s point of view, the X-ray source moves along a spiral or helical locus. The spiral scan enables continuous data acquisition and improves the image quality significantly. Almost all the modern CT devices allow spiral scanning.(a) (b)Fig. 1. Parallel-beam (a) and fan-beam (b) scan geometries.The X-ray photons emitted from the radiation source naturally form a cone traveling away from the source focal spot. A collimator is used in the parallel-beam and fan-beam scanners to restrict the X-ray beam to one single line or a set of lines on a plane, respectively. The first multi-slice spiral CT (MSCT) came into the market in 1998 employing a four-row detector. MSCT was a breakthrough in the CT technology that made possible the sub-minute whole body CT scan a routine clinical exam [13]. To reduce the scan time and improve the X-ray energy efficiency even further a larger area detector is desirable. In MSCT, the cone angle is small, usually of several degrees. Algorithms for fan-beam reconstruction can still be adapted for MSCT. When a cone angle is as large as tens of degrees extended by an area detector, a new challenge must be met in the algorithm development. Cone-beam CT reconstruction has been an active research area for the past decade. Many algorithms have been proposed [14-19]. They can be generally grouped into either analytic or iterative algorithms. For a full review of these algorithms, please see [20]. Recently, Katsevich developed an efficient exact cone-beamreconstruction algorithm [21, 22]. However the truly cone-beam medical X-ray CT scanner has not been popular in the market yet.1.1. Analytic reconstructionThe filtered back-projection (FBP) method is predominantly used with most commercial X-ray CT or PET CT scanners. With FBP the projection data are first filtered then the filtered data is linearly smeared back along ray paths to form image pixels. For example, in the parallel-beam geometry the relationship between the projection data and the object can be described as: [23]∫=),(),(),(t l dl y x f t P θθ , (1.1)where ),(t P θ is the projection data measured at a projection angle θ, and t the detector position in the beam.Using the δ-function to define the line integral, we have∫∫+∞∞−+∞∞−−+=dxdy t y x y x f t P )sin cos (),(),(θθδθ. (1.2)The CT image reconstruction problem is to compute given ),(y x f ),(t P θ. The filtered backprojection (FBP) method can be formulated as:∫=πθθ20),(),(d t D y x f p , θθsin cos y x t +=, (1.3)⎟⎠⎞⎜⎝⎛=∫+∞∞−dr e r t P t D rt j p πθθ221||*),(),( (1.4){}{||),(1r t P F F θ−=},(1.5) where and denote the forward and inverse Fourier transforms, respectively. Equation (1.4) is the filtering step in the form of convolution, which can be carried out in the frequency domain using equation (1.5). The FBP formula for fan-beam CT can be similarly obtained.{}F {}1−FIn the cone-beam geometry, a 2D area detector is used. In this case, there are two types of analytic reconstruction methods: exact and approximate algorithms. The exactalgorithm is able to reconstruct the original object using cone-beam data accurately given sufficiently fine detector resolution and large number of projections. The approximate algorithm often has simpler formulation and faster speed. Not all scanning trajectories permit exact reconstruction. As stated by the Smith-Tuy completeness condition: “If on every plane that intersects the object there lies a vertex (the source point), then one has complete information about the object ” [24]. The completeness condition can be heuristically interpreted as that a cone-beam trajectory must fill out data in 3D Radon space completely. In other words, on any plane that intersects the object being scanned there exists at least one radiation source. For example, a spiral scan curve is complete (even data are longitudinally truncated) for exact cone-beam reconstruction but a circular source trajectory is not.Among many approximate algorithms the Feldkamp-type algorithms are the most popular. The original Feldkamp formula extends the circular fan-beam reconstruction to the circular cone-beam geometry by compensating for the cone angle effect appropriately. Suppose a planar detector is used. The detector position is determined by the Cartesian coordinate. The Feldkamp cone-beam reconstruction formula can be written as:[20]),(v u ∫′′=πβββ202),,(),,(),,(d y x L v u D z y x f cone , (1.12)where the interpolated detector position is determined byβββββsin cos sin cos ),,(y x d x y d y x u so so ++−=′, (1.13) βββsin cos ),,,(y x d d z z y x v so so ++=′,(1.14) andso so d y x d y x L βββsin cos ),,(++=. (1.15) The filtering step is expressed as:⎟⎠⎞⎜⎝⎛⎟⎟⎠⎞⎜⎜⎝⎛++=∫∞+∞−⋅dr e r v u d d v u P v u D u r j soso cone cone πββ221222||*),,(),,(, (1.16) where ),,(v u P cone β are cone-beam projections.1.2. Iterative reconstructionThe iterative reconstruction (IR) methods include statistical reconstruction (SR) algorithms and algebraic reconstruction techniques (ART), but they all compute the final image iteratively through the same top-level loop (Fig. 2). There are many IR algorithms available. Representative algorithms are the maximum likelihood (ML) expectation maximization (EM) formula, [25-27] the simultaneous algebraic reconstruction technique (SART), [28-30] and the Convex algorithm [31, 32]. With ML-EM the image is obtained iteratively as an optimal estimate that maximizes the likelihood of the detection of the actual measured photons based on a statistical model of the imaging system. The EM method can also be deterministically interpreted as the process of minimizing the I-divergence between the estimated and measured projection data in the nonnegative space [33]. The SART algorithm iteratively minimizes the mean square error between the estimated and measured projections in the real space. The Convex algorithm is a statistical reconstruction algorithm for transmission CT, which also aims at maximizing the Poisson likelihood. The IR algorithms are superior to the analytic methods in terms of image quality (contrast and resolution) with noisy and/or incomplete projections.Fig. 2. Iterative reconstruction process.2. Why a single PC insufficientThe combination of the cone-beam geometry and the spiral scan made possible the data acquisition for 3D reconstruction ever faster and easier. However, along with the ease of obtaining large 3D projection data sets is the intense computation needed to reconstruct the volume. A micro-CT scanner with a 4000 by 2000 detector cells will need, depending on resolution settings, more than three hours to reconstruct a decent image volume using the FBP method. The time for whole body CT reconstruction on a medical scanner may take long time too.The situation is even more challenging when the IR technique is used because themajor disadvantage of the IR is its high demand on computation. Besides updating an intermediate image, a single iteration of the IR needs one forward-projection and one backprojection whereas the FBP method only needs one backprojection. The forward-projection step is of the same computational complexity as the backprojection step. As an example, the EM method often requires at least 30 iterations. Therefore, it takes approximately 60 times longer time for a comparable quality reconstruction than the FBP method.Image volumeProjection data Projection data(a) (b) (c)Fig. 3. Memory allocations needed from a scan (a) for FBP (b) and IR (c).Another problem with IR is that it is associated large memory need. For example, to reconstruct a volume of 5123 voxels it may use data from 384 projection views with 256by 512 detectors. The memory buffers for one volume and the projection data will be512MB and 192MB, respectively (Assuming that each data point is stored as 4-byte floating type). Without other intermediate parameters caching we need another two memory allocations of the size of image volume for error update of every voxel and the associated weighting factors, respectively. (Fig. 3) In total at least 1,728 MB (512x3 + 192) or about 1.7GB free memory are needed for IR. It has already exceeded the capacityof most current 32-bit PCs. As a result, even if the IR is highly desirable for 3D image reconstruction, it is impossible to implement such an IR algorithm on a single PC.3. How to accelerate?To speed up the reconstruction, generally three approaches are viable. (Fig. 4) These approaches are not exclusive to each other and actually can be combined. The first approach is to improve the algorithm itself both for a faster convergence rate and a more efficient backprojection implementation. In the FBP case, approximate algorithms are often preferred over exact algorithms for sake of efficiency. In the IR case, the ordered subset (OS) scheme is the most popular trick, which is able to accelerate the IR by anorder of magnitude. [34-36]Katsevich 2002 [21] Fast exact FBP Wang. 1993 [16] Feldkamp-type. Various scan loci Hudson&Larkin 1994 [34] Ordered SubsetsImproved AlgorithmsParallel computing Mainframe parallel computers Guerrini 1989 [42] Vector computer Chen 1990 [46] Hypercube McCarty 1991 [44] Mesh-parallel Atkins 1991 [45] Transputers Hardware Acceleration Proprietary chips on scanners Commodity graphics cards Cabral 1994 [38] Feldkamp-type Shattuck 2002 [48] Fast Internet-connected cluster for EM 3D PETLi 2004 [56] OS-SART on cluster for X-ray CTDistributed computing Li&Ni, 2004 [57, 58] P2P-enhanced network Fig. 4. Development of acceleration methods for CT image reconstruction.The second approach uses the graphics processing unit (GPU) to conduct the expensive processing steps such as forward-projection and backprojection. The modern GPU is able to handle 32-bit floating-point data and may have as many as four channels for color and eight channels for texture rendering that work in parallel [37]. It incorporates common arithmetic operations on the chip. Unlike those specially designed chips that are installed on most commercial CT scanners, the GPU also features a fully programmable interface, and some of them enable high-level programming languages such as C . The commodity texture mapping hardware was reported to support Feldkamp, SART, EM and OS-EM algorithms [38-40].The third approach is for parallelization of the computation. The most time-consuming part in the reconstruction process is backprojection in FBP, and bothforward-projection and backprojection in IR. Various parallelization schemes have been proposed to distribute the workload of backprojection and forward-projection among parallel computing units. They can be categorized into two types according to the hardware used.The first type of parallel computing systems utilizes a centralized parallel system such as a VLSI architecture computer, a vector computer, a large-scale parallel computer, and a shared memory multi-processor computer that is based on the same instruction multiple data (SIMD) structure [41-47]. Early efforts on parallel CT image reconstruction were often based on these mainframe computers. They were expensive, less scaleable and inefficient in terms of hardware usage. The architecture of these parallel machines was usually fixed upon manufacture. Users have little control over the design. This results in incompatible parallel implementations unique to the specific parallel computer. Algorithms would have to be re-developed should the hardware change.Fig. 5. Data partition scheme in the projection data domain. Bluearrow lines and straps illustrate the portion of the projection dataand the corresponding voxels being updated on individual workernodes, respectively.The second type of parallel computing systems employs general-purpose computers that are connected by a fast local area network (LAN). They are built based on the multiple instructions multiple data (MIMD) architecture. As the computing technology advances dramatically most recent parallel implementations become this type [48-52]. The PC cluster comprises several workstation nodes. One of them is the master node and the rest are worker nodes. The master node is usually a control and coordination unit. It distributes data to worker nodes and receives results from them for integration. (Fig. 5) PC clusters have been built on either WinNT or Unix/Linux platforms and could be a hybrid of SIMD computers and MIMD systems. In fact, dual processor nodes are not uncommon for many clusters. Some even installed four processors on each node [50][52].These systems normally call a set of cross-platform message passing libraries as the communication interface among individual nodes. Earlier implementations used the parallel virtual machine (PVM) protocol [53]. Latest ones often take advantages of the message passing interface (MPI) protocol [54]. The multi-threading and the message passing can be combined to fully exploit the computation power of the multi-processor cluster system. In [52], the inter-node communication is coordinated by MPI, while the inner-node data partition is controlled by the OpenMP multi-threading protocol [55]. (Fig.6) Most of these studied have targeted EM and OS-EM algorithms for PET and SPECT. In [56] parallel OS-SART was implemented on a Linux cluster for X-ray CT.Worker nodesFig. 6. Hybrid architecture of SIMD and MIMD. Each node has more than oneprocessor. They achieve multithreading via the OpenMP interface. The inter-nodeinterface is MPI. Note that there can be MPI links among worker nodes. Here onlylinks between master and worker nodes are shown.The parallel performance using a small LINUX cluster is illustrated in Fig. 7 and Fig. 8. The maximum speedup is achieved when number of processors is about 9.Fig. 7. Total iteration time and total computation time vs.number of computational node (servers).Fig. 8. Performance in terms of speedup vs. the number ofcomputing nodes used.The data utilized for the case study is 3-D image of Shepp-Logan phantom is shown in Fig. 9 (a) and (b).(a)(b)Fig. 9. Reconstructed image results: (a) selected sections of reconstructed images; (b) Thereconstructed profile image of 3D Shepp-Logan phantom (The blue blur represents the reconstructedimage cube).4. Distributed parallel reconstruction with InternetDecentralized parallel computing has many desirable features. Among these areflexibility, reliability, and cost-effectiveness. Recent studies favor drastically distributed parallel reconstruction with the modern Internet technologies. For example, in [48], [50], and [52] a java-applet enabled web-interface has been built to submit projection data and start the reconstruction. The remote cluster reconstructs the image and sends it back for analysis. Upon receiving the data the master node may either dispatch them directly to worker nodes or maintain a job queue waiting for regular checks from worker nodes. From the point of view of the distributed network model these tasks were based on the client-server (C/S) topology. The trend is that the client side becomes thinner and simpler leaving only the data submission and job request functionalities. The server completes heavy-duty tasks such as searching a database, calculation, information integration, and image reconstruction in the CT case. If a PC cluster is used instead of a mainframe computer, the master node of the cluster is also the server and is connected with the clients. Other worker nodes only connect to the master node via a high-speed local area network (LAN) (Fig. 9 [57]).(a)(b)(c) Fig. 10. Centralization of C/S model. (a) Clients are connected to a mainframe supercomputer; (b) Clients are connected to a remote gateway first through Internet; (c) Clients are connected to a remote cluster server.5. Peer-to-peer Distributed NetworkThe peer-to-peer (P2P) technology is able to further improve the robustness and performance of a parallel reconstruction system by facilitating distributed computation in a fundamentally different fashion. The P2P model is a fully decentralized and distributed topology. With P2P the client would be directly connected with all other computing peers seamlessly, forming a virtual computer [57, 58]. There are multiple Internet connectionsbetween the client peer and other computing peers (Fig. 8). Because there is no longer there such a thing as the server, one node’s failure will not cause the entire system fail. Each peer is the server and the client at the same time. P2P is not a new concept. The idea of a drastically distributed and decentralized application was tried with USENET as early as in 1980’s. It becomes popular recently, however, largely because of its application in some popular software products such as the Napster (software for MP3 music file sharing) and Gnutella (software for file sharing). [59] As far as Napster is concerned, some may argue that it is actually not a P2P software program in the strict sense because there is still a central server that maintains a database of all peers. Users first install a Napster peer program and then share music files by logging into the server. The server will search the database and match peers who show the file share interests. Without this server, peers would have no idea where other peers are. Therefore, here comes the question what is a P2P network. Generally speaking, in a P2P network a peer should be able to: 1) Find other peers and form or join a small network of those peers who have common interests; 2) Broadcast the resource it owns and discover new resources from other peers. 3) Control its own resource with privilege and share the resource with other peers in collaboration. 4) Perform all the above tasks regardless the operating system and network topology. Of course, there are issues of security, policy, and copyright that are inherent to the P2P network model. Any P2P system that means to be a real application must address these issues properly.(a)(b)Fig. 11. Client/Server (a) and peer-to-peer (b) decentralized topology models.The nice thing comes with P2P in medical imaging is its great flexibility, accountable robustness and scalability. Imagine that a patient at a clinic just receives a CT scan. The doctor wants to see the reconstructed images for diagnosis. The clinic is connected to the hospital’s imaging center via a broadband LAN. The center utilizes a high performance cluster for fast reconstruction. Now suppose the master node has crashed. What will happen next is that the whole cluster becomes unavailable to the outside world even if the rest computing nodes are still functioning. In fact we haveexperienced similar situations when the remote cluster server for parallel computing failed, causing the entire research disrupted. Recently we have built a P2P-enhanced network for iterative image reconstruction. The performance result (shown in Fig. 11) is comparable to a cluster of a similar setup, [57, 58]. It indicates that in a clinic, one can integrate PCs in each room to virtually construction a parallel system for high performance medical image reconstruction. Therefore in the future, P2P technology is very promising, since the personal desktop machines become cheaper and cheaper.Speedup vs. Iteration No.12.00000 10.00000 speedup 8.00000 6.00000 4.00000 2.00000 0.00000 1 2 3 4 5 6 7 8 9 10 no of computational nodes iter=2 iter=3 iter=4 iter=5 iter=6 iter=7 iter=8 iter=9 iter=10Fig. 12. Performance in terms of speedup vs. the number of computing nodes and iterations used, in a P2P distributed computing environment.6. Future WorkDecentralized parallel computing has many merits. Among these are the reliability, low-cost, and economical efficiency. From super-computers to PC clusters; from a local cluster to a remote cluster server, over the decades the computer and the Internet technology have evolved into this trend. However, the computation itself is still centralized regardless of parallelism. Now, the peer-to-peer technology further improves the system reliability and robustness. The concept applies to all the representative iterative and non-iterative algorithms, including EM-type, ART-type, Feldkamp-type and Katsevich-type algorithms for cone-beam image reconstruction, as long as the forward/backward projection is involved. Also, parallel reconstruction techniques are not restricted to transmission and emission CT. Many medical imaging and analysis processes can be accelerated with parallel computing techniques, including optical tomography, image visualization and so on. Recently, the Grid becomes an emerging infrastructure that will change the Internet and the way people think of science and computing. [60-62] In general, the Grid technology views the computing resources, storage, CPU cycle, memory and so forth as a network similar to a electrical network. Any resource, once are grid-enabled, will beready for “plug in-and-use”. Thus, many hospitals may be connected via the Grid and forming a virtual organization (VO). In such a VO, one medical image reconstruction site’s need for parallel computing will be shared across every facility in the Grid. To take advantage of such a highly distributed computing model, however, two things are necessary. One is the high-speed connection such as the TerraByte network. The other is a dedicated reconstruction algorithm that must be sophisticatedly tamed for this diversified distributed parallelism.Reference[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] Seeram E. Computed tomography, physical principles – clinical applications, and quality control. 2nd edition, WB Saunders Co. 2000; 1-8. Burnham CA, Brownell GL. A multi-crystal positron camera. IEEE Transactions on Nuclear Science 1972; 201-205, Chesler DA. Three-dimensional activity distribution from multiple positron scintigraphs. Journal of Nuclear Medicine, 1971; 12: 347-8. Brownell GL. A history of positron imaging. /~glb/alb.html. Kuhl DE, Edwards RQ. Image separation radioisotope scanning. Radiol. 1963; 80: 653-662. Labbe J. SPECT/CT Emerges from the shadow of PET/CT. Biophotonics International. 2003; 50-57. Scarfone C. Single photon emission computed tomography (SPECT). /bae/courses/bae590f/1995/scarfone/index.html. Seeram E. Computed tomography, physical principles – clinical applications, and quality control. 2nd edition, WB Saunders Co. 2000; 23. Goldman LW. Principles of CT and the evolution of CT technology. RSNA Categorical course in diagnostic radiology physics: CT and US cross-sectional imaging. 2000; 33-52. Kalender WA, Seissler W, Vock P. Single-breath-hold spiral volumetric CT by continuous patient translation and scanner rotation. Radiology. 1989; 173: 14. Crawford CR, King KF. Computed tomography scanning with simultaneous patient translation. Med. Phys. 1990; 17: 967-982. Wang G, Lin TH, Cheng PC, Shinozaki DM, Kim HG. Scanning cone-beam reconstruction algorithms for x-ray microtomography. Proceedings of SPIE. 1992; 1556, 99-112. Taguchi K, Aradate H. Algorithm for image reconstruction in multi-slice helical CT. Med. Phys. 1998; 25: 550-561. Feldkamp L, Davis L, and Kress J. Practical cone-beam algorithm. Journal of the Optical Society of America. 1984; 612-619. Kudo H, Saito T. Helical-scan computed tomography using cone-beam projections. IEEE Medical Imaging. 1991; 1958-1962. Wang G, Lin TH, Cheng PC, Shinozaki TM. A general cone-beam reconstruction algorithm. IEEE Trans. Med. Imaging. 1993; 12: 486-496.[17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]Schaller S, Flohr T, Steffen P. A new approximate algorithm for image reconstruction in cone-beam spiral CT at small cone-angles. IEEE Medical Imaging. 1996; 1703-1709. Wang G., Crawford CR, Kalender WA. Multirow detector and cone-beam spiral/helical CT. IEEE Trans. Med. Imaging. 2000; 19: 817-821. Turebell H, Danilsson PE. Helical cone beam tomography. Int. J. Imaging System and Technology. 2000; 11: 91-100. Turbell H. Cone-beam reconstruction using filtered backprojection. Ph.D. thesis, Linköpings Universitet, Sweden 2001. Katsevich A. Theoretically exact filtered backprojection-type inversion algorithm for spiral CT. SIAM J. Appl. Math. 2002; 62: 2012-2026. Katsevich A. An improved exact filtered backprojection algorithm for spiral computed tomography. Advances in Applied Mathematics. 2004; 32 (4): 681–697. Kak AC, Slaney M. Principles of computerized tomographic imaging. SIAM 2001; 50. Tuy HK. An inversion formula for cone beam reconstruction. SIAM J. of Applied Mathematics. 1983; 43: 546-551. Rockmore AJ, Macovski A. A maximum likelihood approach to image reconstruction. IEEE Trans. Nucl. Sci. 1976; NS-23: 1428-1432. Shepp LA, Valdi Y. Maximum likelihood reconstruction for emission tomography. IEEE, Trans. Med. Imag. 1982; MI-1: 113-122. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J. Comput. Assist. Tomog. April 1984; 8(2): 302-316. Andersen AH. Algebraic reconstruction in CT from limited views. IEEE Trans. Med. Imag. 1989; 8: 50-55. Andersen AH, Kak AC. Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm. Ultrasonic Imaging. 1984; 6: 81-94. Jiang M, Wang G. Convergence of the simultaneous algebraic reconstruction technique (SART). IEEE Trans. Image Processing. 2003; 12: 957-961. Lange K. Convergence of EM image reconstruction algorithms with Gibbs smoothing. IEEE Trans. Med. Imaging 1990; 9 (4): 439-446. Lange K, Fessler JA. Globally convergent algorithms for maximum a posteriori transmission tomography. IEEE Trans. Image Processing. 1995; 4 (10): 1430-1438. Snyder DL, Schulz TJ, O'Sullivan JA. Deblurring subject to nonnegativity constraints. IEEE Trans. Signal Processing. 1992; 40: 1143-1150. Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans. Med. Imag. 1994; 13: 601-609. Kamphuis C, Beekman FJ. Accelerated iterative transmission CT reconstruction using an ordered subsets convex algorithm. IEEE Trans. Med. Imaging. December 1998; 17 (6). Erdogan H, Fessler JA. Ordered subsets algorithms for transmission tomography. Phys. Med. Biol. 1999; 44(11). Xu F. Tomographic Reconstruction using graphics hardware. Nov. 2003. Cabral B, Cam N, Foran J. Accelerated volume rendering and tomographic reconstruction using。