Design of a Low-Dose, Stationary, Tomographic Molecular Breast Imaging System Using 3-D Position Sensitive CZT Detectors

Molecular breast imaging (MBI) has been shown to have high sensitivity for lesion detection, particularly in patients with dense breasts where conventional mammography is limited. However, relatively high radiation dose and long imaging time are limiting factors. Most current MBI systems are based on planar imaging. Improved performance can be achieved using tomographic techniques, which normally involve detector motion. Our goal is to develop a low-dose stationary tomographic MBI system with similar or better performance in terms of lesion detection compared to planar MBI. The proposed system utilizes two opposing cadmium zinc telluride detectors with high intrinsic resolution and depth of interaction (DOI) capability, combined with densely packed multipinhole collimators. This leads to improved efficiency and adequate angular sampling, but also to significant multiplexing (MX), which can result in artefacts. We have developed de-MX algorithms that take advantage of the DOI information. We have performed both analytic and Monte Carlo simulations to demonstrate the feasibility, optimize the design and investigate the expected performance of the proposed system. Lesion detectability was preserved with reduction of acquisition time (or radiation dose) by a factor of 2 compared to planar images at the lowest reported dose. The first prototype is under evaluation at Kromek.

breasts [1]. In the USA, between a third and half of screeningaged women are estimated to have breast tissue that is categorized as mammographically dense, which reduces the sensitivity of mammography for the detection of breast cancer [1]. Molecular breast imaging (MBI) has been shown to have much better sensitivity in this patient group. In recent clinical trials, sensitivity for lesion detection was shown to increase by more than a factor of three when combining mammography with MBI [2]. The widespread adoption of MBI has been limited due to a relatively high radiation dose and long imaging time [2]. With a 300-MBq injection of [ 99m Tc]sestamibi, the effective whole-body radiation dose from MBI is ∼five times that of mammography [2], while the dose to the breast itself is lower [3]. Various dose reduction techniques have been developed over the past decade, both in terms of hardware and software [4], [5], [6], [7].
Current MBI systems are equipped with parallel-hole collimators, producing planar (2-D) images, which can be directly interpreted without the need for post-processing. On the downside, they have relatively low contrast and, although some depth information can in principle be obtained [8], this is limited.
With a tomographic system, on the other hand, 3-D images can be generated. Dedicated tomographic MBI systems have been designed involving complex mechanics [9], [10]. The use of standard SPECT scanners equipped with specially designed multislit-slat collimators has been investigated [11]. Simpler, so called tomosynthesis systems, have also been proposed, in which data are acquired at a few different angles either by detector rotation [12] or using a variable angle slant-hole collimator [13], [14]. An alternative for tomographic data acquisition is to use multipinhole (MPH) collimators. Such systems have been suggested, but detector motion was needed to cover the required field-ofview (FOV) [15], [16], [17]. One downside of tomographic systems is that image reconstruction is required. However, once reconstructed, the images provide depth information, as well as higher contrast compared to planar images.
In a relatively recent review on organ specific imaging devices including breast imaging systems, González et al. [18] mentioned that the advantages of such devices include higher detection efficiency, improved spatial resolution, better image contrast recovery, and lower cost. Our goal was to design a low-dose stationary tomographic system which would be suitable for complementary breast screening. The system should have a detection performance similar to or better than current planar systems, but with significant dose reduction and reduced scanning time. Our proposed system is based on dual opposing pixelated cadmium zinc telluride (CZT) detectors with high-intrinsic spatial resolution and capability for depth-of-interaction (DOI) estimation. The system utilizes densely packed MPH collimators which result in significant multiplexing (MX), normally avoided in MPH systems. The 3-D reconstruction relies on use of de-MX algorithms, aided by the DOI information. As the system is stationary, complex mechanics are avoided, and it is suitable for breast cancer screening.
In this article, we describe the overall design concept, optimization, and performance compared with conventional planar MBI. We also present initial results obtained with a prototype system. Parts of this work have been presented at recent conferences [19], [20], [21], [22], [23].

