Radiation dose reduction for chest radiography of infants in intensive care units using a high peak kilovoltage-technique.

BACKGROUND: Chest radiography is an important tool in the care of infants in intensive care units. Image optimization must be monitored to minimize radiation exposure in this susceptible population. OBJECTIVE: To examine the use of a high tube peak kilovoltage technique to achieve radiation dose reduction while maintaining adequate image quality. MATERIALS AND METHODS: A retrospective study was conducted. Radiation doses of chest radiographs performed in the pediatric intensive care units in our institution were calculated. The radiographs were divided into two groups based on the value of the peak kilovoltage used: above and below 60 kilovolts (kV). Image quality was blindly assessed by two fellowship-trained pediatric radiologists. Air kerma, effective dose and quality score for the high versus the low peak kilovoltage group were compared and analyzed. RESULTS: The study included 376 radiographs. One hundred and seven radiographs were performed using peak kilovoltage values equal to or above 60 kV and 269 radiographs were performed using values under 60 kV. The average air kerma for the lower peak kilovoltage group was 56.6 microgray (µGy) (30.7-81.9) vs. 22.9 µGy (11.8-34.4) for the higher peak kilovoltage group (P<0.0001). The mean difference in effective dose between the groups was 11.68 (P<0.0001). The mean difference for the quality score was 0.06 (±0.03, P=0.10), not statistically significant. CONCLUSION: A high peak kilovoltage technique may enable a statistically significant radiation dose reduction without compromising the diagnostic value of the image.

Performance optimization of a tri-hybrid method for estimation of patient scatter into the EPID.

On-treatment EPID images are contaminated with patient-generated scattered photons. If this component can be accurately estimated, its effect can be removed, and therefore a corresponding in vivo patient dose estimate will be more accurate. Our group previously developed a "tri-hybrid" (TH) algorithm to provide fast but accurate estimates of patient-generated photon scatter. The algorithm uses an analytical method to solve for singly-scattered photon fluence, a modified Monte Carlo hybrid method to solve for multiply-scattered photon fluence, and a pencil beam scatter kernel method to solve for electron interaction generated scattered photon fluence. However, for efficient clinical implementation, spatial and energy sampling must be optimized for speed while maintaining overall accuracy. In this work, the most significant sampling issues were examined, including spatial sampling settings for the patient voxel size, the number of Monte Carlo histories used in the modified hybrid MC method, scatter order sampling for the hybrid method, and also a range of energy spectrum sampling (i.e., energy bin sizes). The total predicted patient-scattered photon fluence entering the EPID was compared with full MC simulation (EGSnrc) for validation. Three phantoms were tested with 6 and 18 MV beam energies, field sizes of 4 × 4, 10 × 10, and 20 × 20 cm2 , and source-to-imager distance of 140 cm to develop a set of optimal sampling settings. With the recommended sampling, accuracy and precision of the total-scattered energy fluence of the TH patient scatter prediction method are within 0.9% and 1.2%, respectively, for all test cases compared with full MC simulation results. For the mean energy spectrum across the imaging plane, comparison of TH with full MC simulation showed 95% overlap. This study has optimized sampling settings so that they have minimal impact on patient scatter prediction accuracy while maintaining maximum execution speed, a critical step for future clinical implementation.

A tri-hybrid method to estimate the patient-generated scattered photon fluence components to the EPID image plane.

In vivo dosimetry methods can verify the prescription dose is delivered to the patient during treatment. Unfortunately, in exit dosimetry, the megavoltage image is contaminated with patient-generated scattered photons. However, estimation and removal of the effect of this fluence improves accuracy of in vivo dosimetry methods. This work develops a 'tri-hybrid' algorithm combining analytical, Monte Carlo (MC) and pencil-beam scatter kernel methods to provide accurate estimates of the total patient-generated scattered photon fluence entering the MV imager. For the multiply-scattered photon fluence, a modified MC simulation method was applied, using only a few histories. From each second- and higher-order interaction site in the simulation, energy fluence entering all pixels of the imager was calculated using analytical methods. For photon fluence generated by electron interactions in the patient (i.e. bremsstrahlung and positron annihilation), a convolution/superposition approach was employed using pencil-beam scatter fluence kernels as a function of patient thickness and air gap distance, superposed on the incident fluence distribution. The total patient-scattered photon fluence entering the imager was compared with a corresponding full MC simulation (EGSnrc) for several test cases. These included three geometric phantoms (water, half-water/half-lung, computed tomography thorax) using monoenergetic (1.5, 5.5 and 12.5 MeV) and polyenergetic (6 and 18 MV) photon beams, 10 × 10 cm2 field, source-to-surface distance 100 cm, source-to-imager distance 150 cm and 40 × 40 cm2 imager. The proposed tri-hybrid method is demonstrated to agree well with full MC simulation, with the average fluence differences and standard deviations found to be within 0.5% and 1%, respectively, for test cases examined here. The method, as implemented here with a single CPU (non-parallelized), takes ∼80 s, which is considerably shorter compared to full MC simulation (∼30 h). This is a promising method for fast yet accurate calculation of patient-scattered fluence at the imaging plane for in vivo dosimetry applications.

