Isotope independent determination of PET / CT modulation transfer functions from phantom measurements on spheres

PURPOSE
A PET/CT system's imaging capabilities are best described by its point spread function (PSF) in the spatial domain or equivalently by its modulation transfer function (MTF) in the spatial frequency domain. Knowing PSFs or MTFs is a prerequisite for many numerical methods attempting to improve resolution and to reduce the partial volume effect. In PET/CT, the observed PSF is a convolution of the system's intrinsic imaging capabilities including image reconstruction (PSF0) and the positron range function (PRF) of the imaged β+ emitting isotope. A PRF describes the non-Gaussian distribution of β+ annihilation events around a hypothetical point source. The main aim was to introduce a new method for determining a PET/CT system's intrinsic MTF (MTF0) from phantom measurements of hot spheres independently of the β+ emitting isotope used for image acquisition. Secondary aim was to examine non-Gaussian and nonlinear MTFs of a modern iterative reconstruction algorithm.


METHODS
PET/CT images of seven phantom spheres with volumes ranging from 0.25 to 16 ml and filled either with 18F or with 68Ga were acquired and reconstructed using filtered back projection (FBP). MTFs were modeled with linear splines. The spline fit iteratively minimized the mean squared error between the acquired PET/CT image and a convolution of the thereof derived PSF with a numerical representation of the imaged hot phantom sphere. For determining MTF0, the numerical sphere representations were convolved with a PRF, simulating a fill with either 18F or 68Ga. The MTFs determined by this so-called MTF fit method were compared with MTFs derived from point source measurements and also compared with MTFs derived with a previously published PSF fit method. The MTF fit method was additionally applied to images reconstructed by a vendor iterative algorithm with PSF recovery (Siemens TrueX).


RESULTS
The MTF fit method was able to determine 18F and 68Ga dependent MTFs and MTF0 from FBP reconstructed images. Root-mean-square deviation between fit determined MTFs and point source determined MTFs ranged from 0.023 to 0.039. MTFs from Siemens TrueX reconstructions varied with size of the imaged sphere.


CONCLUSIONS
MTF0 can be determined regardless of the imaged isotope, when using existing PRF models for the MTF fit method presented. The method proves that modern iterative PET/CT reconstruction algorithms have nonlinear imaging properties. This behaviour is not accessible by point source measurements. MTFs resulting from these clinically applied algorithms need to be estimated from objects of similar geometry to those intended for clinical imaging.


