Technical Note: spektr 3.0—A computational tool for x-ray spectrum modeling and analysis

Purpose: A computational toolkit (spektr 3.0) has been developed to calculate x-ray spectra based on the tungsten anode spectral model using interpolating cubic splines (TASMICS) algorithm, updating previous work based on the tungsten anode spectral model using interpolating polynomials (TASMIP) spectral model. The toolkit includes a matlab (The Mathworks, Natick, MA) function library and improved user interface (UI) along with an optimization algorithm to match calculated beam quality with measurements. Methods: The spektr code generates x-ray spectra (photons/mm2/mAs at 100 cm from the source) using TASMICS as default (with TASMIP as an option) in 1 keV energy bins over beam energies 20–150 kV, extensible to 640 kV using the TASMICS spectra. An optimization tool was implemented to compute the added filtration (Al and W) that provides a best match between calculated and measured x-ray tube output (mGy/mAs or mR/mAs) for individual x-ray tubes that may differ from that assumed in TASMICS or TASMIP and to account for factors such as anode angle. Results: The median percent difference in photon counts for a TASMICS and TASMIP spectrum was 4.15% for tube potentials in the range 30–140 kV with the largest percentage difference arising in the low and high energy bins due to measurement errors in the empirically based TASMIP model and inaccurate polynomial fitting. The optimization tool reported a close agreement between measured and calculated spectra with a Pearson coefficient of 0.98. Conclusions: The computational toolkit, spektr, has been updated to version 3.0, validated against measurements and existing models, and made available as open source code. Video tutorials for the spektr function library, UI, and optimization tool are available.

In light of this improved spectral model, this Technical Note reports an update to the  toolkit with TASMICS as the default method for spectral calculation. The code (referred to as  3.0) also includes a new optimization tool to assist with a common issue faced by TASMIP/TASMICS/ users: how to model and match the exposure characteristics of a particular x-ray tube that differs from that in the underlying TASMICS simulation. The code was developed and validated using  release 2013b. The  3.0/TASMICS implementation is detailed below, and the code is freely available for download at: http://istar.jhu.edu/downloads/. Video tutorials for the  function library, GUI, and optimization tool are available at the same link and at the following YouTube links:

2.A. Implementation
The x-ray fluence (photons/mm 2 /mAs at 100 cm from the source) in each energy bin was drawn from the data of Hernandez and Boone 56 and stored in spektrTASMICSdata.mat for beam energies 20-150 kV. The energy-dependent attenuation coefficients in spektrMuRhoElements.mat and spektrMuRho-Compounds.mat were updated to correspond to the average energy in each energy bin for a TASMICS spectrum (1.5, 2.5 keV, etc.) using values from the NIST XCom database and a cubic interpolation of NIST attenuation coefficients for selected compounds, respectively. 1,58 All Microsoft Excel dependencies from the original  release 1 have been removed in  3.0. The basic spectrum calculation is computed with the function call, where kV is the tube potential (kV), mmAl is the added Al filtration (mm), and kV_ripple is the kV ripple (%). For backward compatibility, spectral_model can be set to "TASMIP" (default = "TASMICS") to generate the original  2.0 spectrum. The default normalize flag (=1) normalizes the calculated tube output (mGy/mAs at 100 cm from the source) to match that of a spectrum calculated by  2.0 (TASMIP) at the same kV. Setting the normalize flag to 0 generates a TASMICS spectrum that is not normalized to match the  2.0 tube output. For example, a function call without optional arguments (i.e., kV argument only) generates a 70 kV spectrum using the TASMICS model, with 1.6 mm Al inherent filtration, 0% kV ripple, and normalization of mGy/mAs to match that of the TASMIP/Fewell spectrum. As described below, the choice of 1.6 mm Al filtration matches the inherent filtration of TASMIP. To generate a spectrum equivalent to Eq. (1b) using the TASMIP model, spectral_model is set to TASMIP as in Note that the output of spektrSpectrum() when spectral_ model is set to 'TASMIP' is independent of the state of the normalize parameter. Alternatively, the function call computes a 70 kV TASMICS spectrum with no added filtration, 0% kV ripple, and no normalization of tube output. Although TASMICS is defined for potentials across the orthovoltage range,  3.0 calculations are currently capped at 150 kV for backward compatibility. The code can be extended to 640 kV by adjusting the spektrTASMICSdata.mat file and modifying the functions in the  library to operate on spectra of length [1:640]. A glossary of the main  functions is given in Table I along with new functions introduced in version 3.0.