Technical note: development and validation of a Monte Carlo tool for analysis of patient-generated photon scatter.

Scattered radiation unavoidably generated in the patient will negatively impact both kilovoltage (KV) and megavoltage (MV) imaging applications. Recently, 'hybrid' methods (i.e. combining analytical and Monte Carlo (MC) techniques) are being investigated as a solution to accurately yet quickly calculate the scattered contribution for both KV and MV images. We have developed a customized MC simulation user code for investigating the individual components of patient-scattered photon fluence, which serves as a valuable tool in this area of research. The MC tool is based on the EGSnrc/DOSXYZnrc user code. The IAUSFL flag options associated with subroutine AUSGAB, combined with LATCH tracking, are used to classify the various interactions of particles with the media. Photons are grouped into six different categories: primary, 1st Compton scatter, 1st Rayleigh scatter, multiple scatter, bremsstrahlung, and positron annihilation. We take advantage of the geometric boundary check in DOSXYZnrc, to write exiting photon particle information to a phase-space file. The tool is validated using homogeneous and heterogeneous phantom configurations with monoenergetic and polyenergetic beams under parallel and divergent beam geometry, comparing MC-simulated exit primary fluence and singly-scattered fluence to corresponding analytical calculations. This MC tool has been validated to separately score the primary and scatter fluence components of the KV and MV imaging applications in the field of radiation therapy. The results are acceptable for the various configurations and beam energies tested here. Overall, the mean percentage differences are less than 0.2% and standard deviations less than 1.6%. This will be a critical test instrument for research in photon scatter applications and particularly for the development of hybrid methods, and is freely available from the authors for research purposes.5.

COMP Report: A survey of radiation safety regulations for medical imaging x-ray equipment in Canada.

X-ray regulations and room design methodology vary widely across Canada. The Canadian Organization of Medical Physicists (COMP) conducted a survey in 2016/2017 to provide a useful snapshot of existing variations in rules and methodologies for human patient medical imaging facilities. Some jurisdictions no longer have radiation safety regulatory requirements and COMP is concerned that lack of regulatory oversight might erode safe practices. Harmonized standards will facilitate oversight that will ensure continued attention is given to public safety and to control workplace exposure. COMP encourages all Canadian jurisdictions to adopt the dose limits and constraints outlined in Health Canada Safety Code 35 with the codicil that the design standards be updated to those outlined in NCRP 147 and BIR 2012.

Objective performance assessment of five computed tomography iterative reconstruction algorithms.

OBJECTIVE: Iterative algorithms are gaining clinical acceptance in CT. We performed objective phantom-based image quality evaluation of five commercial iterative reconstruction algorithms available on four different multi-detector CT (MDCT) scanners at different dose levels as well as the conventional filtered back-projection (FBP) reconstruction. METHODS: Using the Catphan500 phantom, we evaluated image noise, contrast-to-noise ratio (CNR), modulation transfer function (MTF) and noise-power spectrum (NPS). The algorithms were evaluated over a CTDIvol range of 0.75-18.7 mGy on four major MDCT scanners: GE DiscoveryCT750HD (algorithms: ASIR™ and VEO™); Siemens Somatom Definition AS+ (algorithm: SAFIRE™); Toshiba Aquilion64 (algorithm: AIDR3D™); and Philips Ingenuity iCT256 (algorithm: iDose4™). Images were reconstructed using FBP and the respective iterative algorithms on the four scanners. RESULTS: Use of iterative algorithms decreased image noise and increased CNR, relative to FBP. In the dose range of 1.3-1.5 mGy, noise reduction using iterative algorithms was in the range of 11%-51% on GE DiscoveryCT750HD, 10%-52% on Siemens Somatom Definition AS+, 49%-62% on Toshiba Aquilion64, and 13%-44% on Philips Ingenuity iCT256. The corresponding CNR increase was in the range 11%-105% on GE, 11%-106% on Siemens, 85%-145% on Toshiba and 13%-77% on Philips respectively. Most algorithms did not affect the MTF, except for VEO™ which produced an increase in the limiting resolution of up to 30%. A shift in the peak of the NPS curve towards lower frequencies and a decrease in NPS amplitude were obtained with all iterative algorithms. VEO™ required long reconstruction times, while all other algorithms produced reconstructions in real time. Compared to FBP, iterative algorithms reduced image noise and increased CNR. CONCLUSIONS: The iterative algorithms available on different scanners achieved different levels of noise reduction and CNR increase while spatial resolution improvements were obtained only with VEO™. This study is useful in that it provides performance assessment of the iterative algorithms available from several mainstream CT manufacturers.

