Towards continuous-to-continuous 3D imaging in the real world

Imaging systems are often modeled as continuous-to-discrete mappings that map the object (i.e. a function of continuous variables such as space, time, energy, wavelength, etc) to a finite set of measurements. When it comes to reconstruction, some discretized version of the object is almost always assumed, leading to a discrete-to-discrete representation of the imaging system. In this paper, we discuss a method for single-photon emission computed tomography (SPECT) imaging that avoids discrete representations of the object or the imaging system, thus allowing reconstruction on an arbitrarily fine set of points.


Introduction
Mathematical tomography centers around integral transforms in which both the object and the image are treated as functions of continuous variables; we refer to an integral transform as a continuous-to-continuous (CC) operator.
The ubiquitous practice of representing real-life objects (such as organs, bones, tissues, etc in the case of medical imaging) with a finite set of intensities over a 2D or 3D grid of pixels or voxels often leads to inaccuracies when the object itself presents many features that cannot be represented using a grid of pixels or voxels. Similarly, sets of discrete data (such as projection images in the case of emission tomography) are typically used as the starting point to perform reconstruction. These discretizations of the object and the data produced by the imaging system are almost always assumed, leading to matrix representations or discrete-to-discrete (DD) operators.
In this paper, we discuss our approach to approximating CC operators for real single-photon emission computed tomography (SPECT) imaging systems in which we use measured calibration data from real detectors and real multi-pinhole imaging systems, thus potentially avoiding discrete representations of the object or the raw data. We use photon-processing detectors to collect data with no binning involved, and show how such data sets can be reconstructed using iterative maximum-likelihood (ML) estimation algorithms. Our approach does not introduce any error due to discretization of the measurement, and it allows reconstructions over an arbitrary sets of points, not necessarily arranged as a 3D uniform grid of voxels.
Prior art in list-mode reconstruction can be found in the work by Levkovitz et al (2001). In it, the authors derive an algorithm for direct reconstruction of list-mode data for positron emission tomography (PET) by starting from an expression of the maximum-likelihood expectation maximization algorithm for binned data, and then they consider the limit of no more than one count per bin. An alternative formulation, also for PET imaging, is provided in Parra and Barrett (1998). Similar to Barrett et al (1997) and Parra and Barrett (1998), our algorithm assumes that the list-mode data are positions of interaction as estimated (Moore et al 2007) from data collected with the gamma-ray cameras in the SPECT scanner. These estimates are allowed to vary over a continuous domain, and probability density functions are evaluated on-the-fly for each item of the list. This paper is organized as follows. Section 2 presents the mathematical notation used in this paper; section 3 discusses photon-processing detectors and introduces a point process we use to mathematically represent list-mode data. In section 4 we discuss the CC imaging operator as a mapping from the object space to data space. This treatment allows us to derive an explicit expression for the kernel of this CC operator. In section 5 we present our reconstruction approach, which is based on the list-mode maximum-likelihood expectationmaximization (MLEM) algorithm. In section 6 we discuss reconstruction results with real data. Future work and conclusions are provided in sections 7 and 8, respectively. This paper is an extended version of the conference proceeding (Caucci et al 2015a) presented by the authors at the 13th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine ('Fully 3D 2015') held in Newport, Rhode Island, USA.

