ABSTRACT

The Radon transform fits badly Single Photon Emission Tomography (SPECT). However, Thin Holes Collimator (THC) and Radon model are widely used. The CACAO project has been proposed to enhance the quality of SPECT images. CACAO is a short hand notation for computer aided collimation tomography. The main idea of this project is to use collimators with much larger holes to increase the sensitivity, and slightly longer holes to increase the spatial resolution. The acquisition sequence includes a translation. The Radon projection is replaced by a 2D sum. A dedicated reconstruction algorithm has been developed. If the physical advantage of the project in terms of sensitivity and spatial resolution is generally admitted, a question remains unanswered: Would the ill-posedness of the inverse problem ruin the quality of the reconstructed images? In this article, a representation of the 2D direct problem matrix is derived. This allows us to compare the two inverse problems (CACAO versus THC). The condition number was used for this comparison. We studied the variation of these condition numbers with several parameters. For a proper set of parameters, the CACAO inverse problem may appear easier to solve and more accurately than the THC one.

Keywords:

Radon Transform, Tomography, SPECT, Large and Long Hole Collimators, Inverse Problems

1. Introduction

Since its invention by H. A. Anger [1] , the gamma camera has been used extensively in medical practice, mostly with a thin-parallel-hole-collimator (THC). It should be noted that, in gamma ray imaging, tomography has taken some time to be widely used in clinical practice. In contrast in X ray imaging, computer tomography (CT) was a breakthrough, due to the sudden and huge improvement in image quality that it brought. The mathematical model currently used in THC-SPECT is based on the Radon transform [2] . However, if this model is well suited to Xray imaging (CT), we think, it is not well adapted for THC-SPECT. In effect, to perfectly fit the Radon transform the collimator would need to have infinitely long and infinitely narrow holes. Such an “ideal” collimator would not allow any photons to reach the detector, and therefore no image could be reconstructed. Needless to say these are the photons which carry the information from the patient to the detector. In practice engineers have made a trade off. They have tried to deviate slightly from this “ideal” collimator but not too excessively. The results is a poor sensitivity: 10,000 photons lost for each 1 which reaches the detector and a poor spatial resolution of 8 to 10 mm (compared to the detector intrinsic resolution ≈ 3 mm). Furthermore, this poor spatial resolution worsens with the source to collimator distance (15 mm in the center of a patient). This compromise is the bottleneck of the conventional THC-SPECT. To visualize the problem Figure 1 shows the motion followed by a gamma camera in a standard THC-SPECT acquisition. Figure 2 shows an example of an acquisition made with this system for an object composed of 2 point sources.

Several authors have tried to increase the sensitivity of the gamma camera. Two-headed cameras or multi-pinholes [3] can increase this sensitivity with the number of heads or pinholes. Zhang et al. [4] , in their interesting study posited the hypothesis that high resolution (very-thin-hole) collimators may not be optimal for SPECT. Furthermore, they demonstrated that the tomographic

reconstruction could help in restoring part of the “lost” resolution due to the use of slightly larger holed collimators. Lodge [5] proposed a further increase in sensitivity through the use of very large holes in one dimension but still thin in the other. To acquire sufficient information, a double rotation motion is applied to the collimator-detector head during the acquisition. The reconstruction program needs to use 2 inversions of the radon transform.

We stopped the introduction here temporarily to depict the plan of the article. A complete understanding of the problem cannot be attained without some knowledge of the alternative model. This is described in Chapter 2. This section presents the definition of the project (2.1) followed by the mathematical analytical model (2.2) and then the physical advantages of the project are briefly described (2.3). After this description, the problem discussed in this article can be described in Chapter 3. This chapter completes this long but necessary intro- duction. The method used for the comparison of the classical THC-Radon inverse problem versus the alternative CACAO project is detailed in Chapter 4. The results are described in Chapter 5. The results are followed by a discussion 6 and a conclusion 7. The methods, Chapter 4 are divided as follow: Subchapter (4.1) gives a reminder of the condition number definition. Then we derive the impulse response function of the CACAO project (4.2) in a total and a partially attenuating collimator model (4.3). Section (4.4) explains the calculation of the transfert matrix. The THC model used is described in Section (4.5). The results begin with a study of the condition number in relation to the geometry of the collimator (5.1). This study is first presented with very small matrix dimensions in order to get a better grasp of the intervening parameters (5.2, 5.3). The following 2 chapters of the results (5.4 to 5.6) show a comparison of the condition number for the two problems. This comparison is studied with more realistic matrix sizes (up to 64 × 64). Finally noiseless and noisy image simulations are depicted (5.7, 5.8).

2. Purpose, the CACAO Project

To improve the quality of SPECT images, and avoid the drawbacks of the Radon modelling, we proposed the CACAO project. Subchapters 2.1 to 2.3 give the necessary background to understand the problem described in Subchapters 2.1.

2.1. Definition

The CACAO project [6] was proposed to avoid the bottleneck caused by conventional thin-hole-collimators. For the sake of simplicity, the CACAO project is presented here in a 2D reduction. This project shares with THC-SPECT a gamma ray detector position and energy sensitive. It differs by having a collimator with larger and longer holes, an added linear motion in the acquisition sequence, and a mathematical model which takes full account of the geometry of the collimator. It must be understood that larger holes means dramatically larger holes. Figure 3 may help to provide an understanding of the CACAO project. The hole has been chosen as large as 2/3rds the length of the object.

The exact response of the system will be calculated in the following section with the help of Figure 5. An initial examination of this Figure, may allow a better understanding of the subject.