A. System Design
The overall design concept is to use dual opposing highspatial resolution detectors with MPH collimators, applying mild breast compression in the same way as in conventional planar MBI systems (Fig. 1). By using detectors with high-intrinsic spatial resolution, detection efficiency can be increased using multiple pinholes with minification (i.e., magnification <1); a principle originally introduced for brain imaging [24].
For the purpose of tomographic imaging, appropriate angular sampling is required. Traditionally, this is obtained via detector motion. However, as we intended to develop a stationary system, we opted for the use of a large number of closely spaced pinholes instead (Fig. 1). This can provide adequate angular sampling, and also leads to increased detection efficiency. The MPH collimators are positioned close to the high-spatial resolution CZT detectors. This results in significant MX, due to overlap of projections from different pinholes, which could lead to artefacts due to ambiguity in the direction of incidence of the detected photons. For this reason, MX is usually avoided, either by increased pinhole spacing or decreased opening angle, leading to lower-detection efficiency, limited angular sampling, reduced FOV and/or the requirement for detector or subject motion.
Various approaches for dealing with MX have been proposed in the past, including combined mixed MX and non-MX data [25], [26], [27], [28], which could be achieved using timed shutters [29], or alternatively using stacked detectors [30], or multiple detector positions [31]. Our approach is to use DOI information to achieve varying degrees of MX, enabling the use of densely packed pinholes with wide opening angle. Fig. 1(b) illustrates how different DOI layers contain different amounts of MX.

B. Simulations
The simulation experiments presented below were performed for a 160 × 160 mm detector size with a pixel size of 1 × 1 mm and five 1-mm thick DOI layers. The MPH collimator consisted of a 5-mm thick tungsten plate with a square grid of square-shaped pinholes. The pinhole apertures were in the range 1-3 mm, and the opening angles in the range 70 • -110 • . Penetration of γ -photons through the collimator material was included in the simulations. The separation between the collimator and the detector was 3 mm.
A number of different phantoms were used, designed for different purposes. We have summarized the properties of all the phantoms in Table I. All the phantoms were 60-mm thick, corresponding to an average breast thickness under mild compression [4]. Fig. 2 shows sections from different phantoms.

C. Image Reconstruction
We have developed algorithms to de-multiplex the data, either during or before reconstruction [20]. In the first algorithm, MX is included in the system matrix. This is referred to as direct reconstruction or the 1-step approach. In the second algorithm, a de-MX procedure is applied to the projection data in order to estimate projection data for each pinhole without MX using an ML-EM algorithm. This procedure is independent of the tomographic reconstruction, which is performed subsequently by a standard algorithm for non-MX data. We refer to this as the 2-step approach. We also implemented a hybrid approach, in which an average of two correction factors is used, based on raw projection data and de-MX projection data, respectively. Our de-MX procedure differs from that presented in [32], where a nested algorithm was used, alternating between de-MX and reconstruction. Tomographic reconstruction was performed using an ordered subsets implementation of a MAP algorithm [33], with a prior designed to equalize the resolution in planes parallel to the detectors. Subsets were chosen based on DOI-layers as follows: {1}, {2 + 5}, and {3 + 4}, resulting in a similar number of counts in the three subsets. DOI layer {1} is closest to the collimator (Fig. 1). Data from both detectors were included in each subset. Noise correlation was not taken into account in the reconstruction. The de-MX and reconstruction methods are described in more detail in the Appendix.
In order to compare the different reconstruction approaches, we performed analytical simulations with a source distribution corresponding to a phantom containing one layer of spherical lesions in four segments (left, right, top, and bottom) with 36 spheres each. The sphere diameter was 6 mm, the smallest centre-to-centre spacing 12 mm, and the tumor-to-background ratios (TBRs) were 5, 10, 15, and 20 in the four segments, respectively. The 144 lesions were placed at a depth of 11 mm in the 60-mm thick phantom. The simulations were done for the different numbers of pinholes, separated by 10, 11, 12, 13, 14 15, 16, 18, and 20 mm, resulting in different amounts of MX. The pinhole apertures were 1.5 mm and the opening angles 90 • . In addition to the three approaches mentioned above, we also reconstructed images from ideal, non-MX data, representing the best possible theoretical scenario. In this case, the data corresponding to different pinholes were stored separately, so there was no ambiguity in the projection space.

