Open Access
Add to BibSonomyAbstract
We investigate in this paper a new reconstruction method in order to perform 3D Terahertz (THz) tomography using a continuous wave acquisition setup in transmission mode. This method is based on the Maximum Likelihood for TRansmission tomography (ML-TR) first developed for X-ray imaging. We optimize the Ordered Subsets Convex (OSC) implementation of the ML-TR by including the Gaussian propagation model of THz waves and take into account the intensity distributions of both blank calibration scan and dark-field measured on THz detectors. THz ML-TR reconstruction quality and accuracy are discussed and compared to other tomographic reconstructions.
© 2014 Optical Society of America
1. Introduction
Terahertz technology (between 0.1 and 4 THz) is now a well-established tool to achieve contact-free and non-destructive testing (NDT) [1–7]. In the field of 3D imaging, X-Ray computed tomography (CT) is an omnipresent technique which provides 3D visualization of dense materials such as human bodies, biological tissues or industrial samples. Nevertheless, this powerful technique cannot be easily applied to soft materials owing to the low absorption of the X-Ray radiation. Alternatively, THz radiation offers a good penetration depth through various soft materials such as plastics, papers or paintings, allowing direct applications in non-destructive inspection. Most demonstrations to date have been performed in 2D, however 3D THz CT is an emerging technique which has been investigated during the last decade [8–17].
Low NEP (Noise Equivalent Power) and fast detectors have been developed and demonstrated in laboratory but transfer to industry is quite limited. Thus detectors are usually mono-to few-pixels large, resulting in a long raster scanning process to acquire 2D images. Additionally, in order to perform 3D THz CT, several radiographs have to be acquired around the sample at different viewing angles.
Tomographic algorithms allowing accurate reconstruction from a limited number of projections are essential in order to overcome these current technical limitations. An important consideration in 3D THz CT concerns the choice of the reconstruction method to be able to visualize the different cross-sectional images and the final 3D volume of the sample. Iterative reconstructions developed for X-ray CT, such as the Simultaneous Algebraic Reconstruction Technique (SART) [18] and the Ordered Subsets Expectation Maximization (OSEM) [19, 20], have been investigated for 3D THz CT in [14, 15]. These methods are especially known for their ability to provide a high quality reconstruction even if the acquired data consists of a few number of projections. However they do not take into account a realistic transmission model.
Thus in this paper we propose to investigate the statistical reconstruction method denoted Maximum Likelihood for TRansmission (ML-TR) tomography[21, 22]. This method is based on Poisson distribution model of transmitted radiation and allows the introduction of some a priori knowledge about the imaged object [23–25]. Specifically, we focus on the implementation denoted the Ordered Subsets Convex (OSC) algorithm since it has an efficient convergence rate despite the noise level and sparsity of acquired data [26, 27].
The paper is organized as follows: in the next section, we describe the OSC algorithm and specific optimizations to take into account THz radiation properties. In particular we include the Gaussian beam propagation model proposed in [15, 16]. Section 3 describes the experimental setup used to perform 3D THz acquisitions and measurements needed by the OSC algorithm. Before concluding, we discuss the reconstruction quality in section 4, by comparing the new algorithm with other standard reconstructions.
2. Ordered subsets convex algorithm for 3D THz tomography
Expectation maximization (EM) refers to a wide class of iterative reconstructions now well established for X-ray transmission tomography [23–25]. These methods reconstruct the 3D structural volume of the sample from a set of radiographs acquired around the object. Assuming Poisson distribution of the acquired data, the probability to observe the measurements R according to the sample μ, denoted p(R|μ), is defined by:
where i is a projection line index and R̂ is the expected intensity measurement, proportional to the number photon counts, at detector position i. R̂ is modeled by the Beer-Lambert law as follows: where I0 corresponds to the maximal photon counts (blank calibration scan), bg is the background noise (dark-field) and wij is a weight coefficient defining the voxel j contribution on the detector measurement i. wij is usually proportional to the distance traveled by the line i in the voxel j. I0 and bg are determined experimentally (cf. section 3.3).Based on the likelihood p(R|μ) and the transmission model R̂, the method consists of max-imazing the log likelihood log p(R|μ), ie. finding the solution μ that minimizes the partial derivative . The solution can not be determined directly since the partial derivative leads to a transcendantal equation owing to the exponential part in the transmission model. Thus several approximations have been proposed based on series expansion or Newton numerical analysis [22, 25]. The latter in particular leads to the ordered subsets convex (OSC) algorithm [26]. We define and adapt the OSC algorithm for THz tomography in the following.
2.1. Ordered subsets convex (OSC) algorithm
The ordered subsets convex (OSC) algorithm [26] consists of iterating in t and subsets s + 1 in order to update each voxel j of the volume μ until convergence of the solution. The volume obtained by update using the subset s is used as starting volume of the next subset s + 1. A main iteration t is completed when all subsets have been processed. The OSC algorithm updates each voxel as follows:
where is the expected photon counts computed from using Eq. (2), αt is a relaxation parameter and S(s) are the radiographs in the subset s. αt = 1 for all t in the following.2.2. THz transmission including Gaussian beam model
In X-ray CT transmission tomography, the beam shape can be considered constant since the wavelength is small with respect to the object size so that there is no diffraction effect. For the same reason, the intensity distribution is uniform over the beam cross-section. This property allows transmission model Eq. (2) to be used directly in all reconstruction methods.
In THz CT imaging, the propagation beam is close to a Gaussian distribution (cf. Fig. 1) which is both determined by the THz wave properties and the lens used to focus the beam. The radius of the beam (from the beam axis) has its minimum value w0 at beam waist. According to the wavelength λ, the radius at the position z from the beam waist is:
where is the Rayleigh range. Moreover, the intensity distribution over the cross-section is given by: where j is the 3D position vector related to the distance r from the beam axis and z (cf. Fig 1(a)), and Imax is the intensity at the center of the beam waist. Such a propagation model is included in Eq. (2) by defining: where ★ denotes the convolution operator. It results that the volume reconstructed by the OSC is no more μ but f, leading to the following update step: where μ is now defined by Eq. (6) in the denominator as well as in the computation of . The overall algorithm, denoted THz ML-TR in the following, processes as follows: 1) compute the projection values from current estimate of f for all i in the subset s using Eq. (2), where μ is given by Eq. (6); 2) update volume f using Eq. (7); 3) repeat steps 1–2 for all subsets; 4) iterate steps 1 to 3 until convergence of the solution. In our reconstructions, we have stopped the iterations if, i) the projection quadratic error is less than 0.5% of the initial backprojected error ∑i [R(i)]2, or if, ii) a maximum of 10 iterations have been performed (to avoid long reconstruction time if i) is hard to achieve). We have usually observed that the algorithm stopped between 7 and 9 iterations, with a residual quadratic error less than 0.5 %.
Fig. 1 (a) Simulated propagation of the beam for a 287 GHz source (position in meters), w0 = 2.3 mm at beam waist according to previously measured beam. (b) Gaussian intensity distribution at beam waist (Intensity (a.u)).
3. Experimental acquisition setup
We first describe in this section the experimental setup used in order to perform THz tomography. Since the beam propagation follows a Gaussian distribution, we detail how to model this distribution from measurements of the beam waist. Then the tomographic acquisition of the sample is described. Finally, we measure and define the blank scan, I0, and dark-field, bg, models used in the THz ML-TR reconstruction.
3.1. Description of the experimental setup
The experimental setup of the 3D millimeter wave tomographic scanner is a Gunn diode and a tripler coupled with a horn antenna (cf. Fig. 2). The output power is 12 mW at 287 GHz. The THz beam is then collimated and focused with a pair of PTFE (Polytetrafluoroethylene) lenses (f = 100 mm focal length and D = 50.8 mm). The sample is positioned on a three-axes motorized stage comprising the X,Y and θ movements, respectively. The detection is performed with a Schottky diode and the beam is modulated at 350 Hz by a mechanical chopper. The amplitude of the transmitted THz signal is acquired with a lock-in amplifier (time-constant: 30ms). Then we checked the 2D transversal profile of the THz beam at the beam waist and outside the Rayleigh zone. At the sample position, the beam profile is homogeneous with a Gaussian circular shape (2.3 mm beam diameter, measured at full width half maximum (FWHM)) in agreement with the theoretical values obtained from the propagation of Gaussian beam models.
Fig. 2 Experimental setup: THz beam, delivered by a Gunn diode and a tripler (power 12 mW at 287 GHz) is collimated and focussed with a pair of PTFE lenses. Acquired sample is positionated on a three-axes motorized stage comprising the X,Y and θ movements.
Using Eq. (5), a simulation of the 3D ray profile can be performed according to the wavelength and the measured w0 ≈ FWHM at beam waist. Fig. 1(a) shows the 3D profile of simulated 287 GHz source focused with w0 = 2.3 mm and Fig. 1(b) illustrates the corresponding Gaussian intensity distribution at beam waist. This model is used in the THz ML-TR algorithm to take the Gaussian propagation Eq. (5) into account during reconstruction. According to the wavelength and beam waist, the Rayleigh range of this setup, positionated between z1 and z2 along the propagation axis such that w(z1) = w(z2) = 2w0 is about 55.13mm.
3.2. Sample acquisition
The sample studied in this paper is a head spray shown on Fig. 3(a). Several regions of this object are out of the Rayleigh range since its largest dimension is greated than 90mm. A 2D transmission radiograph of the sample is obtained by moving the object in the X and Y directions, with a scan step of 1mm in both directions. For 3D tomography, the sample is then rotated in order to provide different radiographs of the object. The overall acquisition of the head spray is processed for 36 projections uniformly distributed around 180° with an angle step of 5° and each projection size is 156 × 100 pixels (acquisition time is about 8min per projection). The measured projections are shown on Fig. 3(b) and Media 1.
Fig. 3 (a) Picture of the acquired head spray. (b) 36 projections of the sample acquired uniformly in the angle range [0°, 180°] with an angle step of 5° between two consecutive projections ( Media 1).
3.3. Measurements of blank scan and dark-field
In order to process the line acquisition Eq. (2), the blank scan, I0, and dark-field, bg, have to be determined. The blank scan corresponds to the maximal measured intensity when the source is lighting the detector. The background corresponds to the residual ambient noise received by the detector when the source is inactive.
First, we have measured the blank scan for 5 projections (the same size as for sample acquisition) by illuminating the detector with the THz source (without a sample). Thus according to the projection size, we have measured 5 × 156 × 100 = 78000 values, in order to estimate the source intensity distribution. The histogram provided in Fig. 4(a) depicts the blank-scan distribution (Y-Axis) according to the measured intensity (in mA). This distribution can be modeled by a normal distribution 𝒩 ( , σI0), where and σI0 = 0.02056 (cf. green curve on Fig 4(a)). From this first approximation, a fitting algorithm is applied to refine the blank scan distribution parameters. We obtain finally and σI0 = 0.0165 (cf. blue curve on Fig. 4(a)). Second, we have processed in a similar fashion the dark-field distribution. The measurements (cf. red histogram on Fig. 4(b)) follow a normal distribution approximated by and σbg = 0.00044. We obtain and σbg = 0.00035 after refinement.
Fig. 4 (a) blank scan distribution : measured distribution (red), first approximate (green) and final model obtained by a Gaussian fitting algorithm (blue). (b) dark-field distribution : measured distribution (red), first approximate (green) and final model obtained by a Gaussian fitting algorithm (blue).
Thus efficient models of blank scan and dark-field in Eq. 2 can be the distributions 𝒩 ( , σI0) and 𝒩 ( , σbg) determined experimentally. In that case, each value I0(i) and bg(i) is chosen such that a value sets (for all i in one projection) follow these distributions, respectively. However, due to the very low observed standard deviations, we simplify these models by defining and for all i.
4. Discussion
In this section, we discuss the quality of THz ML-TR reconstruction compared to backprojection of filtered projections (BFP) and the simultaneous algebraic reconstruction technique developed for THz tomography (THz-SART) in [15,16]. BFP reconstructs each voxel by averaging projection values crossing it. To assume that a voxel is a linear combination of projection values, these values have to represent a linear combination of volume data. Similarly, SART is based on the Karczmarz algorithm to solve iteratively the linear system of equations between voxels and projections. Then, contrary to the THz ML-TR, these methods are based on the absorbance to recover the solution f from all A(i) = ∑j wij [f(j) ★ G(j)]. Projection preprocessing before BFP or THz-SART is done using the same and values than those used in the THz ML-TR method.
4.1. Quality and accuracy observed on reconstructed 2D slices
Fig. 5 shows a set of 6 cross-sectional images along the Y-Axis obtained by the standard BFP, the THz-SART reconstruction and the new THz ML-TR method, respectively. The SART result is obtained after 2 iterations. THz ML-TR convergence is obtained after 10 iterations using 2 subsets. Even if reconstruction time is quite larger, it is negligible compared to the acquisition time.
Fig. 5 Few cross-sections along Y-Axis of the imaged object obtained respectively by BFP, THz-SART and THz ML-TR methods (from the first to the last line, respectively). From left to right: y = 15, 30, 45, 60, 75, and 90.
First we can notice that the reconstructed object is surrounded by many artifacts on BFP results. This method is particularly sensitive to the noise and it is unable to reconstruct accurately from a few number of projections. Indeed, it is common to acquire at least N projections sized N × M to reconstruct accurately a volume of M slices sized N2. Thus in our study case, a correct BFP reconstruction of the head spray could be achieved from 156 acquired radiographs. Even if SART reconstruction quality seems better, one can remark the presence of few streak artifacts (high intensity lines). Moreover, several parts such as the base (left image) and outer surface of the acquired object are not reconstructed whereas the intensity of the other ones is very high.
The artifacts owing to both the BFP and SART methods are clearly identified compared to the THz ML-TR reconstruction in the difference images given in Fig. 6(a–b). The transmission model used in the THz ML-TR method leads to a better uniformity of reconstructed intensities: both the inner and outer parts of the object are recovered with high quality. Since transmission value rapidly approaches zero when the beam is going through a high density material, it is close to a nil measurement (nothing detected on the sensor) experienced when the beam is completely absorbed by the sample. Thus THz ML-TR is less sensitive to full absorption effect. Conversely, full absorption provides very high and locally non-proportional values on the absorbance projections. Thus it leads to streak artifacts in the reconstructions, as visible on the BFP, and, to a lesser extent, on the SART results.
Fig. 6 (a) Difference image BFP at THz ML-TR, and (b) SART and THz ML-TR at position y = 80. (c) Slice at position y = 15 obtained by ML-TR without Gaussian beam model. (d) Difference image between ML-TR and THz ML-TR of cross-sectional images at y = 15.
Moreover, a better accuracy is obtained using the Gaussian beam model in the reconstruction. As an illustration of the improvements, Fig. 6(c) shows ML-TR without Gaussian model, and Fig. 6(d) the difference image between ML-TR and THz ML-TR. We can notice a blur in ML-TR result, even at the centre of the object, owing to the beam propagation considered as linear in the Rayleigh length. Thus optimized reconstruction corrects for blur effect whatever the object position with respect to the Rayleigh range so that Rayleigh length is no longer a limiting factor for 3D THz tomography. As an illustration we can notice an equivalent accuracy of the object shapes which are inside and outside the Rayleigh range in images Fig. 5(c).
4.2. Quality and accuracy discussed from 3D rendering images
In order to confirm the new method efficiency compared to the others, we discuss now the feasability of a 3D segmentation in order to provide a 3D rendering of the sample. This study is realized using a new software developed in our group which performs processing sequences to analyze 3D THz images. We will describe the software capabilities in a future paper.
Usually, one applies prefiltering to reduce noise and preprocessing to enhance contrast before 3D rendering. Even if these processings/filterings lead to a better 3D visualisation, they can reduce or erase relevant data about the sample, leading to a biaised analysis of the object. In other words, the inspection performance is directly linked to the tomographic reconstruction quality. Thus in our study, segmentation and rendering are performed directly from reconstructed images (without any preprocessing and prefiltering), allowing us to discuss tomographic reconstruction quality objectively.
The 3D rendering corresponding to the outer sample surface is shown on Fig. 7. Noise of BFP result makes 3D visualisation absolutely wrong and impractical for any outer shape analysis or measures (cf. Fig. 7(a)). 3D rendering of THz-SART volume (cf. Fig. 7(b)) is slightly better since we can distinguish partially the object structure. However, as we have already noticed in Fig. 5, outer surface disappears on the bottom of the object and few streak artifacts penalize 3D visualisation. Conversely, 3D rendering of the THz ML-TR volume ((cf. Fig. 7(c)) allows a 3D visualisation without any noise or streak artifact. From such a 3D rendered volume, object analysis can thus be performed more accurately. For instance, video recordings (cf. Fig. 8 and Media 2) show 3D rendering view with 2D slice overprints. These visualisations reveal both the inner and outer surfaces composing the object, allowing volume shape analysis. Moreover, specific parts of the object such as the inner tube and the spray mechanism, are clearly distinguishable. Due to the noise and artifacts, similar views are unavailable from equivalent unpreprocessed and unprefiltered BFP and SART reconstructions.
Fig. 8 3D visualisation of the object reconstructed by THz ML-TR. On the left: 3D meshing view, allowing both surface and volume analysis. On the right: 3D dotted view + 2D slice overprints, revealing the inner spray mechanism ( Media 2).
However, despite the superiority of THz ML-TR over SART and BFP, one can remark that holes, or void artifacts appear on the reconstructed volume (denoted by the arrows on Fig. 7(c)). These holes are in very low absorbent regions of the object on several viewing angles. Indeed, on projection images, Fig. 3(b), one can notice that the transmission values in these regions are very close to the maximum transmitted signal. Since this phenomenon is experienced on several projections, it leads to very low reconstructed voxel intensities, which are confused with background by the 3D segmentation algorithm. Thus these void artifacts indicate more an issue of the experimental setup (i.e. not enough projections and/or beam energy not adjusted on the absorption provided by the sample on these regions) rather than a limitation of the reconstruction. Since THz ML-TR provides accurate results yet, it could also be used to determine the experimental setup quality using phantom acquisitions, and to adjust scanner parameters for the sample.
5. Conclusion
In this paper, the OSC algorithm developed to perform ML-TR reconstruction has been optimized for 3D THz tomography. Optimisations consist of i) including the Gaussian beam model in the computation of expected photon counts, ii) modeling the blank scan and dark-field from measurements and, iii) taking into account these models in the algorithm. The new method, called THz ML-TR, has been compared to BFP and THz-SART reconstructions.
Comparisons have demonstrated that the THz ML-TR method outperforms other reconstructions. Including a more realistic transmission model, it provides a better accuracy of outer and inner structures of the object. 3D rendering, usually processed to perform internal structure inspection, has also demonstrated the quality of THz ML-TR reconstruction. Conversely, other methods are hampered by noise or streak artifacts, making structure inspection hard to achieve.
Future works will focus on the development of THz ML-TR for reconstructing from multi-energy and/or multi-intensity acquisitions. We expect that such an approach could solve the hole artifacts observed on the THz ML-TR results. Thus the simplification made here to include dark-field and blank scan models into the reconstruction process will have to be reconsidered and a dual reconstruction update step will be investigated.
Acknowledgments
The authors acknowledge the Agence Nationale de la Recherche (ANR) for their financial support in the InPoSec project (www.inposec.org). Thanks are also due to the referees for their insightful comments.
References and links
1. K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue, “Non-destructive terahertz imaging of illicit drugs using spectral fingerprints,” Opt. Express 11, 2549–2554 (2003). [CrossRef] [PubMed]
2. Y. C. Shen, T. Lo, P. F. Taday, B. E. Cole, W. R. Tribe, and M. C. Kemp, “Detection and identification of explosives using terahertz pulsed spectroscopic imaging,” Appl. Phys. Lett. 86, 241116 (2005). [CrossRef]
3. W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325–1379 (2007). [CrossRef]
4. S. Y. Huang, Y. X. J. Wang, D. K. W. Yeung, A. T. Ahuja, Y. T. Zhang, and E. Pickwell-MacPherson, “Tissue characterization using terahertz pulsed imaging in reflection geometry,” Phys. Med. Biol. 54, 149–160 (2009). [CrossRef]
5. K. Fukunaga and M. Picollo, “Terahertz spectroscopy applied to the analysis of artists materials,” Appl. Phys. A 100, 591–597 (2010). [CrossRef]
6. J.-P. Caumes, A. Younus, S. Salort, B. Chassagne, B. Recur, A. Ziéglé, A. Dautant, and E. Abraham, “Terahertz tomographic imaging of xviiith dynasty egyptian sealed pottery,” Appl. Optics 50, 3604–3608 (2011). [CrossRef]
7. F. Ospald, W. Zouaghi, R. Beigang, C. Matheis, J. Jonuscheit, B. Recur, J.-P. Guillet, P. Mounaix, W. Vleugels, P. V. Bosom, L. V. González, I. López, R. Martínez Edo, Y. Sternberg, and M. Vandewal, “Aeronautics composite material inspection with a terahertz time-domain spectroscopy system,” Opt. Eng. 53(3), 031208 (2014). [CrossRef]
8. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]
9. S. Wang, B. Ferguson, D. Abbott, and X. C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29, 247–256 (2003). [CrossRef] [PubMed]
10. S. Wang and X. C. Zhang, “Pulsed terahertz tomography,” J. Phys. D: Appl. Phys. 37, R1–R36 (2004). [CrossRef]
11. M. M. Awad and R. A. Cheville, “Transmission terahertz waveguide-based imaging below the diffraction limit,” Appl. Phys. Lett. 86, 221107 (2005). [CrossRef]
12. X. Yin, B. W. H. Ng, B. Ferguson, and D. Abbott, “Wavelet based local tomographic image using terahertz techniques,” Digital Signal Process. 19, 750–763 (2009). [CrossRef]
13. A. Brahm, M. Kunz, S. Riehemann, G. Notni, and A. Tünnermann, “Volumetric spectral analysis of materials using terahertz-tomography techniques,” Appl. Phys. B 100, 151–158 (2010). [CrossRef]
14. B. Recur, A. Younus, S. Salort, P. Mounaix, B. Chassagne, P. Desbarats, J-P. Caumes, and E. Abraham, “Investigation on reconstruction methods applied to 3D terahertz computed tomography,” Opt. Express 19, 5105–5117 (2011). [CrossRef] [PubMed]
15. B. Recur, J.-P. Guillet, I. Manek-Hönninger, J.-C. Delagnes, W. Benharbone, P. Desbarats, J.-P. Domenger, L. Canioni, and P. Mounaix, “Propagation beam consideration for 3d thz computed tomography,” Opt. Express 20, 5817–5829 (2012). [CrossRef] [PubMed]
16. B. Recur, J.-P. Guillet, L. Bassel, C. Fragnol, I. Manek-Hönninger, J.-C. Delagnes, W. Benharbone, P. Desbarats, J.-P. Domenger, and P. Mounaix, “THz radiations for tomographic inspection,” Opt. Eng. 51(9), 091609 (2012). [CrossRef]
17. JP Guillet, B. Recur, L. Frederique, B. Bousquet, L. Canioni, I. Manek-Hönninger, P. Desbarats, and P. Mounaix, “Review of Terahertz tomography techniques,” J. Infrared Milli Terahz 35(4), 382–411 (2014). [CrossRef]
18. A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART) : A superior implementation of the ART algorithm,” Ultrasonic Imaging 6, 81–94 (1984). [CrossRef]
19. L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging 1, 113–122 (1982). [CrossRef] [PubMed]
20. H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging 13, 601–609 (1994). [CrossRef] [PubMed]
21. J. A. Fessler, “Hybrid poisson/polynomial objective functions for tomographic image reconstruction from transmission scans,” IEEE Trans. on Im. Processing 4, 1439–1450 (1995). [CrossRef]
22. B. De Man, S. Basu, J.B. Thibault, J. Hsieh, J. Fessler, C. Bouman, and K. Sauer, “A study of four minimization approaches for iterative reconstruction in X-ray CT,” IEEE Nuclear Science Symposium Conference Record M11-3395(2005).
23. K. Lange and R. Carson, “EM reconstruction algorithms for emission and transmission tomography,” J Comput Assist Tomogr 8, 306–316 (1984). [PubMed]
24. K. Lange, M. Bahn, and R. Little, “A theoretical study of some maximum likelihood algorithms for emission and transmission tomography,” IEEE Trans. on Med. Imaging 6, 106–114 (1987). [CrossRef]
25. K. Lange and J. A. Fessler, “Globally convergent algorithms for maximum a posteriori transmission tomography,” IEEE Trans. on Im. Processing 4, 1430–1438 (1995). [CrossRef]
26. C. Kamphuis and F. J. Beekman, “Accelerated iterative transmission ct reconstruction using an ordered subsets convex algorithm,” IEEE Trans. on Med. Imaging 17, 1101–1105 (1998). [CrossRef]
27. H. Erdogan and J. A. Fessler, “Ordered subsets algorithms for transmission tomography,” Phys. Med. Biol. 44(11), R2835 (1999). [CrossRef]