During the linear scan the set collimator + detector travels from the left to the right of Figure 5. This Figure is divided vertically into 3 rows and horizontally into 5 columns. The 5 sectors or columns (from 1 to 5) are limited by 4 special positions of the set that have been labeled with the letters a, b, c, d. The bottom part of the figure shows a point source (S) and the set collimator + detector in the 4 special position limits (detector is shaded blue, collimator is pinstriped). The upper part of the figure is a drawing representing the acquired data after the linear scan. And the text between the upper and the lower parts represents the algebra which will be explained later. Let us start with the left part of the figure,

when the detector is further left than position a, no photons (emitted by S) can reach the detector. At the right of position a and before position b the left part of the detector is illuminated. It is important to remember that the hole is a lot larger than the detector’s spatial resolution. In this example we have chosen a collimator hole with a width of 21 pixels. So if we represent the position of the set detector + collimator in the horizontal axis and the signal recorded by the 21 pixels in the vertical axis, (left part up, right part down) the progressive illumination of the detector during the scan can be represented by a triangle (zone 2). Between positions b and c the detector is completely illuminated. The drawing of the set detector + collimator at position c has been represented behind position b and in a lighter shade of gray. In this area 3 the illumination can be represented by a rectangle. Using symetry we have a triangle in area 4. And nothing in area 5. Finally the data resembles a parallelogram. Due to geometrical considerations the greater the distance between the source and the collimator, the larger and the more slanted the parallelogram. Figure 4 shows a 2D example of the 2 point sources CACAO acquisition from Figure 3. The response of each point source has been colored accordingly in blue and red. One can observe that when the source is near the entrance to the collimator the parallelogram is in an upright position. The center of the parallelogram evolves along sinusoids similar to Figure 2 but limited to only 4 different angles.

2.2. Analytical Description

The analytical description of the CACAO transform is described here in a 2D reduction. With the Radon transform the integration is made on a line (the projection line). With CACAO the integration is made over a trapezoidal surface.

Figure 6 depicts the variables used.

A transverse section of the patient is represented by the density of radio-active sources. The collimator has been simplified to a unique hole (width D, depth P). This hole can rotate around the patient and can also perform a linear motion along the u axis. The coordinate system x, y is considered at rest while the rotating coordinates u, w follow the rotation of the collimator-detector, measured by.

Photons entering the detector are localized by their position in the hole, the position of the hole in the u coordinate measured by, and the angle of rotation. Therefore, the description of the direct problem aims to define the acquisition function, knowing the distribution of the radioactive sources. Elementary geometry considering only photons traveling in the air (no attenuation in the patient), leads to an integral equation for the direct problem. This equation gives in terms of or equivalently. It is obtained by considering the surface integration in a polygon bounded by the 2 extreme rays touching the 2 corners of the collimator entrance hole. This polygon is also limited by the largest camera orbit () and the entrance face of the collimator (nearest distance of observation) [6] .

(1)

where

- represents the intensity of the detected photons with measu- ring the position of the camera head, the coordinate on the detector surface and measures the inclination of the camera head in the coordinates.

- is the distribution of the radioactive source in the rotating coordinates (relative to the detector, u is parallel to the detector, w is perpen- dicular).

- is the angle between the gamma ray direction and the perpendicular line to the surface of the detector (direction w).

- is the source to detector distance.

- represents the law of illumination.

- are the limits of integration in the direction w, P is the depth of the collimator, and represents the farthest limit of integration. (limited by the camera orbit radius).

- is the width of the collimator hole.

2.3. A Brief Reminder of the Physical Advantages of the CACAO Project

These potential advantages have already been published [6] [7] , and are com- mented here for the sake of clarity.

2.3.1. Sensitivity

The sensitivity of a thin parallel hole collimator can be estimated by the following formula (valid in 3D). The sensitivy is the percentage of the number of photons emitted by the patient to the number of photons passing through the collimator.

(2)

where D is the diameter of the collimator hole, P is the length of the collimator hole, K is a coefficient depending on the shape of the holes (rectangular, hexagonal…) and t is the thickness of the septa. This formula is valid for conventional thin-hole-collimators (THC), and it would stand for the CACAO project for a multiholes collimator. Compared to THC, the CACAO sensitivity

could be superior, by a factor of 100 or more by increasing in. To give

some order of magnitude, conventional THCs work with geometrical sensitivities

of for a. In this article, for the CACAO system we have chosen. Such a value allows us to hope for an increase of the sensitivity of the

system by 625 times in 3D. For the sake of simplicity, the description in this paper involves only one hole. But a real implementation of the system will have many holes arranged in a compact pattern. This multi-hole realization will justify the formula (2), and will reduce the extent of the linear scanning motion required during the acquisition, at the price of a larger detector. It is to be noted that the third term of the formula (2) represents the ratio of free surface of the detector to the hidden surface (commonly covered by lead). It will be discussed later.

2.3.2. Spatial Resolution

In the conventional THC case, increasing the hole diameter deteriorates the spatial resolution. In contrast, in the CACAO project this fact is no longer true. This point may be difficult to understand due to the commonly held belief that:

To better understand why this fact is only true for THC, let us draw a parallel with Computerized Tomography. In this system, the signal emanating from a line of voxels is recorded in an acquisition element. The volume corresponding to this line of voxels is thin but quite long (generally 20 cm at the level of the patient’s head and 70 cm at the belly). However a huge number of acquisitions combined with the reconstruction process allow to achieve a milimetre resolution in this long volume.

The following figure examines the physical and geometrical factor governing the spatial resolution of both systems.

The left part of Figure 7 represents a THC collimator with a detector of intrinsic spatial resolution. The spatial resolution deteriorates with

If is the diameter of the thin holes, is the length of these thin holes L the distance source-entrance face of the collimator, the following formula is largely accepted:

(3)

This formula is based on a gaussian model, where the FWHM: Full width at Half Maximum of the convolution of 2 gaussians, is given by a mean square of the corresponding FWHM’s.