INTRODUCTION
In previous work, 1 we had introduced the transconvolution (TC) method for PET/CT image normalization and to improve comparability of quantitative PET measurements performed on different PET/CT systems.TC offers a comprehensive way of mathematically transforming images acquired on different PET/CT systems to images with uniform imaging characteristics, as if they had been acquired by one single virtual PET/CT system (vPET).A prerequisite for the application of the TC method is the knowledge of the point spread functions (PSFs) of all participating PET/CT systems or their associated modulation transfer functions (MTFs).After measuring system-specific PSFs or MTFs, a standardized vPET-PSF and/or vPET-MTF is deconvolved with these PSFs, resulting in specific transconvolution functions (TCFs) for every participating PET/CT system.The TCFs itself are then used to recast (transconvolve) images acquired on these systems into vPET images, allowing for comparable quantification of PET/CT measurements performed on multiple PET/CT systems.
Even though frequently used as such, 2 PSFs in PET/CT tend not to be just isotropic Gaussians, [3][4][5] and simplifying a PSF as such does not meet the requirements of PSFs in the context of TC. 6 One process giving rise to non-Gaussian PSFs lies in the very foundations of how positron ( β + ) decay is detected: After emission, positrons have a nuclide-specific reach in tissue before they thermalize and annihilate together with an electron. 7The energy spectrum of β + emissions and the positron transport up to its maximal range introduce an uncertainty about a positron's point of origin.The result is a source-centered, radial annihilation distribution (positron range function, PRF) with a non-Gaussian cusplike shape. 8he PRF broadens spatial resolution in PET imaging, 8 and different PRFs will therefore give rise to different spatial resolutions when imaging different β + emitters with the same PET/CT system.Any PET/CT system's PSF or MTF is therefore a combination of its intrinsic imaging capabilities (PSF 0 or MTF 0 ) and the PRF of that particular isotope.Both, PRF and PSF 0 , are origin of the so-called partial volume effect (PVE) in PET/CT imaging.The PVE results from signal spillover between adjacent image elements, and hampers quantitative measurement.
For reasons of good reproducibility and ease of handling, measurements of a PET/CT-system's PSF are ideally conducted with a long-lived 68 Ge/ 68 Ga solid state phantom.However, the most abundantly used isotope in clinical PET/CT imaging is 18 F with β + energy different to 68 Ga. 18F emits its majority of β + with an average energy of 249.8 keV and an endpoint energy of 633.5 keV, 9 whereas the β + of 68 Ga leaves the nucleus with a considerably higher average energy of 836.02 keV and an endpoint energy of 1899.1 keV. 10 The resultant different PRFs of 68 Ga and 18 F give rise to different PSFs, even if the measurements are carried out on the same PET/CT system.The problem is intensified when imaging inhomogeneous objects with varying electron density, where spatially variant PRFs need to be taken into account.A method to determine PSF 0 without any associated positron range effects will thus allow the use of 68 Ge/ 68 Ga solid state phantoms for a TC based normalization of 18 F measurements.][13] Our previously suggested heuristic approach for modeling PSFs in the spatial domain, with parameters for width and "wings," 1 was introduced for independent handling of a PSF's full width at half maximum (FWHM) and, more importantly, its full width at tenth maximum (FWTM).This PSF model is hereinafter referred to as the "wing fit" function.It accounts for a PSF's non-Gaussian reach in 3D space, i.e., the wings, by convolving a Gaussian with an exponential decay distribution.The Gaussian component describes the central part of the modeled PSF, while the exponential density function describes its wings.Designed for modeling anisotropic and non-Gaussian PSFs, the wing fit function might still not have enough degrees of freedom for capturing the contributions of the positron annihilation distribution to a final PSF.It is certainly incapable of describing PSFs of modern nonlinear reconstruction methods with their potentially nonmonotonic modulation transfer.A more versatile PSF determination method is required, which is able to describe PSFs regardless of the isotope and image reconstruction method used.This would enable TC to normalize images, where aforementioned PSF determination was done with a different β + emitter than the images itself.
By the virtue of Fourier transform, a MTF is a PSF's representation in spatial frequency space: The MTF's tail, describing high spatial frequencies, is responsible for the imaging system's spatial resolution, whereas the wings of a PSF represent the PVE experienced by that system.Instead of modeling a PSF in position space as previously described, 1 it is therefore possible to model the respective MTF in spatial frequency space.This can be done numerically using splines with multiple interpolation points or constrains to capture a MTF's shape.
Determining PSF 0 becomes possible when the numerical representations of objects used during the iterative fitting process are modeled together with their respective PRFs.As such, the intrinsic PSF 0 can be calculated regardless of the isotope used for imaging and a direct deconvolution with the PRF can be omitted.
Main aim of the current work was to introduce a new method for determining a PET/CT system's PSF 0 or MTF 0 from measurements of hot phantom spheres, independently of the imaged β + emitting isotope.The new method models MTFs iteratively in the spatial frequency domain using linear splines, while the thereof derived PSF is being convolved with numerical representations of hot phantom objects.Secondary aim was to examine the non-Gaussian and nonlinear MTFs of a modern nonlinear iterative reconstruction algorithm.
For this, we compared the new method of MTF determination with the wing fit method, 1,6 using phantom measurements of hot spherical inserts filled alternatively with 18 F or 68 Ga.Our PET/CT system's PSF 0 was successfully determined by utilizing a literature derived model for 18 F and 68 Ga PRFs.The new method's ability of iteratively deconvolving MTFs derived from a nonlinear reconstruction method (Siemens  TrueX) was examined, too.

2.A. Background and theory
Usually, PSFs of a PET/CT system are determined with a given isotope.Common isotopes include, for example, 18 F (PSF F ) or 68 Ga (PSF Ga ).Omitting noise considerations, the image can be understood as the convolution (•) of an object with this PSF, 14 img In PET/CT PSF F or PSF Ga are measured by placing line or point sources in the FOV 15 or are estimated from the image of a known object. 1,16In case the radioactive sources are embedded in a positron absorbing medium, the measured PSF will always be a convolution of the imaging systems intrinsic PSF (PSF 0 ) and the PRF of the utilized isotope, An imaging system's "true" PSF 0 is assumed to emanate from a combination of detector geometry, detector Compton or, in the case of 68 Ga, we get Above, obj 0 stands for the spatial distribution of β + emitter in an object, without its specific annihilation distribution.9][20][21] PRFs, being isotropic, are usually provided as functions of radius r.In this work, the PRF is separated from the PET system's PSF in order to arrive at PSF 0 .The PRF therefore becomes a characteristic of the object, which is assumed to be filled with a specific isotope (obj F or obj Ga ) and possesses the respective annihilation distribution.In the following, we used 18 F with its PRF (PRF F ) and 68 Ga with its PRF (PRF Ga ) to simulate objects filled with these β + emitters, Giving detailed information about the two isotopes of interest, 18 F and 68 Ga, the two most useful PRF models were taken from the wok of Jødal et al. 19 and from Cal-González et al. 20 From those we chose to incorporate the latter in our MTF seeking method with its analytical expression for the radial annihilation distribution function, The parameters a, r 0 , ε, and n and their differing values for 18 F and 68 Ga were obtained from the original work of Cal-González et al. 20 Parameter C was used for scaling.Beyond r 0 , the maximum positron range, values of Eq. ( 5) were set to zero.In contrast to Gaussian PRF approximations, 17 this constitutes a closer equivalent to the expected positron behavior and eliminates nonphysical contribution of distant PRF image volume elements (voxels) to the imaged signal.
The radial annihilation distribution function g 3D (r) bins annihilation events in spherical shells of radius r.To arrive from g 3D (r) at a 3D PRF 3D (r) the following relationship between the two functions was employed: The divergence at the origin of Eqs. ( 5) and ( 6) was handled by defining PRF 3D (0) ≡ 0. For image power conservation, the volume integral of Eq. ( 6) was set to unity.

