Analytical and experimental FWHM of a gamma camera: theoretical and practical issues

Introduction. It is well known that resolution on a gamma camera varies as a function of distance, scatter and the camera’s characteristics (collimator type, crystal thickness, intrinsic resolution etc). Manufacturers frequently provide only a few pre-calculated resolution values (using a line source in air, 10–15 cm from the collimator surface and without scattering). However, these are typically not obtained in situations resembling a clinical setting. From a diagnostic point of view, it is useful to know the expected resolution of a gamma camera at a given distance from the collimator surface for a particular setting in order to decide whether it is worth scanning patients with “small lesion” or not. When dealing with absolute quantification it is also mandatory to know precisely the expected resolution and its uncertainty in order to make appropriate corrections. Aim. Our aims are: to test a novel mathematical approach, the cubic spline interpolation, for the extraction of the full width at half maximum (FWHM) from the acquisition of a line source (experimental resolution) also considering measurement uncertainty; to compare it with the usually adopted methods such as the gaussian approach; to compare it with the theoretical resolution (analytical resolution) of a gamma camera at different distances; to create a web-based educational program with which to test these theories. Methods. Three mathematical methods (direct calculation, global interpolation using gaussian and local interpolation using splines) for calculating FWHM from a line source (planar scintigraphy) were tested and compared. A NEMA Triple Line Source Phantom was used to obtain static images both in air and with different scattering levels. An advanced, open-source software (MATLAB/Octave and PHP based) was created “ad hoc” to obtain and compare FWHM values and relative uncertainty. Results and Conclusion. Local interpolation using splines proved faster and more reliable than the usually-adopted Gaussian interpolation. The proposed freely available software proved effective in assessing both FWHM and its uncertainty.


INTRODUCTION
The spatial resolution of a gamma camera is a measure of its ability to resolve small objects in the field of view. Spatial resolution can also be defined as the minimum distance between two points such that they can be pictured separately. This means that objects placed at a distance smaller than the resolution limit are imaged as a single blurred one.
Minute variations (even of 0.3 mm) in the system's resolution could affect image quality (Dendy, Barber & Bayliss, 1988). Therefore, it is important to know precisely "a priori" what the gamma camera's limits are before scanning a patient. Moreover, when dealing with the absolute quantification of the tracer in SPECT, the measure of the mean radioactivity in a volume of interest (VOI) can be affected by an error proportional to the resolution (Kojima et al., 1989). Resolution is therefore a crucial parameter that measures the reliability of the gamma camera in a specific setting (Soreson & Phelps, 1987). A precise measure of resolution, together with the uncertainty of such measure, can lead to an appropriate qualitative reading of images and to a correct quantitative evaluation.
The overall spatial resolution of a gamma camera system (R s ) depends on different factors, both geometrical and physical, and it is usually expressed in terms of collimator resolution (R c ) and intrinsic resolution (R i ). The R s is typically assessed from the full width at half maximum (FWHM) of the profile of a point-like, or line-like, radiation source.
FWHM can be expressed as a function of the gamma camera's characteristics and the distance between the object and the collimator (Soreson & Phelps, 1987;Cherry, Sorenson & Phelps, 2012;Zaidi, 2006) (the so-called analytical resolution) or computed from the experimental data obtained from the image of a line source (the so-called experimental resolution).
Different methods have been proposed to calculate the experimental FWHM from a point spread function (PSF) or line spread function (LSF) (Hander et al., 1997;Hander et al., 2000;Wasserman, 1998;Metz, Atkins & Beck, 1980), but none of these methods provide the uncertainty of the measure of FWHM.
The aim of this work is to introduce a method for computing FWHM (using a matematical method known as splines) in the case of a parallel-hole collimator and the relative uncertainty from a LSF and compare it to the usually adopted methods. The most reliable one will be chosen using a cost function.
Every algorithm described in this paper was implemented and tested on the Phantom's data and is part of the freely-available package (Resolution Calculator 0.1.zip) developped by our group for educational purposes: http://www.rad.unipd.it/fwhm/.

