
Joint optimization of collimator and reconstruction parameters in x-ray fluorescence computed tomography using analytical point spread function and model observer
Hsin Wu Tseng
Srinivasan Vedantham
Sang Hyun Cho
Andrew Karellas
(andrewkarellas@radiology.arizona.edu)
Issue date 2020 Sep.
Abstract
Objective:
To jointly optimize collimator design and image reconstruction algorithm in X-ray Fluorescence Computed Tomography (XFCT) for imaging low concentrations of high atomic number (Z) elements in small animal models.
Methods:
Single pinhole (SPH) collimator and three types of multi-pinhole (MPH) collimators were evaluated. MPH collimators with 5, 7, and 9 pinholes using lead, stainless steel and brass were considered. A digital cylindrical phantom (0.5 mm3 voxels) of 25 mm diameter and 25 mm height with a central 5 mm diameter and 12.5 mm height cylindrical insert containing gold nanoparticles (2:1 insert: background concentration) was modeled. A 125 kVp, 2 mm Sn filtered x-ray spectrum (0.5 cGy/projection) for gold K-shell XFCT was considered. System matrices were implemented using analytical point spread functions (PSF) for each pinhole collimator. Poisson noise was added to the projection data (16 equiangular views) before image reconstruction using Maximum-Likelihood Expectation-Maximization (ML-EM) algorithm. Signal-present and signal-absent images were generated for the detection task performed by a channelized Hotelling observer (CHO) with 10 Dense Difference-of-Gaussian channels. The area under the curve (AUC) in receiver operating characteristic was used as the image quality metric.
Results:
A stainless steel focusing type MPH with 7 pinholes and 20 iterations of ML-EM provided the highest AUC.
Conclusion:
MPH collimators outperformed SPH collimators for XFCT and consistently high AUCs were observed with focusing type MPH designs with 7 pinholes.
Significance:
The combinations of collimator design and image reconstruction parameters that maximized AUC were identified, which could improve the performance of XFCT.
Index Terms—: X-ray fluorescence (XRF), X-ray fluorescence computed tomography (XFCT), collimator, multi-pinhole collimator, model observer, channelized Hotelling observer, system optimization
I. INTRODUCTION
Computed tomography (CT) is a common medical imaging modality which provides three-dimensional (3D) anatomical information and is routinely used for diagnosis, treatment planning and for monitoring of therapy response. CT also plays a major role in preclinical and small animal imaging. CT relies on measuring the differences in x-ray linear attenuation coefficients, but can be challenging to detect very low concentrations of high atomic number (Z) elements [1], [2]. In contrast, X-ray Fluorescence Computed Tomography (XFCT) relies upon detecting characteristic x-rays, such as those from K-shell or L-shell [3]–[6], enabling imaging and quantification of very low concentrations of high Z elements. Among imaging probes based on high-Z elements, gold nanoparticles (GNPs) conjugated with exogenous moieties targeting specific disease processes are of particular interest and can provide the molecular characteristics of such processes when properly imaged. Thus, the potential of XFCT forin vivo imaging of targeted high Z probes during pharmacologic and biodistribution studies using small animal models is being actively investigated.
Traditionally, synchrotron x-ray sources were used for XFCT [7], [8]. About a decade ago, Cheonget al (2010) built a benchtop XFCT system for imaging of GNPs using a polychromatic x-ray source [6]. Subsequently, Joneset al [5] experimentally demonstrated the feasibility to detect low concentrations of GNPs (0.5 wt. %) using a benchtop XFCT system with a cone-beam x-ray source. Kuanget al [9], [10] also built a benchtop XFCT system with a pencil-beam x-ray source, showing the feasibility of imaging GNPs as well as other high-Z probes [1]. In these systems, x-ray fluorescence (XRF) photons, emitted upon irradiation of the object with x-rays, can be detected by compact spectroscopic x-ray photon counting detectors. For each projection angle, the acquired x-ray spectrum is post-processed to extract the XRF signal, followed by tomographic image reconstruction [11].
XFCT simulation studies and experiments have traditionally used parallel hole collimators [5], [11]. However, miniaturizing parallel hole collimators for small animal imaging is difficult and degradation of the spatial resolution is a problem as the collimator entrance-to-object increases. Therefore, recent studies have focused on the use of pinhole collimators that promise improved image quality [12], [13]. Sasayaet al further proposed multi-pinhole (MPH) collimators for XFCT to increase the sensitivity and validated its performance [14]. The choice of collimator and the design of the pinhole influence the sensitivity and resolution of the system. Typically, collimators are optimized in the projection domain. However, task performance such as the detection task is operated and visualized by the human observer in the reconstructed image domain. Therefore, joint optimization, which optimizes collimator and reconstruction parameters simultaneously, has been proposed by McQuaid et al [15]. They have shown that the image reconstruction algorithm parameters should be considered in the optimization loop to avoid sub-optimal performance. Zhou et al also showed that joint optimization, as opposed to the conventional independent sequential optimization processes, yields improved performance on lesion-detection tasks [16]. For image reconstructions, Feldkamp-Davis-Kress (FDK) is commonly used due to its speed [17]. However, the FDK algorithm does not factor the noise in projection data and the incomplete tomographic information [18]. This could result in sub-optimal performance. Alternatively, good image quality can be obtained by maximum-likelihood (ML) reconstruction. Maximum likelihood expectation maximization (ML-EM) [19] and its accelerated version - ordered subset expectation maximization (OS-EM) [20] are widely used for image reconstruction in nuclear medicine. Since the image reconstruction for XFCT is analogous to tomographic nuclear medicine modalities, ML-EM and OS-EM are appropriate for XFCT. The system matrix information for both the forward projection and the back-projection is required for such iterative algorithms. System matrix elements can be directly measured [21]–[23] or simulated using Monte Carlo methods [24]–[26]. Direct measurement requires fabrication of the system inclusive of the collimator and hence this is beyond the scope of this study. Although the data from Monte Carlo simulations can be compressed due to their sparsity; it still is time consuming and requires storing the large system matrix, which is often not practical for optimization.
High resolution XFCT, primarily due to the availability of high-resolution, pixelated, spectroscopic 2-D detector arrays, is relatively new and enables improved resolution and sensitivity, and reduced acquisition time hitherto unavailable with older XRF imaging. However, with the use of 2D detector arrays, it is important to have an appropriately designed collimator. This study investigates the collimator configurations that are likely to improve performance. Considering the numerous potential configurations, it is important to gain insight as to which configuration(s) are likely to maximize performance, prior to building such a system. This could reduce the costs and the time associated with fabricating and testing each configuration. In this study, a joint optimization procedure was implemented for task-based optimization by using an analytical point spread function (PSF) to reduce the time for generating the system matrix and to avoid storing the large system matrix, and by using a model observer to investigate the sensitivity-resolution tradeoff with task-based optimization. Various pinhole (one single pinhole and three multi-pinhole with 5, 7 and 9 pinholes each) collimators, three pinhole septal materials (stainless steel, brass and lead) and three ML-EM reconstruction parameters (10, 15, and 20 iterations) were simulated. Results for single pinhole and multi-pinhole imaging under different conditions were quantitatively evaluated by the Channelized Hotelling Observer (CHO) with channels previously shown to be correlated to the human visual behavior [36], which also greatly reduced the expense and time associated with reader studies.
II. MATERIALS AND METHODS
A. System Matrix Model
Single or multi-pinhole XFCT systems typically use two opposing collimators coupled to the detectors to provide tomographic information, as shown in (Fig 1). For ML-EM (or OS-EM), discrete models are used to represent the unknown 3-D distribution of the GNPs within the object and the acquired XRF projection data. The relation between the 3-D distribution and the projections can be expressed as
(1) |
wherepi is theith element of the detector,fj is thejth component of the object. The system matrix,Mij comprises the number of XRF photons emitted at the locationj and detected by theith detector element. The system matrix,Mij must be accurate to avoid artifacts during the forward projection and the back-projection processes of the iterative loop during image reconstruction.
Fig. 1.
XFCT system has one detector for transmission imaging and two opposing collimators coupled to the detectors for stimulated emission imaging.
B. Point Spread Function
The pinhole collimator can be described as a double cone (Fig 2). The fraction of photons transmitting through the aperture contributes to the ‘geometric’ sensitivity of the collimator. Ideally, photons incident on the outside of the collimator aperture should be attenuated completely if the collimator is fabricated with highly-absorbing, high Z material. However, depending on the collimator material and the XRF photon energy, partial transmission could occur, and this contributes to the transmission sensitivity of the collimator. For the XRF photons emitted from a point in 3D space, the two-dimensional (2D) PSF of a pinhole collimator contains both the geometric collection and the transmission components. It should be noted that the total sensitivity including the geometric and the transmission components is obtained by the 2D integral of the PSF.
Fig. 2.
Three types of collimators and coordinate systems. (a) Red: NFRCDC, (b) green: FRCDC, and (c) blue: FOCDC based pinhole.γ is the focusing angle pointing to the center of VOI andα[NFRCDC],α[FRCDC], andαmin[FOCDC] are the acceptance angle for NFRCDC, FRCDC, and FOCDC, respectively. The origin of the (ua,va,ta) is the intersection of the uv-plane and the path of the penetrating photon.ρ is the radial distance between the origin of the (u, v, t) and the origin of the (ua,va,ta),θ is the azimuthal angle,φ is the polar angle of the point source, and h is the height of the point source in the (u, v, t) coordinate.θa andφa are the azimuthal and the polar angle of the point source in the (ua,va,ta) coordinate, respectively.
The geometric component of the PSF can be modeled from the shape and the dimensions of the pinhole aperture. However, modeling the penetrative part of the PSF is not trivial because the path length of photons through the collimator needs to be accurately calculated. This path length can be determined by the shape of the pinhole and the relative location between the point source (location of XRF emission) and the center of the pinhole. Pinholes are usually built in the shape of double cone with or without a tilt with respect to the central axis (Fig. 2a). The perpendicular line from the base of the right-circular cone to the vertex (the origin of the coordinate u, v, and t inFig. 2) is defined as ‘focusing axis’. The cone angle (2α) is called the ‘acceptance angle’ and the tilt angleγ is called the ‘focusing angle’. Bal and Acton [27] have derived closed form solutions for the PSFs of three different shapes of pinholes – non-focusing right-circular double-cone (NFRCDC), focusing right-circular double-cone (FRCDC), and focusing oblique-circular double-cone (FOCDC) based pinhole collimator. Among them, the NFRCDC is the special case of the FRCDC withγ = 0°. The equation describing the PSF of NFRCDC in Bal and Acton is consistent with the derivation of Metzler et al [28]. In this study, these three types of pinholes were selected for multi-pinhole collimator designs. The readers are referred to the derivations of the PSFs in Bal and Acton for further details. In order to be consistent with prior literature, we have used the same notation [25]. The final equation of the PSF of FRCDC and FOCDC based pinhole collimators are reproduced below for completeness [27]:
(2) |
whereθa is the polar angle of the point source in the (ua,va,ta) coordinate (Fig. 2),μ is the linear attenuation coefficient of the material atKα1 emission energy of gold (68.8 keV), and ΔL is the path length of the penetrating photon through the material medium. ΔL is zero when the photon passes through the aperture. The ΔL differs for FRCDC and FOCDC [27]. The ΔL for FRCDC is given by:
(3) |
wheredf is the equivalent diameter and is defined as. The other parameters inEqn (3) are:
(4) |
The ΔL of FOCDC is given by [27]:
(5) |
The parameters inEqn (5) are:
(6) |
Here the parameterαγ =αmin +γ andαmin used in FOCDC system is the same acceptance angleα used in FRCDC system (Fig. 2a).
The x-ray scatter and attenuation from the air was neglected because the point source is placed in the air. It has been shown that the scatter within the collimator is approximately 3–5.4% of the total detected photons for 500μm diameter lead pinhole [29]. The scattering within the collimator decreases further when the pinhole diameter increases (1.5 mm used in this study). Hence, the scatter within the collimator was ignored. This analytical form has been shown to be in agreement with the ray-tracing based Monte Carlo simulations with a mean square error of less than 1% [27]. Hence the system matrix elements can be determined with high accuracy using an analytical PSF by moving a point source within the volume of interest (VOI).
C. Pinhole Design and collimator materials
The collimator design parameters are summarized inTable I. Based on their x-ray attenuation characteristics, wide availability and cost, three different materials were studied – lead, brass, and stainless steel. The diameter,d, of NFRCDC and FOCDC-based pinholes were 1.5 mm. For FRCDC (Fig 3),d varied according to Eqn (13) in Bal and Acton [27]. In this study, one single pinhole (NFRCDC-based) and three types (NFRCDC, FRCDC, and FOCDC) of multi-pinhole with 5, 7, and 9 holes were considered. The central pinhole is of NFRCDC style for all multi-pinhole collimators, while the surrounding pinholes can be one of the three types (NFRCDC, FRCDC, and FOCDC). To simplify this study, all surrounding pinholes were considered to be of the same type for each MPH collimator design.
Table I.
Parameters for different pinhole types and pinhole numbers. It should be noted that the type of central hole is NFRCDC for both single and multi-pinhole systems. The unit of d is mm. The unit ofγ andα is degree.
Pinhole Type | NFRCDC | FRCDC | FOCDC | ||||||
---|---|---|---|---|---|---|---|---|---|
Pinhole number╲Parameters | d | γ | α | d | γ | α | d | γ | α |
Central hole | 1.5 | 0 | 64 | None | None | None | None | None | None |
Hole #1 - #4 | 1.5 | 0 | 64 | 1.49 | 5 | 59 | 1.5 | 5 | 59 |
Hole #5 - #8 | 1.5 | 0 | 64 | 1.48 | 10 | 54 | 1.5 | 10 | 54 |
Fig. 3.
Single pinhole and multi-pinhole systems. (a) Single pinhole. (b) Multi-pinhole: 5 holes. (c) Multi-pinhole: 7 holes. (d) Multi-pinhole: 9 holes. For all collimators, the central hole is NFRCDC with 1.5 mm diameter. For multi-pinhole collimators (b-d), the type of the surrounding holes can be NFRCDC, FRCDC, or FOCDC. For NFRCDC and FOCDC based pinholes, the diameter was 1.5 mm. The diameter of FRCDC based pinhole is 1.49 mm and 1.48 mm for holes #1 through #4 and for holes #5 through #8, respectively. The radial distance between the central hole and holes #1 through #4 is mm. The radial distance between the central hole and holes #5 through #8 is 5 mm.
The acceptance angle of NFRCDC-based pinhole is 64°. The focusing angleγ of FRCDC and FOCDC-based pinhole is 5° for pinhole number 1 to 4 (Fig. 3b; 5 pinholes) and is 10° for pinhole number 5 to 8 (Fig. 3(c–d); 7 and 9 holes). Hence, the acceptance angle of both FRCDC and FOCDC is 64° − 5° = 59° for all four surrounding pinholes of the 5-pinholes collimator and is 64° − 10° = 54° for both the 7- and the 9-pinholes collimator. The horizontal and vertical distances between the central pinhole and surrounding four pinholes (number 1 to 4) are 2.5 mm and are 5 mm for pinhole numbers 5 through 8. Finally, the number of projection views is 16 in full rotation (360°). Ideally, these parameters and even the arrangements of holes should be considered for system optimization. However, only the number of holes, the type of pinholes, the collimator materials, and the number of iterations in image reconstruction were considered in this study for simplicity. Examples of the PSF (5 holes) projected on the detector are shown onFig 4.
Fig. 4.
Examples of PSF (5 holes) projected on the detector. Left: PSF generated by the point source located at the corner of the volume of interest. Right: PSF generated by the point source located at the center of the volume of interest.
D. Phantom and system geometry
In this study, a cylindrical phantom of 25 mm diameter and 25 mm height with 0.5 mm voxel pitch (on each side) was simulated (Fig 5). A cylinder of 5 mm diameter and 12.5 mm height inserted centrally within this phantom was used as a signal for the task-based image quality assessment. The signal and background were considered to be a uniform mixture of GNP and water. The GNP concentration for the signal was varied from 0.01% to 0.2% weight fraction, with corresponding variation in the background GNP concentration to maintain a constant 2:1 signal-to-background GNP concentration for the entire study.
Fig. 5.
Dimensions of the cylindrical phantom. The height and diameter of the phantom are both 25 mm. The inserted cylindrical rod represents the signal. The height and the diameter of the rod is 12.5 and 5 mm, respectively. The background is uniform and the signal-to-background GNP concentration is 2:1. For simulation, the phantom is discretized to voxel pitch of 0.5 mm.
The simulated XFCT system comprises the transmission imaging provided by the cone-beam CT component and the tomographic XRF emission imaging system provided by the XFCT system. Theoretical and empirical studies for the transmission cone-beam CT imaging have been previously reported [30], [31]. Hence, the system geometry pertinent to XFCT alone is described. It is assumed that the phantom is co-aligned with the axis of rotation and the distances from the x-ray source to the axis of rotation is 150 mm. The distances from the center of the cylindrical phantom to the central pinhole and the XFCT detector are 30 mm and 90 mm, respectively.
E. Simulation of XRF emission
The number of XRF photons emitted over 4π,Nemitted can be approximated as:
(7) |
InEqn. (7), ω = 0.95 is the probability ofK-shell interaction with the GNP [32], η=0.762 [33] is theKα1 fluorescence yield of GNP,Nincident,j is the number of photons from the cone-beam x-ray source incident on voxelj of the phantom and is dependent on the dose per projection,J is the total number of voxels in the phantom,E is the energy of the incident photon,Kab is theK-edge absorption energy of gold (80.7 keV),EkV is the applied tube voltage (125 kV) ,qnorm is the incident x-ray spectrum normalized to unit area,μGNP is the linear attenuation coefficient of the mixture of gold and water for the specified GNP concentration, andt is the x-ray path length through voxelj that is dependent on the fan angle,βfan, and cone angle,γcone.
A 125 kVp x-ray spectrum from a tungsten anode x-ray tube with 0.8 mm beryllium inherent filtration and 2 mm of added tin filtration (Figure 6) was modeled usingSPEKTR 3.0 [34]. A relatively high energy spectrum is used with photon energies above theK-shell binding energy of gold. Previously reported data from Monte Carlo simulations for water-equivalent dose [30] was used to determineNincident,j for the radiation dose of 0.5cGy/projection. This computation provided the number of XRF photons emitted over 4π. The number of XRF photons reaching the surface of the collimator was determined by multiplying the solid angle (Ω) subtended by the collimator with respect to the phantom.
Fig. 6.
The 125 kV, tungsten anode x-ray spectrum with 0.8 mm inherent beryllium filtration and 2 mm added tin filter. The x-ray spectrum is normalized to unit area.
For each condition studied, projection data were generated by using the analytical PSF-based system matrix. Poisson noise with noise-level based on the selected GNPs concentration was added to the sinogram before image reconstruction processes. All simulations were implemented on MatLab (Mathworks, Natick, MA) using in-house developed CUDA (Nvidia Corporation, Santa Clara, CA) codes and processed on a workstation (T7810, Dell Technologies, Round Rock, TX) with a Graphic Processing Unit (GPU) card (GeForce GTX 1080, Nvidia Corporation, Santa Clara, CA).
F. Image Reconstruction
Although ML-EM and OS-EM were both appropriate for XFCT, ML-EM was selected for all studies. ML-EM can be described as [19]:
(8) |
wheref is the image,M is the system matrix, andp is the projection data. For acceleration of the image reconstruction processes, GPU-accelerated ML-EM algorithm was implemented.
G. Channelized Hotelling Observer
The goal of the system optimization is to provide the best possible task-based performance rather than optimizing the system for maximizing sensitivity or resolution. For this study, the binary detection task of detecting the signal (central rod in the phantom inFig 5) was selected. A numerical model observer was used for quantitative image quality assessment as it enables rapid computation of task-based figure of merit and avoids the need for time-consuming human observer studies. Observers that perform binary tasks such as the detection of a signal can be considered as mapping the image data to a decision variable. If the observer models are linear, then this mapping can be described mathematically as:
(9) |
wheret(r) is the observer’s discriminant function,f is the reconstruction image, superscript indexT means the transpose operator, andw is the linear image template. In general,t(r) andw(r) both depends upon the signal positon. In this study, the signal position is fixed at the center of the phantom; thus,t(r) andw(r) can be simplified ast andw, respectively. The Hotelling observer (HO) is a linear observer [35] that maximizes the signal-to-noise ratio (SNR) of the observer’s test statistics. The HO is difficult to implement due to the large dimensionality of the image data. The covariance matrices required for HO need to be estimated and inverted, and this task is computationally nontrivial. The Channelized Hotelling Observer (CHO) [36][37] is easy to implement because the dimensions of the covariance matrices are reduced by the selected channels. The CHO operates on the image with linear channels and then uses the channel outputs to generate test statistics. Usually, the dimension of the designed channel is much smaller than the dimension of the image data. Hence, the CHO mitigates the computational difficulties inherent in HO. In this study, we chose ten dense difference of Gaussian channels (D-DOG) proposed by Abbey and Barrett [38] with the same channel parameters. The CHO template can be calculated as:
(10) |
wherev1 andv0 stands for signal-present and signal-absent data sets, and represent the mean ofv1 andv0, respectively,Kv1 andKv2 are the covariance matrices estimated byv1 andv0, respectively. It has been shown that D-DOG channels have similar image template as the one estimated by the human observer without the assumptions pertaining to the noise texture in images [39]. Since the goal is to evaluate the trends in the AUC and not to provide a numerical match to the human observer AUC, the internal noise was not included. The signal-known-exactly/background-known-exactly (SKE/BKE) task considered here is appropriate as the location of the signal can be determined from the transmission CT.
H. Receiver operating characteristic and area under the curve
The reconstructed images were all independent and were divided into two groups – training and testing group. There were 400 reconstructed images in total, which were equally divided into two groups – 200 for training and 200 for testing group. Within each training and testing group of 200 images, they were separated into 100 signal present and 100 signal absent images (examples are shown onfigure 7.). Training images after the channelization step were used to calculate the image template based onEqn. 10 followed by applying this template on the channelized testing image to generate the test statistics usingEqn. 9. Two test statistics were generated; one from the signal present images and another from the signal absent images. Moving a decision threshold, a receiver-operating-characteristic (ROC) curve was generated.
Fig. 7.
Example of signal present (left) and signal absent (right) images. To display the signal clearly, lower noise levels were used in this example.
The area under the ROC curve (AUC) was used as the figure of merit. There are several ways to estimate the AUC variance. Bootstrap [40], jackknife [41], and shuffle [41] [42] are popular methods and their performance, such as bias, has been studied [43][44][45][46][47]. In this work, the AUC variance was estimated by a completely nonparametric and unbiased approach, the one-shot method [48]. The one-shot method is a variant of the common Dorfman, Berbaum, Metz [49][50][51] method with the main difference in that it is not a resampling dependent technique.
I. Analysis
The AUC is dependent upon the number of XRF photons detected. This implies that increasing the radiation dose to the tissue, increasing the GNP concentration and improving the collimator efficiency can all contribute to higher AUC values. Since we fixed the radiation dose per projection and the objective is to optimize the collimator type/design, the GNP concentration needs to be chosen from a range so that the AUC values for all cases can be clearly compared when the collimator type/design is varied. An initial set of simulations was conducted with SPH and 5-pinhole MPH lead collimators with the GNP signal concentration varied from 0.01% to 2% weight fraction (while maintaining a constant signal-to-background GNP ratio of 2:1).Figure 8 shows the results from this initial study. While the entire range of GNP concentrations were appropriate for SPH, GNP concentrations in the range of 0.01% to 0.05% were appropriate for MPH collimators. Since there were 30 different collimators (3 SPH; 9 NFRCDC; 9 FRCDC; 9 FOCDC) based on the combinations of the types of pinholes, the number of pinholes and collimator material, in conjunction with 3 settings for ML-EM iterations, there were a total of 90 conditions. Hence, we chose 0.03% GNP weight fraction for the signal and 0.015% weight fraction for the background to maintain a signal-to-background GNP concentration ratio of 2:1 in subsequent studies.
Fig. 8.
An initial study of AUC as a function of GNP concentrations were conducted with SPH (cross) and 5-pinhole MPH (circle) Pb collimators to determine the appropriate GNP concentration for subsequent studies. While the entire range of GNP concentrations were appropriate for SPH, GNP concentrations in the range of 0.01% to 0.05% were appropriate for MPH collimators.
All 90 conditions were investigated using the fixed GNP concentration for the signal and the background. The highest AUC among all conditions investigated is denoted asAUCmax. The uncertainty in theAUCmax estimate was obtained from the one-shot method mentioned earlier. For the sample size (200 images per condition with 100 each of signal present and signal absent images) and an assumed statistical power of 0.9, the lowest AUC value, denoted asAUCLB, that was statistically non-inferior (α = 0.05) toAUCmax was determined using exact binomial proportions (SAS version 9.4, SAS Institute Inc., Cary, NC). All conditions withAUC ≥AUCLB were analyzed for the frequency of collimator material, pinhole type, number of pinholes, and number of ML-EM iterations. These conditions provide guidance as to the combination of system design and image reconstruction parameters that are likely to provide consistently high AUC.
III. RESULTS
For all conditions studied, the MPH-based XFCT systems outperformed the SPH-based XFCT system, irrespective of the number of ML-EM iterations. The results are presented in terms of the collimator material as it is often preselected based on cost and ease of manufacturing. For each collimator material, the AUCs with 10 ML-EM iterations for different collimator designs are shown first. Then, the best performing collimator designs for that collimator material are selected to show the AUCs with varying number of ML-EM iterations for brevity.
A. Lead
For collimators made of lead,Fig. 9a shows the AUC for SPH and for MPH collimators with varying number of pinholes after 10 iterations of ML-EM for image reconstruction. For NFRCDC and FOCDC based MPH designs, the highest AUC values were obtained with 7 pinholes. The AUC for FRCDC based MPH increased with the number of pinholes. Overall, FOCDC MPH with 7 pinholes and FRCDC with 9 pinholes provided the highest AUC values (Fig. 9a) and hence their AUC with varying number of ML-EM iterations are shown inFig. 9b. The AUC for FOCDC with 7 pinholes showed a decreasing trend with an increasing number of ML-EM iterations, whereas the AUC for FRCDC with 9 pinholes did not show a noticeable variation (0.9791 ± 0.0082 for 10 iterations to 0.9858 ± 0.0080 for 20 iterations).
Fig. 9.
AUC for lead collimator material. Error bars represent ±1 standard deviation. (a) AUC vs. number of pinholes with 10 ML-EM iterations. Red: SPH; Cyan: NFRCDC based MPH; Green: FRCDC based MPH; Blue: FOCDC based MPH. Overall, the FOCDC MPH with 7 pinholes provided the highest AUC with comparable performance for FRCDC MPH with 9 pinholes. (b) AUC vs. number of ML-EM iterations. Blue: FOCDC MPH with 7 pinholes; Green: FRCDC MPH with 9 pinholes. FOCDC MPH with 7 pinholes showed a decreasing trend for AUC with an increased number of ML-EM iterations, whereas FRCDC with 9 pinholes did not show an appreciable trend in AUC.
B. Stainless Steel
For stainless steel material, FRCDC with 7 pinholes and NFRCDC with 7 pinholes provided similarly high AUCs (0.9770 ± 0.0084 and 0.9767 ± 0.0088, respectively) after 10 ML-EM iterations (Fig. 10a). For NFRCDC with 7 pinholes, there is a decreasing trend in AUC with an increasing number of ML-EM iterations (Fig 10b). For FRCDC with 7 pinholes, there was a marginal increase in AUC when the number of ML-EM iterations is increased. Among the conditions studied for stainless steel material, the highest AUC of 0.9909 ± 0.0043 is achieved with 20 ML-EM iterations for FRCDC with 7 pinholes (Fig. 10b).
Fig. 10.
AUC for stainless steel collimator material. Error bars represent ±1 standard deviation. (a) AUC vs. number of pinholes with 10 ML-EM iterations. Red: SPH; Cyan: NFRCDC based MPH; Green: FRCDC based MPH; Blue: FOCDC based MPH. (b) AUC vs. number of ML-EM iterations. Green: FRCDC with 7 pinholes; Cyan: NFRCDC with 7 pinholes. NRFCDC showed a decreasing trend in AUC with an increase in number of ML-EM iterations, whereas FRCDC with 7 pinholes showed a marginal increase and the highest AUC is achieved with 20 ML-EM iterations.
C. Brass
For brass collimator material, FRCDC with 9 pinholes offered the highest AUC (0.9791 ± 0.0082) with 10 ML-EM iterations (Fig. 11a) followed by NFRCDC with 7 pinholes (0.9436 ± 0.0155) and FOCDC with 7 pinholes (0.9435 ± 0.0156) (Fig. 11a). In this specific material, the highest AUC value of the 7 pinholes NFRCDC (0.9706 ± 0.0092) and the 7 pinholes FOCDC (0.9705 ± 0.0090) can be reached when with 15 iterations of ML-EM. The AUC value of the 9 pinholes FRCDC decreases greatly from 10 to 15 iterations and reaches the lowest value at 20 iterations (0.7012 ± 0.0370). The best possible image quality can be reached using the 9 pinholes FRCDC type MPH with 10 ML-EM iterations.
Fig. 11.
AUC for brass collimator material. Error bars represent ±1 standard deviation. (a) AUC vs. number of pinholes with 10 ML-EM iterations. Red: SPH; Cyan: NFRCDC based MPH; Green: FRCDC based MPH; Blue: FOCDC based MPH. (b) AUC vs. number of ML-EM iterations. Green: FRCDC with 9 pinholes; Cyan: NFRCDC with 7 pinholes. Blue: FOCDC with 7 pinholes. The 9-holes-FRCDC with 10 iterations provides the highest AUC value.
The AUCs using the model observer were obtained with 100 signal present cases and 100 signal absent cases, for a total of 200 cases. Among the 90 conditions investigated, the maximum AUC (AUCmax) ± standard deviation (SD) of 0.9909 +/− 0.0043 was achieved with stainless steel MPH collimator of FRCDC type with 7 pinholes and with 20 iterations of ML-EM. For the observedAUCmax, theAUCLB was 0.967. The MPH collimator designs along with number of ML-EM iterations that were statistically non-inferior toAUCmax are summarized inTable II. Analyzing the frequency of collimator material, pinhole type, number of pinholes and number of ML-EM iterations for all conditions withAUC ≥AUCLB, indicated that stainless steel material, FRCDC type, 7 pinholes MPH and 15 iterations of ML-EM, occurred most frequently. The aforementioned combination in terms of collimator parameters is identical to that withAUCmax.
Table II.
MPH designs (with image reconstruction parameters) that were statistically non-inferior to theAUCmax (i.e. serial number 1).
Material | Type | # of pinholes | # of ML-EM iterations | AUC Value | Standard Deviation | Serial Number |
---|---|---|---|---|---|---|
Stainless Steel | FRCDC | 7 | 20 | 0.9909 | 0.0043 | 1 |
Stainless Steel | FOCDC | 7 | 15 | 0.9906 | 0.0046 | 2 |
Lead | FOCDC | 7 | 10 | 0.9898 | 0.0044 | 3 |
Stainless Steel | FRCDC | 7 | 15 | 0.9877 | 0.0055 | 4 |
Stainless Steel | FOCDC | 7 | 20 | 0.9872 | 0.0055 | 5 |
Stainless Steel | FRCDC | 9 | 20 | 0.9863 | 0.0058 | 6 |
Lead | FRCDC | 9 | 20 | 0.9858 | 0.008 | 7 |
Brass | FRCDC | 7 | 20 | 0.9854 | 0.0056 | 8 |
Brass | FRCDC | 7 | 15 | 0.9833 | 0.0062 | 9 |
Lead | FRCDC | 9 | 10 | 0.9791 | 0.0082 | 10 |
Brass | FRCDC | 9 | 10 | 0.9791 | 0.0082 | 11 |
Stainless Steel | FRCDC | 9 | 15 | 0.9774 | 0.0076 | 12 |
Stainless Steel | FRCDC | 7 | 10 | 0.977 | 0.0084 | 13 |
Stainless Steel | NFRCDC | 7 | 10 | 0.9767 | 0.0088 | 14 |
Lead | FRCDC | 9 | 15 | 0.9765 | 0.0104 | 15 |
Brass | NFRCDC | 7 | 15 | 0.9706 | 0.0092 | 16 |
Brass | FOCDC | 7 | 15 | 0.9705 | 0.009 | 17 |
IV. DISCUSSION
Our study showed that XFCT system optimization is necessary for detecting low GNP concentration (Fig 9). For the detection task considered, MPH collimators outperformed SPH collimators for all conditions studied. Results indicated that the focusing-type pinholes - FRCDC and FOCDC, in general, provided higher AUCs for most conditions. In general, for the FOCDC type MPH, consistently higher AUCs were observed with 7 pinholes. For the FRCDC type MPH, higher AUCs were observed with 9 pinholes for lead and brass collimators and with 7 pinholes for stainless steel. In the case of the NFRCDC type MPH, higher AUCs were obtained when using 7 pinholes for all materials studied here.
Although lead is a common material for fabricating collimators, we demonstrated that for the x-ray energies appropriate forK-shell GNP imaging, stainless steel and brass are also suitable for XFCT. Depending on the manufacturing process (machining, casting, layering, vapor deposition) the suitability of alternatives to lead can potentially expedite manufacturing and reduce cost.
The framework provided in this study allows for optimizing the XFCT system parameters and image reconstruction parameters. Although the study was focused on collimator designs and ML-EM reconstruction, it can be extended for other aspects of the system design such as the detector and other reconstruction algorithms [52]. Importantly, this framework speeds up the process by using analytical PSF, GPU acceleration, and model observer for task-specific optimization.
Our study had limitations. X-ray scatter from the object was not modeled in this study. Although adding Poisson noise to the detected XRF emission images is appropriate due to the discrete characteristic x-ray energies, it would be of interest to study if similar results and trends can be observed when x-ray scatter is included with complex objects/phantom (i.e. varying and multiple signals and heterogeneous backgrounds), and for discrimination tasks. Also, our study assumed an ideal energy -resolved, spectroscopic, photon-counting detector; whereas in practice, pulse-height and charge-sharing contribute to reduction in both energy and spatial resolution. Arrangement and size of pinholes, modeling of x-ray scatter using Monte Carlo simulations and incorporating photon-counting detector characteristics are all subjects of future investigations.
V. CONCLUSIONS
A joint optimization process using analytical PSF and ML-EM reconstruction in conjunction with CHO was developed for XFCT system optimization. This method simplifies the data storage requirements for the large system matrix with the use of analytical PSF, reduces the computational difficulty associated with inverting the large covariance matrix, and reduces the time required to conduct human observer studies by using a CHO numerical observer model. For the detection task, the task-based imaging performance was quantified using the well-known AUC metric for both SPH and MPH collimators of varying design and material. The framework developed in this study can be extended to include other aspects of XFCT such as the arrangement and the size of pinholes and other reconstruction algorithms. Our study, for the phantom and the task considered, showed that MPH collimators outperformed SPH collimators for XFCT and consistently high AUCs were observed with the focusing type MPH designs with 7 pinholes.
ACKNOWLEDGMENT
This work was supported in part by grant R01 EB020658 from the National Institute of Biomedical Imaging and Bioengineering (NIBIB) of the National Institutes of Health (NIH). The contents are solely the responsibility of the authors and do not necessarily represent the official views of the NIH or the NIBIB.
Contributor Information
Hsin Wu Tseng, Medical Imaging, University of Arizona, Tucson, Arizona United States.
Srinivasan Vedantham, University of Arizona Arizona Health Sciences Center, Tucson, Arizona United States.
Sang Hyun Cho, Department of Radiation Physics, The University of Texas MD Anderson Cancer Center, Houston, Texas United States 77030.
Andrew Karellas, University of Arizona Arizona Health Sciences Center, Tucson, Arizona United States.
REFERENCES
- [1].Bazalova Met al. , “Investigation of X-ray fluorescence computed tomography (XFCT) and K-edge imaging,” IEEE Trans. Med. Imaging, vol. 31, no. 8, pp. 1620–1627, 2012. [DOI] [PubMed] [Google Scholar]
- [2].Feng Pet al. , “Analytic Comparison Between X-Ray Fluorescence CT and K-Edge CT,” IEEE Trans. Biomed. Eng, vol. 61, no. 3, pp. 975–985, 2014. [DOI] [PubMed] [Google Scholar]
- [3].Bazalova-Carter Met al. , “Experimental validation of L-shell x-ray fluorescence computed tomography imaging: phantom study,” J. Med. Imaging, vol. 2, no. 4, p. 043501, 2015. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [4].Ahmad Met al. , “X-ray luminescence and X-ray fluorescence computed tomography: New molecular imaging modalities,” IEEE Access, vol. 2, pp. 1051–1061, 2014. [Google Scholar]
- [5].Jones BLet al. , “Experimental demonstration of benchtop x-ray fluorescence computed tomography (XFCT) of gold nanoparticle-loaded objects using lead- and tin-filtered polychromatic cone-beams,” Phys. Med. Biol, vol. 57, no. 23, pp. 457–467, 2012. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [6].Cheong SKet al. , “X-ray fluorescence computed tomography (XFCT) imaging of gold nanoparticle-loaded objects using 110 kVp x-rays,” Phys. Med. Biol, vol. 55, no. 3, pp. 647–662, 2010. [DOI] [PubMed] [Google Scholar]
- [7].Cesareo R and Mascarenhas S, “A new tomographic device based on the detection of fluorescent x-rays,” Nucl. Instrum. Methods, vol. 277, pp. 669–672, 1989. [Google Scholar]
- [8].Rust GF and Weigelt J, “X-Ray Fluorescent Computer Tomography with Synchrotron Radiation,” IEEE Trans. Nucl. Sci, vol. 45, no. 1, pp. 75–88, 1998. [Google Scholar]
- [9].Kuang Yet al. , “First demonstration of multiplexed X-Ray Fluorescence Computed Tomography(XFCT) imaging,” IEEE Trans. Med. Imaging, vol. 32, no. 2, pp. 262–267, 2013. [DOI] [PubMed] [Google Scholar]
- [10].Kuang Yet al. , “Development of XFCT imaging strategy for monitoring the spatial distribution of platinum-based chemodrugs: Instrumentation and phantom validation,” Med. Phys, vol. 40, no. 3, pp. 1–7, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [11].Jones BL and Cho SH, “The feasibility of polychromatic cone-beam x-ray fluorescence computed tomography (XFCT) imaging of gold nanoparticle-loaded objects: A Monte Carlo study,” Phys. Med. Biol, vol. 56, no. 12, pp. 3719–3730, 2011. [DOI] [PubMed] [Google Scholar]
- [12].Li Let al. , “Full-field fan-beam x-ray fluorescence computed tomography with a conventional x-ray tube and photon-counting detectors for fast nanoparticle bioimaging,” Opt. Eng, vol. 56, no. 4, p. 043106, 2017. [Google Scholar]
- [13].Sasaya Tet al. , “Dual-energy fluorescent x-ray computed tomography system with a pinhole design: Use of K-edge discontinuity for scatter correction,” Sci. Rep, vol. 7, no. 1, p. 44143, 2017. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [14].Sasaya Tet al. , “Multi-pinhole fluorescent x-ray computed tomography for molecular imaging,” Sci. Rep, vol. 7, no. 1, p. 5742, 2017. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [15].McQuaid SJet al. , “Joint optimization of collimator and reconstruction parameters in SPECT imaging for lesion quantification,” Phys. Med. Biol, vol. 56, no. 21, pp. 6983–7000, 2011. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [16].Zhou Let al. , “Strategies to jointly optimize SPECT collimator and reconstruction parameters for a detection task,” IEEE Itnl. Symposium on Biomedical Imaging, pp. 394 – 39, 2009. [Google Scholar]
- [17].Feldkamp LA, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A, vol. 1, no. 6, pp. 612–619, 1984. [Google Scholar]
- [18].Tuy HK, “An inversion formula for cone beam reconstruction,” SIAM. J. Appl. Math, vol. 43(3), pp. 546–552, 1983. [Google Scholar]
- [19].Shepp LA and Vardi Y, “Maximum Likelihood Reconstruction for Emission Tomography,” IEEE Trans. Med. Imaging, vol. 1, no. 2, pp. 113–122, 1982. [DOI] [PubMed] [Google Scholar]
- [20].Hudson HM and Larkin RS, “Accelerated Image Reconstruction Using Ordered Subsets of Projection Data,” IEEE Trans. Med. Imaging, vol. 13, no. 4, pp. 601–609, 1994. [DOI] [PubMed] [Google Scholar]
- [21].Furenlid LRet al. , “FastSPECT II: A second-generation high-resolution dynamic SPECT imager,” IEEE Trans. Nucl. Sci, vol. 51, no. 3, pp. 631–635, 2004. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [22].Rowe RKet al. , “A Stationary Hemispherical SPECT Imager for three-dimensional brain imaging,” J. Nuclear Medicine, vol. 34, no. 3, pp. 474–480, 1993. [PubMed] [Google Scholar]
- [23].van der Have Fet al. , “System Calibration and Statistical Image Reconstruction for the U-SPECT-I Stationary Pinhole System,” IEEE Trans. Med. Imaging, vol. 27, no. 7, pp. 960–971, 2008. [DOI] [PubMed] [Google Scholar]
- [24].Rafecas Met al. , “Use of a monte carlo-based probability matrix for 3-D iterative reconstruction of MADPET-II data,” IEEE Trans. Nucl. Sci, vol. 51, no. 5, pp. 2597–2605, 2004. [Google Scholar]
- [25].Li X and Furenlid LR, “ A SPECT system simulator built on the SolidWorks TM 3D design package ,” Proc. SPIE, vol. 9214, pp. 92140G, 2014. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [26].Lazaro Det al. , “Fully 3D Monte Carlo reconstruction in SPECT: A feasibility study,” Phys. Med. Biol, vol. 50, no. 16, pp. 3739–3754, 2005. [DOI] [PubMed] [Google Scholar]
- [27].Bal G and Acton PD, “Analytical derivation of the point spread function for pinhole collimators,” Phys. Med. Biol, vol. 51, no. 19, pp. 4923–4950, 2006. [DOI] [PubMed] [Google Scholar]
- [28].Metzler SDet al. , “Analytic determination of pinhole collimator sensitivity with penetration,” IEEE Trans. Med. Imaging, vol. 20, no. 8, pp. 730–741, 2001. [DOI] [PubMed] [Google Scholar]
- [29].van der Have F and Beekman FJ, “Photon penetration and scatter in micro-pinhole imaging: A Monte Carlo investigation,” Phys. Med. Biol, vol. 49, no. 8, pp. 1369–1386, 2004. [DOI] [PubMed] [Google Scholar]
- [30].Vedantham Set al. , “A framework for optimizing micro-CT in dual-modality micro-CT/XFCT small-animal imaging system,” Proc. SPIE, vol. 10393, pp. 103930R1–7, 2017. [Google Scholar]
- [31].Deng Let al. , “Investigation of transmission computed tomography (CT) image quality and x-ray dose achievable from an experimental dual-mode benchtop x-ray fluorescence CT and transmission CT system,” J. Xray. Sci. Technol, vol. 27, no. 3, pp. 431–442, 2019. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [32].Thompson Aet al. , “X-ray Data Booklet: Center for X-ray Optics and Advanced Light Source,” 2001. [Online] Available:https://xdb.lbl.gov/xdb.pdf.
- [33].Chantler CT, “Detailed tabulation of atomic form factors, photoelectric absorption and scattering cross section, and mass attenuation coefficients in the vicinity of absorption edges in the soft x-ray (Z= 30–36, Z = 60–89, E = 0.1–10 keV)—addressing convergence issues of earlier work,” J. Synchrotron Radiat, vol. 8, pp. 1124 , 2001. [DOI] [PubMed] [Google Scholar]
- [34].Punnoose Jet al. , “Technical Note: spektr 3.0-A computational tool for x-ray spectrum modeling and analysis,” Med. Phys vol. 43, no. 8, pp. 4711–4717, 2016. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [35].Barrett HHet al. , “Stabilized estimates of Hotelling-observer detection performance in patient-structured noise,” Proc. SPIE, vol. 3340, pp. 27–43, 1998. [Google Scholar]
- [36].Myers KJ and Barrett HH, “Addition of a channel mechanism to the ideal-observer model,” J. Opt. Soc. Am. A, vol. 4, no. 12, p. 2447, 1987. [DOI] [PubMed] [Google Scholar]
- [37].Gallas BD and Barrett HH, “Validating the use of channels to estimate the ideal linear observer,” J. Opt. Soc. Am. A, vol. 20, no. 9, p. 1725, 2003. [DOI] [PubMed] [Google Scholar]
- [38].Abbey CK and Barrett HH, “Human- and model-observer performance in ramp-spectrum noise: effects of regularization and object variability,” J. Opt. Soc. Am. A, vol. 18, no. 3, p. 473, 2001. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [39].Abbey CKet al. , Jr., “Human-observer templates for detection of a simulated lesion in mammographic images,” Proc. SPIE, Med. Imag, vol. 4686, pp. 25–36, 2002. [Google Scholar]
- [40].Efron B and Tibshirani RJ, “The accuracy of a sample mean,”,in An introduction to the bootstrap: Mono-graphs on statistics and applied probability, Chapman and Hall, New York, NY, 1993. [Google Scholar]
- [41].Dorfman DDet al. , “Receiver operating characteristic rating analysis. Generalization to the population of readers and patients with the jackknife method.,” Invest. Radiol, vol. 27, no. 9, pp. 723–31, 1992. [PubMed] [Google Scholar]
- [42].Fukunaga K and Hayes RR, “Effects of sample size in classifier design,” IEEE Trans. Pattern Anal. Mach. Intell, vol. 11, no. 8, pp. 873–885, 1989. [Google Scholar]
- [43].Wagner RFet al. , “Finite-sample effects and resampling plans: applications to linear classifiers in computer-aided diagnosis,” Proc. SPIE, vol. 3034, pp. 467–477, 1997. [Google Scholar]
- [44].Wagner RFet al. , “Components of variance in ROC analysis of CADx classifier performance,” Proc. SPIE, vol. 3338, pp. 859–875, 1998. [Google Scholar]
- [45].Beiden SVet al. , “Components-of-variance models and multiple-bootstrap experiments: An alternative method for random-effects, receiver operating characteristic analysis,” Acad. Radiol, vol. 7, no. 5, pp. 341–349, 2000. [DOI] [PubMed] [Google Scholar]
- [46].Obuchowski NA and Rockette HE, “Hypothesis testing of diagnostic accuracy for multiple readers and multiple tests an anova approach with dependent observations,” Commun. Stat. - Simul. Comput, vol. 24, no. 2, pp. 285–308, 1995. [Google Scholar]
- [47].Obuchowski NA, “Multireader receiver operating characteristic studies: A comparison of study designs,” Acad. Radiol, vol. 2, no. 8, pp. 709–716, 1995. [DOI] [PubMed] [Google Scholar]
- [48].Gallas BD, “One-shot estimate of MRMC variance: AUC.,” Acad. Radiol, vol. 13, no. 3, pp. 353–62, 2006. [DOI] [PubMed] [Google Scholar]
- [49].Dorfman DDet al. , “Monte Carlo validation of a multireader method for receiver operating characteristic discrete rating data: Factorial experimental design,” Acad. Radiol, vol. 5, no. 9, pp. 591–602, 1998. [DOI] [PubMed] [Google Scholar]
- [50].Roe CA and Metz CE, “Dorfman-Berbaum-Metz method for statistical analysis of multireader, multimodality receiver operating characteristic data: Validation with computer simulation,” Acad. Radiol, vol. 4, no. 4, pp. 298–303, 1997. [DOI] [PubMed] [Google Scholar]
- [51].Roe CA and Metz CE, “Variance-component modeling in the analysis of receiver operating characteristic index estimates,” Acad. Radiol, vol. 4, no. 8, pp. 587–600, 1997. [DOI] [PubMed] [Google Scholar]
- [52].La Rivière PJ and Vargas PA, “Monotonic penalized-likelihood image reconstruction for X-ray fluorescence computed tomography,” IEEE Nucl. Sci. Symp. Conf. Rec, vol. 4, no. 9, pp. 1991–1995, 2005. [DOI] [PubMed] [Google Scholar]