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

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 obtain $R_s^2=R_c^2+R_i^2$, which gives us (1)${\displaystyle 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 object-to-collimator 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{△}{HIB}$) and triangles $\stackrel{△}{HIB}$ and $\stackrel{△}{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{△}{PAB}$ and $\stackrel{△}{P^{\prime }{A}^{\prime }B}$ and the fact that $t\ll D\left(=\overline{A^{\prime }{B}^{\prime }}\right)$ we infer that ${\displaystyle \overline{P^{\prime }{A}^{\prime }}:\overline{A^{\prime }{B}^{\prime }}=\overline{PA}:\overline{AB}}$ which is equivalent to: (2)${\displaystyle 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)${\displaystyle L_{eff}=L-\frac{2}{\mu }.}$ The constant μ is the linear attenuation coefficient of the material of the collimators (usually lead, μ = 2.49 mm⁻¹ at 140 KeV). Thus from (2) and (3): (4)${\displaystyle 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 ^99mTc 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 ^99mTc) 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.

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 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 ${\displaystyle 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 $\tilde{x}=x_K$ were found. Two points $z_1<\tilde{x}$ and $z_2>\tilde{x}$, which are the closest to $\frac{h}{2}$, were used and their distance ${\displaystyle FWH{M}_j=|z_1-z_2|}$ was ascertained.

For this case the following cost was defined: (5)${\displaystyle 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 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_{\bar{a}}\left(x\right)=\sum _{i=1}^n{a}_i{\varphi }_i\left(x\right)$ depending on some parameters $\bar{a}=\left(a_1,\dots ,a_n\right)$, 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_{\bar{a}}$ that is a nonlinear combination of some known functions.

The algorithm looks for: ${\displaystyle {\bar{a}}^{\ast }=arg\underset{\bar{a}\in \stackrel{n}{\mathbb{R}}}{min}J\left(\bar{a}\right)=arg\underset{\bar{a}\in \stackrel{n}{\mathbb{R}}}{min}\Vert y_i-f_{\bar{a}}\left(x_i\right){\Vert }^2.}$ Since $f_{\bar{a}}$ is not always linearly dependent on the parameters $\bar{a}$, an iterative optimization algorithm was used to estimate ${\bar{a}}^{\ast }$. 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 (6)${\displaystyle C_2\left(\bar{a}\right)=\frac{J\left({\bar{a}}^{\ast }\right)}{N}=\frac{\Vert y-f_{{\bar{a}}^{\ast }}\left(x\right){\Vert }^2}{N}.}$ The Gaussian function was used (Zaidi, 2006): ${\displaystyle f_{\bar{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⁴ 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):

1. 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.

2. s(x_i) = y_i, i = 1, …, N, i.e., s interpolates the data.

3. $s_i^{\left(k\right)}\left(x_i\right)=s_{i+1}^{\left(k\right)}\left(x_i\right),i=2,\dots ,N-1,k\le m-1,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 well-known approximation properties (de Boor, 2001), and for the ability to provide a model-independent 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₁ and z₂ a minimal distance away from the half of the maximum: ${\displaystyle z_i=arg\underset{x\in \underset{i}{J}}{min}|s\left(x\right)-h/2|,i=1,2}$ where J₁ and J₂ are sets of 10⁴ equidistant points of the intervals $\left(x_1,\tilde{x}\right)$ and $\left(\tilde{x},x_N\right)$ respectively.

The distance between these two points gives a good estimate of the FWHM_j ${\displaystyle FWH{M}_j=|z_1-z_2|.}$ The cost is defined as in Method 1: (7)${\displaystyle 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}.}$

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.¹

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 zoomed-in 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 source-to-collimator distance (radius).

• all the methods used revealed a cost that decreased in proportion to source-to-collimator 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/.