2.B. Validation
Spectra computed using  3.0 (TASMICS) were compared to  2.0 (TASMIP) calculations. As shown in Fig. 1(a), spectra computed with normalize set to 1 are in close agreement. Similar agreement was demonstrated for spectra computed over the range 30-140 kV (not shown; Pearson's R 2 coefficient >0.93 for all cases). The level of agreement is more fully quantified in Fig. 1(b), which plots the difference in photon fluence (x-rays/mm 2 /mAs at 100 cm from the source) in each energy bin for spektrSpectrum(kV) computed using  2.0 and 3.0, pooling calculations over 30-140 kV in 5 kV intervals. A median discrepancy of 4.15% was observed over the range 10-150 keV, slightly higher than that reported by Hernandez and Boone. 56 The slight discrepancy between  2.0 and  3.0 calculations is attributed to differences in fitting within each energy bin: TASMIP ( 2.0) uses a polynomial to fit photon fluence within a given bin, whereas TASMICS ( 3.0) uses a localized, piecewise third-order fit to the photon fluence in an energy bin.
Larger variations in the low-energy bins (10-15 keV) are likely due to measurement errors (presumably due to challenges in measurement of x-ray spectra, such as electronic noise, charge pile up, and x-ray scatter) and/or the resulting polynomial coefficients on which  2.0 is based.  3.0 is based on the MC simulation 57 underlying TASMICS and does not depend on such factors. Variation in the 60 keV bin is due to the increased energy resolution in  3.0. Specifically, the 2 keV energy resolution Fewell et al. data, upon which  2.0 is based, sums both the K α1 (57.98 keV) and K α2 (59.32 keV) tungsten edges into the 60 keV bin; 55  3.0, on the other hand, better resolves the characteristic radiation in 1 keV bins. 20 Similarly for the highenergy bins, while  2.0 likely suffers from inaccurate polynomial fitting in the 135 keV energy bin, 54  3.0 uses a spectral model defined up to 640 kV (Ref. 20) and produces a more accurate fit. Although the discrepancy in the high-energy bin appears large in terms of percent difference, the error is small in terms of the absolute fluence.

 function Description
spektr Launch  graphical UI airKerma = spektrAirKerma(q) Calculate the mGy/mAs for spectrum q at 100 cm from the focal spot fluencePerAirKerma = spektrFluencePerAirKerma(q) Calculate the fluence per air kerma for spectrum q at 100 cm from the focal spot X = spektrExposure(q) Compute mR/mAs for spectrum q at 100 cm from the focal spot Compute the fluence per exposure for spectrum q at 100 cm from the focal spot [mu, rho] = spektrMuRhoElement(Z) Return µ(E) and ρ for element Z in units of mm −1 and g/cm 3  Compute the filtration which provides best match of calculated and measured x-ray tube output and TASMIP/ 2.0 differ in their basic spectra and tube output characteristics, e.g., the mGy/mAs at 100 cm from the source. To account for such differences,  3.0 includes optional inherent filtration by 1.6 mm Al to match the inherent filtration of TASMIP. Normalizing the tube output (mGy/mAs at 100 cm from the source) of  3.0 calculations to those of  2.0 provides a close match between  2.0 and 3.0 spectra for beam energies 30-140 kV as shown in Fig. 1

3.B. Optimization ("tuning") of SPEKTR input parameters
An optimization tool called spektrTuner() was developed to assist in matching spectral calculations to measurements for a particular x-ray tube in terms of the output (mGy/mAs or mR/mAs). Similar to the process used by Sisniega, Desco, and Vaquero, 23 the algorithm performs an optimization to match calculated and measured outputs through variation of the thickness of Al and W filtration assumed in the calculation. Larger anode angles typically require higher filtration thickness due to the larger effective path of xrays produced at depth within the anode. The optimization algorithm and an illustration of the 2D search space are illustrated in Fig. 2 for in-air exposure measurements.
There are several means by which  calculations could be tuned to match the output of a particular x-ray tube, e.g., adjusting added filtration to match the HVL and/or mGy/mAs measured at various kV. The spektrTuner() function performs a best-fit between measured and calculated tube output as follows. The user first measures either in-air exposure (mR) or air kerma (mGy) at a particular kV and source-detector distance (SDD) at M settings of mAs, storing measurements in a [ As illustrated in Fig. 2(a), the function first approximates the measured beam via spectrum, a [150 × 1] vector with units of photons/mm 2 /mAs at 100 cm from the source generated using spektrSpectrum() at the specified measurement kV. Added filtration in the measured beam is accounted for via addFilt, where each row contains a filter element (atomic number) and its corresponding thickness (mm). The density of the filter material is assumed to be that reported by NIST for the element at standard temperature and pressure and can be retrieved via the function spektrMuRhoElement(). Measurements are normalized by mAs and scaled to SDD = 100 cm by the inverse-square law. The user provides an estimate of the inherent filtration, estimAl and estimW, as a [2 × 2] vector representing two material types [viz., Al (Z = 13) and W (Z = 74) in the first column and their respective thicknesses in the second (e.g., 2 and 0.01 mm)]. In the presence of multiple minima, the minimum closest to the estimates is returned. The tube output (mGy/mAs or mR/mAs) is computed and compared to the quotient of the measurements and mAs arguments. The sum-of-squared-difference is taken as a cost function that is minimized by adjusting the inherent Al and W filtration estimation via the simplex algorithm [computed using the  fminsearch() function]. Figure 2(b) shows an illustrative sweep across the search space of Al thickness (0.2-0.75 mm) and W thickness (0-0.0035 mm) and the optimal combination that minimizes the objective function. The spektrTuner() function returns a [2 × 2] vector containing the atomic number and thickness (mm) of Al and W that provides a best match to the measured tube output. The measurements in Fig. 2 were obtained on an xray imaging bench 59 incorporating an x-ray tube (RAD13, 16 • anode angle, and 0.4 mm focal spot size, Varian, Salt Lake City, UT) and silicon diode (Diagnostic Dose Diode, RadCal Corporation, Monrovia, CA) placed in air at SDD = 745 mm. Measurements were collected over the range 60-140 kV at six mAs stations each. A separate tuning of Al and W filtration was computed for each kV, showing close agreement between measurement and calculation (Pearson coefficient R 2 = 0.98) as shown in Fig. 2(c).

OTHER ENHANCEMENTS TO SPEKTR FUNCTIONS AND UI
Aside from the underlying TASMICS parameterization, a variety of enhancements and bug fixes have been implemented in  3.0. For example, the UI window (shown in Fig. 3) permits automatic resizing, the Spatial Filter interface has been removed, and all Excel file dependencies from the original release were eliminated. As illustrated in Fig. 3(A), the plotting tool includes standard tools for automatic axis scaling, pan, zoom, and a data cursor.
As illustrated in Fig. 3(B), the underlying spectral parameterization can be selected via radio button between TASMICS and TASMIP. The "Tube Select" drop-down menu has also been modified to automatically display the inherent filtration for any of the tube presets stored therein. The parameters for each tube can be adjusted within a new function, tubeSettings(), which keeps a "library" of tube presets with various filtration and kV ripple settings. A field was also added to the x-ray tube settings frame to allow addition of Cu filtration in the basic spectrum calculation.
The Added Filtration tool, shown in Fig. 3(C), was improved by modifying the underlying function, spektr-Beers(), to accept either the element atomic number (Z) or chemical symbol (as a string). Similarly for filtration by compounds, the function spektrBeersCompoundsNIST() was modified to accept either the compound index (C) listed in spektrCompoundList.m or the compound name (as a string). Also, a bug was corrected involving an error in the density of GaAs and Gd 2 O 2 S (5.31 and 7.44 gm/cm 3 , respectively) which were reversed in the original  release.
As shown in Fig. 3(D),  3.0 calculates a variety of basic metrics associated with a given spectrum, including exposure (mR/mAs at 100 cm from the focal spot), air kerma (mGy/mAs at 100 cm from the focal spot), HVL (as well as 2nd HVL, 3rd HVL for any element), conversion of the absolute spectrum (photons/mm 2 /mAs at 100 cm from the source) to a normalized probability density spectrum, fluence per unit exposure (x-rays/mm 2 /mR), fluence per unit air kerma (x-rays/mm 2 /mGy), and mean energy (keV).
Finally, as shown in Fig. 3(E), the ability to load previously computed spectra was updated with a simple Load button with folder browsing, and similarly for saving a spectrum via a simple Save button. A pushbutton was created to reset all fields in  to default values and clear variables from memory.

SUMMARY
With the development of the TASMICS algorithm 56 offering higher spectral resolution, broader energy range, and improved overall spectral characteristics with respect to modern x-ray tubes, this work presents an updated implementation of the  function library and UI for research in medical physics and x-ray imaging. A key improvement in this model is to avoid the errors associated with energy bin interpolation in the previous TASMIP and  tools.
Despite this improvement, slight differences can be expected between  3.0 calculations and measurements of x-ray spectra or tube output characteristics due to differences in anode angle. To help mitigate such differences, the  code scales fluence calculations for beam energies between 20 and 150 kV with the option to match the tube output (mGy/mAs at 100 cm from the x-ray source) computed by TASMICS to that in the previous  2.0 implementation (which in turn matches the measurements by Fewell et al. 55 ). The resulting  3.0 calculations match within a ∼4% error in mGy/mAs over the range 10-150 keV. This value is similar to the mean percent difference calculated by Hernandez and Boone (2.7%). 20 Moreover,  3.0 includes a new utility to help users match spectral calculations to measurements with a particular x-ray tube using the spektrTuner() optimization. This utility computes the thickness of Al and W filtration (positive or negative thickness) that minimizes the sum-of-squared difference between measured and calculated tube output (mGy/mAs or mR/mAs) using a simplex optimization. Measurements on an x-ray imaging bench demonstrated agreement with tuned  calculations with a Pearson Coefficient R 2 = 0.98 for beam energies ranging 60-140 kV. Alternative forms of spek-trTuner() could be developed to perform best match to other objective functions, for example, HVL. Other modifications improved the functionality of the UI, including better display of inherent filtration for various x-ray tube presets, removal of Excel dependencies, simplification of the input parameters in the spektrBeers() and spektrBeersCompoundNIST(), and automatic resizing of the UI.
These enhancements update  to provide a  interface to the TASMICS x-ray spectrum parameterization and will hopefully be of use to researchers in x-ray spectral analysis, image quality modeling, MC simulations, polyenergetic image reconstruction algorithms, and other areas of research for systems in the diagnostic x-ray energy range.

ACKNOWLEDGMENTS
This work was supported by NIH Grant No. R01-CA-112163. Thanks to Dr. J. W. Stayman (Biomedical Engineering, Johns Hopkins University) for assistance with x-ray bench measurements and useful discussion. Thanks also to Dr. J. M. Boone (Radiology, University of California -Davis), Dr. K. Taguchi and Dr. M. Mahesh (Radiology, Johns Hopkins University), and Dr. E Samei (Radiology, Duke University) for valuable discussion.