Fast analytical scatter estimation using graphics processing units.

PURPOSE: To develop a fast patient-specific analytical estimator of first-order Compton and Rayleigh scatter in cone-beam computed tomography, implemented using graphics processing units. METHODS: The authors developed an analytical estimator for first-order Compton and Rayleigh scatter in a cone-beam computed tomography geometry. The estimator was coded using NVIDIA's CUDA environment for execution on an NVIDIA graphics processing unit. Performance of the analytical estimator was validated by comparison with high-count Monte Carlo simulations for two different numerical phantoms. Monoenergetic analytical simulations were compared with monoenergetic and polyenergetic Monte Carlo simulations. Analytical and Monte Carlo scatter estimates were compared both qualitatively, from visual inspection of images and profiles, and quantitatively, using a scaled root-mean-square difference metric. Reconstruction of simulated cone-beam projection data of an anthropomorphic breast phantom illustrated the potential of this method as a component of a scatter correction algorithm. RESULTS: The monoenergetic analytical and Monte Carlo scatter estimates showed very good agreement. The monoenergetic analytical estimates showed good agreement for Compton single scatter and reasonable agreement for Rayleigh single scatter when compared with polyenergetic Monte Carlo estimates. For a voxelized phantom with dimensions 128 × 128 × 128 voxels and a detector with 256 × 256 pixels, the analytical estimator required 669 seconds for a single projection, using a single NVIDIA 9800 GX2 video card. Accounting for first order scatter in cone-beam image reconstruction improves the contrast to noise ratio of the reconstructed images. CONCLUSION: The analytical scatter estimator, implemented using graphics processing units, provides rapid and accurate estimates of single scatter and with further acceleration and a method to account for multiple scatter may be useful for practical scatter correction schemes.

Estimation of dose-area product-to-effective dose conversion factors for neonatal radiography using PCXMC.

Dose-area product-to-effective dose (E) conversion factors for chest, abdomen and abdomen-chest neonatal radiographs were computed. Seven patient models in the Monte Carlo software, PCXMC, were defined, representing neonates ranging in weight from 0.5 to 6.0 kg. Conversion factors for a tube potential range of 50-80 kVp at two beam filtrations (3.0 mm Al and 3.0 mm Al+0.1 mm Cu) were calculated. For 133 neonatal radiographs, effective dose values determined using these conversion factors were compared with those obtained from PCXMC simulations customised for each radiograph. For a 3.0-kg newborn irradiated at 60 kVp/3.0 mm Al beam filtration, the conversion factors were 2.58, 1.90 and 1.91 µSv (mGy cm(2))(-1) for chest, chest-abdomen and abdomen radiographs, respectively. Average dose difference between the conversion factors and customised dose calculations was 16 %. Disagreement in effective dose was most strongly correlated with under-collimation in the lateral direction.

Survey of clinical doses from computed tomography examinations in the Canadian province of Manitoba.