D. Initial Feasibility Study
Monte Carlo (MC) simulations were performed using the GATE software [34] with the new camera design and also with that of an existing parallel-hole collimation system [35] for comparison. The simulation model [36] was comprised of three parts. The first one uses GATE to model photon propagation and interactions in CZT detectors which create the initial charge cloud distribution. The second part is based on COMSOL multiphysics software (COMSOL AB, Stockholm Sweden) to model the drift of the charge cloud taking into account Einstein diffusion, Coulomb repulsion, and charge trapping to calculate charge induction efficiency in all surrounding pixels using the method described by Prettyman [37]. The last part combines the previous two and adds Fano noise and energy resolution factor to account for the electronics noise. A wide energy window of 110-150 keV was applied to the simulated energy spectra to include the escape peaks due to Cd (23.2 keV) and Te (27.5 keV) x-ray fluorescence.
These simulations allowed us to evaluate the feasibility of the proposed design and to compare its performance with planar MBI characterization described in [35]. Two detector configurations were simulated for that purpose. The first one corresponded to the design of the commercial CZT MBI camera which has been characterized in this article. The camera has the field of view of 22 × 16 cm, with 5-mm thick detectors with 1.6-mm square pixels. The second one corresponded to our initial conceptual design. It was chosen to be close to the described planar MBI camera, but also to reflect initial detector configuration and performance estimates based on the existing detector data and geometry. As a result, the simulated detector array was comprised of 7.3-mm thick CZT detectors with 2-mm pixels. The pinhole apertures were 1.0 mm and the opening angles 75 • . In order to allow for direct comparison with published values, the simulated activity distribution corresponded to a physical phantom used previously for evaluation of planar MBI [35]. The phantom was 60-mm thick, containing 48 cylindrical lesions, 3-10 mm in diameter, with six different planar TBRs in the range 2.0-5.3. The lesion detection sensitivity was calculated as the number of lesions with a contrast-to-noise ratio (CNR) > 3 divided by the total number. The threshold for lesion visibility in terms of CNR is usually considered to be in the range 3-5. This is known as the Rose criterion [38].

E. Design Optimization
The system design parameter optimization was performed by means of analytical simulations. The parameters investigated here were pinhole size, opening angle and hole separation. Collimator-to-detector separation was fixed to a value obtained at an earlier stage (3 mm). Two different approaches were used in the optimization.
The first approach was based on maximization of the CNR. We define contrast as C = (A-B) A hot-lesion phantom was designed based on the planar phantom described in the previous section. The thickness was the same (60 mm). However, here we used spherical instead of cylindrical lesions, as our purpose was 3-D imaging. These were placed at a depth of 15 mm. The TBRs were calculated so that the contrasts in a planar projection would be similar, resulting in TBRs in the range 10-130. We excluded a number of lesions with large diameter and/or TBR, as well as three lesions with small diameter and/or TBR, leaving a total of 18 spheres. About half of the image plane was left empty for estimation of the background noise. Projection data were generated by analytical forward projection and Poisson noise was added corresponding to the number of counts acquired over 5 min, after injection of 150 MBq of [ 99m Tc]-MIBI. Based on planar background count-density and detector efficiency values presented by Long et al. [35], we estimated that the background activity concentration corresponding to a 150-MBq injection would be 760 Bq/mL. Images were reconstructed with 1-5 iterations and three subsets. Noise estimates were obtained as the COV in a uniform region of the phantom, and the average CNR over all 18 lesions was used as the outcome metric. For details see [23].
The second optimization approach used was based on a numerical observer study using the scan statistics model [39]. The algorithm searched for lesions of a particular size, defined a priori. We therefore simulated data for a phantom containing 20 spherical lesions, all with a diameter of 5 mm and a TBR of 10, placed at random locations in the phantom. Data were generated and reconstructed as described above. The area under the localization receiver operating characteristics curve (LROC-AUC) was the outcome measure used. Optimization was performed in terms of the parameters mentioned above as well as the number of iterations.

F. Final Performance Evaluation
For the purpose of evaluating the performance of the final system design, compared to a conventional planar system, analytical simulations were performed using the optimized geometry and the hybrid reconstruction approach.
Here, we used a phantom with a more realistic range of spherical lesion sizes (4-8 mm) and TBR values (5-25). The phantom thickness was 60 mm. Simulations were performed for an injected activity of 300 MBq of [ 99m Tc]-MIBI and an acquisition time of 10 min, which are typical values for a standard clinical MBI study [2]. Further simulations were performed with acquisition times of 5, 2.5, and 1 min. The CNR was estimated for all lesions, and the number of lesions that exceeded the Rose criteria (CNR = 3-5) determined. Results for planar images, with the Gaussian smoothing (FWHM = 5 mm), were compared with tomographic reconstructions.