The right part of Figure 7 represents a CACAO collimator with the same detector intrinsic spatial resolution and the same source collimator distance. By using a similar integration as classical tomography the dimensions of the volume studied by an element of the detector are not important here. The spatial resolution will be governed by the steepness of the slope deliminating this volume. Mathematically a derivative of this slope, will give a gaussian, from which a FWHM can be derived. Clearly this steepness will be in the order of:

(4)

where P is the depth of the collimator hole (CACAO) and and L are the same intrinsic spatial resolution and object to collimator distance as previously described. Two points need to be emphasized here: First D which greatly influences the sensivity in the THC case is no longer present. Hence the sensitivity-resolution trade-off is different to the THC case. Second by choosing a longer P, longer than L we can obtain a better spatial resolution than the intrinsic one. To make the things clearer let us fix some dimensions: In the comparison chapter, we will use pixels. With a 3 mm pixel width (comparable to the resolution of the actual gamma camera) we will have

mm. This value is slightly greater than conventional system currently in use (order of magnitude: 35 mm or less). This figure gives potentially a spatial resolution of the acquisitions of 4.5 mm at 100 mm from the entrance of the collimator, if not degraded by other factors. It must be emphasized that P also governs the speed with which the spatial resolution worsens with the source to collimator distance. With the parameters given, the resolution of the THC case will degrade roughly 2 time as fast as the CACAO case. Also this effect will not be clearly apparent in this article (due to the very small size of the images presented). Finally the step of the linear scanning motion could be another way to improve the spatial resolution in the CACAO case. According to Shannon criteria a step of half a pixel size is optimum. Our calculations and simulations use a step equal to.

2.3.3. Quality of Collimation

Another advantage of the CACAO project is the possibility to use larger septa (the walls between the holes of the collimator). The width of the septa has already appeared in the formula (2). Currently THCs use very thin septa (t = 0.015 mm for example) which let 5% of the photons to go through the septa and be detected. Due to the huge number of holes (40,000), a tiny increase of the thickness of the septa will dramatically reduce the sensitivity of conventional systems. On the contrary, for CACAO with only 10 holes or so, in the 600 mm field of view of the camera, a thickness of 2 mm or more will have a very limited effect in the sensitity. By the same token, it can greatly reduces the septal penetration.

2.3.4. Pixelated Detectors

Finally, it has been published that a conjoint use of the CACAO project with pixelated detectors presents a synergistic improvement [8] . These semiconductor detectors can greatly improve the spatial resolution of the currently used detector. For a NaI detector the intrinsic resolution () is of the order of 3 mm. With a pixelated detector of 1mm the intrinsic resolution is divided by 3. The diameter of the hole of a thin-hole-collimator would have to be divided by at least the same value to benefit from this improvement. The sensitivity would be reduced by due to the same token. With CACAO there is no need to reduce the surface of the hole. The improvement of the intrinsic spatial resolution will automatically improve the sharpness of the edge of the parallelogram, without impairing the sensitivity.

3. The Problem

The physical advantages of the CACAO project in term of increased sensitivity are clear and have possibly been accepted by many readers. With larger holes, one can collect more photons. It is also clear that the information is carried by these photons. However two arguments can be raised against this project.

The first is the following: Do the photons collected with CACAO carry the same information quality as the photons collected with THC SPECT in terms of location?

The second question is: Would the gain obtained by the increased number of collected photons be hampered by the difficulty of solving the inverse problem? In other words: which factor will be the greatest: 1) the noise amplification inherent in the mathematical calculation of the reconstructed image or 2) the improvement of the signal to noise ratio of the acquisition brought about by the increased number of photons?

Both these questions are addressed in this article. In fact a positive answer to these questions based on information theory argument has already been published [9] . In this article we studied the multiplex volume (the volume of the object sending photons to a single detector pixel), and we showed that: one can obtain better reconstructed image with a large multiplex volume with precise boundaries (CACAO), than with a small multiplex volume with blurred boundaries (THC).

But we offer here another mathematical argument, based on linear algebra theory. This article studies the difficulty in solving the CACAO problem and compares it to the conventional THC-SPECT problem. To derive this comparison we first calculated the direct matrix of both the CACAO project and the THC-SPECT. The physical considerations for both systems were chosen so as to be equivalent. Then we calculated the condition number of these matrices. The condition number of a matrix measures the accuracy of the solution of the linear equation system and its noise sensitivity. Therefore, with equivalent noise in the data (second member of the equation), a better (lower) condition number of the transfer matrix will give a better solution, e.g. better images.

4. Methods

4.1. Condition Number, Inverse Problem

The condition number [10] , here noted has the following interest:

Considering the linear system, with noisy:

If then:

(5)

For practical reasons the tomographic reconstruction is nowadays calculated by computer and the reconstructed image has a finite number of pixels:

, which represents the acquisition is also discrete. For a rectangular matrix m by n (m rows and n columns), the estimation is done by a singular value decomposition that is: Where and are orthogonal matrices (respectively mxm and nxn) and is a matrix whose only non zero elements are along the diagonal and they are called singular values. In this article, the condition numbers are calculated by the ratio of the greatest singular value over the smallest non zero singular value. In this article, we use very overdetermined systems () which are not rank deficient, the number of non zero singular values is equal to the number of columns: n. Every comparison between the two systems (CACAO and THC) will be done with the same n.

In SPECT the main source of noise is due to photon counting (shot noise) and is well modelled by a Poisson distribution. For a large number of collected

photons, depends on the inverse square root of the number of collected

photons.

Here is an illustration of the meaning of the condition number. Let us take two extreme cases. We will begin by the easiest problem to inverse. For example, a problem whose direct matrix is identity. In this case, whatever the matrix dimension, the condition number will be one. One is the smallest condition number possible. In this simple case, the inverse problem is obviously very easy because no calculations are needed. The data measured is the results, without any further calculations.