Mathematical formulation
In the absence of detector noise, imaging systems are often modeled as linear continuous-todiscrete (CD) operators ℋ that map the object function f (r) to a set of numbers g. In abstract form, this mapping is written as (Barrett and Myers 2004): in which f corresponds to the object function f (r) and belongs to the Hilbert space of squareintegrable functions of the variable r. In our treatment, we will consider object functions f (r) that depend only on the 3D continuous variable r; more general cases of object functions depending on other variables are possible. We will further restrict our attention to SPECT imaging, and will assume that f (r) Δr Δt is the mean number of gamma-ray photons emitted isotropically from a small volume Δr (centered around point r) and during time interval Δt. Hence, the units of f (r) are photons/(s · m 3 ). Finally and unless otherwise stated, we will ignore changes in nuclear activity due to radioactive decay.
If we assume that g is a vector with M components denoted as g 1 , …, g M , then the expression in (1) can be written in component form as in which ∫ S … denotes integration over the support S of the object f (r), and the function h m (r) represents the system response at point r in the object and for the mth element of g.
In the development of reconstruction algorithms, it is customary to characterize the imaging system with a system matrix H, so that a discrete representation f of the object function f (r) can be related to the noise-free data g according to the following equation (Barrett and Myers 2004): (3) In the expression above, f is a vector with N components f 1 , …, f N and H is an M × N matrix with components h m,n . In component form, the expression in (3) becomes g m = ∑ n = 1 N ℎ m, n f n .
(4) Because both f and g are vectors of finite size, this representation is often referred to as discrete-to-discrete (DD) representation (Barrett and Myers 2004).
One can relate the function f (r) to the discrete vector f, and the kernel h m (r) of ℋ to the elements h m,n of the matrix H by first assuming that there exist functions ϕ 1 (r), …, ϕ N (r) so that f(r) = ∑ n = 1 N f n ϕ n (r), (5) in which f 1 , …, f N are real numbers. Upon substitution, we obtain: which shows ℎ m, n = ∫ S ℎ m (r)ϕ n (r)d 3 r .
If functions ϕ 1 (r), …, ϕ N (r) are orthonormal, the numbers f 1 , …, f N are calculated from f (r) according to f n = ∫ S f(r)ϕ n (r)d 3 r .
It is worth noting that the condition in (5) is very restrictive. Real-world objects are often complicated and cannot, in general, be described as a finite sum as in (5). Hence, the formalism of (3) is fundamentally flawed and the usual way to circumvent this fallacy is to consider approximated versions of many of the expressions above.
In this work, we use a more accurate representation of the imaging system, which avoids introducing approximations. Mathematically, we consider a mapping ℒ between the object f and a new function u(A): u = ℒf, (10) in which u belongs to a Hilbert space and corresponds to the function u(A) of the continuous variable A (to be discussed below). Because variable A is allowed to take on continuous values, the mapping above, characterized by the operator ℒ, is a continuous-to-continuous (CC) mapping.