G. Prototype System
A small proof-of-principle prototype system was built based on the existing Kromek's DMatrix gamma imager. It was comprised of a 2 × 2 array of 7.3-mm thick CZT detectors with DOI capability, each one containing an 11 × 11 array of 2-mm pixels. Two spectral corrections have been implemented. Contrast versus noise curves for the 1-step, 2-step, and hybrid reconstruction approaches, compared to ideal reconstruction. The pinhole separations are indicated at the beginning and at the end of the lines. 1) Charge sharing correction to add back the energy split between the neighboring pixels. 2) DOI correction to correct for the depth dependence of the detected (anode) energy. The depth of the interaction was calculated using the cathode-over-anode (C/A) ratio after the application of charge sharing correction. This method was initially suggested in [40] and [41] for coplanar grid CZT detectors and was later shown to work on pixelated detectors [42]. The subpixelisation into 1mm virtual pixels ("digital" subpixelisation) was implemented by using energies of the neighboring pixels surrounding the pixel with the biggest energy fraction.
The prototype was equipped with a small 3-D printed tungsten MPH collimator mounted on a positioning jig. The collimator was 5-mm thick and had an area of 45 × 45 mm 2 , comprising 49 pinholes with 1-mm aperture and 75 • opening angle. The collimator parameters were optimized using the previously described CNR-based method. The imaging performance of the prototype was characterized using a 1.1 MBq, 1-mm diameter 57 Co point-source mounted on a 3axis motorized stage assembly (adapted from [43]). The stages provide continuous movement of the source in space at a variable speed, thus creating a virtual 3-D activity distribution. Images were reconstructed with 50 ML-EM iterations. Fig. 3 illustrates the mean of the estimated TBR values in all lesions versus noise (COV) in the 144-lesion phantom for varying degrees of MX. The 1-step and 2-step approaches demonstrate opposite characteristics with increased MX. The hybrid approach provides a more stable reconstruction with a gain in performance even with significant MX. We call this the "bow and arrow" graph. The arrow corresponds to the ideal, non-MX data. The data-points at the right-hand side of the graph correspond to the MX-free case, with 20-mm pinhole separation and no overlap between pinhole projections. The various points along the front of the bow correspond to  Fig. 4 presents the results of the GATE simulated comparison of lesion detection sensitivity between a parallel hole CZT scanner [35] and our proposed design for the planar contrast phantom. Our system demonstrates a higher sensitivity in general, which also decreases more slowly with decreasing acquisition time.

C. Design Optimization
With the CNR-based optimization we found that the best choice, based on a large number of simulations, was defined by the following parameters: pinhole aperture = 1.75 mm, pinhole separation = 9 mm, and pinhole opening angle (including penumbra) = 85 • . This configuration resulted in 64% multiplexed data.
In the AUC-based optimization, the overall best result was obtained for the following parameters: pinhole aperture = 2 mm, pinhole separation = 9 mm, and pinhole opening angle = 80 • , with a MX-fraction of 54%. These values are close to those of the CNR-based optimization. Fig. 5 shows one slice from the reconstructed phantom used in the CNR-optimization with both suboptimal (MX-free) and optimized parameters; demonstrating a clear improvement in contrast.

D. Final Performance Evaluation
Comparison of reconstructed versus planar data is presented in Fig. 6. The sensitivity, in terms of the number of lesions with CNR exceeding the Rose criteria of 3, 4, or 5 is presented in Fig. 7 for a range of acquisition times. Comparing the sensitivity at 5 and 2.5 min, it can be seen that reduction in acquisition time (or reduced activity/dose) by a factor of ∼2 can be achieved without loss in terms of performance. System optimization results; reconstructed plane of simulated phantom containing spheres of different sizes and contrasts using nonoptimized (left) and optimized parameters (right).  Fig. 8 shows photographs of the prototype system and Fig. 9 shows the measured 57 Co spectrum before and after charge summing and DOI-based energy correction. The average uncorrected energy resolution (FWHM) of all four modules is 2.1%, which improves to 2.0% after both corrections. The depth resolution was estimated around the middle of the detector using a 250-μm tungsten collimator. The FWHM was estimated as ∼1.2 mm. Fig. 10 shows some initial experimental results with an array of point sources in air with a total number of counts of 1.5 m. The coronal sections show the improvement in depthresolution obtained with subpixelisation and DOI estimation. Fig. 11 shows an image obtained with the prototype system from an acquisition with three layers of point sources with a total number of counts of 2.8 m. The in-plane spatial resolution (FWHM) in the top and bottom planes was 1.67 mm and the resolution perpendicular to the plane was 2.38 mm. In the middle plane, the corresponding values were 2.92 mm and 5.75 mm, respectively.

