research papers
Dualenergy Xray analysis using synchrotron computed tomography at 35 and 60 keV for the estimation of photon interaction coefficients describing attenuation and energy absorption
^{a}School of Physics, Monash University, Clayton, VIC 3080, Australia, and ^{b}Department of Radiation Therapy, University of Otago, Wellington, New Zealand
^{*}Correspondence email: stewart.midgley@monash.edu
A novel method for dualenergy Xray analysis (DEXA) is tested using measurements of the Xray μ. The key is a mathematical model that describes elemental cross sections using a polynomial in The model is combined with the mixture rule to describe μ for materials, using the same polynomial coefficients. Materials are characterized by their electron density N_{e} and statistical moments R_{k} describing their distribution of elements, analogous to the concept of effective In an experiment with materials of known density and composition, measurements of μ are written as a system of linear simultaneous equations, which is solved for the polynomial coefficients. DEXA itself involves computed tomography (CT) scans at two energies to provide a system of nonlinear simultaneous equations that are solved for N_{e} and the fourth statistical moment R_{4}. Results are presented for phantoms containing dilute salt solutions and for a biological specimen. The experiment identifies 1% systematic errors in the CT measurements, arising from thirdharmonic radiation, and 20–30% noise, which is reduced to 3–5% by preprocessing with the median filter and careful choice of reconstruction parameters. DEXA accuracy is quantified for the phantom as the mean absolute differences for N_{e} and R_{4}: 0.8% and 1.0% for soft tissue and 1.2% and 0.8% for bonelike samples, respectively. The DEXA results for the biological specimen are combined with model coefficients obtained from the tabulations to predict μ and the mass energy at energies of 10 keV to 20 MeV.
Keywords: linear attenuation coefficient; mass energy absorption coefficient; dualenergy Xray analysis.
1. Introduction
Synchrotron light sources (Boldeman & Einfeld, 2004) combined with wavelengthdispersive Xray optics, a samplerotation stage and an area detector enable monoenergetic Xray computed tomography (CT) studies with excellent spatial resolution (Stevenson et al., 2012). The ability to rapidly change beam energy allows materials analysis based upon the energy and compositional dependence of the Xray μ. Absorptionedge discontinuities can be exploited to quantify elemental concentration via methods such as Kedge subtraction (KES) (Zhu et al., 2014). However, the Kedges for elements in biological samples are below 5 keV, so KES is restricted to very small samples (e.g. excised specimens) or introduced contrast agents. At higher energies, μ varies smoothly (i.e. with continuous first derivative), decreasing at higher energies E and increasing for higher Z materials, according to the powerlaw relationship Z^{m}/E^{n} where the photoelectric effect dominates (with m = 3–4 and n = 2–3) and the Klein–Nishina formula where is important. Traditional methods of analysis utilize basis materials decomposition, which assumes that the sample is composed of only two or more materials. Measurements at different energies are written as linear simultaneous equations and solved for the contribution from each basis material.
We use a mathematical model describing the N_{e} and a nonlinear function of elemental composition (Jackson & Hawkes, 1981; Torikoshi et al., 2003; Midgley, 2004). For mixtures and compounds, the distribution of elements is described using statistical moments representing the mean its variance, skewness and kurtosis. These moments account for the compositional dependence for each interaction process; atomic form factors (Hubbell & Överbö, 1979) and factors (Hubbell et al., 1975) have weak contributions from all moments, whilst the photoelectric effect is dominated by the fourth statistical moment R_{4} which is related to the concept of effective (White, 1977; Jackson & Hawkes, 1981). §2.1 reviews the theoretical model describing μ for mixtures. Dual and multipleenergy Xray analysis (DEXA and MEXA) obtain measurements at two or more photon energies, which are written as simultaneous equations and then solved for N_{e} and R_{4}.
for materials proportional to electron densityThe energy ) is described by the same model. Both parameters μ and are required for radiation dosimetry calculations where μ accounts for attenuation and an estimate for the transferred to the medium (kerma) by the primary beam (Greening, 1981) is obtained via
(Seltzer, 1993In this expression, I_{t} is the transmitted beam fluence (number of photons per unit area) at depth in the medium with mass energy , and the summation is over the spectrum of beam energies. Calculation of absorbed dose must also account for scattered radiation arising from interactions at other locations and depths, whereby kerma is multiplied by the backscatter factor (BSF) estimated from measurements (Harrison, 1982) and via Monte Carlo methods. For diagnostic Xray energies (approximately 20–150 keV) interacting with tissues, the BSF is unity for a pencil beam and approximately 1.10 to 1.35 for fanbeam to widebeam conditions.
The aims of this study were to acquire CT data at diagnostic Xray energies (i.e. suitable for animal or human imaging), against which to test our DEXA model and its applications for radiation dosimetry, namely predicting μ and over a broad range of photon energies. §2 describes the DEXA methodology, phantom design and methods for characterizing attenuation measured at the Australian Synchrotron Imaging and Medical Beamline (IMBL). §2.6 reviews the necessary CT data processing steps and quality control procedures for beamline characterization. The CT results for the phantoms and a biological specimen are presented in §3, where we examine DEXA accuracy and predict interaction coefficients at photon energies of 10 keV to 20 MeV.
2. Materials and methods
The following sections review DEXA theory, describe the instrumentation for monoenergetic CT, the design of a phantom for characterizing measured attenuation, necessary quality control procedures for the experiment, CT data acquisition and processing methods.
2.1. Theory and mathematical methods for DEXA
The energy and compositional dependence of atomic cross sections σ can be described by nonlinear and nontrivial mathematical expressions that account for each interaction process. Jackson & Hawkes (1981) used a Taylor expansion to third order to rewrite these expressions as separate functions of Z and E. Torikoshi et al. (2003) recast the parameterization of Jackson & Hawkes as a function of electron density, with composition expressed as a function of an effective (which is denoted R_{4} below), plus further coefficients that represent photoelectric and scattering interactions. An alternative approach (Midgley, 2004) describes the compositional dependence of cross sections expressed per electron using a polynomial,
where coefficients S_{k} are a function of photon energy. The attenuation coefficients for mixtures and compounds are obtained by combining this expression with the mixture rule. The constituent elements have atomic volume density n_{a}(Z), and the material is characterized by the parameters
where N_{e} is the electron density. The compositional parameters R_{k} have the same units as analogous with the concept of effective (White, 1977), representing the statistical moments for variance, skewness and kurtosis. The model is written either as
or as the nonlinear expression
requiring four parameters for biological tissues at energies 30–150 keV, five coefficients at lower energies or for a broader compositional range, and fewer coefficients at higher energies. The parameters R_{k} are correlated (Midgley, 2011), so it is not possible to use equation (4) for multienergy analysis (MEXA) and recover all compositional parameters (Midgley, 2005). Instead we rewrite the model as a function of R_{4} attributed to the strong compositional contribution from photoelectric interactions, by introducing the ratios
which can be approximated by a constant or, for tissue, by a slowly varying function of R_{4} (Midgley, 2011). With materials now characterized by just two parameters, DEXA uses the following nonlinear model:
The parameterization given by equations (4), (5) or (7) can be written in matrix form . Prediction is the forward calculation of based upon knowledge of matrix and the vector of mixture parameters . The model coefficients S_{k} can be obtained by fitting the polynomial of equation (2) to the tabulations. We have used the Lawrence Livermore National Laboratory (LLNL) (Cullen et al., 1989; Boone & Chavez, 1996) and the National Institute of Standards (NIST) tabulations (Hubbell & Seltzer, 1995), with numerical results presented in the appendix of Midgley (2004). A complementary method uses μ measurements for different mixtures of known density and composition, to write linear simultaneous equations based upon equation (4). We use SVDFIT.C (Press et al., 1992) to solve for the coefficients S_{k}. DEXA and MEXA use knowledge of and measured to estimate the material parameters N_{e} and R_{4} by solving nonlinear equations based upon equation (7). For this purpose we use an algorithm (Midgley, 2013) based upon the Levenberg–Marquardt method (Press et al., 1992).
2.2. Calculation of the mass energy absorption coefficients for tissues
The compositional dependence of the energy N_{e} and R_{4} allowing to be calculated via equation (7). An estimate for the mass energy is obtained by dividing by the ρ related to N_{e} via
can also be described using a polynomial. DEXA deliverswhere the ratio of average ^{−1} for body tissues, and can be represented as a function of R_{4}:
to molecular weight MW is mol gwhere a_{0} = 0.8029, a_{1} = 0.01433, a_{2} = 1.567 and a_{3} = 3.645 for body tissues. Thus the mass energy is evaluated as
2.3. IMBL instrumentation for monoenergetic CT
The Australian Synchrotron storage ring is operated at 3 GeV and utilizes topup mode to compensate for the beam halflife (approximately 30 h) to maintain a current near 200 mA. The IMBL source is superconducting multipole wiggler radiation with a broad spectrum of energies extending beyond 100 keV and with a total radiative power of several tens of kW (Stevenson et al., 2010, 2012).
The experimental configuration is illustrated in Fig. 1. The beam was filtered by the invacuum attenuators located in hutch 1A (thickness 5 mm graphite and 0.5–3.5 mm Al; oriented at 45° to give longer path lengths), plus air and beryllium windows in hutches 2A and 3B. Beam energy selection used a dualcrystal bent Laue monochromator (DCBLM) where the mechanical bending serves to refocus the divergent beam. This configuration cannot maintain focus whilst one crystal is detuned, for example, to reject harmonic radiation. Beam energy is calculated via the Bragg equation, with the energy angle calibration being confirmed prior to the experiment by scanning across Kabsorption edges in a series of metal attenuators (Zr to Sb). The available decreases at higher energies, and the useable energy range is approximately 20–80 keV. The beam is a broad fan that was collimated to a beam area of approximately 100 mm wide by 5 mm high. At the time of the experiment, a mechanical arm controlling bending of the second crystal had become detached, thereby compromising the ability to focus the beam.
The experiment was conducted in hutch 3B at 146 m from the source using the Ruby detection system (Hall et al., 2013), which comprises a luminescent screen, optical coupling via mirror and lens (Nikon MicroNikkor 105 mm f2.8 macro) and a CMOS camera (PCO Edge, 2560 × 2160 pixels with native size 6.5 µm) that is located outside of the primary beam. We used the thickest available screen, a Gd_{2}O_{2}S:Tb phosphor taken from a mammography film cassette (Kodak MinR), with screen coating weight 31.7 mg cm^{2} or 43 µm thickness (Liaparinos et al., 2006). The detector was operated at nearmaximum pixel size of 45 µm and field of view 108 mm by 36 mm, and placed 100 cm behind the rotation stage, whereby the coarse angular resolution minimizes the influence of refraction contrast (Midgley, 2007) in the raw CT projections.
Data acquisition used the Australian Synchrotron Experimental Physics and Industrial Control System (EPICS) to rotate the sample, control the frame grabber and camera via a plugin for ImageJ, and save data as 16 bit TIFF files. Data acquisition rates are 50 fps (20 ms frame^{−1}) for the camera, 35 fps (29 ms frame^{−1}) for the frame grabber and 3 fps (330 ms frame^{−1}) for EPICS. At the time of this experiment, frame averaging was not supported. Thus each frame was 330 ms involving 20 ms (6%) camera time, 29 ms (9%) data transfer and 270 ms (83%) of unused beam time. CT data acquisition used continuous rotation over 180° with the beam fully exposed for the duration.
The incident beam was monitored by an ADC IC105 free air et al. (2013) with 60 mm path length, 11 mm electrode gap and approximately 150 mm width. In addition, a collimated silicon diode detector was monitoring scattered radiation arising from primary beam interactions with air in the experiment hutch.
(Advanced Design Consulting USA Inc., Lansing NY) as described by Crosbie2.4. Sample details
Two phantoms were assembled containing liquids of known density and composition. The phantoms comprised a polypropylene cylindrical container (internal diameter 62 mm, length 100 mm) with Perspex insert to support a hexagonal array of 19 polypropylene shell vials (outer diameter 8 mm, length 40 mm) containing aqueous ethanol and salt solutions. Sample details are listed in Table 1, with R_{4} spanning that of tissues from fat to adult cortical bone.

The aqueous solutions of varying concentrations were prepared by mixing a known mass of solute and solvent, and measuring their density with a 10 ml laboratory pycnometer. Uncertainties for the hydration state of the solute can introduce systematic errors to estimates for concentration. These were reduced by consulting published tabulations of ; Perry & Green, 2007). The phantoms are water based with similar dimensions to the biological sample to ensure similar scattering properties. Their electron density and compositional parameters were calculated as follows. Weight fractions were converted to an approximate chemical formula cf(Z) with the contribution from the least abundant element rounded to unity. The molecular weight is
against concentration (Nikolski, 1964The atomic density of each element in units of mol cm^{3} is
and the mixture parameters are calculated via equation (3).
An ex vivo biological specimen (Wistar rat) was also subject to CT for DEXA. The animal was bred according to local regulations and codes of practice governing biological research facilities. The carcass was retrieved from a batch of unallocated and thus surplus animals recently euthanized by submersion in carbon dioxide, sealed in a plastic ziplock bag, placed in a cylindrical container, stored in a freezer held at 193 K, and returned after the experiment.
2.5. Experimental methods
The samples were placed on a Huber rotation stage. The ex vivo specimen was mounted in a PVC tube glued to a plastic board, which was joined with screws to the base plate of the rotation stage. During data acquisition, a liquid nitrogen filled dewar directed a stream of cold dry gas towards the body, thereby preventing the buildup of water ice and slowing the process of thawing to room temperature.
Prior to each CT scan, the presence of beam harmonics was investigated using the Ruby detector to measure transmission through a copper step wedge. This important quality control check investigates whether the exponential attenuation law holds for a range of sample thicknesses (Creagh & Hubbell, 1990). At energies below 35 keV, we identified harmonics that led to significant beam hardening with measured μ decreasing with thickness.
Projection data sets for CT were acquired for the phantoms and three regions of the frozen Wistar rat. Before and after each CT scan, series of darkfield images were acquired with the beamline shutters closed, followed by series of flatfield images with the shutters opened and no sample in the beam. Data acquisition used a beam collimated to 74 mm by 2.4 mm and the camera operated to acquire just one frame per view, acquiring 1100 projections over 180°, in 360 s at two photon energies 35 and 60 keV.
2.6. Preprocessing corrections and CT reconstruction methods
The transmitted intensity signal recorded by the camera was subject to the following corrections: subtraction of the dark field, division by the flat field to remove spatial nonuniformities, and taking the logarithm to express the result as the raysum, representing the sum of μ and thickness t along the projection line through the sample. As the beam profile can change over time, flatfield images were collected before and after each CT scan, and their weighted average was used to estimate the flat field for each view.
Since the beam divergence across the detector was less than 0.1°, CT reconstruction utilized parallelbeam geometry filtered back projection (FBP) with the filtration step implemented as a realspace convolution. A variety of apodization functions (Ramachandran & Lakshminarayanan, 1971; Shepp & Logan, 1974; Webb, 1982, 1988) were available to attenuate the highfrequency noise that is otherwise amplified by the ramp filter during FBP. The backprojection step provided three interpolation options: nearest neighbour (NINT), linear (LINT) and spline (SPLINT) interpolation (Press et al., 1992). In §3.4 we compare results for SPLINT against LINT. Preprocessing and CT reconstruction were coded using the C programming language. Quantitative analysis of quality control measurements with the camera and the CT results was conducted using ImageJ (Rasband, 1997; Schneider et al., 2012).
3. Results
DEXA requires quality control measures to minimize the influence of systematic and random errors. To this end, we investigated the influence of harmonic radiation, and assessed camera performance and the temporal stability of the incident beam. We also estimated the radiation dose delivered during each CT scan. The present detection system requires data preprocessing to reduce noise. The model coefficients S_{k} in equations (2) to (7) are estimated from measurements of μ with the phantoms of known density and composition. DEXA accuracy was assessed using the same data set. In addition, the biological sample was subject to DEXA and the results used to predict μ and at other photon energies.
3.1. Estimation of harmonics content
Transmission measurements with the copper step wedge were expressed as raysums and fitted to a quadratic function with zero offset,
where monoenergetic μ for copper (Hubbell & Seltzer, 1995) is 63.0 cm^{1} (35 keV) and 14.2 cm^{1} (60 keV). Figs. 2(a) and 2(b) show the difference between measured raysums and the linear term in this model representing the relative error in forming the raysums. The systematic error is larger for thicker and higherZ materials.
The harmonic content was estimated as follows. We assumed a beam with the principal and third harmonic only, with relative intensities I_{0}(E_{0}) and I_{0}(3E_{0}). This model was fitted with NIST attenuation coefficients for copper against our raysum measurements. The results show that measured harmonic contamination relative to the principal is 1.0% at 35 keV and 1.5% at 60 keV. This information was used to calculate the relative error in estimating raysums (Midgley, 2006) for tissues as shown in Figs. 2(c) and 2(d).
3.2. Analysis of the camera dark signal and flat field
The camera dark signal was uniform and without structure, stable over time and independent of beam energy, with a mean pixel value of 100 and a noisetosignal ratio (NSR) of 3.2%, evaluated as the ratio of standard deviation to mean.
The incident beam was collimated to 1618×38 pixels (72.8 mm by 1.7 mm). Flatfield images did not contain bright or dark spots that can arise from diffraction by the beamline optics (Black & Long, 2004). The spatial profile across the flat field was brightest in the centre, falling by 10–30% at the edges, due to a mechanical fault that compromised the ability to focus the beam and the absence of a shading correction for the camera lens. For a region of interest (ROI) occupying the centre third of the flat field, the NSR was 7–12% at 35 keV and 4–8% at 60 keV.
The temporal stability of the flat field was quantified by subtracting images summed over 30 frames before and after each CT scan (6–10 min apart), and dividing by the mean. As the storage ring was operated in topup mode, it was not necessary to apply a decay correction. Subtracted flatfield images did not contain structure, indicating good spatial stability over the duration of a CT scan. Results are summarized in Table 2, where the average difference ranges from 0.5% to a few percent. Noise in the projections propagates into the CT reconstruction and is considered in detail in §3.4.

3.3. Beam dosimetry
The I_{ic} = dQ/dt in the irradiated volume containing the mass of air m_{air}. Assuming that the corrections for electron energy loss and ion recombination are unity, the air kerma rate for the incident primary beam is estimated via
measures the rate of charge productionwhere w_{air} = 34 eV is the mean for air (Greening, 1981). We use equation (1) expressed as the ratio of air kerma to beam fluence (averaged over the total beam area) to convert the measured air kerma rate to an estimate for the beam fluence rate.
The irradiated volume has a length of 60 mm, and the camera images show that the width and depth are 1635 ×55 rectangular pixels of size 45 µm. It contains a mass of air of 1.28×10^{5} kg based upon a dry air density of 1.205×10^{6} kg cm^{3} and assuming that the chamber correction factor for departures from standard temperature and pressure is unity. Table 3 summarizes the ionchamber measurements and their conversion to beam fluence. The air kerma per unit fluence was calculated for a monoenergetic beam via equation (1) using the for air from the NIST tabulation (Hubbell & Seltzer, 1995). The air kerma rate was estimated via equation (13) leading to an estimate for the incident fluence rate in SI units and expressed per pixel per frame (for 494 pixels per mm^{2} and 0.33 s per frame).

Each CT scan acquired 1100 views in 6 min, delivering an air kerma at the sample of 15 Gy at 35 keV and 3.5 Gy at 60 keV. For a uniform cylindrical sample of diameter 2R, the entrance skin air kerma (ESAK) is reduced due to sample rotation by a factor of , which is approximately 0.5 for the samples considered in this study. The beam is a narrow fan, with a (Harrison, 1982), so the absorbed dose to the skin is approximately 7 Gy at 35 keV and 2 Gy at 60 keV.
3.4. CT noisereduction strategies
The ratio of random errors for the raysum to those for the incidentbeam fluence I_{o} is given by propagation of errors analysis (Rose & Shapiro, 1948; Nördfors, 1960; Midgley, 2006) leading to the expression
where the detector records a background (or dark) signal I_{B} = B I_{o}. No background signal (B = 0) leads to a broad minimum, where the optimum thickness for CT is . Predicted results are shown in Fig. 3 with crosses denoting the average raysums for our samples. Our measurements were for radiologically thin samples, where the mean raysum is 0.6–1.3 and the error amplification factor is 1.8–3.6. The mean pixel value for the flat field is 28000 at 35 keV and 6400 at 60 keV, whilst the NSR (Table 4) for an average of 60 flat fields is approximately 10% at 35 keV and 5% at 60 keV. For FBP reconstruction with a tophat apodization filter and LINT, the reconstructed NSR is 20–30%. The following preprocessing strategies were employed to improve the reconstructed image quality.

The FBP reconstruction replaced the tophat apodization factor (Webb, 1988) with the Hamming window to suppress the influence of noise amplification by the ramp filter, improving the reconstructed NSR by a factor of 0.6. The interpolation process for back projection used SPLINT over LINT, further improving the reconstructed NSR by a factor of 0.8. We investigated two preprocessing strategies for filtering the raw data: pixel averaging by rebinning, and median filtering, with results presented in Table 4.
Rebinning to (N ×N) pixels (results not shown) improves the NSR by 1/N, but at the cost of reducing spatial resolution. The median filter was investigated using the algorithm SELECT.C of Press et al. (1992) with mask size (2M+1)^{2}. Results presented in Table 4 show that the NSR is reduced (i.e. improved) as 1/M. Mask sizes greater than M = 2 involve longer processing times with degraded spatial resolution. Therefore we used median filtering with M = 2 to preprocess all CT data sets, and improve the NSR by a factor of 0.5.
3.5. Beamline characterization for DEXA
Measured μ for materials in both phantoms illustrated in Fig. 4 were written as linear simultaneous equations according to equation (4) and solved for the model coefficients S_{k} using the methods outlined in §2.1. Predicted cross sections should increase with Z. This was not the case for , which is outside of the sampled compositional range. Therefore the measured μ data were supplemented with NIST values for hydrogen and helium with of unity. With this addition, the predicted cross sections were physically meaningful for the compositional range of interest hydrogen to calcium. Model coefficients S_{k} are presented in Table 5 and expressed in Fig. 5 as the atomic per electron. Fig. 6 shows the differences between measurements and μ predicted by the model, where the full range is % at 35 keV and % at 60 keV.

3.6. DEXA results
Measurements for the liquid samples were subject to DEXA using equation (7) to write nonlinear simultaneous equations that were solved for N_{e} and R_{4} using the model coefficients S_{k} obtained in §3.5 and the methods outlined in §2.1. The DEXA algorithm used an estimate for coefficients f_{k} based upon a parameterization of values for mixtures of lipid and water and for mixtures of water and compact bone (Midgley, 2013). Fig. 7 compares DEXA results for the phantoms against the known values for N_{e} and R_{4}.
The volumetric CT data sets for the biological sample were spatially coregistered using the Advanced Normalization Tools (ANTs) software package (Avantes et al., 2011), and underwent the DEXA on a pixelbypixel basis. The scanned volume was too narrow in the axial direction to provide enough landmarks for full threedimensional coregistration, so we used similar axial planes and a twodimensional affine transformation. The DEXA results were rejected when outside the range and , by assigning the floatingpoint value `not a number' (NAN). Results for the pelvis region are presented in Fig. 8.
3.7. Predicted interaction coefficients at other photon energies
In order to predict μ and for other beam energies, we used material parameters estimated via DEXA, the linear models of equations (4) and (10), and parameterization coefficients S_{k} obtained by fitting a polynomial to the tabulations (Midgley, 2004). The results for μ are presented in Fig. 9, and similar images are produced for with an example shown in Fig. 10(c). The images underwent quantitative analysis using circular ROI with diameters 2.5–8 mm. The DEXA results were expressed as the coordinates (N_{e}, R_{4}) as per Fig. 7 and are as follows: body of the spinal vertebrae (, ), muscle (, ) and peritoneal fat (, ). The measured tissue parameters are within the expected ranges for similar human tissues (for numerical values see the appendices of Midgley, 2011). Fig. 10 compares results for both μ and against values for similar human tissues evaluated using the mixture rule and NIST tabulation (Berger et al., 1990).
4. Discussion
The raw projection data contained a significant amount of noise, some minor artefacts and systematic errors arising from harmonic radiation. The camera design emphasizes robustness and cost over ) is reduced from 10–30% to approximately 3–5%.
with the electronics protected from the onset of radiation damage by being located outside of the primary beam; this comes with the cost of weak optical coupling efficiency. The random error arising from noise in the detection system is 20–30% per frame, and propagates into the CT reconstruction. We suppressed this noise contribution by preprocessing with a median filter, using FBP reconstruction with a Hamming apodization window to attenuate highfrequency noise that would otherwise be amplified by the ramp filter, and by means of cubic spline interpolation during back projection. The NSR in each reconstructed slice (Table 4CT reconstructions for the phantom exhibited weak ring artefacts and linear streaks (see the reconstructions at 35 keV in Fig. 4). The rings arise from small changes in the beam profile that were not fully removed by the flatfield correction. The causes were small temporal variations in the electronbeam current between topup injections, and spatial variations arising from mechanical vibration of the monochromator. Differences between flatfield images acquired before and after each scan, listed in Table 2, indicate the magnitude of these influences is of the order of 1%. Data acquisition involves continuous sample rotation with asynchronous data transfer over a network, with a small and variable time lag before writing to disk. The lag only becomes apparent after many scans and for periods of increased network traffic, leading to the angular sampling range being marginally less than the required 180° (Webb, 1988), which produces linear streaks in the reconstructions.
Transmission measurements with the copper step wedge (Figs. 2a and 2b) identified 1.0% (35 keV) and 1.5% (60 keV) contribution from thirdharmonic radiation. The camera falls at higher energies, reducing the influence of the third harmonic by 34% (105 keV/35 keV) and 9% (180 keV/60 keV), so the ratio for third harmonic to the fundamental wavelength is 5–10%. Figs. 2(c) and 2(d) quantify the systematic errors for tissues due to beam hardening by harmonic radiation. Propagation of errors analysis for random errors in raysum measurements is illustrated in Fig. 3 as the ratio of errors for the raysum to those for the incident beam. Both are small because our samples are radiographically thin with mean raysums of 0.7–1.4 at 35 keV and 0.5–0.9 at 60 keV. Thirdharmonic radiation leads to measured raysums being underestimated, and for our samples the systematic errors are approximately −0.6% for soft tissues and −1.2% for bonelike materials.
CT measurements of the phantoms were used to obtain coefficients S_{k} describing the compositional dependence of atomic cross sections (Fig. 5), including some forwardscattered radiation as measured by the instrumentation. Synthetic data for hydrogen and helium with of unity extend the compositional sampling below R_{4} = 6 to avoid nonphysical cross sections and thus ensure that the predicted values increase with Z. The goodness of fit is shown in Fig. 6, where the difference between measured and modelled μ is less than 1.2–2.4% with standard deviation 0.7–1.0%.
The same data was used to test the DEXA algorithm, with results presented in Fig. 7. Systematic and random errors in the μ measurements are shared amongst N_{e} and R_{4} in the DEXA solution such that one is underestimated and the other overestimated. Propagation of errors analysis (Midgley, 2013) shows how these errors are shared in the DEXA solution. The error for R_{4} is weighted by the fractional compositional crossproduct (representing the relative contribution of all R_{k} coefficients to μ) which is a function of beam energy and composition. At 35–60 keV the fractional compositional crossproduct is 0.20–0.05 for soft tissue and 0.7–0.4 for bone. In the present case the DEXA accuracy (, ) spans the range (%, −1.9% to +3.2%) for soft tissues and (%, −1.9% to +0.9%) for bonelike materials. The absolute fractional differences for (N_{e}, R_{4}) have average values (0.8%, 1.0%) for soft tissue and (1.2%, 0.8%) for the bonelike materials.
Fig. 8 presents CT results for the pelvis CT scan of the ex vivo sample showing the lower intestine, fat, muscle, the spine, skin and the ziplock bag. The intestines are filled with gas pockets and partially digested food, in particular Ridley AgriProducts mouse cubes, which consist of a mixture of grain, roughage and flakes of a mineral supplement containing 0.9% potassium and 1.2% calcium by weight. The DEXA maps for N_{e} and R_{4} are similar due to the strong correlations between these quantities for biological tissues. This was exploited via the model parameters f_{k}, allowing the linear model of equation (5) with five or more material parameters to be rewritten as the nonlinear model of equation (7), now a function of two variables. The DEXA results for N_{e} and R_{4} span the expected range for each parameter. Regions of black denote NAN, which indicates that the DEXA algorithm has failed to deliver meaningful results for regions of gas (where μ is small), at boundaries (due to partial volume effects), and where the artefacts (e.g. ring artefacts) are strong. The failure at boundaries is accentuated by using twodimensional spatial coregistration instead of full threedimensional methods.
Calculations for μ and are presented in Fig. 9 and in Fig. 10 against predictions based upon the NIST and LLNL tabulations for similar human tissues. The results show reasonable agreement between predictions for rodent tissue and expectations for nonadult human tissues. The comparison supports the methodology outlined in this paper in using DEXA to characterize the electron density and composition for tissues, and to predict tissue interaction coefficients at arbitrary photon energies.
Dosimetry estimates (§3.2) are summarized in Table 3. For a typical CT scan with 6 min of CT data acquisition, the radiation absorbed dose delivered to the sample is substantial. Achieving lowerdose CT would involve reducing the (e.g. by lowering the wiggler magnetic field and/or operating the monochromators on the shoulder region of the rocking curve) and utilizing a detection system with thicker phosphor and more efficient optical coupling. In addition, frame averaging may improve image quality, and the radiation dose might be further reduced by employing a fast shutter to block the beam while the system is busy with data transfer.
5. Conclusions
Necessary quality control steps for DEXA were outlined, and involve minimizing systematic and random errors in the raw data. It is important to test whether the exponential attenuation law holds (Creagh & Hubbell, 1990) by measuring transmission through a step wedge, whereby systematic errors become apparent when μ decreases with increasing attenuator thickness. Random errors arise from a variety of noise sources.
The IMBL with Ruby detection system was capable of CT at 20–80 keV with a nearparallel beam of size 110×2.4 mm. The experiment used 45 µm pixel size and acquired 1100 parallel projections over 180°, which were reconstructed via FBP to 90 µm voxels. Systematic errors were identified in the projection data, which arose from the thirdharmonic radiation with relative intensities of 5–10%. Their contribution was reduced to 1–1.5% of the recorded signal by the decrease in the at higher energies for the Gd_{2}O_{2}S phosphor. Our samples are radiographically thin with mean raysums of 0.5–1.5; thus the systematic errors in the measured raysums are low, from −0.6% for soft tissues to −1.2% for bonelike tissues. Weak optical coupling is a feature of the present detector design. As a consequence, noise in the flatfield images is approximately 10% at 35 keV and 5% at 60 keV, with temporal variations (Table 2) of a few per cent over the course of each CT scan. The radiation dose for each CT scan was estimated in §3.4 based upon air kerma rates measured by the inbeam Strategies were suggested for reducing the beam intensity and dose rate for diagnostic imaging, and improving the detection system to capture more of the incident signal and to deliver better image quality.
Reconstruction by FBP to 90 µm voxels, using a tophat apodization window (i.e. a pure ramp filter) with back projection based on linear interpolation delivered significant noise, approximately 30–20% at 35–60 keV. The situation was improved by preprocessing with a (5×5) median filter, replacing the top hat with a Hamming filter, and using spline interpolation. The respective noise improvements are 0.2 ×0.6 ×0.8 = 0.1, reducing the reconstructed noise to 3–5%.
The system was characterized, using measurements with liquid samples of known density and composition to determine model parameters that describe the attenuation measured by the beamline. Here it is important to measure the and the results expressed as atomic cross sections at each beam energy. DEXA was tested against the same measurements, written as nonlinear simultaneous equations based upon equation (7), which were solved using a modified Levenberg–Marquardt algorithm. Results were summarized as the difference from true values (, ), approximately (0.8%, 1.0%) for softtissuelike samples and (1.2%, 0.8%) for bonelike samples. Torikoshi et al. (2003) conducted a similar synchrotron CT experiment with tissue substitute materials and liquids, to deliver N_{e} and R_{4} within a few per cent of expected values. The important difference between these complementary approaches is the methods for obtaining the model coefficients: from the tabulations (Torikoshi et al., 2003; Midgley, 2004) versus via a `calibration' experiment as used here.
of each sample, and to choose compositions that fully span and evenly sample those of tissues characterized by . For the salt solutions, the amount of hydration water in the solute can be unknown, so we used tabulations of density and concentration to check our concentration estimates. A suitable algorithm was identified for solving the linear simultaneous equations based upon equation (4)The applications for DEXA are tissue characterization and the calculation of photon interaction coefficients at other beam energies. The latter were explored using CT scans at the same energies for a frozen biological sample and parameterization coefficients S_{k} obtained from published tabulations. Results were presented as volumetric maps representing N_{e} and R_{4}, and furthermore these results were successfully used to predict maps of μ and at energies from 10 keV to 20 MeV. This information is required for attenuation corrections, and is also used to calculate absorbed dose via equation (1) for diagnostic imaging across all modalities and for radiation therapy.
Acknowledgements
This research was undertaken on the Imaging and Medical Beamline at the Australian Synchrotron, Victoria, Australia. Travel funding was provided by the New Zealand Synchrotron Group. Sample preparation used the laboratory facilities of Professor Peter Rogers within the Department of Obstetrics and Gynaecology, The University of Melbourne. The assistance of Andrew Stevenson and Chris Hall in driving the beamline optics and the CT data acquisition systems (during experiment AS133/IM/6962 19–21 September 2013) is gratefully acknowledged. Spatial coregistration was facilitated by Brad Moffat and Chris Steward of the Brain Imaging Laboratory, Department of Radiology, Melbourne University.
References
Avantes, B. B., Tustison, T. & Song, G. (2011). Advanced normalisations tools (ANTS). Report. Penn Image Computing and Science Laboratory, University of Pennsylvania, PA, USA. Google Scholar
Berger, M. et al. (1990). NIST Report NBSIR 87–3597. National Institute of Standards and Technology, Bethesda, MD, USA. Google Scholar
Black, D. & Long, G. (2004). NIST Report SP 960–10. National Institute of Standards and Technology, Bethesda, MD, USA. Google Scholar
Boldeman, J. & Einfeld, D. (2004). Nucl. Instrum. Methods Phys. Res. A, 521, 306–317. CrossRef CAS Google Scholar
Boone, J. & Chavez, A. (1996). Med. Phys. 23, 1997–2005. CrossRef CAS PubMed Google Scholar
Creagh, D. C. & Hubbell, J. H. (1990). Acta Cryst. A46, 402–408. CrossRef CAS Web of Science IUCr Journals Google Scholar
Crosbie, J., Rogers, P. A. W., Stevenson, A. W., Hall, C. J., Lye, J. E., Nordström, T., Midgley, S. M. & Lewis, R. A. (2013). Med. Phys. 40, 062103. CrossRef PubMed Google Scholar
Cullen, D. et al. (1989). Report UCRL 50400. US Department of Commerce, Springfield, VA, USA. Google Scholar
Greening, J. (1981). Fundamentals of Radiation Dosimetry. Bristol: Hilger. Google Scholar
Hall, C. et al. (2013). J. Instrum. 8, C06011. Google Scholar
Harrison, R. (1982). Phys. Med. Biol. 27, 1465–1474. CrossRef CAS PubMed Google Scholar
Hubbell, J. et al. (1975). J. Phys. Chem. Ref. Data, 3, 417–538. Google Scholar
Hubbell, J. & Överbö, I. (1979). J. Phys. Chem. Ref. Data, 8, 69–105. CrossRef CAS Google Scholar
Hubbell, J. & Seltzer, S. (1995). Report NISTIR 5632. National Institute of Standards and Technology, Gaithersburg, MD, USA. Google Scholar
Jackson, D. & Hawkes, D. (1981). Phys. Rep. 70, 169–233. CrossRef CAS Google Scholar
Liaparinos, P. F., Kandarakis, I. S., Cavouras, D. A., Delis, H. B. & Panayiotakis, G. S. (2006). Med. Phys. 33, 4502–4514. Web of Science CrossRef PubMed CAS Google Scholar
Midgley, S. (2004). Phys. Med. Biol. 49, 307–325. CrossRef PubMed CAS Google Scholar
Midgley, S. (2005). Phys. Med. Biol. 50, 4139–4157. CrossRef PubMed CAS Google Scholar
Midgley, S. (2006). Radiat. Phys. Chem. 75, 936–944. CrossRef CAS Google Scholar
Midgley, S. (2007). Phys. Med. Biol. 52, 5173–5186. CrossRef PubMed CAS Google Scholar
Midgley, S. (2011). Phys. Med. Biol. 56, 2943–2962. CrossRef CAS PubMed Google Scholar
Midgley, S. (2013). Phys. Med. Biol. 58, 1185–1205. CrossRef CAS PubMed Google Scholar
Nikolski, B. (1964). Chemists' Reference Book, Vol. III, 2nd ed. Moscow: Leningrad. Google Scholar
Nördfors, B. (1960). Ark. Fys. 18, 37–47. Google Scholar
Perry, R. & Green, D. (2007). Perry Chemical Engineers Handbook, 8th ed. New York: McGraw Hill. Google Scholar
Press, W. S. A. T. et al. (1992). Numerical Recipes in C. The Art of Scientific Computing. 2nd ed. Cambridge University Press. Google Scholar
Ramachandran, G. & Lakshminarayanan, A. (1971). Proc. Natl. Acad. Sci. USA, 68, 2236–2240. CrossRef CAS PubMed Google Scholar
Rasband, W. (1997). ImageJ. Report. US National Institute of Health, Bethesda, MD, USA. Google Scholar
Rose, R. & Shapiro, M. (1948). Phys. Rev. 74, 1853–1864. CrossRef Google Scholar
Schneider, C. A., Rasband, W. & Eliceiri, K. (2012). Nat. Methods, 9, 671–675. CrossRef CAS PubMed Google Scholar
Seltzer, S. (1993). Radiat. Res. 136, 147–170. CrossRef CAS PubMed Google Scholar
Shepp, L. & Logan, B. (1974). IEEE Trans. Nucl. Sci. 21, 21–43. CrossRef Web of Science Google Scholar
Stevenson, A. W., Hall, C. J., Mayo, S. C., Häusermann, D., Maksimenko, A., Gureyev, T. E., Nesterets, Y. I., Wilkins, S. W. & Lewis, R. A. (2012). J. Synchrotron Rad. 19, 728–750. Web of Science CrossRef CAS IUCr Journals Google Scholar
Stevenson, A. W., Mayo, S. C., Häusermann, D., Maksimenko, A., Garrett, R. F., Hall, C. J., Wilkins, S. W., Lewis, R. A. & Myers, D. E. (2010). J. Synchrotron Rad. 17, 75–80. Web of Science CrossRef CAS IUCr Journals Google Scholar
Torikoshi, M., Tsunoo, T., Sasaki, M., Endo, M., Noda, Y., Ohno, Y., Kohno, T., Hyodo, K., Uesugi, K. & Yagi, N. (2003). Phys. Med. Biol. 48, 673–685. Web of Science CrossRef PubMed Google Scholar
Webb, S. (1982). Phys. Med. Biol. 27, 419–423. CrossRef CAS PubMed Google Scholar
Webb, S. (1988). The Physics of Medical Imaging. Bristol: Adam Hilger. Google Scholar
White, D. (1977). Phys. Med. Biol. 22, 219–228. CrossRef CAS PubMed Google Scholar
Zhu, Y., Samadi, N., Martinson, M., Bassey, B., Wei, Z., Belev, G. & Chapman, D. (2014). Phys. Med. Biol. 59, 2485–2503. Web of Science CrossRef CAS PubMed Google Scholar
© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.