2.B.1. Phantom measurements
A water filled flask made of polytetrafluoroethylene (PTFE) with a cylindrical section was used to house a single "hot" insert.The total volume of the flask for filling water was 500 ml. Figure 1 shows the sheer plan of this phantom.The water in the flask contained no activity ("cold background").The hot insert was one out of seven spheres on rods with different volumes of 0.25, 0.5, 1, 2, 4, 8, and 16 ml.Sphere wall thickness was 1 mm.These spherical inserts were produced using rapid prototyping techniques from epoxy resin and were filled with 18 F or 68 Ga solution in water (200-400 kBq/ml).After filling, filler holes were heat sealed with thermoplastic Parafilm (American National Cam, Neenah, WI, USA).Single hot inserts were centered onto the cylindrical section of the flask phantom for individual measurements.The phantom was placed longitudinally in the PET/CT system with the spheres centered in the tomograph's field of view (FOV).
Point sources were fabricated by drenching polymeric adsorbent beads (Dowex Optipore L493, Dow Chemical Company, Midland MI, USA) of 0.30-0.85mm diameter with concentrated (200-400 MBq/ml) 18 F or 68 Ga isotope solution.Beads were left to dry and were than individually centered in 14 ml Falcon tubes filled with molten paraffin.
Differences in positron ranges that might have arisen from the use of paraffin instead of water as a β + absorbing material were found to be below the resolution limits of clinical PET/CT systems.Cal-González et al. state a r 0 of 2.30 mm in adipose tissue vs a r 0 of 2.00 mm in water. 20fter curing of the paraffin, the tubes were inserted one by one in the same water filled flask as used for the spherical inserts.Measurements were always performed with a single bead centered in the FOV following the same image acquisition protocol as for spherical inserts.
PET/CT images were acquired and reconstructed on a Biograph mCT 128 True-V (Siemens Medical Solutions USA, Knoxville, TN) in 3D mode, with CT based attenuation correction and using the time of flight (TOF) option.Except for the two smallest spheres, acquisitions were run until 10 6 counts/ml was reached.The two smallest spheres were imaged until the acquisition of 10 6 counts, and the point sources until the acquisition of 2 × 10 5 counts per source.
For linearity, 14,22 images were reconstructed with filtered back projection (FBP), using the Fourier rebinning (FORE) algorithm provided by the PET/CT manufacturer. 23,24TrueX was the nonlinear high resolution reconstruction method of choice and was run with one iteration and 21 subsets.Image data were reconstructed into a matrix of 512 × 512 volume elements (voxels) in 222 slices of size 1.59 × 1.59 × 1 mm.The anisotropic voxel sizing aided in validating the new method's ability to determine spatially anisotropic MTFs.A postreconstruction Gaussian filter with a full width at half maximum (FWHM) of 2 mm was applied onto the FBP data and a FWHM of 1 mm onto the TrueX data.Data were stored and handled in the digital imaging and communication in medicine (DICOM) standard. 25

2.B.2. Numerical analysis
Data analysis and all numerical calculations were performed on a HP Z620 Workstation (Palo Alto, CA) running Microsoft Windows 7 as the operating system.All data were analyzed using either our in-house developed Java based multiparadigm software framework or IGOR Pro, 6.22 (Wavemetrics, Portland, OR).
Measured point sources were integral normalized to unity and subsequently Fourier transformed to obtain the corresponding MTFs.Profiles at z = 0 were binned into a common radial MTF.Axial MTF profiles were taken from values at x = 0 and y = 0.
All acquired images were cropped into 64 × 64 × 64 voxel subimages around the centroid of the represented hot insert for further numerical processing.The size of the cropping defines, besides the maximal spatial frequency, number and size of available MTF frequency steps in Fourier space and thus defines available fit precision.