Analytical resolution
The system resolution R s depends on the collimator resolution R c and on the intrinsic resolution R i (Soreson & Phelps, 1987;Cherry, Sorenson & Phelps, 2012;Zaidi, 2006). Using the convolution mathematical theory (Cherry, Sorenson & Phelps, 2012),we due to the fact that R s will be positive. The intrinsic resolution R i is linked to the properties of the detector and electronics. For the given energy of a photon, R i could be considered independent of the object-tocollimator distance, whereas the collimator's resolution depends largely on the geometrical layout and can be expressed as a function of a number of parameters: • x: distance between the object and the collimator's surface; • L: collimator's hole length; • D: collimator's hole size; • c: crystal's thickness, including an estimate of the gap between collimator and crystal and between crystal and image plane. An estimate of the average depth of interaction in the crystal has also been considered; • t: thickness of the septa where L, D, c, t are declared by the manufacturer as well as R i . Figure 1 schematically shows the geometrical layout of a point source.
To calculate R c , consider gamma rays coming from a point source P (as in Fig. 1) and particularly rays PA (parallel to the septa) and PB (angular limit). Now, since the radiation profile is similar in shape to an isosceles triangle ( which is equivalent to: Instead of L, L eff is usually used which is a length that is weighted to take the septal penetration into account, The constant µ is the linear attenuation coefficient of the material of the collimators (usually lead, µ = 2.49 mm −1 at 140 KeV). Thus from (2) and (3): (4) Figure 1 Geometrical layout of a point source acquisition by a gamma camera equipped with a parallel hole collimator (D is the diameter of the holes of the collimator and t is the septal thickness).

Image acquisition
A planar static image of a line source, filled with about 200 MBq of 99m Tc activity and inserted in the center of a NEMA SPECT Triple Line Source Phantom (as in Fig. 2) was acquired using a Thriple-Head Irix Marconi-Philips gamma-camera (256 × 256 matrix, 180 s) equipped with a parallel-hole, ultra-high resolution collimator. The line source was placed in air, water and a radioactive solution (about 30 kBq/ml of 99m Tc) to reproduce different background conditions. A loss of resolution (Cherry, Sorenson & Phelps, 2012) was expected as a consequence of an increasing scattering effect. Planar images were acquired with the Phantom at increasing source-to-collimator distances (134, 164, 194, 224, 254 and 284 mm) and exported in DICOM format. Figure 3 shows an image derived from the acquisition of a line source. It is an N × N data matrix with the number of radioactive counts in N points at N different heights. For each image, an N × J submatrix was visually selected (on the middle third of the line) so as to obtain near-constant data profiles For each j-th row of the submatrix, FWHM j was calculated from the data (x i ,y i ) i=1,...,N using the three methods described below. The FWHM value was assessed as the average of FWHM j (j = 1,...,J). The standard deviation σ and the variation coefficient Cv were calculated to estimate the absolute and relative uncertainties respectively where

From data to experimental resolution
Another way to quantify the uncertainty of FWHM is the use of a quadratic cost defined case by case. The maximum error in FWHM is expected to be proportional to the square root of such cost (Walter & Pronzato, 1997).

Method 1: direct calculation
This first method was intentionally a basic one to prove that a naive approach will lead to an unreliable FWHM value. The maximum pixel value h = max(y i ) = y K for a proper index K and its argument x = x K were found. Two points z 1 <x and z 2 >x, which are the closest to h 2 , were used and their distance FWHM j = |z 1 − z 2 | was ascertained. For this case the following cost was defined:

Method 2: Gaussian-global interpolation
Data (x i ,y i ) were modelled as a deterministic function with a small level of noise. The process called least-squares is reliable for choosing a function close to the data. Mathematically speaking, the least-squares approximation of a given data set looks for the best fit minimizing a suitable functional. Usually the functional is the sum of the squares of all deviations of a function chosen from the data. The linear least-squares approximation consists in finding a function fā(x) =  n i=1 a i φ i (x) depending on some parametersā = (a 1 ,...,a n ), with φ i ,i = 1,...,n a set of known (basis) functions. Nonlinear least-squares approximation can also be constructed (see below), providing a function fā that is a nonlinear combination of some known functions.
The algorithm looks for: a * = arg min a∈R n J(ā) = arg min a∈R n ∥y i − fā(x i )∥ 2 .
Since fā is not always linearly dependent on the parametersā, an iterative optimization algorithm was used to estimateā * . If the Optimization Toolbox/optim package is installed in MATLAB/Octave, the proposed software will use the well-known Levenberg-Marquardt algorithm; if not, it will use the Gauss-Newton algorithm (Nocedal & Wright, 2006). The cost used in this method, which is a mean square deviation (as the cost used in direct calculation) was The Gaussian function was used (Zaidi, 2006): which has a resolution FWHM j = 2 √ log(2) |a 2 | .

Method 3: splines-local interpolation
The third method proposes the use of splines of degree 1 (linear), 2 (quadratic) or 3 (cubic) calculated on a huge number of points (we used 10 4 interpolation points in the experiments).
A spline of order m is a function s(x) defined by the following (de Boor, 2001;Lancaster &Šalkauskas, 1986): This means that the polynomial pieces are continuous up to order m − 1 in each inner point.
This approach is called local interpolation. Cubic splines were chosen for their well-known approximation properties (de Boor, 2001), and for the ability to provide a modelindependent interpolation. In other words, the cubic spline approach is able to accurately follow the shape suggested by discrete data on a pointwise basis instead of searching a global fitting function.
As in Method 1 the algorithm searches for two points z 1 and z 2 a minimal distance away from the half of the maximum: where J 1 and J 2 are sets of 10 4 equidistant points of the intervals (x 1 ,x) and (x,x N ) respectively.
The distance between these two points gives a good estimate of the FWHM j The cost is defined as in Method 1:

Computation of analytical and experimental resolution curves
On the basis of formulas (1) and (4), R s can be expressed as follows: Using parameters (L,D,c,t,R i ) declared by the manufacturer (Table 1) the following values were calculated for the analytical resolution: The formula (8) was also used to fit the experimental FWHM data obtaining an experimental fitting curve.
The parameters that describe the experimental fitting curve were computed using the weighted least-squares method, as defined in Lancaster &Šalkauskas (1986).
To caclulate parameters an expression equivalent to (8) was used: where y is the vector of FWHM s and y 2 is the pointwise square of y. The vector p = (p 1 ,p 2 ,p 3 ) can be obtained by solving the normal equations of the weighted least squares method where V is the Vandermonde matrix, f = y 2 and Σ = diag(1/σ 2 1 ,...,1/σ 2 n ) the weight matrix, and σ i is the standard deviation of the i th FWHM value.
Usually the manufacturer provides the FWHM values at 0 and 100 mm. These two values (not measurable in our experimental setting), which weigh 90% less than the experimental FWHM ones, were used to regularize the fitting curve only in a scatterless condition (in air).

RESULTS
The experimental results, obtained by scanning the phantom in air, are given in Table 2 where the absolute (σ ) and relative (Cv) deviations between the experimental and analytical FWHM are also shown.    It should be noted that: -the direct method, as expected, demonstrated higher σ , Cv and mean cost values than the other two methods.
-the FWHM values obtained using global and local interpolation methods were nearly identical, and σ and Cv were also very similar for local and global interpolation methods.
-the mean cost calculated using splines is significantly lower than when the other methods are used (at least 50% less than the cost of global interpolation).
-all previous observations were also valid in the case of different scattering conditions.
-as expected, the FWHM increases with source-to-collimator distance (radius).
-all the methods used revealed a cost that decreased in proportion to source-tocollimator distance.

DISCUSSION
The direct method is very quick and easy but it proves more "costly" than the other methods. The resulting FWHM does not differ significantly from the others because it is obtained as the mean of a number of profiles. It is important to use the means of multiple profiles if the direct method is the only one available, but the high Cv should discourage the use of this method. The results obtained using local and global interpolation were nearly identical in terms of FWHM, σ and Cv but different in term of mean cost. Local interpolation is therefore more reliable than global interpolation when a small sample of profiles is chosen.
It should be noted that the widely used gaussian interpolation is much more time-consuming (about 25 times) than spline calculation. Given its lower cost and higher computational speed, the splines method is a very good choice for calculating FWHM from a static image. Where images are used qualitatively, reported differences in FWHM (between local and global approaches) are irrelevant. If images are also used for a quantitative approach, it is mandatory to have a FWHM value as reliable as possible. Therefore, in this latter situation, the spline method seems to be a better choice.
When the FWHM obtained using analytical calculations is compared with the results of the splines method, the difference, in air, range from about 3% to 5% (up to 0.5 mm for the largest radius); the greater the difference, the larger the radius. As expected, this trend becomes worse with a greater degree of scattering (even considering, as in our data, a low background activity concentration). For example the difference between the experimental FWHM with background and analytical resolution ranges from 13% to 15% (up to 1.7 mm for the largest radius). The software presented in this work is able to quantify this uncertainty effectively in terms of σ , Cv and cost.

CONCLUSIONS
Three mathematical methods for assessing the experimental resolution obtained from static data were inputed and tested on Phantom-derived data. Local interpolation using splines proved more reliable and faster than the usually adopted gaussian interpolation (the global method).
An open source package for calculating analytical and experimental FWHM was developed in MATLAB/Octave and proved effective in assessing both FWHM and its uncertainty. A similar PHP web-based application was also developed for open access. Both tools enable a graphical and numerical comparison of experimental and analytical FWHM. These programs are freely available at: http://www.rad.unipd.it/fwhm/.