Analytical and experimental FWHM of a gamma camera: theoretical and practical issues
 Published
 Accepted
 Received
 Academic Editor
 Noriko Hiroi
 Subject Areas
 Bioengineering, Biophysics, Radiology and Medical Imaging
 Keywords
 Gaussian, Splines, Gamma camera, FWHM, Resolution, Local interpolation, Global interpolation
 Copyright
 © 2015 Cecchin et al.
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
 Cite this article
 2015) Analytical and experimental FWHM of a gamma camera: theoretical and practical issues. PeerJ 3:e722 https://doi.org/10.7717/peerj.722 (
Abstract
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 precalculated 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 webbased 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, opensource 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 usuallyadopted 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 pointlike, or linelike, 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 socalled analytical resolution) or computed from the experimental data obtained from the image of a line source (the socalled 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 parallelhole 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 freelyavailable 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 obtain ${R}_{s}^{2}={R}_{c}^{2}+{R}_{i}^{2}$, which gives us (1)$R}_{s}=\sqrt{{R}_{c}^{2}+{R}_{i}^{2}}\hspace{0.17em$ 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 objecttocollimator 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 $\overline{PA}$ (parallel to the septa) and $\overline{PB}$ (angular limit). Now, since the radiation profile is similar in shape to an isosceles triangle ($\stackrel{\u25b3}{HIB}$) and triangles $\stackrel{\u25b3}{HIB}$ and $\stackrel{\u25b3}{HFG}$ are similar, then the FWHM is about half of the base $\left(\right.{R}_{c}\simeq \overline{AB}$).
Because of the similitude of the triangles $\stackrel{\u25b3}{PAB}$ and $\stackrel{\u25b3}{{P}^{\prime}{A}^{\prime}B}$ and the fact that $t\ll D\left(=\overline{{A}^{\prime}{B}^{\prime}}\right)$ we infer that $\overline{{P}^{\prime}{A}^{\prime}}:\overline{{A}^{\prime}{B}^{\prime}}=\overline{PA}:\overline{AB}$ which is equivalent to: (2)$L:D=\left(c+L+x\right):{R}_{c}.$ Instead of L, L_{eff} is usually used which is a length that is weighted to take the septal penetration into account, (3)${L}_{eff}=L\frac{2}{\mu}.$ 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)${R}_{c}=D\left(1+\frac{x+c}{{L}_{eff}}\right).$
Materials and Methods
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 ThripleHead Irix MarconiPhilips gammacamera (256 × 256 matrix, 180 s) equipped with a parallelhole, ultrahigh 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 sourcetocollimator distances (134, 164, 194, 224, 254 and 284 mm) and exported in DICOM format.
From data to experimental resolution
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 nearconstant data profiles For each jth 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 $Cv=\frac{100\cdot \sigma}{FWHM}.$ 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 $\stackrel{\u0303}{x}={x}_{K}$ were found. Two points ${z}_{1}<\stackrel{\u0303}{x}$ and ${z}_{2}>\stackrel{\u0303}{x}$, which are the closest to $\frac{h}{2}$, were used and their distance $FWH{M}_{j}={z}_{1}{z}_{2}$ was ascertained.
For this case the following cost was defined: (5)${C}_{1}\left({z}_{1},{z}_{2}\right)=\frac{{\left(y\left({z}_{1}\right)h/2\right)}^{2}+{\left(y\left({z}_{2}\right)h/2\right)}^{2}}{2}.$
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 leastsquares is reliable for choosing a function close to the data. Mathematically speaking, the leastsquares 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 leastsquares approximation consists in finding a function ${f}_{\stackrel{\u0304}{a}}\left(x\right)=\sum _{i=1}^{n}{a}_{i}{\varphi}_{i}\left(x\right)$ depending on some parameters $\stackrel{\u0304}{a}=\left({a}_{1},\dots ,{a}_{n}\right)$, with ϕ_{i}, i = 1, …, n a set of known (basis) functions. Nonlinear leastsquares approximation can also be constructed (see below), providing a function ${f}_{\stackrel{\u0304}{a}}$ that is a nonlinear combination of some known functions.
The algorithm looks for: ${\stackrel{\u0304}{a}}^{\ast}=arg\underset{\stackrel{\u0304}{a}\in \stackrel{n}{\mathbb{R}}}{min}J\left(\stackrel{\u0304}{a}\right)=arg\underset{\stackrel{\u0304}{a}\in \stackrel{n}{\mathbb{R}}}{min}\Vert {y}_{i}{f}_{\stackrel{\u0304}{a}}\left({x}_{i}\right){\Vert}^{2}.$ Since ${f}_{\stackrel{\u0304}{a}}$ is not always linearly dependent on the parameters $\stackrel{\u0304}{a}$, an iterative optimization algorithm was used to estimate ${\stackrel{\u0304}{a}}^{\ast}$. If the Optimization Toolbox/optim package is installed in MATLAB/Octave, the proposed software will use the wellknown Levenberg–Marquardt algorithm; if not, it will use the GaussNewton algorithm (Nocedal & Wright, 2006).
The cost used in this method, which is a mean square deviation (as the cost used in direct calculation) was (6)${C}_{2}\left(\stackrel{\u0304}{a}\right)=\frac{J\left({\stackrel{\u0304}{a}}^{\ast}\right)}{N}=\frac{\Vert y{f}_{{\stackrel{\u0304}{a}}^{\ast}}\left(x\right){\Vert}^{2}}{N}.$ The Gaussian function was used (Zaidi, 2006): $f}_{\stackrel{\u0304}{a}}\left(x\right)={a}_{1}{e}^{{a}_{2}^{2}{x}^{2}$ which has a resolution $FWH{M}_{j}=2\frac{\sqrt{log\left(2\right)}}{\left{a}_{2}\right}$.
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):