The purpose of this study was to document CT doses for common CT examinations performed throughout the province of Manitoba. Survey forms were sent out to all provincial CT sites. Thirteen out of sixteen (81 %) sites participated. The authors assessed scans of the brain, routine abdomen-pelvis, routine chest, sinuses, lumbar spine, low-dose lung nodule studies, CT pulmonary angiograms, CT KUBs, CT colonographies and combination chest-abdomen-pelvis exams. Sites recorded scanner model, protocol techniques and patient and dose data for 100 consecutive patients who were scanned with any of the aforementioned examinations. Mean effective doses and standard deviations for the province and for individual scanners were computed. The Kruskal-Wallis test was used to compare the variability of effective doses amongst scanners. The t test was used to compare doses and their provincial ranges between newer and older scanners and scanners that used dose saving tools and those that did not. Abdomen-pelvis, chest and brain scans accounted for over 70 % of scans. Their mean effective doses were 18.0 ± 6.7, 13.2 ± 6.4 and 3.0 ± 1.0 mSv, respectively. Variations in doses amongst scanners were statistically significant. Most examinations were performed at 120 kVp, and no lower kVp was used. Dose variations due to scanner age and use of dose saving tools were not statistically significant. Clinical CT doses in Manitoba are broadly similar to but higher than those reported in other Canadian provinces. Results suggest that further dose reduction can be achieved by modifying scanning techniques, such as using lower kVp. Wide variation in doses amongst different scanners suggests that standardisation of scanning protocols can reduce patient dose. New technological advances, such as dose-reduction software algorithms, can be adopted to reduce patient dose.

Estimation of organ and effective doses from newborn radiography of the chest and abdomen.

Neonatal intensive care patients undergo frequent chest and abdomen radiographic imaging. In this study, the organ doses and the effective dose resulting from combined chest-abdomen radiography of the newborn child are determined. These values are calculated using the Monte Carlo simulation software PCXCM 2.0 and compared with direct dose measurements obtained from thermoluminescent detectors (TLDs) in a physical phantom. The effective dose obtained from PCXMC is 21.2 ± 0.7 µSv and that obtained from TLD measurements is 22.0 ± 0.5 µSv. While the two methods are in close agreement with regard to the effective dose, there is a wide range of variation in organ doses, ranging from 85 % difference for the testes to 1.4 % for the lungs. Large organ dose variations are attributed to organs at the edge of the field of view, or organs with large experimental error or simulation uncertainty. This study suggests that PCXMC can be used to estimate organ and effective doses for newborn patients.

Antiscatter grid use in pediatric digital tomosynthesis imaging.

The objective of this study was to assess the effect of antiscatter grid use on tomosynthesis image quality. We performed an observer study that rated the image quality of digital tomosynthesis scout radiographs and slice images of a Leeds TO.20 contrast-detail test object embedded in acrylic with and without a grid. We considered 10, 15, 20 and 25 cm of acrylic to represent the wide range of patient thicknesses encountered in pediatric imaging. We also acquired and rated images without a grid at an increased patient dose. The readers counted the total number of visible details in each image as a measure of relative image quality. We observed that the antiscatter grid improves tomosynthesis image quality compared to the grid-out case, which received image quality scores similar to grid-in radiography. Our results suggest that, in order to achieve the best image quality in exchange for the increase in patient dose, it may often be appropriate to include an antiscatter grid for pediatric tomosynthesis imaging, particularly if the patient thickness is greater than 10 cm.

A GPU implementation of EGSnrc's Monte Carlo photon transport for imaging applications.

EGSnrc is a well-known Monte Carlo simulation package for coupled electron-photon transport that is widely used in medical physics application. This paper proposes a parallel implementation of the photon transport mechanism of EGSnrc for graphics processing units (GPUs) using NVIDIA's Compute Unified Device Architecture (CUDA). The implementation is specifically designed for imaging applications in the diagnostic energy range and does not model electrons. No approximations or simplifications of the original EGSnrc code were made other than using single floating-point precision instead of double precision and a different random number generator. To avoid performance penalties due to the random nature of the Monte Carlo method, the simulation was divided into smaller steps that could easily be performed in a parallel fashion suitable for GPUs. Speedups of 20 to 40 times for 64(3) to 256(3) voxels were observed while the accuracy of the simulation was preserved. A detailed analysis of the differences between the CUDA simulation and the original EGSnrc was conducted. The two simulations were found to produce equivalent results for scattered photons and an overall systematic deviation of less than 0.08% was observed for primary photons.

Epp: A C++ EGSnrc user code for x-ray imaging and scattering simulations.

