A new virtual ring-based system matrix generator for iterative image reconstruction in high resolution small volume PET systems

A common approach to improving the spatial resolution of small animal PET scanners is to reduce the size of scintillation crystals and/or employ high resolution pixellated semiconductor detectors. The large number of detector elements results in the system matrix—an essential part of statistical iterative reconstruction algorithms—becoming impractically large. In this paper, we propose a methodology for system matrix modelling which utilises a virtual single-layer detector ring to greatly reduce the size of the system matrix without sacrificing precision. Two methods for populating the system matrix are compared; the first utilises a geometrically-derived system matrix based on Siddon’s ray tracer method with the addition of an accurate detector response function, while the second uses Monte Carlo simulation to populate the system matrix. The effectiveness of both variations of the proposed technique is demonstrated via simulations of PETiPIX, an ultra high spatial resolution small animal PET scanner featuring high-resolution DoI capabilities, which has previously been simulated and characterised using classical image reconstruction methods. Compression factors of 5×107 ?> and 2.5×107 ?> are achieved using this methodology for the system matrices produced using the geometric and Monte Carlo-based approaches, respectively, requiring a total of 0.5–1.2 GB of memory-resident storage. Images reconstructed from Monte Carlo simulations of various point source and phantom models, produced using system matrices generated via both geometric and simulation methods, are used to evaluate the quality of the resulting system matrix in terms of achievable spatial resolution and the CRC, CoV and CW-SSIM index image quality metrics. The Monte Carlo-based system matrix is shown to provide the best image quality at the cost of substantial one-off computational effort and a lower (but still practical) compression factor. Finally, a straightforward extension of the virtual ring method to a three dimensional virtual cylinder is demonstrated using a 3D DoI PET scanner.

A common approach to improving the spatial resolution of small animal PET scanners is to reduce the size of scintillation crystals and/or employ high resolution pixellated semiconductor detectors. The large number of detector elements results in the system matrix-an essential part of statistical iterative reconstruction algorithms-becoming impractically large. In this paper, we propose a methodology for system matrix modelling which utilises a virtual single-layer detector ring to greatly reduce the size of the system matrix without sacrificing precision. Two methods for populating the system matrix are compared; the first utilises a geometrically-derived system matrix based on Siddon's ray tracer method with the addition of an accurate detector response function, while the second uses Monte Carlo simulation to populate the system matrix. The effectiveness of both variations of the proposed technique is demonstrated via simulations of PETiPIX, an ultra high spatial resolution small animal PET scanner featuring high-resolution DoI capabilities, which has previously been simulated and characterised using classical image reconstruction methods. Compression factors of × 5 10 7 and × 2.5 10 7 are achieved using this methodology for the system matrices produced using the geometric and Monte Carlo-based approaches, respectively, requiring a total