Photon-processing detectors
The CC systems discussed here do not collect data as sets of pixel counts. Instead, they make use of photon-processing detectors (Caucci et al 2013, Caucci et al 2018 to collect data. We define a photon-processing detector as any detector that (Caucci et al 2015b): • uses a gain mechanism and multiple sensors to obtain multiple measurements for a single absorbed photon; • collects data from all sensors that respond to each event at full precision; • uses the sensor data and maximum-likelihood estimation to estimate a vector of parameters (or 'attributes'); • stores all ML estimates of attributes in a list, at full precision.
For example, in the case of SPECT imaging with Anger cameras (Anger 1958, Peterson andFurenlid 2011), a photon-processing detector utilizes all available photomultiplier tube (PMT) signals that have been properly sampled and converted to digital numbers to estimate attributes (such as 2D or 3D position within the camera's crystal, photon energy, and so on) of gamma-ray photons interacting with the camera. These maximum-likelihood estimates are not restricted to belong to finite sets of values (as in the case of DD systems), but they can take on any-potentially continuous-value.
The vector of attributes for the jth detected photon is denoted as A j , for j = 1, …, J, where J is the total number of photon detected during some time interval of length T. Vector A j is a vector with K components, and it can consist of any attribute we can estimate from the detector outputs, such as the 2D or 3D position of interaction of the photon with the crystal, the energy deposited by the photon in the crystal, the time of arrival of the photon and the direction of propagation. We use the 'hat' symbol to mean that A j is calculated by means of an estimation algorithm (such as an ML estimation algorithm (Dempster et al 1977)) from noisy detector outputs. The noise in the detector outputs is random, which makes A j a random quantity as well.
The list of attributes A = A 1 , …, A J can mathematically be represented as a random point process (Lehovich 2005, Caucci and Barrett 2012, Caucci 2012, Jha 2013) u: where δ(…) is a K-dimensional Dirac delta function. Any imaging system that stores the list of attribute vectors A = A 1 , …, A J (as opposed to pixel counts g 1 , …, g M ) is said to be storing the data in list-mode format. A schematic representation of a SPECT photonprocessing imaging system is shown in figure 1.
The formalism we defined here is not restricted to SPECT imaging with Anger cameras. In fact, applications of photon-processing detectors to imaging with charged particles (such as electrons or alpha and beta particles) as well as imaging in the visible range have recently appeared in the literature (Ding et al 2014a, 2014b, Barrett et al 2017.

The imaging operator ℒ
Though both ℋ and ℒ operators are linear and operate on the same object, the CC operator ℒ is fundamentally different from the CD operator ℋ. A CD operator necessarily has an infinite-dimensional null space (Barrett and Myers 2004); the system maps a vector in an infinite-dimensional Hilbert space to a finite set of numbers, so there is an infinite set of vectors f null that yield no data at all: ℋf null = 0. Because real objects almost always have null components for any CD system, the object itself cannot be recovered from the discrete data, and even simple linear functionals such as integrals of the object over voxels cannot be estimated.
On the other hand, a CC operator such as ℒ may have a null space, but for the setup of section 3 and if K ⩾ 3, the null space is not demanded by dimensionality considerations.
Thus, it is possible that some-conceivably all-of the estimability problems inherent in the world of conventional digital imaging may disappear . Notice that replacing the operator ℋ with ℒ does not, in general, ensure that an inverse problem that was ill-posed (Natterer 1986, Bertero andBoccacci 1998) in the CD case becomes wellposed when a CC operator is used. As an example, although the Radon transform maps a continuous function to another continuous function, the problem its inverse solves is illposed.
To derive an explicit expression for the kernel of the CC operator ℒ, we begin by considering the probability density function for attribute vector A for a given object f; we will denote such a quantity as pr(A | f), in which the vertical line '\vert' denotes conditioning with respect to a known quantity. In principle, one could use the Boltzmann transport equation to relate the source distribution f to the spectral photon radiance (Barrett and Myers 2004) at the detector's face, and then express pr(A | f) as a blurred version of the properly scaled and marginalized spectral photon radiance (Caucci et al 2015b), with the blurring accounting for the inaccuracies in estimating A from the sensor data. In our treatment, however, we follow a derivation that explicitly takes into account the construction of the random point process u in (11) from list-mode data A. By the properties of conditional probability, we have (Caucci and Barrett 2012): in which pr(A | r, c) is the probability density function for measuring A under the assumption that a gamma-ray photon was emitted from point r ∈ S. The probability density function pr(r | f) is given by in which s(r) is the system's sensitivity function, which is defined here as the probability for a photon emitted at point r to produce any signal. In other words: Notice that pr(r | f) as defined above does not describe the probability density function for the emission point for an object f. In fact, pr(r | f) is the probability density function for an emission from point r and an object f to produce a signal in the detector. A photon emitted from point r can just be blocked by the pinhole or, if it does not, it can go through the scintillator in the camera and produce no signal. The probability density function pr(r | f) as defined in (13) not only takes into account the distribution of the emission points (defined by f) but also the probability of the emitted gamma-ray photon to produce a signal (defined by the system's geometry).
In principle, s(r) can be calculated (either analytically or using numerical methods, such as Monte Carlo estimation (Fishman 1996)), at any point from knowledge of the system geometry, or it can be estimated by moving a point-like source across the system's field of view and then by taking s(r) = αJ r in which J r is the number of attribute vectors collected during a sufficiently long time interval and when the point-like source was located at point r, and α is a constant that makes s(r) the probability for a photon emitted at point r to produce any signal.
The probability density function pr(A | r) can be further rewritten to explicitly take into account estimation of A from noisy signals (e.g. discretized PMT outputs produced by the camera in the case of SPECT) and propagation of light from point r ∈ S to the camera's entrance face: For any given object f, the average u of u defined in (11) is in which Pr(J | f, T) is the probability of measuring exactly J attribute vectors while imaging the object f for time T, and J(f, T ) is the average of J under the same circumstances. It turns out (Lehovich 2005, Caucci andBarrett 2012): Upon substitution in the last expression in (16), we obtain

Author Manuscript
Author Manuscript Author Manuscript Author Manuscript which shows that the kernel of ℒ is: [ℒ](A, r) = T pr(A | r)s(r) . (19)

System description
The imaging system we used for our study, called FastSPECT II, was developed at the Center for Gamma-Ray Imaging, University of Arizona. FastSPECT II is a small-animal SPECT imager built with modular scintillation cameras (Milster 1987, Milster et al 1990 and list-mode data-acquisition electronics. The modular gamma-ray scintillation camera designed for FastSPECT II comprises a 5 mm thick NaI(Tl) scintillation crystal, a 15 mm thick quartz light guide, and a 3 × 3 array of 1.5 inch diameter end-on photomultiplier tubes. Each camera has an input face measuring about 120 × 120 mm 2 (Furenlid et al 2004, Chen 2006. Each entry in the data list collected by a gamma-ray camera corresponds to a detected scintillation event and consists of a camera identifier, the nine signal values present in the 3 × 3 array of photomultipliers, and a time stamp. The cameras are stationary and arranged as two rings of eight on opposite sides of a pair of central plates (Furenlid et al 2004, Chen 2006. As shown in figure 2, FastSPECT II attains good sensitivity over a large volumeabout 57 cm 3 -well suited for small-animal imaging. As reported in Chen (2006), the sensitivity of FastSPECT II at the center of the FOV is about 267 counts per second (cps) per each million gamma rays emitted per second. Hence, s(r center ) ≈ 2.67 · 10 −4 .

Estimation of photon attributes
We used PMT signals p 1 , …, p I (I = 9 for FastSPECT II) from each camera to perform ML estimation of photon attributes. The estimated attributes consisted of the 2D position of interaction between a gamma-ray photon and the camera's entrance face. Hence, A = (X, Y ). This ML estimate is defined as in which Pr(p 1 , …, p I | A) is the probability of measuring PMT signals p 1 , …, p I when a gamma-ray interaction occurs at point A = (X, Y) and the 'arg max A ' notation in (20) denotes the value of A that maximizes Pr(p 1 , …, p I | A). The quantity L(A; p 1 , …, p I ) = Pr(p 1 , …, p I | A), in which p 1 , …, p I are assumed fixed is called 'likelihood' (Barrett and Myers 2004). Its logarithm, l A; p 1 , …, p I = log L A; p 1 , …, p I , is often referred to as 'loglikelihood'.
Calculating the probability Pr(p 1 , …, p I | A) for a given A requires knowledge of the forward model and its statistical properties. A calibration procedure is often used to obtain information about the statistics of p 1 , …, p I for points A on a regular grid (Peterson andFurenlid 2011, Chen 2006). Other approaches, including Monte Carlo estimation, have been considered (Hunter et al 2013). We will assume that p 1 , …, p I conditioned on A are independent and follow Poisson statistics (Hunter 2007; under such hypotheses, we have: in which p i (A) is often called mean detector response function (MDRF). For any i = 1, …, I, the MDRF p i (A) is the mean of p i under the assumption that a gamma-ray interaction has occurred at point A. Although p i (A) is usually estimated or measured for a finite set of points (typically arranged in a uniform pattern), in practice it is a smooth function of A and it is therefore possible to interpolate these samples so that p i (A) can be evaluated for any A.
After an ML estimate A is calculated as in (20), the likelihood L A; p 1 , …, p I is compared to a threshold L 0 (A). If L A; p 1 , …, p I < L 0 (A), then estimate A is discarded. This technique, originally introduced in Milster et al (1990), is called 'likelihood thresholding' and it is used to ignore estimates A that poorly match the statistical model of (21). This is often due to scattering of gamma-ray photons and subsequent loss of photon energy (Chen 1997). Calculation of the position-dependent likelihood threshold L 0 (A) has been detailed in Chen (2006). Briefly, for any fixed A 0 , raw PMT data p 1 , …, p I collected in the calibration step are reprocessed to create histograms of log-likelihood l A 0 ; p 1 , …, p I . Sample mean and variance of the log-likelihood values are calculated. The algorithm iteratively discards values of log-likelihood that are outside two standard deviations around the sample mean, and recomputes the sample mean and variance until the mean converges to a fixed value with an error less than one decimal place (Chen 2006). A value l 0 A 0 is set at four standard deviations below the sample mean. Finally, the likelihood threshold L 0 (A 0 ) is set as Reduced spatial resolution has been observed near the edges and corners of the crystal , Moore 2011) with a larger-than-expected number of position estimates clustering near the edges of the crystal; an example is shown in figure 4. This is likely due to a combination of effects, including non-monotonic signal variation as a function of position in these areas. Methods to cope with this problem have been proposed; they generally rely on the estimation of the total energy the gamma-ray photon has deposited in the crystal. In this work, we simply decided to ignore all estimates that are less than 4 mm away from any of the crystal edges.
Graphics processing unit (GPU) code was developed (Caucci and Furenlid 2015) to process PMT data and perform the 2D ML estimation of position of interaction in (20) for each detected gamma-ray photon. Our algorithm for ML estimation of position of interaction used an iterative contracting-grid approach (Furenlid et al 2005). We took full advantage of the texture mapping unit of the GPU devices to interpolate the PMT mean detector response function (MDRF) of each camera. The CUDA cubic B-spline interpolation library (Ruijters et al 2008) was used to perform on-the-fly interpolation of MDRF samples with low-degree spline functions (Cox 1972, Schumaker 1981). The main advantage of this contracting-grid approach is that 2D estimates can be performed with arbitrary precision by simply running more iterations of the contracting-grid search algorithm. Other customizations, such as contracting factor after each iterations and the size of the grid, are possible (Caucci andFurenlid 2015, Caucci et al 2018). Moreover, the same algorithm is well suited for the estimation of other photon attributes provided that the proper statistical model Pr(p 1 , …, p I | A) is available.

The list-mode MLEM reconstruction algorithm
A list-mode variant Barrett 1998, Caucci 2012) of the maximum-likelihood expectation-maximization (MLEM) (Shepp and Vardi 1982) algorithm has been implemented on GPUs to perform reconstruction from list-mode data A in which c = 1, …, C denotes the camera and J c is the number of recorded events for camera c. In the formalism discussed in this paper, the list-mode MLEM algorithm is written as in which T is the exposure time, f n (l) is the estimate of f n after l iterations and pr A (c) | r, c is the quantity calculated in (15) in which we have extended our notation to make it clear that attribute vector A (c) corresponds to a recorded event in camera c. The iterative expression above requires an initial guess f n (0) for n = 1, …, N, which we will assume to be a nonnegative constant; in other words, f n (0) = f 0 > 0 for all n. Derivations of the list-mode MLEM algorithm are available in Levkovitz et al (2001) and Caucci (2012). Some stopping criteria can be found in Gaitanis et al (2010).
Critical to our treatment is the evaluation of pr(A | r, c). We followed an approach similar to the one delineated in Furenlid et al (2004) and Chen (2006), in which calibration data was collected and a Gaussian-shaped function was fit to the calibration events. Details can be found in appendix. With our approach, evaluation of pr(A | r, c) for fixed camera c and r only requires six coefficients.
System calibration for FastSPECT II entails moving a point-like radioactive source across the field of view to collect list-mode data from which estimate fitting coefficients (Chen 2006). As this process is often time-consuming, source decay must usually be taken into account. Moreover, one can only consider a finite set of positions at which collect calibration data. The alternative is to derive a complete analytical model for pr(A | r, c), but this would require precise knowledge of camera responses and their positions, as well as knowledge of position of all pinholes and their shapes. For this reason, we opted for estimating fitting coefficients for points r belonging to a 3D grid of points Γ, and whenever fitting coefficients at r not in Γ are needed, they are interpolated (as discussed in appendix), so that pr(A | r, c) can be evaluated quickly for any attribute vector A and point r and for any cam era c.
To summarize, we use photomultiplier tube signals to estimate the components of attribute vectors A with arbitrary precision. Our ML algorithm implemented on GPUs produces estimates A in real-time during the acquisition scan. Event attributes are collected into lists , …, A J c (c) (one list for each camera c = 1, …, C), which are then used in the iterative expression in (23) to obtain estimates f (l) . By properly fitting measured calibration data, it is possible to evaluate pr(A | r, c) for any A and r and for all the cameras c = 1, …, C.
By (19), this amounts to directly evaluate the kernel of the continuous-to-continuous operator ℒ, which we use to obtain an estimate f to f that satisfies the forward model of (18).

Phantom studies
In a first set of reconstructions, we imaged a Jaszczak-like phantom (shown in figure 3) whose bores were filled up with an aqueous solution of 99m Tc-based sodium pertechnetate. The metastable nuclear isomer 99m Tc decays to 99 Tc with the emission of a gamma ray with a photon energy of 140 keV.
A total of about 13.55 · 10 6 photons were collected while imaging the phantom with FastSPECT II. Of these, around 32.41% were discarded during likelihood thresholding (see section 5.2) or because estimated positions were less than 4 mm away from any of the crystal edges (see figure 4 for examples). This overall process left around 9.16 · 10 6 counts available for reconstruction.
List-mode data were reconstructed with the algorithm in (23) and at points r n arranged on a non-constant spacing grid. More specifically, we considered a coarse grid with a specific region of 2× finer spacing. Post-reconstruction nearest-neighbor resampling of the coarse regions was used to ease visualization without altering the process of image reconstruction. Figure 5(a) shows a cross-section of the reconstructed data in which the region of finer spacing was emphasized with a dashed box. (Notice that during preparation of the phantom, some of the bores were not properly filled with the radiotracer.) Some line profiles through the post-processed reconstructed data are shown as figure 5(b).

Reconstructions of KPC mouse
Imaging studies of FastSPECT II were performed in a genetically-modified K-ras LSL.G12D/+ ; p53 R172H/+ ; PdxCre mouse model (Hingorani et al 2005, Westphalen and  ultimately progress over months to overt carcinoma with extensive stromal desmoplasia, similar to the most common morphology of PDAC observed in humans. 99m Tc-labeled 3P 4 -RGD 2 dimer ( 99m Tc-3P 4 -RGD 2 ), which binds to α v β 3 integrin and indicates angiogenic activity in tumor stroma, was selected to image PDAC angiogenesis (Wang et al 2009). Integrin α v β 3 is a cell-surface receptor with an exposed arginine-glycineaspartate (RGD) binding site for a variety of extracellular matrix (ECM) proteins (Davis 1992).
Three hours after intravenous tracer injection (1.0 mCi of 99m Tc-3P 4 -RGD 2 ) into the mouse, FastSPECT II imaging was performed. A total of 1.71 × 10 6 events across the 16 cameras were collected, 34.21% of which were classified as scattered by likelihood thresholding and, therefore, rejected. This left about 1.13 × 10 6 events available for reconstruction. Highresolution reconstructions (figure 6) showed two lesions with high radioactive uptake (arrows) in the upper right abdomen (left panel), which were consistent with the tumors observed in the mouse pancreas by postmortem examination (right panel).
All the reconstructions were obtained after 50 iterations of the method in (23). In other words, we assumed that the final reconstruction was f (50) . Renderings of f (50) were obtained with AMIDE 5 .
To test the main point of this paper-that reconstruction methods based on continuous-tocontinuous operators provide an advantage over algorithms that assume a discrete-to-discrete model for the imaging system-we considered and compared different sets of reconstructions in which: i. the number of iterations of the contracting-grid algorithm for the estimation of attribute vectors A was increased, thus allowing estimates A to take on values approximating a continuous domain; or, ii. finer and finer system interpolation (see appendix A.2) was performed, thus allowing reconstructions on arbitrary sets of points (which are taken as approximations of continuous functions).
Reconstruction Visual assessment of figure 7 shows that as the size of the voxel array used in the reconstruction increases, resolution improves and tumors are well resolved. Similarly, as the event estimates A = X, Y are allowed to take value on a continuous domain, image resolution improves and the tumors (see left panel in figure 6) are better resolved. To better 5 http://amide.sourceforge.net/index.html appreciate this fact, we report in the top row of figure 8 the reconstructions we obtained for a voxel grid of size 417 × 321 × 321 and different number of contracting-grid iterations. In the bottom row, we reported the absolute value of the difference between pairs of reconstructions, which were displayed as images. Notice the different color scales used in the images at the bottom of each figure.
We ran our codes on a server equipped with four Intel® Xeon® CPU E5-2698 processors, 256 GB of RAM and eight NVIDIA® Tesla® P100-SXM2 GPU accelerator cards. The software configuration included Linux openSUSE Leap 15.0, NVIDIA® CUDA SDK release 10.0, and GNU C compiler 7.3.1. Our codes heavily relied on GPU processing for the estimation of photon attributes (see section 5.2) and the calculation in (23). We designed our implementations to fully utilize all the GPU accelerators available. Estimation of photon attributes for all the 1.71 × 10 6 events took less than half of a second, no matter how many iterations of the contracting-grid algorithm were performed. In fact, the algorithm just needs a few iterations to 'localize' an event to a small region of the search space. Subsequent iterations are much faster, as they rely on cached data for the evaluation of the MDRF fitting splines. On the other hand, reconstruction time is highly affected by the image-space granularity and the number of iterations during reconstruction. Table 1 reports running time for 50 iterations and different sizes of the voxel array. These results show that running time is proportional to the number of reconstruction points.

Future work
One possible follow-up research in this field might consider objective assessment of image quality  for continuous-to-continuous imaging systems. For example, a task-based figure-of-merit can be used to compare a continuous-to-continuous imaging system with an imaging system that use the traditional discrete-to-discrete model. In the case of a detection problem, such as detection of tumor necrosis, a meaningful figure-of-merit is the area under the receiver operating characteristic curve (Van Trees 1968). For an estimation problem, one could derive an expression for the list-mode Fisher information matrix (Fisher 1925), which also provides a lower bound on the variance any unbiased estimator (Van Trees 1968). Some preliminary work has been done in the context of estimating activity within a certain region of interest, and it was shown that continuous-tocontinuous systems provided more accurate tracer uptake estimates compared to continuousto-discrete systems .
Although this paper has been limited to imaging with gamma-ray photons, the same methodology can be applied to imaging with charged particles, such as alpha and beta particles. Some work  has already been done in this direction, and it was shown that estimating additional parameters (such as direction and/or particle's residual energy) besides positions, does provide enormous benefits in terms of reduced null space.
Other methods for reconstruction of a continuous functions have been introduced in the past decade. Among them, is the Backus-Gilbert method (Backus andGilbert 1968, Kirsch et al 1988), originally introduced for the estimation of geological models from a finite set of noisy measurements. A more recent method  for reconstruction of continuous functions uses the singular value decomposition of the imaging operator. A task-based comparison between these methods and the reconstruction approach developed in this paper will further shed some light on the advantages of image reconstruction algorithms based on continuous-to-continuous model.

Conclusions
We presented a new approach to image reconstruction for single-photon emission computed tomography (SPECT) that does not use a discrete-to-discrete representation of the imaging system. In fact, our approach models the imaging system via a continuous-to-continuous operator that maps the object (i.e. a function of continuous variables) to a function that takes as input an event attribute, which is assumed to take values on a continuous set. The reconstruction algorithm we propose is based on maximum-likelihood. It uses measured calibration data, and we showed how these calibration data can be interpolated to arbitrary precision.
We tested our approach on real data, which we collected with FastSPECT II. We used 99m Tc 3P 4 -RGD 2 peptide to image a genetically-modified KPC mouse that developed pancreatic ductal adenocarcinoma in two different locations. By varying key estimation/interpolation parameters in our code, we were able to attain reconstructions on finer and finer grids, thus approximating reconstructions on a continuous domain. We did not perform a formal assessment of image quality. Instead, we relied on visual comparison of reconstructed volume renderings, which showed increased details as the limit of a continuous-tocontinuous system was approached.
Our algorithms were implemented for a GPU architecture, which offers enormous parallelcomputing capabilities as well as large amounts of fast memory. In particular, we took advantage of texture interpolation to rapidly retrieve calibration data. The algorithm we proposed to interpolate the system response avoids conditional branching, which are very time consuming on a GPU architecture. Moreover, our interpolation and reconstruction algorithm provides greater flexibility as they allow reconstructions on an arbitrary domain, not necessarily restricted to a uniform grid of points. Although GPUs were shown to be well suited for implementation of our parallel codes for maximum-likelihood estimation, other parallel architecture (e. g., field-programmable gate arrays, or FPGAs) might be considered as well.
Calculation of these fitting coefficients at any point r and for any camera c requires two separate steps: i. system calibration to estimate fitting coefficients for r belonging to a uniform grid; ii. interpolation of the fitting coefficients obtained in the previous step.
The Gaussian model of (A.1) was chosen for its simplicity. Depending on the system geometry, other choices are possible. This includes sums of Gaussian functions (for the case of a multi-pinhole imaging system), or a more general model for p(X, Y ; r, c) in (A.1) based on, fore example, the 2D Pearson type VII distribution (Kotz 1975).

A.1. System calibration
In a typical system calibration procedure, a point source is moved across the field of view by means of computer-controlled motorized stages. The position of the point source will be denoted as r n , in which we defined n = (n x , n y , n z ) for n x = 1, …,N x , n y = 1, …, N y and n z = 1, …, N z . Points r n are arranged on a uniform grid, which we will denote as Γ, and adjacent points are separated by a known distance δ.
During calibration, list-mode data A n (c) are collected at each of the N x N y N z points positions r n ∈ Γ and for each of the C cameras. We will assume that A n (c) contains J n (c) attribute vectors . Notice that A j (and, consequently, X j and Y j ) do depend on n and c, but we have discarded these dependencies for notational convenience.
We used our GPU-based maximum-likelihood code to estimate attribute vectors Plots of a r n , c , μ X r n , c , μ Y r n , c , σ X, X 2 r n , c , σ Y , Y 2 r n , c and ρ r n , c for one of the cameras and calculated according to (A.4)-(A.11). These plots show that estimated coefficients vary smoothly as the point r n is moved over the field of view (position along the axes are in units of millimeters).

A.2. Interpolation
The system calibration process outlined above is often performed on a coarse grid Γ of points r n . Time constraints (due to radioactive decay in the point source used for calibration) along with the necessity of collecting a few thousand events at each point r n limit the separation δ between two adjacent points to around 2 mm (Chen 2006). With a typical shortlived isotope used in medical imaging (e.g. 99m Tc), collecting enough calibration data at each point of a finer grid while covering the entire field of view (see figure 2) is infeasible. For this reason, we have developed a method for fast calculation of the fitting coefficients μ X r, c , μ Y r, c , σ X, X 2 r, c , σ Y , Y 2 r, c and σ X, Y 2 r, c any point r and for any camera c from quantities estimated in (A.4)-(A.11). Figure A1 reports sample plots of a r n , c , μ X r n , c , μ Y r n , c , σ X, X 2 r n , c , σ Y , Y 2 r n , c and ρ r n , c as estimated according to (A.4)-(A.11) and along perpendicular planes intersection at the center of the field of view. Interpolation of v n x , n y , n z .
The first step in our method is to diagonalize the following matrix: K r n , c = σ X, X 2 r n , c σ X, Y 2 r n , c σ X, Y 2 r n , c σ Y , Y 2 r n , c (A.12) so that: K r n , c = R ϕ r n , c ΛR ϕ r n , c R ϕ r n , c = cos ϕ r n , c −sin ϕ r n , c sin ϕ r n , c cos ϕ r n , c (A.14) and Λ = λ X r n , c 0 0 λ Y r n , c . (A.15) Quantities ϕ(r n , c), λ X (r n , c) and λ Y (r n , c) are calculated as follows: ϕ r n , c = 1 2 arctan 2σ X, Y 2 r n , c σ X, X 2 r n , c − σ Y , Y 2 r n , c , (A.16) λ X r n , c = 1 2 σ X, X 2 r n , c + σ Y , Y 2 r n , c + Δ , (A.17) λ Y r n , c = 1 2 σ X, X 2 r n , c + σ Y , Y 2 r n , c − Δ , (A.18) in which, Δ = σ X, X 2 r n , c − σ Y , Y 2 r n , c 2 + 4 ρ r n , c σ X, X r n , c σ Y , Y r n , c 2 1/2 . (A.19) Interpolation was performed according to the diagram show in figure A2, in which we used the placeholder 'v' to denote one among a, μ X , μ Y , ϕ, λ X or λ Y .
For simplicity, we assumed that any interpolated value was indexed by a subscript vector (n x , n y , n z ) in which at least one in n x , n y or n z is non-integer. In figure A2 we used the 'floor' (⌊…⌋) and 'ceiling' (⌈…⌉) notation to calculate subscripts for the non-interpolated values of v that 'surrounds' v n x , n y , n z and are used in the interpolation. These noninterpolated values, which by construction corresponds to points in Γ, are marked in figure  A2 with squares. Weighted two-value average was used to interpolate two values at a time as follows: v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n x − n x , v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n x − n x , v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n x − n x , v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n x − n x , v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n y − n y , v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n y − n y , v n x , n y , n z = W v n x , n y , n z , v n x , n y , n z , n z − n z , (A.20) in which W(v 1 , v 2 , γ) is the weighting two-value average function defined as: W v 1 , v 2 , γ = (1 − γ)v 1 −1 + γv 2 −1 −1 if v = λ X or v = λ Y , (1 − γ)v 1 + γv 2 otherwise. (A.21) In other words, W(v 1 , v 2 , γ) is the two-value weighted arithmetic mean between v 1 and v 2 (with weights 1 -γ and γ, respectively) if v = a, v = μ X , v = μ Y or v = ϕ; while W(v 1 , v 2 , γ) is the two-value weighted harmonic v = λ X or v = γ Y .
Notice that the scheme in (A.20) operates on one dimension at a time; this is emphasized in figure A2 by solid lines connecting the pairs of points used in (A.20). Furthermore, for any integer number n, we have ⌊n⌋ = ⌈n⌉ = n, so that the scheme in (A.20) along with the definition of W(v 1 , v 2 , γ) in (A.21) effectively performs no interpolation along the dimension(s) for which the corresponding index(es) in (n x , n y , n z ) are integer. Although this does perform unnecessary computation when one of more indeces in (n x , n y , n z ) are integer, in practice, it is well suited for parallel implementation due to the absence of conditional branching. Example of a photon-processing imaging system (adapted from Caucci et al (2013)).  Plot of the sensitivity across three planes, showing almost uniform sensitivity over a large portion of the system's field of view. The volume for which the sensitivity exceeds 5% its maximum value is about 57 cm 3 .   Scatter plots of raw estimated positions (after likelihood thresholding but before discarding estimates less than 4 mm away from any of the crystal edges) on the camera's face obtained while imaging the phantom. Notice how some of the estimates cluster near the edges and corners of the crystal.