PURPOSE: Easy particle propagation (Epp) is a user code for the EGSnrc code package based on the c+ + class library egspp. A main feature of egspp (and Epp) is the ability to use analytical objects to construct simulation geometries. The authors developed Epp to facilitate the simulation of x-ray imaging geometries, especially in the case of scatter studies. While direct use of egspp requires knowledge of c+ +, Epp requires no programming experience. METHODS: Epp's features include calculation of dose deposited in a voxelized phantom and photon propagation to a user-defined imaging plane. Projection images of primary, single Rayleigh scattered, single Compton scattered, and multiple scattered photons may be generated. Epp input files can be nested, allowing for the construction of complex simulation geometries from more basic components. To demonstrate the imaging features of Epp, the authors simulate 38 keV x rays from a point source propagating through a water cylinder 12 cm in diameter, using both analytical and voxelized representations of the cylinder. The simulation generates projection images of primary and scattered photons at a user-defined imaging plane. The authors also simulate dose scoring in the voxelized version of the phantom in both Epp and DOSXYZnrc and examine the accuracy of Epp using the Kawrakow-Fippel test. RESULTS: The results of the imaging simulations with Epp using voxelized and analytical descriptions of the water cylinder agree within 1%. The results of the Kawrakow-Fippel test suggest good agreement between Epp and DOSXYZnrc. CONCLUSIONS: Epp provides the user with useful features, including the ability to build complex geometries from simpler ones and the ability to generate images of scattered and primary photons. There is no inherent computational time saving arising from Epp, except for those arising from egspp's ability to use analytical representations of simulation geometries. Epp agrees with DOSXYZnrc in dose calculation, since they are both based on the well-validated standard EGSnrc radiation transport physics model.

A novel hybrid reconstruction algorithm for first generation incoherent scatter CT (ISCT) of large objects with potential medical imaging applications.

This work presents a first generation incoherent scatter CT (ISCT) hybrid (analytic-iterative) reconstruction algorithm for accurate ρ{e}imaging of objects with clinically relevant sizes. The algorithm reconstructs quantitative images of ρ{e} within a few iterations, avoiding the challenges of optimization based reconstruction algorithms while addressing the limitations of current analytical algorithms. A 4π detector is conceptualized in order to address the issue of directional dependency and is then replaced with a ring of detectors which detect a constant fraction of the scattered photons. The ISCT algorithm corrects for the attenuation of photons using a limited number of iterations and filtered back projection (FBP) for image reconstruction. This results in a hybrid reconstruction algorithm that was tested with sinograms generated by Monte Carlo (MC) and analytical (AN) simulations. Results show that the ISCT algorithm is weakly dependent on the ρ{e} initial estimate. Simulation results show that the proposed algorithm reconstruct ρ{e} images with a mean error of -1% ± 3% for the AN model and from -6% to -8% for the MC model. Finally, the algorithm is capable of reconstructing qualitatively good images even in the presence of multiple scatter. The proposed algorithm would be suitable for in-vivo medical imaging as long as practical limitations can be addressed.

A 1st generation scatter CT algorithm for electron density breast imaging which accounts for bound incoherent, coherent and multiple scatter: A Monte Carlo study.

Breast CT is an emerging modality that reconstructs 3D linear attenuation coefficient (µ) images of the breast. Its tomographic nature reduces the overlap of structures and may improve tissue visualization. Current prototype systems produce large levels of scatter that could be used to reconstruct electron density (ρ _{e}) images. This could potentially enhance diagnosis. We are developing a first generation bench top CT system to investigate the benefits of simultaneous imaging µ and ρ _{e} of the intact breast. The system uses an algorithm capable of reconstructing ρ _{e} images from single Klein-Nishina scatter. It has been suggested that this algorithm may be impractical since measurements include coherent, bound incoherent and multiple scatter. To investigate this, the EGSnrc Monte Carlo (MC) code was used to simulate scans using a first generation system. These simulations were used to quantify the dose per scan, to provide raw data for the ρ _{e} reconstructions and to investigate corrections for multiple and coherent scatter since these can not be directly related to ρ _{e}. MC simulations show that the dose coefficients are similar to those of cone beam breast CT. Coherent scatter is only ∼9% concentrated in scattering angles < 8°. Electron binding reduces the number of incoherently scattered photons but this reduction can be included in the quantification of scatter measured by the system. Multiple scatter was found to be the major source of errors and, if not corrected for, can result in an overestimation of ρ _{e} by more than a factor of two. Empirical corrections, based on breast thickness or radiological path, can be used to reconstruct images where the variance in ρ _{e} error is half of that found in images derived from primary photons only. Although some practical challenges remain in creating a laboratory system, this work has shown that it is possible to reconstruct scatter images of the breast with a 4 mGy dose and further experimental evaluation of this technique is warranted.