To show an example fo a difficult (badly conditioned) problem: Imagine a matrix full of 1 with on the diagonal. A small leads to a very difficult problem. In fact, the smaller the value, the bigger the condition number of the matrix is. Increasing the dimensions of the matrix will make the problem worse. This example is sufficient to obtain a bad condition number, but not necessary. It is even rather extreme because when the system become rank deficient.

4.2. Impulse Response Function in the Total Attenuation Model (Perfect Collimator)

To shorten the notation, all the distances are expressed as a ratio with the width of the pixel of the image, and of the detector (for the sake of simplicity we consider the image pixel having the same width of the detector pixel). All the dimensions are then given by a pure number without units.

We will develop in this chapter an analytical calculation of the response function of the CACAO project. This calculation will be done for a point source placed at the coordinate and, and for a given angle (Figure 5). In this chapter we consider the walls of the collimator, made in an ideal attenuating medium with an infinite power of attenuation. In the next chapter the correction needed to represent a real attenuating medium will be studied.

The response function for other values of will be calculated with a standard rotation of the coordinates. For a point source situated at other values than a simple translation will give the result. Figure 5 has been introduced in a preceding chapter. It shows the 5 zones studied (① to ⑤) separated by 4 special positions (a, b, c, d). These positions of the set collimator + detector correspond to specific values.

First, the extreme lateral zones (① and ⑤ on Figure 5). They are defined for:

and. In these extreme zones, the detector will not

receive any photon (full attenuation model), so. The septa of the collimator do not allow any rays (starting from the point S) to reach the detector.

For the central zone (③ on Figure 5) the whole detector will

light up and the illumination will follow the Lambert’s cosine law and inverse

square law.: where represents a constant value which

depends on the brightness of the source; d is the distance from the source to the detector, and is the angle between the normal to the detector surface and the ray considered. The element of the detector considered is placed at the abcissa. The distance d is given by the Pythagore theorem:.

The angle can be calculated by Replacing these values in the

illumination law:

(6)

This law of illumination is also valid for the two remaining domains (② and ④ on Figure 5).

However, in these domains, the detector is not completely lit up. Notice, for example, the 2nd domain (Figure 8) the detector is lit up under the abcissa and similarly, in the area 4 it is lit up above with:

(7)

The upper part of Figure 5 shows the parallelogram representing the support of the impulse function in the plane for the 3 central zones. The limit of

this parallelogram are: and for the upper and the lower

edges and and (7) for the lateral edges.

4.3. Septal Penetration

To extend the previous calculation in a more realistic model, we study the shadow areas (for example above). This area can be illuminated by rays which have crossed the edge of the collimator in the entrance plane. Let us begin with the case. Figure 9 gives the notations used.

By considering the similar triangles made by the w axis and the ray’s path, elementary geometry leads to the length of the path in the attenuating medium: (). The Equation (6) can be extended in the shadow area by:

(8)

where represents the linear attenuation coefficient in the considered medium. A similar expression for can be calculated for. Using the corresponding for the attenuation of 140 KeV photon by the lead, we observed (Figure 10) a rapid decrease of the intensity of illumination in the dark area. In

consequence, we did not further extend the model to hold the width of the septa. (see section 5.2). We did not pay attention to rays crossing the edge of the collimator close to the detector. With a single hole collimator this event will be irrelevant because the detector and the collimator can be fitted to the same size.

4.4. Direct Matrix Calculation

This section describes the calculation of the direct matrix, in a discrete-discrete form [11] . As we derived a continuous impulse function in the preceding chapter, we have to make it discrete. The emitting source is considered composed of point sources situated in the center of the pixels. This choice is guided by the possibility of high spatial resolution study. At the detector level, we integrate the response function over the width of the detector pixel (model well adapted to semiconductor detectors). This approach is used for both cases.

A program (written in Python) calculates the impulse function. This impulse function, is composed by the response function (6) in the lit up zones and the penetration septal function (8) in the shadow areas. For a given dimension of the image to be reconstructed, the coordinates of all the pixel centers situated in a disc (defined by its radius) are calculated. The pixel size is fixed to 3 mm to match actual gamma detectors. Due to the rotation of the set detector + collimator, it is useless to reconstruct the pixels situated in the corners of the image.

As a simple example, for a 4 × 4 image matrix only 12 pixels are considered (Figure 11). The matrix of the direct problem, will be composed by 12 columns. The columns are sorted as in Figure 11. Let be the matrix which gives the complete projections in digital form. Each row of this matrix will thus correspond to all the possible measurements acquired. These measurements are calculated by integrating the function in the width of the detector pixels. As depends from three variables, we needed to choose an order of these measurements along the column of the matrix. The order chosen is:. That is the inner loop runs on the middle loop runs on and the outer loop on. In a digital form, the Equation (1) becomes:.

We further calculated the extent of the value where is different from zero. Figure 12 exhibits the geometry of the problem.

We introduce the radius of the circular area (Figure 11) and the distance between the entrance plane of the collimator and the center C of rotation (see Figure 12). The range of the possible values is calculated (Figure 12):

(9)

Due to the possible septal penetration, we use a value slightly superior to the one of Equation (9). It is preferable to add rows of zero to the matrix than to work with incomplete data. We verified that rows full of zeros do not affect the condition number. Figure 13 gives a small extract of the matrix represented

in gray scale, with, , , , ,

, ,. Where represents the number of values. The dimensions of the complete matrix is i.e., multiplied by the number of pixels in the circular area. Only 50 rows of the matrix are depicted in Figure 13.

4.5. The THC Direct Problem Matrix