on each subinterval I_{i} = [x_{i}, x_{i+1}], i = 1, …, N−1, s_{Ii} = s_{i}∈𝒫^{m}(ℝ) where 𝒫^{m}(ℝ) is the space of real polynomials of degree ≤ m.

s(x_{i}) = y_{i}, i = 1, …, N, i.e., s interpolates the data.

${s}_{i}^{\left(k\right)}\left({x}_{i}\right)={s}_{i+1}^{\left(k\right)}\left({x}_{i}\right),i=2,\dots ,N1,k\le m1,m\le N$. 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 wellknown 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: ${z}_{i}=arg\underset{x\in \underset{i}{J}}{min}s\left(x\right)h/2,\phantom{\rule{1em}{0ex}}i=1,2$ where J_{1} and J_{2} are sets of 10^{4} equidistant points of the intervals $\left({x}_{1},\stackrel{\u0303}{x}\right)$ and $\left(\stackrel{\u0303}{x},{x}_{N}\right)$ respectively.
The distance between these two points gives a good estimate of the FWHM_{j} $FWH{M}_{j}={z}_{1}{z}_{2}.$ The cost is defined as in Method 1: (7)${C}_{3}\left({z}_{1},{z}_{2}\right)=\frac{{\left(s\left({z}_{1}\right)h/2\right)}^{2}+{\left(s\left({z}_{2}\right)h/2\right)}^{2}}{2}.$
Computation of analytical and experimental resolution curves
On the basis of formulas (1) and (4), R_{s} can be expressed as follows: (8)${R}_{s}\left(x\right)=\sqrt{{p}_{1}{x}^{2}+{p}_{2}x+{p}_{3}}.$ Using parameters (L, D, c, t, R_{i}) declared by the manufacturer (Table 1) the following values were calculated for the analytical resolution: (9)${p}_{1}\cong 0.0010,\phantom{\rule{2em}{0ex}}{p}_{2}\cong 0.1468,\phantom{\rule{2em}{0ex}}{p}_{3}\cong 2.4267$
L  D  c  t  FWHM  FWHM  R_{i} 

(mm)  (mm)  (mm)  (mm)  @ 0 mm  @ 100 mm  (mm) 
58.4  1.78  19  0.152  4.8 mm  6.7 mm  4.1 
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 leastsquares method, as defined in Lancaster & Šalkauskas (1986).
To caclulate parameters an expression equivalent to (8) was used: $y}^{2}={p}_{1}{x}^{2}+{p}_{2}x+{p}_{3$ 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 $\left({V}^{T}\Sigma V\right)p={V}^{T}\Sigma f$ where V is the Vandermonde matrix, f = y^{2} and $\Sigma =diag\left(1/{\sigma}_{1}^{2},\dots ,1/{\sigma}_{n}^{2}\right)$ 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.^{1}
Experimental FWHM  
radius (mm)  134  164  194  224  254  284  
Direct Met.  FWHM (mm)  7.31  8.13  8.96  9.71  10.62  11.37 
σ (mm)  0.80  1.19  0.87  0.87  1.18  0.80  
Cv (%)  10.9  14.6  9.7  9.0  11.1  7.0  
mean cost  147543  137062  88423  66376  44393  34223  
Global Int.  FWHM (mm)  7.53  8.20  8.88  9.65  10.54  11.35 
σ (mm)  0.24  0.22  0.22  0.21  0.33  0.30  
Cv (%)  3.2  2.7  2.5  2.2  3.1  2.7  
mean cost  3.75  3.11  2.33  2.06  2.15  1.99  
Local Int.  FWHM (mm)  7.50  8.23  8.92  9.70  10.54  11.36 
σ (mm)  0.26  0.27  0.26  0.26  0.37  0.31  
Cv (%)  3.4  3.3  2.9  2.6  3.5  2.7  
mean cost  1.55  0.95  0.79  0.49  0.35  0.25  
Analytical FWHM  
FWHM (mm)  7.70  8.50  9.33  10.17  11.03  11.89 
The analytical and experimental (both in air and with different level of background) FWHM curves are shown in Fig. 4.
Figure 5 shows an overall comparison between zoomedin details of line sources acquired in all conditions.
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 sourcetocollimator distance (radius).

all the methods used revealed a cost that decreased in proportion to sourcetocollimator 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 timeconsuming (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 Phantomderived 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 webbased 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/.