2.B.3. MTF fitting in Fourier space
PSFs were determined by iteratively fitting an approximation of the real image (img calc ) to the real image img F or img Ga .The approximated img calc was calculated by multiplying in spatial frequency space F {} a numerical representation of the imaged spheres with the respective MTF.This method is herein after referred to as the "MTF fit" method, For the inverse Fourier transform of Eq. ( 7), the imaginary part of the MTF was kept at zero, which gave an origincentered and symmetric PSF.
In the spatial frequency domain, MTFs can be fitted with any desirable precision up to cutoff frequency ξ c .Beyond ξ c the modulation transfer drops to zero and MTF fitting is not possible.An initial ξ c is given by the voxel size s x y z of the original PET/CT data [Eq.(8)] and constitutes the last interpolation point in the fitting procedure, Here, the original sampling of the acquired images reconstructed into voxels of 1.59 × 1.59 × 1 mm provided an upper limit for ξ c,r of 0.314 cycles/mm in the radial direction and axially an ξ c, z of 0.50 cycles/mm.MTFs were described by a linear spline of 16 points length forming a polygonal curve.Where applicable, the radial and axial MTF components were determined independently by fitting 14 equidistant points of each direction at different spatial frequencies.Keeping the value of the zeroth order of both components fixed to unity, an elementwise multiplication of the two components provided the final 3D MTF.Image power conservation dictated that a scaling factor had to be calculated from the quotient of the measured image power and the simulated object power.This factor was then applied to the MTF.The last value of the MTF spline at ξ c was fixed to zero, while values between the fitting points were linearly interpolated to describe the full MTF.
Before MTF fitting, an initial guess was provided for the first iteration in the form of a Gaussian MTF with a FWHM of 0.168 mm −1 .The FWHM of the corresponding Gaussian PSF was with 5.26 mm close to the FWHM of typical clinical PET/CT systems.Afterwards, the three steps during one MTF fitting iteration were as follows: Step ( Steps 1 to 3 are iteratively repeated until a minimal MSE is reached.
Here, ∆MTF was a 14D vector.As an additional constrain, the MTF was never allowed to fall below zero.
The MTF fitting algorithm was a combination of simulated annealing and gradient descent as described previously. 1 It tried to minimize the MSE between img calc and the actual measured image, by iteratively multiplying the MTF with an numerical representation of the spherical inserts in spatial frequency space (obj 0 , obj F , or obj Ga ).These spheres were calculated in a 256 × 256 × 256 matrix and Fourier-resampled to the 64×64×64 matrix of the subimage.This approximated, with observable Gibbs artifacts [Figs.2, 3(a), and 3(b)], ideal sampling of the simulated objects into the 64×64×64 matrix.
For the first step of the MTF fit algorithm, numerical representations of spheres were modeled accordingly as F. 2. MTF fit method.Diagram depicts the procedure for determining PSF Ga (red boxes and arrows) or alternatively PSF 0 (blue boxes and arrows) by iteratively convolving (⃝) a numerical representation of a known object (obj 0 or obj Ga ) with a MTF guess.Shown are maximum intensity projections of 3D data sets from the 8 ml 68 Ga sphere measurement at iteration i. Steps labeled 1 to 3 are explained in the main text.
without (obj 0 ) or with the respective PRF (obj F and obj Ga ).In the latter case, the simulated spheres were convolved with a PRF 3D as derived from Eqs. ( 5) and ( 6) for 18 F and 68 Ga, respectively.Thus, the PRF became part of the object [Eqs.(4a) and (4b)] used in Eq. (2a) or (2b).As such the true PSF (PSF 0 ) of a PET/CT system could be determined independently of the isotope utilized for imaging.
In the second step of every fitting cycle, MSE was calculated both in the spatial domain (MSE f ) and in the frequency domain (MSE F ); the two values were then used for simulated annealing energy states of the system.Chances for progressing into a chosen direction by ∆MTF were calculated according to the Boltzmann probability distribution.Only if the annealing temperature allowed it in both MSE domains, a successful step was counted.This approach was chosen, because in the spatial domain only one voxel contributes to the PSF tip, whereas the wings are represented by thousands of voxels.In the frequency domain, the situation is reversed.Here, high spatial frequencies (aka the tip) are represented by many voxels in frequency space (froxels), whereas low frequencies F. 3. (a) and (b) Examples of radial profiles from FBP reconstructed 8 ml spheres ( 18 F and 68 Ga), corresponding simulated spheres (obj), and its convolution img calc,F and img calc,Ga with the determined PSF after the last iteration of the MTF fitting procedure.Profiles are normalized to the expected activity in the imaged object.Profiles from calculated images img calc have both been determined by using obj F or obj Ga .(c) MTFs determined from the 8 ml sphere measurements for 18  occupy only a comparable small amount of froxels.Combining both domains renders voxel (or froxel) weighting for MSE calculations superfluous, and image noise, which tends to occupy high spatial frequencies, is not conveyed easily into the PSF or MTF.The use of both domains for MSE calculations thereby counteracts overfitting and unwanted noise amplification, making the algorithm more robust.
Annealing temperature was successively reduced by 5% after every unsuccessful step, leading to a reduction in step size |∆MTF|.The starting temperature as well as particle mass was chosen to allow at the beginning a step size |∆MTF| of 0.044 mm −1 .Fit convergence was reached when no success within the limits of numerical precision could be achieved for the last 90 iterations.
Imperfect centering of the image in its frame relative to the simulated object made fitting of a translational offset in all three spatial directions necessary.The measured image was incrementally translated by ∆s in the spatial domain using a complex valued multiplier m(ξ) in frequency space until a best fit was achieved, Six iterations of fitting ∆s in every spatial direction were performed intermittently to every 45 iterations of the main fitting procedure.Together with the 28 free MTF points, thus a total of 31 constrains were handled.Figure 2 provides a schematic overview of the MTF fitting procedure.

2.B.4. Statistical analysis
For the two different isotopes, individually fitted MTFs from all seven spheres were pooled together and averaged into two respective MTFs (MTF avg ).Point source MTFs were also averaged (MTF PS , n = 5).Averaged MTFs are shown with ± standard deviation.Comparison of point source and fit derived averaged MTFs was done by linear regression analysis providing R 2 ± standard deviation and by the rootmean-square deviation (RMSD) between the two different MTFs,

3.A. Comparison of MTF fit with point source measurements
MTFs were determined from 18 F and 68 Ga sphere measurements using the MTF fit.At first, only imaging data reconstructed with FPB were analyzed.Figure 3 illustrates radial profiles of the imaged 8 ml sphere, the corresponding numerical representations, the determined MTFs, and the resultant PSFs.
When using obj 0 for the fitting procedure, MTFs derived from all seven 68 Ga measurements (MTF Ga ) showed expectedly a worse frequency transfer than the 18 F MTFs (MTF F ).
MTFs from 68 Ga measurements thus gave correspondingly broader PSFs than their 18 F counterparts.Averaged isotope specific MTFs and PSFs compared well with the respective point source measurements with respect to FWHM, FWTM, and 50% cutoff frequency (ξ 50% ).The fit derived MTFs followed almost exactly point source measurement derived MTFs for most spatial frequencies [Figs.4(a) and 4(b), solid lines].For 18 F, the RMSD was 0.035 both radially and axially.The respective RMSD for 68 Ga was 0.023 radially and 0.039 axially.Any difference seen in the corresponding PSFs between point source measurements and spherical measurements was below the resolution capabilities of our PET/CT system, i.e., 4.5 mm (NEMA 2007) and 2 mm (Siemens True X + Time of Flight). 24TF 0 was determined from 68 Ga measurements when using obj Ga for the fitting process.This MTF 0 converged with the fit derived MTF F and with the MTF F derived from point source measurements [Figs.4(a) and 4(b), dashed lines].No meaningful differences in MTFs derived either from the use of obj 0 or obj F were observed, but the match between different MTFs degraded towards higher spatial frequencies [Figs.4(a) and 4(b), dashed lines].Due to the voxel anisotropy, the radial components of a MTF showed a lower intrinsic variance at higher frequencies compared to axial components.
There was also a trend visible in all fit determined MTFs towards overemphasizing higher spatial frequencies."Humps" seen at higher spatial frequencies were responsible for some ringing as seen in the thereof derived PSFs.This also explains the negative values seen in these PSFs [Fig.4(c)].When using obj F and obj Ga for fitting, PSFs aligned and effectively determined the imaging systems PSF 0 [Figs.3  and 4(d), dashed lines].Table I compares FWHM, FWTM, and ξ 50% for averaged fit and point source derived PSF or MTFs, respectively.From these values, it can also be seen that the 68 Ga annihilation distribution affects FWTM more and leaves FWHM largely unaffected.
In Figs.4(c) and 4(d), values of PSF voxels were summed and binned according to their radial distance from the PSF peak.Such a representation unveils the stronger contribution of radially distant PSF parts to the image signal.It also allows for a better comparison of PSFs of dissimilar shape than peak normalized PSF profiles.From Figs. 4(c) and 4(d) it can be seen that 18 F and 68 Ga PSFs have their maximally contributing voxels binned around r = 3 mm.The broader 68 Ga PSFs however have a higher contribution of more distant voxels to the image signal compared to 18 F PSFs. Use of obj Ga effectively determines PSF 0 , which exhibits a similar relative signal contribution as PSF F and as point source derived MTFs [Fig.4(d)].

3.B. Comparison of the MTF fit with the wing fit method
For comparison, PSFs were fitted with the wing fit method, using the same FBP phantom sphere measurements as above.In this wing fit, a parameter p allows modeling of a PSFs wings, independently of its width.The width is modeled with a parameter w.The 0.25 ml sphere gave PSFs with repeatedly inconsistent parameters p and w (data not shown) and was thus omitted from analysis.Resulting averaged MTFs and PSFs from the remaining spheres are shown in Fig. 5.The wing fit was able to accurately describe radial and axial components of MTF Ga across the entire spatial frequency range.In 18 F measurements, and when trying to determine PSF 0 by using obj Ga in 68 Ga measurements, the wing fit diverged from point source derived MTFs at higher spatial frequencies.Except for 68 Ga measurements, FWHM, FWTM, and ξ 50% were usually closer to the point source determined values when using the T I. Averaged PSF parameters (FWHM and FWTM) and the respective MTF parameters (ξ 50% ) in radial and axial directions derived by MTF fitting (MTF Fit) and PSF fitting (Wing fit) according to the numerical representation (Object) for their determination.Reference parameters from averaged imaged point sources (PS) are provided as well.MTF fit method (Table I) compared to the wing fit.Figure 6 depicts differences between wing fit determined MTFs and point source derived MTFs at different spatial frequencies.

Radial
Here, the RMSD for 18 F amounted to 0.063 radially and to 0.041 axially.For 68 Ga the RMSD was 0.004 in radial direction and 0.029 in axial direction.Wing fit determined MTFs exhibited its highest differences to point source MTFs in radial direction (9.4%±5.8%modulation transfer) at a spatial frequency of 0.168 cycles/mm when imaging 18 F. MTF fit determined MTFs exhibited its highest differences to point sources in radial direction (8.0% ± 5.7% modulation transfer) at a spatial frequency of 0.133 cycles/mm when imaging 68 Ga.
F. 6. MTF differences between point source derived MTFs and MTFs determined with the MTF fit and wing fit method in radial (a) and axial (b) directions from FBP phantom measurements.Here, both fit methods used Obj 0 as their numerical object representation.Regression analysis revealed that the best conformity with 18 F point source derived MTFs was found for MTFs determined with the MTF fit method when using obj 0 . 68Ga point source derived MTFs were best described with the wing fit method, also using obj 0 .When comparing MTF 0 determined with obj Ga and MTFs from 18 F point source measurements, the intrinsic MTF 0 was best approximated from 68 Ga measurements by using the new MTF fit method.Table II provides slope, R 2 , and p values from the regression analysis of fit derived MTFs with point source derived MTFs.

3.C. MTFs from 18 F TrueX reconstructions
The MTF fit method was also used to determine MTFs from 18 F images reconstructed with the Siemens TrueX algorithm for all seven 18 F filled spheres.The spatial frequency transfer of those MTFs was better than the frequency transfer of MTFs from FBP reconstructions.TrueX MTFs had a frequency transfer close to unity up to a certain cutoff frequency around ξ 50% (0.12 cycles/mm), where afterwards the transfer declined sharply to near zero.This cutoff frequency varied clearly with sphere size and was generally shifted towards higher frequencies for smaller spheres [Fig.7(a)].Comparing the fit determined MTFs with the averaged MTF derived from 18 F point source measurements revealed that TrueX transfers spatial frequencies depending on the imaged object-a result undemonstrable with point source measurements [Fig.7(b)].Similar to MTFs from FBP reconstructions, frequency transfer followed a monotonically decreasing curve when imaging point sources.In contrast thereof, the averaged MTF from TrueX reconstructed images of spheres showed a maximum of around 0.063 cycles/mm [Fig.7  method compared favorably with MTF measurements from point sources.At the same time it held advantages over the wing fit method, which by itself had been deemed accurate enough for TC purposes. 6sing existing models of PRFs and the MTF fitting method, it was possible to determine the PET/CT systems intrinsic PSF 0 , regardless of the two positron emitter used.It could be shown that within the accuracy of the imaging parameters of our PET/CT system PSF F ≈ PSF 0 , and use of obj F was not necessary for its determination due to the short positron range of 18 F.In contrast, a PET system with higher resolution capabilities, such as a small animal PET, will most likely make it necessary to use 18 F PRF models in order to determine its intrinsic PSF 0 . 26Here, a relevant difference between isotopespecific PSF and PSF 0 was only observed for the higher energy β + emitter 68 Ga.Nevertheless it was shown that by using obj Ga it became possible to determine PSF 0 .Both PSF 0 from 18 F and 68 Ga experiments were nearly identical, validating the new approach.While the wing fit method was as good in describing PSF Ga as the MTF fit method, the MTF fit was better suited for PSF 0 determination from 68 Ga measurements.This was most likely due to the imaging parameters chosen for high spatial resolution, giving rise to a cusp shaped PSF 0 .In this case, the resultant PSF 0 had less of a Gaussian shape, which was difficult to describe with the wing fit method.When imaging the high energy β + -emitter, the wing fit method provided better low-noise MTF approximations than the MTF fit method [c.f.Fig. 5(a)].Contrary to the MTF fit method, the wing fit method was not able to determine the PSF/MTF from the smallest sphere.

4.B. MTFs from TrueX reconstructions
The novel MTF fitting method was also able to determine MTFs of nonlinear reconstruction methods.Although erratic, these MTFs experienced a strong dependency on the imaged sphere's size already after just one iteration of the TrueX algorithm.As such, MTFs from TrueX reconstructions did not match MTFs from point source measurements.Even though point source derived MTFs were able to show a better frequency transfer when reconstructed with TrueX than when using FBP [Fig.7(b)], usefulness of these measurements is limited in the context of real object imaging.This would hold true for any case where spatially variant PSFs are super positioned, such as when using line sources or grid patterns with nonlinear image reconstruction methods: PSFs (or MTFs) obtained this way will experience similar limitations as point source derived MTFs. 27For TrueX reconstructed images, the MTF fit method seemed superior to point source derived MTFs, as it reflects true frequency transfer of the imaged object.Therefore, position, object, and background dependent MTFs from nonlinear iterative PET/CT image reconstructions need to be estimated from objects similar to those intended for clinical imaging.It remains to be shown how MTFs of nonlinear reconstruction methods depend on the size, heterogeneity, gradient, texture, background, and other imaging features of the imaged object.Our findings suggest that the standard method for lesion delineation-the application of a fixed relative threshold 28 -will therefore yield inconsistent and incomparable results with clinical TrueX images.

4.C. Cutoff frequency and sampling
At higher frequencies, MTFs fitted in the spatial frequency domain deviated more from the MTFs obtained with point sources.The voxel sizing and postreconstruction filtering of the actual acquisition protocol were chosen to explore these imaging properties at their resolution limits of the employed PET/CT system.The results therefore reflected likely the imaging systems inability to properly convey spatial frequencies above ∼0.25 cycles/mm [Figs.4(a) and 4(b)], rather than limitations of the fitting method itself.Corresponding to a maximal resolution of 2.0 mm in the spatial domain, this compares well with the specification stated for HD PET (True X + TOF) by the PET/CT manufacturer. 24It can be noticed that this true ξ c is self-determining and equals to that spatial frequency, at which the fit becomes erratic.Any spatial frequency above this true ξ c can thus be ignored for PSF determination.The point source derived FBP MTFs also tended towards zero around this ξ c .

4.D. Strengths of the MTF fit method
The MTF fit method shows several strengths compared to the other methods.It can fit non-Gaussian and complex shaped MTFs, as seen for the axial direction of the 18 F point source measurement [Fig.4(b)].Also, MTF/PSF and recovery curve measurements are obtainable in one measurement.The contribution of more voxels to the signal is also advantageous over point source measurements, for which several averaging acquisitions are required.Compared to point source measurements, the MTF fit method is more insensitive to differences in phantom placement and sampling errors.It also makes it easier to obtain imaging data with a well-defined activity concentration.In contrast to a simple deconvolution of the image with the numerical object representation, 29 the MTF fit will not lead to high frequency noise amplification.While this problem was somewhat mitigated by Lodge et al. 29 through the use of a large object with no background activity, the determination of spatially variant MTFs was thereby hampered.In contrast to the method employed by Lodge et al., 29 the MTF fit is able to determine spatially variant MTFs in different places across the FOV by using smaller objects.

4.E. Limitations and pitfalls of the MTF fit method
Several limitations of the presented MTF fit method should be taken into consideration: The number of constrains in the spatial frequency domain limits fit accuracy.Increasing the amount of fit points comes at the expense of increased fit duration and, more importantly, the possibility of overfitting MTFs.
Furthermore, so far only radially centered PSFs can be determined, limiting the method to the center of the FOV.Fitting of the phase would become necessary to provide for directionality of anisotropic or skewed PSFs/MTFs.Doubling the amount of constrains would most likely require unconventional fitting methods for high dimensional fit spaces.
Also, the Nyquist-Shannon sampling theorem has to be obeyed throughout the fitting process for an accurate PSF determination.Here, this made ideal sampling of the numerical object representations necessary, if sinc-function shaped artifacts are to be avoided in the final MTFs.This was seen axially, where the chosen voxel sizing of 1 mm was below the resolution capabilities of our PET/CT system.In such a case, the MTF fit method can introduce artifacts into the MTF.The MTF fit algorithm converged properly as long as the initial FWHM of the initially guessed Gaussian MTF was within ±50% of the PET/CT system's stated FWHM (data not shown).

4.F. Relevance and implications
For the first time, the new MTF fit method can accurately determine non-Gaussian PSFs.These MTFs/PSFs can then be used for the calculation of TCFs of PET/CT systems.In PET/CT system characterization for clinical studies, determination of MTF 0 or PSF 0 could here complement or even supersede the nowadays commonly used recovery curves.
Also, such PSF determination can be done in the context of routinely performed quality assurance measurements.Filling and imaging of hot spheres are usually much easily achieved than the preparation and imaging of point sources.If a long lived 68 Ge-phantom is available, it is possible to determine PSF 0 of 18 F images acquired on the same PET/CT system.Measuring always with the same long lived 68 Ge-phantom will add reproducibility and accuracy in the characterization of PET/CT systems and its different (nonlinear) reconstruction algorithms.For transconvolution of 18 F images from different PET/CT systems to a common virtual PET system, PSF 0 of these individual PET/CT systems can then be determined with the favorable 68 Ge-phantom.
From the result PSF F ≈ PSF 0 , it follows that clinical PET/CT systems, for which a measured PSF 0 is used in iterative image reconstructions algorithms, are mainly optimized for 18 F imaging.Here, high energy β + emitters will most likely yield incomparable images of similarly sized lesions.
In clinical applications, where nonlinear image reconstruction methods are preferred over FBP, estimating a lesion associated local MTF can aid in its quantification by using this MTF/PSF for PVE correction. 30

CONCLUSIONS
With the new fit method, non-Gaussian MTFs can be accurately determined from PET/CT images of known objects.By applying a literature derived PRF model, the new method is suitable for the determination of a PET/CT system's intrinsic MTF 0 or PSF 0, regardless of the β + emitter used for data acquisition.This allows the use of 68 Ge-phantoms for accurate transconvolution of 18 F images in multicenter clinical trials.
Modern iterative PET/CT reconstruction algorithms have nonlinear imaging properties, which depend on the imaged object itself.Their MTFs are not accessible by point source measurements.

F. 1 .
(a) Sheer plan of the PTFE phantom flask with water (light grey) and the hot 8 ml sphere insert (dark grey).(b) FBP PET/CT image of the phantom flask and the 8 ml sphere insert filled with 68 Ga.The filler hole with its paraffin sealing appears at the sphere opposite to the fastening rod.Medical Physics, Vol.43, No. 10, October 2016 F and 68 Ga.(d) PSFs derived from the MTFs shown.Dotted lines actually depict PSF 0 of the PET/CT system.For comparison of shapes, PSFs are normalized to their peaks.Scale bars provided for spatial measure.Medical Physics, Vol.43, No. 10, October 2016

F. 4 .
MTF fitting method.Average fitted radial (a) and axial (b) MTFs determined from FBP images.Simulated objects (spheres) were used for MTF fitting, either with or without their associated PRF.The dotted lines represent de facto the imaging system intrinsic MTF 0 or its PSF 0 respectively.The error bars on the MTFs represent ± standard deviation.(c) and (d) The relative contribution of PSFs to the signal as function of PSF radius.Three dimensional PSFs were reconstructed from the averaged MTFs above and compared with PSFs obtained from point source measurements.Then, values of all PSF voxels with same radius around the most intense central voxel were binned in 1 mm intervals and summed up to the values in the graph.(c) shows the PSFs determined with point sources and object models not using any PRFs, whereas (d) shows PSFs obtained with models incorporating 18 F and 68 Ga PRFs compared with the 18 F point source PSF.For clarity, the line for 68 Ga is not shown anymore.Legend in (a) applies to all.

F. 5 .
PSF wing fit method.Average radial (a) and axial (b) MTFs determined from FBP images.Simulated objects (spheres) were used for MTF fitting, either with or without their associated PRF.The dotted lines represent MTF 0 or its PSF 0 respectively.The error bars on the MTFs represent ±standard deviation.(c) and (d) The relative contribution of PSFs to the signal as function of PSF radius.Three dimensional PSFs were reconstructed from the averaged MTFs above and compared with PSFs obtained from point source measurements.Then, values of all PSF voxels with same radius around the most intense central voxel were binned in 1 mm intervals and summed up to the values in the graph.(c) shows the PSFs determined with point sources and object models not using any PRFs, whereas (d) shows PSFs obtained with models incorporating 18 F and 68 Ga PRFs compared with the 18 F point source PSF.The line for 68 Ga is not shown anymore.Legend applies to all.
(b)].Error bars in Fig. 7(b) hint also a higher uniformity of point source MTFs over fit derived MTFs.
. MTFs from FBP reconstructions Our results confirm the applicability of fitting MTFs in the spatial frequency domain for determining PSFs from linearly reconstructed spherical phantom measurements.The MTF fit F. 7. (a) Individual MTFs from TrueX reconstructed 18 F hot spheres and its average.The 0.25 sphere is excluded from the average trace in red.(b) Averaged MTF from TrueX point source measurements compared with the average trace from (a) and with the averaged FBP 18 F MTF (grey) taken from (a). (See color online version.)Medical Physics, Vol.43, No. 10, October 2016

Medical Physics, Vol. 43, No. 10, October 2016
1) Calculated image: The inverse Fourier transform of the MTF provides a PSF, which is convolved with obj 0 or alternatively, when determining PSF 0 , convolved with obj F or obj Ga to arrive at a calculated image img calc .
18 II.Regression analysis: Slope and R 2 are derived from regression against either the averaged radial18F point source MTF or 68 Ga point source MTF.p-values < 0.001 for all cases.