Several possibilities have been proposed to derive the THC-SPECT transform in a discrete form. To be comparable with the previous calculation, we needed a model which took into account the septal penetration and the variation of the collimator response with the source to collimator distance. A Gaussian model was chosen to fulfill a compromise between simplicity and accuracy. We used the experiment driven measurements of Youngho Seo et al. [12] . These authors proposed the following parameters for the Gaussian FWHM for a VPC-45 LEHR collimator and for a Tc^99m source without scatter or attenuation correction in the sources:

(10)

where all the measurements are expressed in centimeters, and where w is the source to collimator distance. As described earlier, the calculation is limited to the points located in a circular area (Figure 11). To match the former description, we evaluate the coordinates of each pixel center in this circular area. The corresponding Gaussian functions are integrated through the boundaries of the detector pixels (fixed to 3 mm). We then obtain a matrix () giving the acqui- sition values () (projections) for each point of the object ()..

An example of matrix for the radon direct problem is given in Figure 14 for the reconstruction of an image size, with and. The dimensions of the complete matrix are but only a limited number of rows is shown in gray scale in Figure 14.

4.6. Simulation and Reconstruction by the Least Square Method

4.6.1. SNR, SNRG, ppp

To illustrate the effect of the condition number, we have simulated some acquisitions of a resolution pattern (64 × 64). The noise free acquisition was calculated for both systems (THC and CACAO) by a simple product with both direct matrices. A gaussian noise was drawn in both acquisitions. Let be the average number of simulated photons per non zero acquisition pixel. The standard deviation of the simulated gaussian noise is chosen to. To study the response at various level of noise, has been varied from 10² to 10¹⁰ by step of 10². In each set the drawing is performed 10 times and averaged to produce the final measurement. For each sample the signal to noise ratio is calculated using the following formula:

(11)

where means the average calculated on the non zero pixels:

(12)

and represents the deviation between the ideal and the degraded signal. is the number of photon in the pixel i.

is the number of photon in the pixel i in the degraded image.

is the number of photon in the pixel i in the ideal image.

This formulae can be applied either to the acquisition or the reconstructed image. For a reconstructed image obtained from a noise free simulated acquisition the degraded image comes from the roundoff errors, and the ideal image is the initial simulated object.

For noise free acquisition-reconstruction this SNR is expressed in decibels (db) is then replaced by the mean number of photon per pixel in the original image,) is the value of the pixel in the original image pattern (), and is the reconstructed image with roundoff errors ():

(13)

For the noisy data, a snrg is calculated by the ratio of the SNR of the reconstructed image () versus to the SNR of the noisy acquisition (or) versus noise free ones. The snrg will be used in chapter 5.8

(14)

4.6.2. Resolution Pinstripe Pattern

For the whole simulation we have chosen a resolution pinstripe pattern define as follow: If is the intensity of the column i of the image,

if i is even

otherwise.

The square image formerly described is then chopped according to and Figure 11. Then the rows are concatenated in a vector.

4.6.3. Acquisition Simulation

The pinstripe pattern is then simply multiplied by both direct matrix to give the noise free acquisitions.

(15)

For noisy simulation, a Gaussian noise is added to each non-zero pixel of the acquisition. However to vary the SNR, the data are scaled by a simple multi- plication. The numbers of non zero pixels () were slightly different for both acquisition data (8960 for the classical acquisition and 19,232 for the CACAO acquisition). See chapter 5.8. The Gaussian drawing has a standard deviation of (pseudo-Poisson).

(16)

But the SNR’s were nearly the same at acquisition level (depending only on the draw). To illustrate that point, the following Figure 15 and Figure 16 show both acquisitions for a pseudo-Poisson expectancy of 10 photons per pixels..

For each simulated acquisition sample we reconstructed an image by the least square method. This method is easily done by the calculation of a pseudoInverse The following formulas are used to calculated both least square estimate:

(17)

(18)

In the last equation, and represent the acquisitions without noise, and with noise.

5. Results

Three parts can be differentiated in the results: in the first we study the condition number with the parameters of the model, with calculation on small matrix size, subchapters: 5.1 to 5.3.

In the second part we make the condition number comparison: THC versus CACAO, with more realistic matrix sizes (up to 64 × 64), subchapters: 5.4 to 5.6.

The final part is devoted to the simulations and reconstructions of noise free and noisy images. 5.7 to 5.8.

5.1. Variation of the CACAO Condition Number with the Collimator Dimensions

Figure 17 shows an example of the variations of the condition number of the CACAO direct problem with the value of D ranging from 2 to 55 for 2 different P values. The condition number increases for larger values of. But this variation appears very low and follows a steep fall for small D values. (, , , ,).

In the following chapters, we will study the variations of the condition numbers for both systems (THC and CACAO) with the common parameters of the problem. We will first focus on the parameters that are important to understand the validity of this paper.

5.2. Variation of the Condition Numbers with the Cut-Off Value

As we have used a Gaussian model for the THC system and an exponential attenuation for the CACAO system, these two functions extend to infinity. These functions were integrated numerically (Python with Scipy and Numpy package). Obviously, very small values (under for example) are meaningless. We thus introduce a cut-off value to chop the results. This decreases the number of non zero elements in our matrices. In addition, this choice speeds up the calculation of the matrix and the condition number. As the condition number involves the smallest non zero singular value, it was compulsory to study the variation of the condition number with the cut-off value. It is to be noted that the cut-off is applied to the impulse response function and that no cut-off was applied to the singular values of the matrix.

Table 1 shows the condition numbers calculated for the two systems for various cutoff. The condition number levels off in a wide range of cut-offs. For higher values of cut-offs the condition number began to decrease () as depicted in the table for the CACAO row. The condition number of the THC began to decrease for bigger cut-off values (10²). In the following calculations the cut-off was fixed to well inside the plateau range for both case (THC and CACAO).