IV. DISCUSSION
We have proposed a novel detector system for low-dose MBI. Although dose reduction in MBI may seem like  a desirable outcome, it generally leads to reduced lesion detection sensitivity [44]. Furthermore, by simply scanning for longer one can readily increase the sensitivity for the same radiation dose. However, there are practical as well as financial aspects to consider. Therefore, it is necessary to assess all these factors in relation to a smallest detectable lesion size and contrast, taking into account the fact that small lesions with low contrast may not be clinically relevant. However, what the limits should be is not entirely clear.
Our novel system is based on MPH collimators. Usually, MPH geometry is restricted due to the need to avoid MX. The addition of DOI information to determine layers within the detector with variable degrees of MX aids in the solution of the de-MX algorithm. As a result, a relatively large number of apertures can be used with better angular sampling and detection efficiency. The results of both MC and analytical simulations as well as real measured data from the prototype Fig. 9. Three measured 122 keV spectra of 57 Co: raw spectrum, after charge sharing correction, and after DOI-based energy correction.  For geometric reasons, the detection efficiency of the scanner will be lower at the edges of the FOV than in the center due to reduced sampling, which could affect lesion detection in these areas. However, in reality this will be less of an issue as the actual detector array will be wider than the simulated ones (20 × 16 instead of 16 × 16 cm). On the other hand, for lesions near the chest wall, there could be additional problems due to background radiation from organs within the body. This will be investigated.
System design optimization can be done by maximizing CNR, aiming for images with high contrast and low noise level. On the other hand, the use of numerical observers is a way of directly estimating the system's capability in terms of lesion detection. The fact that, in this case, the two different optimization criteria gave similar results increased the credibility of the outcome. Human observer studies are currently in progress to further verify these results.
Direct comparison of results from planar and tomographic studies is not straightforward. Tomographic images may have higher contrast, but they are also affected by noiseamplification during the reconstruction process and the noise has a characteristic appearance. Tao et al. [7] showed that the administered dose could be reduced from 300 to 150 MBq for planar MBI, using the BM3D collaborative filtering algorithm. We did not use such a filter here, but it is likely that application of this type of denoising algorithms to the 3-D images from our proposed system would result in a similar benefit. (This will be further investigated in the future.) In Fig. 7, a 5min acquisition with an administered dose of 300 MBq is equivalent to a 10-min acquisition with a dose of 150 MBq. The graphs suggest that a further reduction by a factor of 2 in scanning time or administered dose is possible with our proposed system without loss of lesion detection sensitivity in comparison to planar imaging.
A tomographic system provides improved depth resolution, and thereby an improved 3-D localization of tumors, which would be helpful in planning and performing biopsies. In planar MBI, two different views are routinely obtained: craniocaudal and mediolateral oblique projections [1]. This may no longer be needed with a tomographic system. It has been shown that an administered activity of ∼300 MBq results in adequate image quality in planar MBI, which can possibly be reduced by a factor of 2 [4], [5]. Our results indicate that with our proposed system the activity can be reduced by another factor of 2, and in addition the scanning time could also be reduced by a factor of 2 by reducing the number of views. Such a significant improvement would overcome the main impediments in the MBI adoption and will position our low-dose MBI technology on par with the mammography dose for complementary screening of patients with dense breast tissue. Further evaluations will be performed using the proof-of-principle prototype.
V. CONCLUSION A new principle for MBI system has been introduced, based on CZT detectors with high-spatial resolution and DOI capability combined with densely packed MPH collimators. A stationary MBI tomographic system has been designed that provides 3-D localization without additional views. Feasibility of the new principle has been verified using MC simulations. The geometry has been optimized using two different approaches, with similar results. Analytical simulations suggest that significant dose reduction can be achieved, opening a pathway to reduce the dose to the mammography level. A feasibility prototype has been designed and built, providing promising initial results.

RECONSTRUCTION ALGORITHMS
The simulated projection data was stored in independent detector blocks, one for each pinhole. These are small 3-D matrices. For the purpose of the de-MX operation of the 2-step reconstruction method, we also utilize virtual 2-D data planes, one for each pinhole. We imagine these planes as being located on the object side of the pinholes. The de-MX algorithm consists of an iterative procedure where data are forward and back-projected between the virtual planes and the detector blocks. Each iteration can be described as follows: where P k i,j is the 2-D virtual plane and V k i,j the 3-D detector block for detector i and pinhole j after k iterations, N d and N p are the numbers of detectors and pinholes per detector, respectively, A is a matrix for transformation from the 2-D to the 3-D data format, the matrix B represents the MX operator, and Q i is the measured data for detector i. The MX process was implemented by simply adding data from different pinholes in detector pixels where there was overlap.
The correction factors corresponding to the three different reconstruction approaches can be described as follows: for the 1-step, 2-step, and hybrid approach, respectively, where Q and Q d are the raw and de-MX projection data, respectively, f k is the current image estimate after k iterations, and FP M {·} and FP{·} represent the forward projection operation with and without MX, respectively (see [20] for more details).