Effect of scatter and an antiscatter grid on the performance of a slot-scanning digital mammography system.

The use of a grid increases perceptibility of low contrast objects in mammography. Slot-scan mammography provides a more dose efficient reduction of the scattered radiation reaching the detector than obtained with an antiscatter grid in screen-film or flat-panel digital mammography. In this paper, the potential of using a grid in a slot-scan system to provide a further reduction of scattered radiation is investigated. The components of the digital signal: primary radiation, off-focus radiation, scattered radiation, and optical fluorescence glare in a CsI(Tl) detector were quantified. Based on these measurements, the primary and scatter transmission factors (Tp, Ts), scatter-to-primary ratio (SPR), signal-difference-to-noise ratio (SDNR), and the SDNR improvement factor (K(SDNR)) were obtained. Our results showed that the SPR ranged from 0.05 to 0.19 for breast thicknesses between 2 and 8 cm, respectively. The values of K(SDNR) ranged from 0.85 to 0.94. Because the slot-scanning system has an inherently low SPR, the increase in dose required when the grid is used outweighs the benefit of the small increase in SDNR. It is possible that greater benefit could be achieved by using a grid with a higher Tp, such as obtained using air-core technology.

Automatic exposure control for a slot scanning full field digital mammography system.

Automatic exposure control (AEC) is an important feature in mammography. It enables consistently optimal image exposure despite variations in tissue density and thickness, and user skill level. Full field digital mammography systems cannot employ conventional AEC methods because digital receptors fully absorb the x-ray beam. In this paper we describe an AEC procedure for slot scanning mammography. With slot scanning detectors, our approach uses a fast low-resolution and low-exposure prescan to acquire an image of the breast. Tube potential depends on breast thickness, and the prescan histogram provides the necessary information to calculate the required tube current. We validate our approach with simulated prescan images and phantom measurements. We achieve accurate exposure tracking with thickness and density, and expect this method of AEC to reduce retakes and improve workflow.

Image quality assessment via segmentation of breast lesion in X-ray and ultrasound phantom images from Fischer's full field digital mammography and ultrasound (FFDMUS) system.

Fischer has been developing a fused full-field digital mammography and ultrasound (FFDMUS) system funded by the National Institute of Health (NIH). In FFDMUS, two sets of acquisitions are performed: 2-D X-ray and 3-D ultrasound. The segmentation of acquired lesions in phantom images is important: (i) to assess the image quality of X-ray and ultrasound images; (ii) to register multi-modality images; and (iii) to establish an automatic lesion detection methodology to assist the radiologist. In this paper we developed lesion segmentation strategies for ultrasound and X-ray images acquired using FFDMUS. For ultrasound lesion segmentation, a signal-to-noise (SNR)-based method was adapted. For X-ray segmentation, we used gradient vector flow (GVF)-based deformable model. The performance of these segmentation algorithms was evaluated. We also performed partial volume correction (PVC) analysis on the segmentation of ultrasound images. For X-ray lesion segmentation, we also studied the effect of PDE smoothing on GVF's ability to segment the lesion. We conclude that ultrasound image qualities from FFDMUS and Hand-Held ultrasound (HHUS) are comparable. The mean percentage error with PVC was 4.56% (4.31%) and 6.63% (5.89%) for 5 mm lesion and 3 mm lesion respectively. The mean average error from the segmented X-ray images with PDE yielded an average error of 9.61%. We also tested our program on synthetic datasets. The system was developed for Linux workstation using C/C++.

Statistical image reconstruction for polyenergetic X-ray computed tomography.

This paper describes a statistical image reconstruction method for X-ray computed tomography (CT) that is based on a physical model that accounts for the polyenergetic X-ray source spectrum and the measurement nonlinearities caused by energy-dependent attenuation. We assume that the object consists of a given number of nonoverlapping materials, such as soft tissue and bone. The attenuation coefficient of each voxel is the product of its unknown density and a known energy-dependent mass attenuation coefficient. We formulate a penalized-likelihood function for this polyenergetic model and develop an ordered-subsets iterative algorithm for estimating the unknown densities in each voxel. The algorithm monotonically decreases the cost function at each iteration when one subset is used. Applying this method to simulated X-ray CT measurements of objects containing both bone and soft tissue yields images with significantly reduced beam hardening artifacts.