5.3. Variation of the Condition Numbers with the Number of Acquisition Angles

The variations of the condition numbers of the system matrices, versus the number of acquisition angles, exhibit a rather similar behaviour for both systems. As the number of angles increases, the condition number decreases until it reaches a plateau (Figure 18).

Values for this figure: Pixel width = 3 mm, THC: VPC-45 LEHR, CACAO:

, , , , , ,. Number of acquisition angles (), ranging from 8 to 36 by step 2. There are obvious differences between the two systems. For THC a small number of angles is meaningless. Even for a very small number of unknowns, a minimal number of acquisition angles is compulsory in order to have more data than unknowns. With the CACAO project, due to the addition of a linear scanning motion, one angle of acquisition may allow an image reconstruction (Figure 3 and Figure 5). The examination of Figure 18 shows us that the plateau is more rapidly reached in the CACAO approach than with the THC approach. This is, of course, counterbalanced by the scanning linear motion required by the CACAO system. For this choice of parameters, a plateau of value 25.1 is reached with the CACAO approach, while the THC system levels off at 62.8. In the following, the comparison was performed with a large number of acquisition angles, to stay in the plateau.

5.4. Variation of the Condition Numbers with the Orbital Radius

The orbital radius R_g determines the source to detector distance. For both acquisitions (THC or CACAO) increasing R_g increases the condition number as shown in Figure 19.

Pixel width = 3 mm, , , , , Orbital radius (), ranging from 15 to 65 pixels. THC: VPC-45 LEHR,. CACAO:, ,.

5.5. Variation of the Condition Numbers with the Size of the Reconstructed Image

The size of the reconstructed image matrix is obviously a very important factor. The caculations performed so far have been made with a dramatically small image size, for the sake of simplicity. For real applications, however, enormous matrix size becomes compulsory. Figure 20 exhibits the influence of the image size, while ranging from to. While the image radius follows the

image matrix size with the formula:.

R_g follows the formula:.

The number of acquisition angles for the THC configuration was chosen high:. Calculations have been made to check that this number was high enough (variations of the result less than between the last 2 steps while increasing the number of angles).

Pixel width = 3 mm, , to, to 3196, to 31.9, THC: VPC-45 LEHR, CACAO:, ,.

With a 64 × 64 matrix size the number of columns of the matrix was 3196 (number of pixels in the circular cache).

Table 2 gives the values corresponding to Figure 20.

5.6. Study of the Spectrum of the Singular Value Decomposition

Figure 21 was calculated from the complete spectrum of the singular values for the classical THC approach and the CACAO project. This figure does not represent the singular values themselves but rather the ratio of the greatest

singular value over the ith.. As the range of variation of these ratios is

very large we have used a logarithmic scale. The singular values are sorted in descending order. Therefore the greatest is the first () and both curves begin with a value of 1 (zero on the logarithmic scale). The condition number of the whole matrix is given by the right hand side extremity of the curves. Both curves have the same number of points (3196). The THC curve is lower or equal to the CACAO curve for the first 2252 values, then climbs to its highest value.

Figure 22 shows the position of the cross point which is where the two curves on Figure 21 meet. Figure 22 shows the cross point position variation depending image size. The position is measured as a percentage of the matrix rank (number of singular values). The curve begins by a small climb, due to the limited number of singular values. For example the point of the curve, only the last singular value ratio for THC is below the CACAO one, and therefore the percentage is limited by. Then, the curve reaches its maximum level and begins to show a clear decreasing trend.

5.7. Reconstruction Results from Noise Free Acquisitions

In the following two sections the parameters of the collimation and of the acqui- sition will remain fixed. These parameters are given in the following Table 3.

The same digital pinstripe resolution pattern (64 × 64) (24) is first multiplied by both direct matrices (THC and CACAO), to give the 2 (noise free) acquisitions. The two estimates are calculated by a simple matrix multiplication with both pseudo Inverses (18). Due to the relatively small size of the image and the absence of noise, both results are excellent. These two estimates are not depicted here, because the human eye cannot detect difference from the original image (24). However, the computer can: both SNR in decibels are given in the last row of the previous Table 3.

5.8. Noisy Simulations Results

For the two chosen systems, (whose parameters are given in Table 3) and for 5 choices of signal to noise ratio at the acquisition level, we simulated noisy acquisitions then we reconstructed the corresponding least square estimates. By the same token, for both systems, the acquisition is different and simulated noises are different. To average these necessary differences, 10 simulations of noise were drawn at each noise level. for both systems. Two very noisy acquisitions are depicted in Figure 15 and Figure 16. For each of these 100 simulations (5 × 10 × 2) we calculated the gain in SNR by formula 14. The different values of snrg are depicted in Figure 23.

Even though some fluctuations are clearly visible in the different simulations, there is a clear gap between the two systems. The following Table 4 gives the

mean of all the snrg for both systems.

Figure 24 depicts 10 reconstructed images next to the original pinstripe pattern.

6. Discussion

6.1. Condition Number in the Literature

The value of the condition is well described in the inverse problem literature and linear algebra textbooks [10] [13] . The SPECT and PET literature has also recognized and used this tool [14] - [20] . Our results, for example, agree with the work of [19] which shows that a sufficient number of acquisitions angles in THC are necessary to reduce the condition number. These authors showed that the system even becomes rank deficient if the number of projection angles () is too low. As depicted in this article the result of the least square reconstruction is well predicted by the condition number of the direct matrix. In addition to this, for a standard regularization procedure by Tikhonov and Phillips [21] [22] the condition number is also the parameter to consider, as claimed by Hansen [23] . In this last article the author emphasized the similarities between the Tikhonov method and the truncated singular value decomposition (TSVD) which can also be used to regularize the inverse problem. Here again the condition number is of paramount importance. The TSVD approach was applied by Jorgensen and Zeng in [20] to multipinhole SPECT. The results, depicted in Figure 5 of their article, comparing the TSVD version to the global condition number do not show dramatic differences between both approaches. In addition, in the second part of this article these authors used only the condition number and advocated the fact that it is easier to calculate. We can thus conclude that even if the condition number is not a perfect and universal way to evaluate the difficulty of an inverse problem, it does exist, it is simple, it is general and it is particularly accurate.