Introduction
Positron emission tomography (PET) is a functional imaging technique used in clinical diagnostic applications and biomedical research. It is used as a non-invasive in vivo imaging modality to observe biochemical processes in small animal models, in particular for the study of diseases and to assist in the development of new treatments (Rowland and Cherry 2008).
A critical factor in PET image quality is the method used to reconstruct the images from photon coincidence data. Statistical iterative image reconstruction techniques such as the maximum-likelihood expectation-maximisation (MLEM) algorithm (Shepp and Vardi 1982) and the Bayesian reconstruction method (Mumcuoglu et al 1996) have now largely supplanted classical analytic image reconstruction algorithms such as filtered backprojection (FBP) (Kak and Slaney 2001). These iterative algorithms have been shown to significantly improve the quality of the reconstructed image at the cost of increased computational complexity and memory requirements in comparison to FBP (Tohme and Qi 2009).
A central element of all iterative statistical image reconstruction algorithms is the system matrix. For a reconstructed image with a total of I pixels generated using data from a scanner with J possible detector pairs or possible lines of response (LoR), each element p ij of the system matrix represents the probability that the annihilation of a positron emitted from the ith image pixel ( ∈ [ … ] i I 1, 2, ) is detected by the jth ( ∈ [ … ] j J 1, 2, ) detector pair. This system matrix models the underlying physics of the data acquisition process and is utilised in both forward and back projection operations during each iteration of the image reconstruction process. Therefore, image quality (measured via metrics such as spatial resolution, contrast recovery capability and background noise) is highly sensitive to the accuracy of the system matrix. Several approaches have been developed for obtaining an accurate estimate of the system matrix for a given PET system, including analytical calculation, experimental measurements and Monte Carlo simulations.
One well-known analytical method is the Siddon's ray-tracer (Siddon 1985). In this technique, the value of each system matrix element is calculated by the length of the intersection between a specific LoR and an image pixel. For 3D reconstruction, an approach that calculates the volume of intersection between a specific tube-of-response (ToR) and an image voxel is demonstrated to achieve a better contrast-noise trade-off than the basic Siddon algorithm (Scheins et al 2006). The geometric approximation of the LoR/ToR can be further improved by using multiple-ray-tracing schemes (Moehrs andDefrise 2008, Aquiar et al 2010). In addition, the orthogonal-distance based ray-tracer (normalised to the full width at half maximum (FWHM) of the system point spread function (PSF)) has also been evaluated in several studies, and shown to offer significant improvement in contrast recovery and spatial resolution compared to Siddon's method (Pratx et al 2009, Aquiar et al 2010. Alternatively, the system matrix can be estimated experimentally by scanning point or line sources in all locations within the region of interest and recording the associated projection profiles (Frese et al 2003, Tohme and Qi 2009, Alessio et al 2010. This approach potentially provides the most accurate system matrix by taking all physical effects into consideration during the photon detection process. However, it requires a complete data set acquired by a fully functional scanner; it cannot be used in feasibility studies for a newly proposed scanner architecture. Monte Carlo (MC) simulations provide an alternative method for obtaining the system matrix, by including a very detailed model of the scanner geometry and data acquisition process (Veklerov et al 1998, Rafecas and Mosler 2004, Lee and Choi 2012. The major drawback of this method is the potentially lengthy simulation time needed to obtain sufficient statistics for an accurate estimate of the system matrix . Specialised tomographyoriented Monte Carlo simulators such as the Geant4 Application for Tomographic Emission (GATE) (Jan et al 2004) or the more system matrix-oriented egs_pet (Kawrakow et al 2008, Zhang et al 2010 provide some optimisation and/or simplification of the simulated physics processes to increase the speed of simulation, providing an accurate and computationally feasible method of estimating the system matrix.
While all of these methods for system matrix construction are effective, for certain scanner geometries-in particular, for advanced scanners with large numbers of detector elementsthe system matrix can become impractically large. For such scanners, a variety of methods have been developed to reduce the computational cost and in-memory storage requirements of the system matrix. On-the-fly system matrix calculation eliminates the storage requirement entirely, however, this approach has an obvious substantial computational cost (which is imposed for each image reconstruction). Matrix factorisation methods may be used to factorise the system matrix into several components (typically consisting of a geometric projection matrix, a sinogram blurring matrix and an image blurring matrix) which are well-suited to storage in a sparse-matrix format (Qi et al 1998, Zhou andQi 2011). The resulting matrices have far fewer non-zero elements, and the approach provides the flexibility of being able to choose the most appropriate method (analytic, experimental or Monte Carlo) for obtaining each of the aforementioned components. Alternatively, by exploiting the geometric symmetry of the scanner (axial and planar symmetries) (Zhang et al 2010) as well as by using a polar rather than cubic voxelisation method (Cabello and Rafecas 2011), the size of the system matrix can be reduced. In addition, Harraiz et al have proposed the concept of quasi-symmetry (Herraiz et al 2006), in which a number of quasi-equivalent LoR classes are constructed to categorise LoRs associated with different rings along the axial direction for a 3D small animal PET scanner, providing further reductions in system matrix storage requirements. In order to compress the system matrix for a 4-layer DoI PET scanner, Yamaya et al proposed a method to combine deep pairs of DoI layers into shallow pairs of DoI layers whose detector response function highly correlate. As a result, a 117 petabyte system matrix was successfully compressed such that it was able to be fully memory-resident (Yamaya et al 2005(Yamaya et al , 2008. Scheins et al proposed a scanner-independent system matrix generation and compression method which incorporates a virtual generic cylinder model (GCM) and rotation-symmetric voxels (Scheins et al 2011). This method generates a system matrix for an equivalent single-layer virtual detector cylinder located just inside the detector ring of the physical scanner; therefore, the expected matrix compression ratio delivered by rotation-symmetric voxels method can be achieved for a wide range of PET scanner configurations. The system matrix is populated using full detector-to-detector volumes of response (VoRs) rather than the one-dimensional lines of response used in Siddon's method. Respective compression factors of 288 and 320 were achieved for two different PET scanners, with excellent reconstructed image quality. However, although this method is applicable to scanners with DoI capabilities, this is not explored in Scheins et al (2011) or, to our knowledge, in subsequent work.
To provide sufficient DoI information and increase intrinsic spatial resolution, successive generations of PET scanners have tended to reduce the size of the scintillation crystals-or, in the extreme case, may eliminate the scintillator altogether, instead using direct gamma photon detection in pixellated semiconductor detector arrays (typically fabricated from relatively thick silicon or CZT) (Domenico et al 2007, Habte et al 2007, Park et al 2007, Li et al 2013. This geometry results in a very large number of potential LoRs, especially for 3D scanners with a multi-ring geometry. In addition, the emerging innovative design of zoom-in PET with additional detector modules placed inside a conventional detector ring makes its system matrix irregular and dependent on specific scan configuration parameters. Therefore, it is necessary to develop a scanner-independent technique that can achieve significant system matrix compression ratio for unconventional PET scanners with extremely large amounts of DoI information.
In this paper, we present a methodology for achieving very high rates of system matrix compression for PET systems featuring detectors with extremely fine granularity in the radial dimension. This is achieved by exploiting the ability of the virtual cylinder model to represent multiple layers of detector elements in the radial dimension as a single detector volume in the virtual detector layer. In this work, we specifically consider PET systems with a single axial layer with extreme DoI; therefore the resulting virtual cylinder is termed a virtual ring. Whereas Scheins et al use a geometric approach to calculate the values of the system matrix based on volumes of response, we compare an enhanced geometric approach based on Siddon's ray-tracer, incorporating an accurate detector response model, with a method based on a Monte Carlo simulation of a flood-filled disc insert.
The proposed techniques are evaluated by applying them to PETiPIX, a representative ultra high resolution DoI PET system (Li et al 2013), and are shown to provide the expected very high system matrix compression ratio. The proposed methods for system matrix generation are shown to yield images exhibiting excellent image quality. Images reconstructed from GATE simulations of various point sources and phantoms demonstrate excellent spatial resolution, contrast recovery coefficient (CRC), coefficient of variance (CoV), and structural similarity index (which is closely correlated with image quality as perceived by human eyes). The Monte Carlo system matrix generation method is shown to produce the best overall image results at the expense of a substantial one-off computational effort. Furthermore, the 2D virtual ring method has been extended to a virtual cylinder for use in 3D PET systems, where it also achieves a very high factor of reduction in system matrix memory requirements. The effectiveness of the Monte Carlo system matrix generator is demonstrated using a virtual cylinder model constructed for a 3D DoI scanner; data resulting from Monte Carlo simulations of the 3D scanner imaging a 3D Ultra Micro Hot Spot Phantom TM 10 model are used to reconstruct a high quality 3D image (Data Spectrum Corporation 2007).
The remainder of the paper is organised as follows: section 2 describes the proposed method and discusses the implementation of the virtual ring model and its 3D extension. Section 3 presents the results of applying the virtual ring model with both system matrix generation methods to the PETiPIX scanner; results from MLEM image reconstruction performed using each system matrix construction method are compared using Monte Carlo simulations of point sources and a variety of standard phantoms, and evaluated using several objective and perceptual figures of merit. A successful demonstration of the 3D extension of virtual ring method is also presented in this section. The implications of the results are discussed in section 4; conclusions and proposed future work are discussed in section 5.

The virtual ring method
For small-bore animal PET systems based on scintillation crystals with DoI capability or pixellated semiconductor detector arrays, scanner geometry is not well-approximated by the conventional cylindrical model with a single layer of detector elements around the ring.
As shown in figure 1(a), detector modules are represented as rectangular detector modules arranged around the FoV, with the inner edges of the modules forming a regular polyhedron. Each detector module consists of a rectangular array of individual semiconductor pixels and/ or scintillation crystals; for high resolution PET systems, the number of elements of these arrays can be very high, particularly where DoI information is obtained through pixellation in the radial direction.
Consider a true coincidence originating from a pixel b somewhere in image space as shown in figure 1(b), with photons recorded by two detector elements S i and S j , located at ( ) x y , i i and ( ) x y , j j respectively. Detector pair d ij is defined by S i and S j , and p(b, d ij ) is the probability that an emission from pixel b is detected in detector pair d ij .
A virtual ring with its centre co-located with the centre of the FoV and a radius of R is created as shown in figure 1(b). The virtual detector ring is comprised of square detector elements of width W. R is chosen such that the ring surrounds the imaging subject, yet is fully contained within the actual scanner FoV, so that the line of response of any true coincidence that can be detected by the actual ring will intersect with the virtual ring at exactly two points. The Cartesian coordinates of these two points of intersection ( ) ′ ′ x y , i i and ( ) ′ ′ x y , j j are calculated using the method proposed by Rhoad et al (1991): (1) , the corresponding virtual detector pixels ′ S i and ′ S j can be identified, and a detector pair v between detector positions on the virtual ring is determined. The activity contributions from the LoRs between all detector pairs in the actual scanner geometry that can be represented as an LoR between virtual detector pair v are added together, as shown in figure 1(c). Thus, *( ) n v is the total number of counts measured for detector pair v in the virtual ring. p (b, v) can then be calculated as *( ) n v divided by the counts emitted from image pixel b. As a result, the system matrix from the actual detector ring is transformed into an equivalent form for the virtual ring, which may be used for image reconstruction in lieu of the original system matrix. Provided the resolution of the virtual ring is sufficient to avoid aliasing, the resulting matrix should produce indistinguishable results when applied to image reconstruction. Appendix A provides a brief comparison of image reconstructions for a simple PET geometry performed with an uncompressed system matrix and the virtual ring method to demonstrate that the impact of the virtual ring on the resulting image quality is minimal.
The virtual ring reduces the dimensions of the system matrix by eliminating all potential LoRs outside of the FoV, and by summarising the DoI information into an equivalent compact form. Consequently, significant compression of the system matrix size is achieved, particularly for scanners with extremely high DoI resolution.

Image reconstruction
The image reconstruction algorithm used throughout this paper is based on MLEM. Following the method proposed by Shepp and Vardi (1982), it is assumed that two annihilation photons are detected in coincidence by a particular pair of detector elements, denoted d, where ∈ [ … ] d D 1, 2, , and D is the total number of possible detector pairs in the physical scanner. Let λ( ) b represent the expected number of detected events originating within pixel b, one of a total of B pixels inside the scanner field of view, ∈ [ … ] b B 1, 2, , . The recorded data for each potential detector pair are then denoted 1, 2, . The MLEM algorithm is then defined as follows: is the next estimated value of pixel b based on the current estimate of λ ( ) b old , and p(b, d) is the probability that an emission from pixel b is detected in detector pair d.
Based on the virtual ring system matrix element, p(b, v), defined in section 2.1, the MLEM algorithm iteratively optimises λ( ) b as: V 1, 2, , and V is the total number of detector pairs in the virtual ring.

Scanner case study: PETiPIX
The performance of the virtual ring method is demonstrated in simulation using the PETiPIX scanner as a case study. PETiPIX is a small animal PET scanner with DoI capabilities, which offers exceptionally high spatial resolution (Li et al 2013). The trade-off for its high spatial resolution is relatively low detection sensitivity, which is a consequence of the use of direct detection of gamma photons and recoil electrons in silicon rather than traditional indirect detection via a scintillator crystal. As shown in figure 2, PETiPIX utilises four high resolution Timepix pixellated silicon detectors in an edge-on configuration, for maximum sensitivity and DoI detection. Each detector consists of an array of × 256 256 pixels, with pixel dimensions of 55 μm × 55 μm × (300-1000) μm (the thickness depending on the specific detector type; timepix detectors can be fabricated on silicon or CZT with various thicknesses). The complete scanner has a 15 mm wide square field of view (FOV), and the detector ring can be translated axially in steps of 0.0125 mm, allowing multiple 2D scans to be performed for different slices across the entire mouse brain. For the highest-quality images, the PETiPIX detector ring is positioned at the desired axial position, and two acquisitions are performed, with the ring rotated in the azimuthal dimension by 45 degrees for the second acquisition. However, in this paper, the simulated PETiPIX scanner is configured to perform a single acquisition from a fixed position. Coincidences between any two detector arrays are permitted, including adjacent detectors, in order to maximise the field of view.
PETiPIX has already been extensively evaluated and characterised via Monte Carlo simulation, using the GATE Monte Carlo simulator with classical image reconstruction techniques (Li et al 2013). An overall spatial resolution of 0.29 mm FWHM was obtained using filtered backprojection, and sensitivity was estimated to be 0.01%. The relatively low sensitivity (compared to that provided by other small animal PET scanners) is due to the low stopping power of silicon, as well as the small solid angle covered by the detectors. Most photons are not detected via photoelectric interactions with the detector; instead, most interactions are detected by tracking the recoil electrons resulting from Compton scattering. The recoil electrons are typically registered by a small number of adjacent detector pixels along their track. Since the total track length is significantly less than mean positron range in water, recoil electron range does not significantly degrade the resulting image. The point of interaction is taken as the pixel with the most energy deposited along the track.
In this paper, PETiPIX is chosen as a case study due to the high resolution of the detector modules in the radial direction. The geometry of PETiPIX is a pathological case for classical system matrix construction as a single Timepix module contains 65 536 small pixel elements. For an image consisting of × 256 256 pixels with the width of image pixels being 55 μm, the dimensions of the uncompressed PETiPIX system matrix are estimated to be 25.8 billion (LoRs) × 65 536 (image pixels). The estimated size of this system matrix is 6.76 PB, assuming a precision of 4 bytes per element. Thus, the original system matrix remains impractically large, necessitating further reductions in system matrix size.

Implementation of the virtual ring method for PETiPIX
For PETiPIX's 15 mm × 15 mm FOV, the positron-emitting source will always be separated from the detectors by at least 1-2 mm. In this paper, the positron sources are contained entirely within the a disc measuring 12 mm in diameter. The virtual ring radius r is set to 6.5 mm, which ensures that all true coincidences detected by the actual scanner will have two points of intersection with the virtual ring.
It is also possible to calculate the maximum virtual detector pixel size which will provide a reconstructed image of equivalent quality to the uncompacted system matrix, by applying equation (2) in Moses (2011). The theoretical spatial resolution of the original scanner is first calculated, neglecting the effects of positron range and other non-idealities, and then the dimensions of a detector element on the virtual detector which will give the same spatial resolution are computed. For PETiPIX, this calculation suggests a pixel size of 70 μm. However, in practice, a larger detector element can be used to exchange a small degradation in image quality for a considerably smaller system matrix; if positron range and other non-idealities are also included in the simulation (or an experimental system is used) then these will almost  certainly be the limiting factors to scanner performance rather than virtual pixel size. For most cases, this is an acceptable compromise between accuracy and performance.
To determine an appropriate compromise virtual pixel size, The virtual ring is initially segmented into 45, 90, 180 and 360 detector elements to observe the effects on image quality with several test phantoms (arrays of point sources). Very little additional improvement is observed when the number of elements is increased from 180 to 360, and improvements beyond this point are not expected to be measurable. Therefore, 360 virtual detector elements are used, each with a width of 110 μm. The reconstructed image consists of × 256 256 pixels, with the width of image pixels being 0.05 mm.
2.4.1. Siddon ray-tracer with detector response function based system matrix. Siddon's raytracer is a well-known method for computing geometrical efficiency in emission tomography systems (Siddon 1985). As shown in figure 3, the probability that an emission from pixel b is detected by detector pair d, p (b, d), is calculated as the length of intersection between detector pair d and image pixel b. This model is most accurate when the scanner has a standard ringtype geometry and the detector elements are relatively small compared to the FoV. However, in unconventional scanner geometries, such as PETiPIX or ring type scanners with small FoV, it oversimplifies the system matrix by ignoring the detector response function.
As the Timepix detector has very high pixel density (over 330 pixels mm −2 ), it can be represented as a 14 mm × 14 mm continuous silicon slab. Each element of the virtual ring system matrix p (b, v) can therefore be calculated as: where p 1 (b, v) is the length of intersection between virtual ring detector pair v and image pixel b, and p 2 (v) is the virtual ring detector response function calculated as: where μ ρ is × − 8.748 10 2 cm 2 g −1 , ρ is 2.3290 g/cm −3 for silicon. As shown in figure 3, L 1 and L 2 are the length of intersections between the LoR and the two Timepix detectors respectively. As a result, a virtual ring detector pair together with its detector response function calculated using the above equation is able to accurately represent all possible detector pairs of the PETiPIX scanner along that particular LoR.
In the remainder of this paper, the virtual ring system matrix derived from Siddon's raytracer method is denoted V-S and the system matrix based on Siddon ray-tracer with virtual ring detector response function is denoted V-SD.
2.4.2. Monte carlo based system matrix. To generate the data for constructing the system matrix based on Monte Carlo simulation, an accurate model of the scanner is simulated in GATE version 7.0. For the PETiPIX scanner, a simulated flood-filled disc measuring 12 mm in diameter containing a solution of 18 F in water is placed at the centre of the scanner FoV. The disk is loaded with an activity of 37 MBq (although the actual activity is irrelevant to the simulation as execution time is only dependent on the total number of decays). To decrease the required simulation time, positron range and photon non-collinearity are not included. Approximately 1.7 billion true coincidences are recorded in the list-mode format. This data set is mapped to the virtual ring and used to generate the Monte Carlo-based system matrix (V-MC).
Two well-known statistical metrics are used to evaluate the compactness and quality of the resulting V-MC matrix: the number of non-zero matrix elements (N ), and the mean relative error (MRE) σ rel , originally introduced by Rafecas and Mosler (2004). The latter of these metrics is calculated as , Mosler 2004, Ortuno et al 2010).
As the flood-filled disc has a diameter of 12 mm and the virtual ring is 13 mm in diameter, it is not possible for true coincidences to occur between virtual ring elements which are separated by less than a particular number of detector elements. For a given geometry, the minimum valid difference D between detector indices can be calculated; for PETiPIX, the detector pair must be separated by a minimum of 45 virtual detector elements for a coincidence to be considered valid. Therefore, the raw data from the Monte Carlo simulation are filtered to eliminate all events with D less than 45, as these coincidences must be due to inter-detector scattering (randoms, i.e. coincidences originating from independent events which occur within the timing window, are already filtered out based on their event ID).

Phantoms for evaluation and figures of merit
Phantom simulations were performed using GATE 7.1, with the same scanner and detector models as used for the system matrix generator. Simulations were performed on a cluster consisting of 90 Intel Xeon E5-2695 v3 CPU cores. The physics processes, including gamma interactions (photoelectric effect, Compton scattering and Rayleigh scattering), electron ionisation, multiple scattering and Bremsstrahlung were modelled using the Standard GEANT4 electromagnetic physics Package. The detector energy window was set to 10-650 keV with a constant energy resolution of 1.2 keV for all energies. The wide energy resolution is required because most photons interact with the detector via Compton scattering rather than the photoelectric effect (due to the low stopping power of silicon) and therefore must be indirectly detected by the detection of recoil electrons (and occasionally by low-energy scattered photons). The timing resolution was set to 10 ns. All phantom simulations included modelling of positron range and annihilation photon non-collinearity.
2.5.1. Point source phantom. Spatial resolution was assessed by simulating a spherical water phantom, 13 mm in diameter, containing eleven 37 MBq 18 F point sources. As shown in figure 4, the point sources are uniformly distributed along two radii (with orientations of θ ϕ ( ) = (°°) , 0, 0 and (°°) 0 , 45 ) within the water phantom, separated by by 1 mm increments. 430 000 coincidences were recorded, and the images reconstructed by 2D-FBP, MLEM with Siddon's ray-tracer system matrix (MLEM-V-S), MLEM with Siddon's ray-tracer with virtual ring element response function (MLEM-V-SD), and MLEM with Monte Carlo-based virtual ring system matrix (MLEM-V-MC). The spatial resolution was measured as the FWHM of the point source intensity profiles measured in the radial direction along each of the two radii. Gaussian curve fitting was not applied as the point source intensity profile is expected to be a cusp-like shape due to the blurring from non-zero positron range.

Contrast phantom.
The contrast phantom, also shown in figure 4, consists of five rods inside a hollow cylinder measuring 10 mm in diameter with an axial length of 0.3 mm. The diameters of the hot rods are 0.5 mm, 1 mm and 1.5 mm respectively; the ratio of activity in the hot rods to the warm background was : 5 1. The largest two rods, measuring 2 mm and 2.5 mm, are cold rods with zero activity in each. A total of 10 6 coincidences are generated including true, random and scatter events. The images are reconstructed using the MLEM-V-S, MLEM-V-SD and MLEM-V-MC algorithms. A total of 300 iterations was used for the iterative reconstruction algorithms.
The contrast recovery coefficient (CRC) is a widely used metric for quantitative evaluation of the reconstructed images. Three regions of interest (RoIs) are selected for this contrast phantom: the cold rods, hot rods and the background respectively. The CRCs of the cold dots and the hot dots are calculated as follows:  (11) where C b is the average number of counts recorded in the warm background region, C c is the average number of counts in the cold RoI, C h is the average number of counts in the hot RoI and R is the true activity concentration ratio between the hot RoI and the warm background (R was set to 5 in the above simulation). In this paper, all cold and hot regions in the contrast phantoms are considered in the calculation of hot and cold RoI CRCs. The normalised noise level is commonly denoted as the coefficient of variation (CoV), which is calculated as where μ s is the mean value and σ s is the standard deviation of the background.
The complex wavelet structural similarity (CW-SSIM) index is an effective method of objectively assessing perceptual image similarity based on structural differences. It is considered to be more consistent with the subjective perceptual response of the human eye perception than other objective image quality metrics, and is widely employed in the field of medical imaging (Ng et al 2012, Casti et al 2015. The detail of this method can be found in Sampat et al (2009).

MOBY phantom.
A digital mouse phantom, MOBY, is used to assess overall image quality (Segars et al 2004). A section through the brain is selected which includes sections of cortex, striatum, amygdala, hypothalamus and thalamus. The relative uptake of FDG for each voxel in each region is based on data published by Welch et al (2013) and is shown in figure 5. A total of 10 7 coincidences are recorded. The images are reconstructed using the MLEM-V-SD and MLEM-V-MC algorithms with 10, 50 and 100 iterations.

Validation with 3D DoI PET
The virtual ring method presented in this paper can be also extended to 3D PET systems, and is especially beneficial when applied to scanners with extensive 3D DoI capability. Instead of using a virtual ring as in the 2D case (discussed in section 2.1), a virtual cylinder is employed to transform the original system matrix into an equivalent single-layer form, such that the extensive 3D DoI information is summarised into a more condensed form by locating the points of intersection between lines of response and the virtual cylinder. This method is implemented and a proof-of-concept validation based on a 3D DoI small animal PET scanner is presented.
Since PETiPIX is not intended for use as a 3D imaging system, but rather for high-resolution examination of a particular slice of interest, the scanner used to demonstrate the 3D validation in this paper is a 3D extension of a specific configuration of the variable-geometry 2D DoI CMRPET scanner, previously introduced and experimentally validated in Safavi-Naeini et al (2011), Safavi-Naeini (2013. CMRPET3D has an inner diameter of 46 mm and an axial FoV of 31 mm consisting of 8 rings, with an inter-ring gap of 1 mm. Each ring contains 12 detector modules, each of which containing a × 4 4 array of LYSO crystals in an edge-on configuration attached to a matching array of SiPMs as the detector and readout system. The crystal size is × × 3 3 3 mm 3 . Depth of interaction information is provided by CMRPET3D's pixellation structure in the radial direction as well as the arrangement of different rings in the axial direction. A 3D model of the scanner is shown in figure 6. To obtain its system matrix via the Monte Carlo method, an accurate model of this scanner is simulated in GATE. A simulated flood-filled cylinder, 30 mm in diameter and 30 mm in height, is placed at the centre of the scanner FoV and filled with 37 MBq 18 F. Positron range and photon non-collinearity are not included. A total of 2 billion true coincidences are recorded and processed to construct the Monte Carlo-based 3D system matrix.
The virtual cylinder diameter and height are set to 35 mm and 30 mm respectively, ensuring that two points of intersection exist between the LoRs of all detected true coincidences and the virtual cylinder. The virtual cylinder is segmented into 20 rings without gaps, and each ring consists of 90 virtual voxels, for a total of 1800 virtual cylinder voxels. The virtual cylinder voxel dimensions are therefore × × 1.22 1.22 1.5 mm 3 . The reconstructed image resolution is × × 128 128 128 voxels, with voxel size of × × 0.25 0.25 0.25 mm 3 . The image volume covers the full range of potential locations for positron annihilation.
To validate the virtual cylinder system matrix, a 3D Ultra-Micro Hot Spot phantom, filled with a total activity of 74 MBq 18 F was simulated in GATE. It consists of a hollow cylinder, 30 mm in diameter and with a height of 30 mm, filled with the 18 F solution. It accommodates six sectors of hot rods, 25 mm in height, with diameters of 2.4 mm, 2 mm, 1.7 mm, 1.3 mm, 1 mm and 0.75 mm respectively. The centre-to-centre spacing within each sector is twice the rod diameter. × 3 10 8 coincidences are recorded and the 3D image is reconstructed by a modified MLEM-V-MC algorithm designed for 3D image reconstruction with 20 iterations. The implementation is fully 3D, and includes all LoRs with oblique angles. Table 1 compares the sizes of the virtual ring-based system matrices resulting from the V-S, V-SD and V-MC generation methods. The number of non-zero elements, compression ratio, compressed size in computer memory and an approximate computational requirements are also shown in this table. In figure 7(a), the number of non-zero elements are plotted versus number of detected coincidences per image pixel for uncorrected V-MC and corrected V-MC (D = 45) respectively. It can be observed that the corrected V-MC has fewer non-zero elements and converges faster than the uncorrected V-MC. The MREs of both corrected and uncorrected V-MC generators are shown in figure 7(b).  Table 2 shows the uniformity of image reconstructions of a simulated flood-filled disc with the same dimensions as those used to populate the Monte Carlo-based system matrix. The uniformity of the Monte Carlo-based system matrix is the best by a significant margin, and in fact the uniformity metrics continue to decline as a larger number of coincidence pairs are used. The V-S method performs very poorly since an accurate model of the detector response function is not included.   Note: The total number of elements for the uncompressed PETiPIX system matrix is estimated to be × 1.69 10 15 . The compression ratio shown here results from applying the virtual ring method first and then removing all zero elements. The compression ratio due to the removal of all zero elements (by using a sparse matrix representation) is estimated to be around 10 2 . Computation time is approximate.  Figure 10 compares the average FWHM of all 11 point source intensity profiles using each matrix construction method as a function of the number of MLEM iterations, ranging from 1 to 100. Figure 11 compares images of the contrast phantom reconstructed using the three variants of the MLEM algorithm under evaluation for 10, 50, 100 and 300 iterations. Figure 12 compares the intensity profiles of the images reconstructed using MLEM-V-SD and MLEM-V-MC. The profiles are measured horizontally across the two hot sources (0.5 mm and 1.5 mm in diameter, respectively), for images resulting from 50 iterations of each algorithm.

Contrast phantom
Quantitative    (2008); for both metrics, a completely flat field should yield a value of zero. A total of × 1.2 10 8 coincidences are used for each method. Figure 16 shows images of the MOBY phantom reconstructed using the MLEM-V-SD and MLEM-V-MC methods with 10, 50 and 100 iterations. The image is shown with the 'jet' colour map.

Validation of 3D DoI PET
Based on the flood-filled cylinder simulation, a virtual cylinder-based system matrix is generated. The number of non-zero elements is approximately × 9.2 10 8 , which corresponds to a compression ratio of × 7.3 10 3 . Approximately 4000 core-hours (a total of around 45 h on 90 Intel Xeon E5-2695 v3 CPU cores) of processor time is used for this simulation. The compressed memory-resident matrix size is 24 GB. No symmetry-based compression methods were utilised in this simulation; a further reduction in memory requirements to around 2 GB can be expected by exploiting the transaxial symmetry of the scanner, since there are 12 modules around the complete ring.   Figure 17 shows the reconstructed images of the hot rod secontion of a simulated Ultra-Micro Hot Spot phantom, generated using the MLEM-V-MC method with 20 iterations (further significant improvements were not observed with more MLEM iterations for this simulation). The hot rods with respective diameters of 2.4 mm, 2 mm and 1.7 mm are clearly resolvable.

Discussion
The virtual ring method achieves an extremely high compression ratio for the PETiPIX system matrix. This is due to the very large amount of DoI information, a result of the ultra small pixellated detector elements of the PETiPIX scanner. By storing the virtual ring system matrix in sparse matrix format, the final memory consumption of matrices generated using the V-S, V-D and V-MC methods are 0.54 GB, 0.54 GB and 1.19 GB respectively; this is much less than the physical memory of a typical (c. 2015) desktop PC. Symmetry-exploiting compression along the transaxial direction of the PETiPIX scanner has not been applied here; however, further improvements in compression ratio can be expected.
As shown in figure 7(a), the number of non-zero elements for the corrected V-MC converges rapidly when more than × 5 10 3 coincidences are recorded per pixel. In contrast, convergence is much slower in the case of the uncorrected system matrix, with the number of non-zero elements increasing almost linearly after × 1 10 4 recorded coincidences per pixel. This is because randoms and scatter events (D < 45) significantly contribute to the whole dataset, occupying new elements of the system matrix with very low counts.
The value of σ rel for the uncorrected V-MC shown in figure 7(b) decreases from 0.94 to 0.7 as the number of detected coincidences per image pixel increases to about × 7 10 3 . Beyond this point, σ rel starts to increase again, since the increasing number of randoms and scatter events (D < 45) degrade the matrix. For the corrected V-MC, the σ rel monotonically decreases from 0.9 to a limiting value of approximately 0.52 once the number of detected coincidences per pixel exceeds × 2 10 4 . At this point, the system matrix has converged; additional simulation time will result in no further significant improvement to the accuracy of the matrix. The reconstructed images shown in figure 8 demonstrate that the FBP algorithm will cause significant image distortion. This is because the PETiPIX scanner is stationary during the data acquisition process, resulting an incomplete sinogram (covering only half of the full set of full angular information). This causes analytical reconstruction methods such as FBP to perform poorly, particularly when some key structure of the images from particular angles (in the case of PETiPIX, angles along diagonals between detector modules) lack sufficient recorded LoRs. As a result, the point sources distributed along the horizontal radius cannot be resolved accurately; however, the point sources along the diagonal radius at θ ϕ ( ) = (°°) , 0, 45 can be resolved well due to the more complete data at this angle. In contrast, the iterative MLEM algorithm with an accurate system matrix suffers very little from the incomplete angular coverage, as demonstrated in figures 8(b)-(d). Point sources along both radii are reconstructed clearly.
As can be seen in figure 9, intensity profiles taken through images reconstructed using the MLEM-V-MC method provide the highest peak-to-valley ratio compared to those using on MLEM-V-S and MLEM-V-SD. The average FWHM of the point sources shown in figure 10 illustrates that the MLEM-V-MC algorithm results in the highest spatial resolution of approximately 0.16 mm when more than 40 MLEM iterations are used. MLEM-V-SD yields slightly (a) better spatial resolution (approximately 0.18 mm) compared to the MLEM-V-S (approximately 0.19 mm) method due to additional detector response modelling. Figure 11 shows that as the number of iterations increases, application of the MLEM-V-S algorithm results in an increasingly distorted image, especially at its geometric centre. This is a consequence of the fact that the MLEM-V-S method does not accurately model the unique geometry the PETiPIX scanner; it assigns a very high value for the system matrix elements of central image pixels, which is valid for a conventional ring-structured PET scanner, but not for unconventional scanner geometries such as PETiPIX. By applying an accurate detector response model, shown in figures 11(e)-(h), or employing accurate MC simulations, as in figures 11(i)-(l), all hot and cold regions can be resolved unambiguously. It should be noted that some distortion of those images is still evident (e.g. the in the region between two cold disks). This is due to the fact that the recorded coincidences only account for half of its angular coverage, even though this fact has been addressed during the construction of the MLEM-V-SD and  MLEM-V-MC system matrices. These distortions will disappear if the scanner is rotated by 45 degrees in order to cover the remaining angles. However, the excellent images obtained with MLEM-V-SD and MLEM-V-MC algorithms demonstrate the potential of these approaches when only partial angular coverage is available. The line profiles shown in figures 12(a) and (b) demonstrate that MLEM-V-MC results in a less noisy image than MLEM-V-SD, as its system matrix is more complete (having more nonzero elements), at the cost of greatly increased computational workload for the generation of the system matrix. However, it should be noted that these pre-computed system matrices can be repeatedly used unless the scanner has some structural changes.
For the cold RoI, as the number of MLEM iterations increases from 1 to 500, images reconstructed using MLEM-V-MC exhibit a higher CRC (nearly 80%) compared to both MLEM-V-S (30%) and MLEM-V-SD (70%). The three methods have similar results for hot RoIs (between 72% and 78%), while MLEM-V-SD is slightly better than MLEM-V-S and MLEM-V-MC.
When taking the image background noise, or CoV, into consideration, it can be seen from figure 14(a) that the image reconstructed using MLEM-V-MC exhibits a higher CRC with a  (d)). Similar trends can be observed for hot ROIs as well ( figure 14(b)), showing that the MLEM-V-MC method yields a lower CoV when a given CRC is required.
As shown in figure 15, the CW-SSIM index of the image reconstructed using MLEM-V-MC is greater than 0.8 when a sufficient number of iterations are applied (between 20 and 200). This is consistently higher than for the images constructed using MLEM-V-SD (0.7-0.8). The CW-SSIM index of the images reconstructed using MLEM-V-S is lower than 0.3 even when 50 or more iterations are applied, again due to the distortion which is visible in figures 11(c) and (d). Figure 15 demonstrates that MLEM-V-MC produces reconstructed images closer to the original image in terms of the human visual perception model.
The MOBY phantom shown in figure 16 demonstrates the potential of PETiPIX to image a realistic low-contrast radiotracer distribution in a mouse brain. The different FDG concentrations in each structure can be resolved well in the images generated using MLEM-V-MC with increasing sensitivity as more iterations are used. For the images reconstructed from MLEM-V-SD, the cortex is not resolved clearly compared to the other parts of the brain. This is a consequence of the better CRC offered by MLEM-V-MC over MLEM-V-SD, allowing the relatively low contrast ratio between the cortex and the whole brain background ( : 1.06 1) in the brain model to be resolved.
The 3D images reconstructed using the MLEM algorithm with the system matrix generated using the virtual cylinder model demonstrate that the method works well when extended to 3D scanner geometries, while providing very high levels of system matrix compression. Further improvements in compression ratio can be obtained when conventional symmetry-exploiting compression methods along the transaxial and axial planes are additionally employed. The spatial resolution obtained with the CMRPET3D scanner can be estimated from the line profile shown in figure 17; the three largest hot rods (with diameters of 2.4, 2.0 and 1.7 mm) are clearly resolvable. The spatial resolution obtained is comparable (slightly superior) to the previously-published results for the original two-dimensional CMRPET system of approximately 2 mm. However, the primary purpose of this demonstration is to show that the technique provides good 3D image reconstruction for 3D scanners while offering significant system matrix compression. In principle, the proposed method will also work with 3D scanners with much smaller crystals (and hence greater spatial resolution and better system matrix compression ratio), with the obvious proviso that generation of the system matrix will require a correspondingly greater quantity of CPU time.

Conclusion
In this paper, an innovative method for constructing a compact yet accurate system matrix for MLEM-based image reconstruction using the concept of a virtual detector ring has been developed. This method is validated using the PETiPIX small animal PET scanner design, which consists of ultra-high-resolution pixellated silicon detectors in an edge-on configuration.  The virtual ring system matrix obtained by two different methods (Siddon's ray-tracer with virtual ring detector response and Monte Carlo simulation) achieves a compression factor of over × 2 10 7 and demonstrates excellent performance in terms of spatial resolution, CRC, CoV and CW-SSIM index for various phantoms. Reconstructed images of the MOBY mouse brain phantom with a realistic radiotracer distribution demonstrate that by applying the proposed methods, the PETiPIX scanner is able to resolve brain structures with sub-millimetre resolution. The MLEM algorithm with the virtual ring method is therefore shown to be an effective and flexible method for PET image reconstruction by a range of well-known image quality metrics, particularly when used with the Monte Carlo-based system matrix generator.
Finally, a 3D extension of the virtual ring technique to a virtual cylinder has also been developed and successfully demonstrated with a small-animal 3D PET scanner design, showing that the proposed methods are both practical and effective for a wide variety of PET system geometries.
CW-SSIM index decreases for both methods slightly, dropping to approximately 0.9 at 500 iterations. This demonstrates that the images reconstructed using both virtual ring and uncompressed system matrix are very close to the original image in terms of human visual perception. In conclusion, the additional procedure of applying a virtual ring to the raw simulation data will not cause a significant decrease in accuracy compared to an uncompressed system matrix, provided that the number of virtual detector elements is chosen properly. Figure A1. Comparison of CRC (cold RoI and hot RoI) and CW-SSIM for MLEM with between 1 and 500 iterations performed with an uncompressed system matrix 'ground truth' and the virtual ring method. K Li et al Phys. Med. Biol. 60 (2015) 6949