6.2. Why Have Not We Used another Performance Index?

The medical imaging literature is full of papers intended to measure the performance or the quality of a tomographic image. Some papers even try to evaluate the usefulness of the images related to a specific task [24] . Some of these performance indices apply to the reconstructed images: we can cite the Standard Deviation [25] or the Contrast Recovery Coefficient (CRC) [26] or the SNR gain [27] . But all these performance measurements do not apply to the same thing. The condition number applies solely to the geometry of the system, meaning, the system direct matrix. With the exception of the last 2 subchapters we did not mention images. The condition number applies whatever the image studied, and in addition applies whatever the reconstruction algorithm used. It is simply an intrinsic property of the transfer matrix. In the last two subchapters we needed to choose an image (pinstripe resolution pattern), a noise distribution and a reconstruction algorithm. In this chapter we use the snrg in a similar manner to [27] . With more hypothesis (the type of lesion to detect in the observer model) we can go further with the specialization, but with a narrower result. This was not the purpose of this article which tries to compare 2 geometries with a minimum of extra hypothesis.

6.3. So, What Does the Results Mean?

In chapter 3 we introduced 2 questions. The first was: Do the photons collected with CACAO carry the same information quality as the photons collected with THC SPECT in terms of location?

If the location were imprecise, the direct matrix would be rank deficient. The fact that the rank of the matrix equals the number of unknowns means that, in the absence of noise, the reconstructed image would be perfect. In linear algebra we say: the matrix has full rank. This result is the consequence of the added linear scan in the acquisition sequence. This ensures an overdetermined linear system. Hence the system is accurate in terms of location (comparable to the THC-SPECT).

We introduced the condition number to evaluate the difficulty of solving an inverse problem, which helped us respond to the second question. The graphs speak for themselves: when the dimension of the system increases and when rotation radius increases, the CACAO project presents a smaller condition number than the THC-SPECT.

6.4. Can the Geometry Explain These Results?

The limited performance of the THC SPECT is probably due to the gaussian shape of the response. In fact it is quite difficult to distinguish the difference between the sum of 2 Gaussians, shoulder to shoulder, from a sole gaussian slightly taller. And when we have to separate a large number of Gaussians of various widths the problem becomes worse. Moreover, these Gaussian responses differ only slightly with the source to collimator distance. On the contrary, the added scanning motion combined with the increased depth of the collimator hole produce responses which are easier to distinguish. This difference leads to the difference in condition number.

People familiar with the difficulty of solving deconvolution problems may find our results strange. Large kernel deconvolutions have a reputation for ill-posedness and large condition number. Several years ago we published a reconstruction algorithm for the CACAO problem with a deconvolution step [6] . A transformed matrix of the system led to a deconvolution problem. Condition number calculations of this transformed matrix led to us abandoning this method. Effectively, the use of the complete and direct response matrix (without a deconvolution step) leads to a better conditioned problem.

The difficulty of the conventional tomographic inverse problem is often underestimated. It is this difficulty which allows us to believe that the approach of the rotating slat [5] may be not optimal. Although this rotating slat can dramatically increase the acceptance angle of the collimator hole (at least in one direction), it needs a double application of the inverse radon transform. Considering the bad conditioning of the inversion of only 1 stage radon transform, we have some doubt as to the quality of the reconstruction. Another quite different approach is the one proposed by Zhang [4] . In this work, the increase of the acceptance angle (referred as cone angle in this article) is very weak: 6˚ for the largest collimator hole considered. In addition, these authors reported that “The blurring of this last collimator hole is too severe to be corrected…”. This may be related to the fact that this article uses a gaussian shaped response, which is hard to deconvolute.

6.5. What Are the Consequences for the Noise Amplification?

Not only, is the condition number related to the geometry of the tomographic model, it is a property of the response matrix of the model. The direct consequence is that for the same, then for CACAO will be lower than for THC. This is assuming a proper choice of the system parameters for CACAO was made, and for the same image dimension, with the same source to collimator distance.

6.6. What Can Be Said about the Simulated Images?

To illustrate this noise amplification, we conducted some simulations. These simulations raise a number of questions. First the choice of the pseudo inverse reconstruction method (unpenalized least square) may be criticized because, as already mentioned, it is not the best method for solving difficult inverse problems. It is known to amplify the noise. It was chosen here for several reasons. First it is a well known method, and intensively theoretically studied. Second it is quick, simple and has great efficacy in exact or near exact data. Everyone can judge it by looking at the simulations with a high number of ppp. Finally every regularized method which would have given better results at low SNR would have needed the choice of a frequency cut-off, a regularization parameter or a stopping criteria for an iterative method. The OSEM (which is the accepted method in emission tomography) not only needs the choice of the number of iterations and the number of subsets, but it has been developed for the THC-SPECT, meaning for a large number of acquisition angles. Such a big number does not hold for CACAO and a comparison based on this algorithm may have been unfair to CACAO. For all these reasons, reinforced by the fact that the condition number has been calculated in the L2 Norm, we chose a least square estimator. A Gaussian noise is then a choice consistent with the least square reconstruction. We are well aware that in emission tomography the noise is better modelled by a Poisson statistic. We chose a value of the standard deviation of the Gaussian noise, because it is known that the Poisson statistic tends towards the Gaussian one as soon as the number of photons per pixel (ppp) exceeds 10⁴. With these choices the result presented in Figure 24 may at first glance appear a little gloomy. Quite the opposite, it must be emphasized that the pinstripe image test with 3 mm width for the line is extremely difficult to observe with a standard gamma camera equipped with a thin hole collimator. So the result can only be attained with a very high signal-to-noise ratio for the acquisition and a well adapted model of the point spread function of the collimator. Table 4 perfectly sums up the concept of condition number. From Equation (5) and Equation (14) we can see the resemblance:

(19)

Figure 23 and Figure 24 show the power of the condition number, predicting a lower snrg for the THC SPECT than for the CACAO project for a large range of SNR at the acquisition level. Table 4 gives a ratio for the 2 snrg of 15.87 compared to 41.84 for the ratio of the condition numbers.

6.7. What Can We Say about Highly Regularized Reconstruction?

Our simple model allows us to go further in the comparison. We can effectively predict the performance of a regularized method. This performance is depicted in Figure 21. The condition number related to a TSVD truncated at the ith

element is given by the ratio:. The explanation of the figure is: if one gets

rid of half the singular values, the THC is preferable. On the contrary, if one wants to collect the maximum amount of detail and precision in the images, and tries to work with more singular values, then the CACAO approach would be preferable. Keeping a large amount of singular values may seem foolish right now, due to the huge uncertainties caused by the use of the thin-hole-collimator. In fact, the CACAO approach implies more collected photons, meaning a better signal to noise ratio at the acquisition level. Moreover, as shown in Figure 22 the point where the curves cross tends to shifts to the left when the size of the image increases. In other words for big data (more realistic) the THC will perform better than CACAO, but only with dramatically strong regularization (e.g. very low spatial resolution).

6.8. And What about Real Systems?

This study was done analytically, with an approximated and simple hypothesis. The attenuation or the scattering of the gamma rays in the body of the patient, or the scattering in the collimator, were not taken into account in this simplified study. We are aware that nowadays advanced methods of reconstruction in SPECT use much more complex models [28] [29] . It should be remembered, however, that the most impairing phenomenon in SPECT is the photon counting noise (sensitivity of a conventional collimator, which means that for each detected photon 10,000 are lost!)

El Fakhri in [30] has shown that the second most important factor influencing the contrast and the spatial resolution of the images is the collimator response. It is to be noted that our work proposes improvements for both of these factors.

The feeling that scattering would impair the CACAO reconstruction more than the THC reconstruction was opposed to our previous works. Obviously, a bad condition number could magnify the noise more than the signal. Smith [31] has shown that attenuation and scattering degrade the condition number of the SPECT problem. But there is no reason why this phenomena would degrade the CACAO problem more than that of the THC. Faced with the improved condition number at the start, the opposite seems more likely. In fact, regarding our condition number results, all the noise sources (in the patient) impairing conventional THC are expected to be reduced in the CACAO project.

Another point that may be controversial in this study is the lack of a complete 3D calculation. We did not try it due to the time consuming burden of the calculation, added to the complexity of the programming. In addition, this simple 2D calculation allows for a better understanding of the involved phenomenon combined with a simpler calculation, and so is simpler to verify. It also allows a simple analysis of the whole spectrum of singular values. Even without a true 3D study, this article presents an interesting message because we can already project an application of this article in a real 2D experiment with a slit hole.

6.9. The Problem of DOI and High Count Rate

There are 2 physical phenomena which can be predicted as being worse for CACAO than THC. One is the depth of interaction in the detector. Obviously with an angle of 45˚ for the extreme rays the intrinsic spatial resolution of the detector may be degraded. At 140 KeV, 5 mm of CdTe absorb 85% of the incoming photons. If we reduce the thickness by 3 mm we can keep a 3mm spatial resolution at the price of a slight reduction in sensitivity (from 85% to 70%). Another possibility is a direct measure of the DOI which is a fashionable subject with semiconductor detectors [32] [33] [34] .

The other physical limitation is the high count rate which comes directly from the dramatic increase in sensitivity provided by the new collimator. Obviously a lot of gamma ray detectors present a decrease in their performance with increasing count rates. Here again the semiconductor pixellated detector is less subject to deterioration than the large NaI scintillator which could only count one gamma at a time.

6.10. Can We Predict the Future?

This article is purely theoretical and it will not replace a fine Monte Carlo simulation or a real experiment. It is of course difficult to predict the future. Even the result of an experiment does not guarantee the future. We should remember that the first electronic microscope did not provide a better spatial resolution than the existing optical microscope. Today the spatial resolution of the electronic microscope is 10⁴ times smaller. It is also worth mentioning that the neutron, the positron, the rounded shape of the earth and some stars had been theoretically predicted before they were discovered and universally accepted. So we can let the reader believe what they want and we just offer the following Table 5 to sum-up the pros and cons.

7. Conclusion

The mathematical inverse problem of SPECT with THC is hard to solve. The CACAO project may not only improve the signal to noise ratio at acquisition level, but it may also facilitate the calculation of a good reconstructed image and the use of a more complete set of data. To put it differently, this article seems to demonstrate that the photons recorded by the CACAO project carry better localisation information than that carried by the THC photons. For a low resolution image, with a small image size and a small gyration radius, the THC system is the best method.

Aknowledgements

This study was supported by the “Pays de Loire” MPIA Grant. The authors wishes to thank all the students who have worked on this subject and the colleagues of the UA who have helped us.

Disclosures

No disclosures.

Cite this paper

Jeanguillaume, C., Simonnet, R., Cavaro-Menard, C., Douiri, A., Bouali, I. and Loeb, J.J. (2017) An Alternative Model to Radon Transform for Gamma Ray Emission Tomography. Advances in Molecular Imaging, 7, 13-47. https://doi.org/10.4236/ami.2017.72002

References