Optimized Bayes variational regularization prior for 3D PET images

Accepted Manuscript Title: Optimized Bayes variational regularization prior for 3D PET images Author: Eugenio Rapisarda Luca Presotto Elisabetta De Bernardi Maria Carla Gilardi Valentino Bettinardi PII: DOI: Reference:

S0895-6111(14)00064-0 http://dx.doi.org/doi:10.1016/j.compmedimag.2014.05.004 CMIG 1264

To appear in:

Computerized Medical Imaging and Graphics

Received date: Revised date: Accepted date:

4-1-2013 21-2-2014 2-5-2014

Please cite this article as: Rapisarda E, Presotto L, De Bernardi E, Gilardi MC, Bettinardi V, Optimized Bayes variational regularization prior for 3D PET images, Computerized Medical Imaging and Graphics (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.05.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Optimized Bayes variational regularization prior for 3D PET images Eugenio Rapisarda a,b,c), Luca Presotto b,c), Elisabetta De Bernardi b,d), Maria Carla Gilardi a,b,d) and Valentino Bettinardi b) a)

IBFM-CNR, Institute for Molecular Bioimaging and Physiology, 20090 Segrate, Italy Department of Nuclear Medicine, Scientific Institute San Raffaele, 20132 Milano, Italy c) Department of Physics, University of Milano-Bicocca, 20126 Milano, Italy d) Department of Health Sciences – Tecnomed Foundation, University of Milano-Bicocca, 20900 Monza, Italy

cr

Eugenio Rapisarda (corresponding author) Tel. +39 (0)2 26435175, fax +39 (0)2 26415202; email: [email protected]

ip t

b)

us

Luca Presotto Tel. +39 (0)2 26435175, fax +39 (0)2 26415202; email: [email protected]

an

Elisabetta De Bernardi Tel. +39 (0)2 26435175, fax +39 (0)2 26415202; email: [email protected] Maria Carla Gilardi Tel. +39 (0)2 26432793, fax +39 (0)2 26415202; email: [email protected]

M

Valentino Bettinardi Tel. +39 (0)2 26433887, fax +39 (0)2 26415202; email: [email protected]

d

Header title: Optimized variational regularization prior for PET images

Ac ce p

te

Abstract A new prior for variational Maximum a Posteriori regularization is proposed to be used in a 3D One-StepLate (OSL) reconstruction algorithm accounting also for the Point Spread Function (PSF) of the PET system. The new regularization prior strongly smoothes background regions, while preserving transitions. A detectability index is proposed to optimize the prior. The new algorithm has been compared with different reconstruction algorithms such as 3D-OSEM+PSF, 3DOSEM+PSF+post-filtering and 3D-OSL with a Gauss-Total Variation (GTV) prior. The proposed regularization allows controlling noise, while maintaining good signal recovery; compared to the other algorithms it demonstrates a very good compromise between an improved quantitation and good image quality. Key words—3-D image reconstruction, point spread function, image regularization, positron emission tomography (PET).

Introduction Positron emission tomography (PET) and PET/computed tomography (CT) imaging systems have been greatly improved in the last few years both from the hardware and software point of view. For example, all the physical effects limiting the spatial resolution of PET can be taken into account by a “global” Point Spread Function (PSF), as representative of the system response, which – if implemented in an iterative reconstruction algorithm – allows taking into account the system resolution loss and recovering it [1-4]. Unfortunately, a common effect of iterative reconstruction techniques is the increase of noise as iterations proceed, due to the ill-posed nature of the reconstruction problem [5]. On the other hand, a high number of iterations is usually needed to recover a substantial percentage of the signal, especially when implementing PSF modelling [6]. A post-filtering applied to the image can limit noisy patterns, but it usually causes a loss Page 1 of 20

Ac ce p

te

d

M

an

us

cr

ip t

of spatial resolution and a reduction in the quantification accuracy, which should obviously be avoided, particularly if a strategy of spatial resolution recovery is employed. Regularization techniques, instead, could help in controlling noise amplification during the reconstruction process (and, thus, they allow increasing the number of iterations), provided that the regularization method is able to preserve the spatial resolution as much as possible. Many different regularization strategies have been proposed (see e.g. [7] for a review). The most widely used ones introduce a spatial interaction between voxels in the image – often by using Markov random field models – and thus associate, to the image, a global energy function whose value is based on local potentials evaluated on some subsets of the entire image. The effects of the regularization process depend on the characteristics of the local potential employed, which ideally should provide a smoothing behaviour inside each region and preserve the transitions between adjacent regions. Among the different strategies to define the energy function, the variational approach [8, 9] allows taking into account not only the mean activity in the considered districts of the image, but also their rapidity of variation between adjacent regions by employing the gradients in the image. Many different penalty functions have been proposed [10]. The most used Bayesian regularization strategy is based on the Gaussian prior, which smoothes all the gradient intensities with equal relative strength, resulting in a strong smoothing effect. On the other hand, a commonly used prior for edge preservation is the Total Variation (TV) [11], which is able to maintain mainly unchanged the transitions between different regions in the image, but tends to produce staircasing effects in regions gradually changing in space. A modification of the Gaussian prior is the Huber prior [12] or, equivalently in the variational framework, the Gauss-Total Variation (GTV) prior [10], which introduces a parameter to discriminate between different regularization behaviours to be applied respectively to background regions – Gaussian (GR) component – for noise suppression and signal regions – Total Variation component – for edge preservation. Another promising strategy is represented by the generalized Gaussian or p-Gaussian prior (PR), as proposed in [13]: this regularization scheme has the capability of strongly smoothing background regions, while maintaining higher detail in signal regions. In this work a new variational regularization prior has been proposed by modifying the PR to reduce the resolution loss introduced by it. The prior has been implemented in a 3D Ordered Subsets Expectation Maximization (OSEM) reconstruction algorithm accounting for the PSF in the image space and validated, on phantom data and clinical images, by comparing it with OSEM+PSF images, post-filtered OSEM+PSF images and the images obtained with the Gauss-Total Variation prior. Moreover, since the presented regularization priors depend on different parameters, an optimization strategy (based on the qualitative and quantitative content of the images) has also been proposed and employed to compare the priors at their best performance. Materials and methods

Bayesian regularization and variational approach Most iterative algorithms seek the image  maximizing the so-called likelihood function, i.e. the probability p y  of obtaining the recorded projections y given image  . Since all images are in principle equally

 

probable, due to the ill-posed nature of the reconstruction problem the process often converges to a noisy image, with the noise level increasing at each iteration. To reduce this effect, it is possible to favour a certain class of images by using Bayes’ theorem and introducing an a-priori probability p   for the image  . The Bayesian formulation usually leads to the computation of a maximum a posteriori (MAP) estimate of the image as the maximizer of the Bayes’ probability

p  y  

p  y   p   p y 

The a-priori probability p   is often chosen with a functional form following Gibbs’ distribution

p(  )    Ne  U  

where  sets the contribution of the regularization prior on the final result and U   is the Gibbs’ energy function. ( N is a normalization coefficient to let the Gibbs’ distribution be interpreted as a probability). Page 2 of 20

Usually the energy function is chosen as a sum of potentials, often defined on pairwise cliques of neighbouring voxels b of the image N

U    

  

b 1 k N b

bk

b

 k 

In this work, instead, a variational approach has been followed. The Gibbs energy function has been chosen of the form

ip t

U        d 

with    0 is the prior function (or, in the following, simply the prior) and  spans the image domain  .

cr

The symbol  indicates the gradient and  the  2 norm, while      is a continuous function

us

describing the activity distribution in  . Taking the logarithm of the a posteriori probability and neglecting constants, the Bayesian formulation is reduced to the maximization of the function

J    L y         d 

 

an

where L y  represents the log-likelihood.

[k ]

λb  1  D λ[ k ]



 

d

λb

[ k 1]

M

Since the log-likelihood functions used in PET image reconstruction (derived from Poisson model) are concave, J   will be concave if U ( ) is convex: this is guaranteed if    0 . In this case, the required maximization can be obtained by using standard nonlinear optimization methods, such as the ExpectationMaximization (EM) algorithm [14] with “one-step-late” (OSL) method [15]. The iterative updating rule, after k iterations, for the OSL-EM with a variational approach may be written as

Bb

b

yd Pd b

[k ]

where b[k] indicates the content of the voxel b in the image  at iteration k, yd are the counts recorded along line of response (LOR) d, the “projection” operator (projector) is defined as Pd () b  () b p bd , the

te



Ac ce p

“backprojection” operator (backprojector) is defined as B b () d 

 ()

b

d

pbd and

d

D   

 U             

   

results from the variational maximization strategy of the a posteriori probability. To increase computational efficiency, in the implementation described below an Ordered Subsets Expectation Maximization (OSEM) algorithm [16] has been used, i.e. the sum in the backprojector is calculated at each iteration only on a subset of the acquired entire dataset of projections. The functional D  has been discretized as described in [10]. In this work a Maximum a Posteriori (MAP) OSL [11] multiplicative approach, as presented in [17], was used to implement a PET image reconstruction algorithm also accounting for different physical corrections:

λb

[ k 1]

[k ]

λb  B b Ad   1  D λ [ k ]  B Ad







b

Bb

Ad y d Ad Pd b

[k ]

 Rd  S d

where Ad is the attenuation correction factor relative to LOR d, yd are the counts recorded along LOR d (precorrected for normalization, decay and dead time), while Rd and Sd are the estimations of random and scattered coincidences relative to LOR d, respectively. The integrated PET/CT system used in this work was a Discovery-690 (D690) by General Electric Medical Systems (GEMS, Milwaukee, WI, USA) [18]. This integrated system has also the capability of acquiring data in Time of Flight (TOF) modality, but for this work this capability was not used, to focus only on the specific effects introduced by the variation regularization. The spatial resolution of D690 was measured (and the projector and backprojector were modified accordingly) following the method specified in [4] to take the PSF of the PET system into account in the reconstruction process. Page 3 of 20

On the final image the effects of a variational regularization strategy depend on the mathematical characteristics of the prior  x  . In particular, convex functions (   x   0 ) smooth the image (with

increasing effect for larger values of   x  ), lowering the noise level but reducing the sharpness of the edges

in the image, while concave functions (   x   0 ) enhance the edges (with increasing effect for larger

values of   x  ), even if leading to possible amplification of noisy textures and creation of “patchy”

us

cr

ip t

artefacts. Both these behaviours are, in principle, contemporarily needed for PET, since a reconstruction algorithm for PET should be able to smooth uniform areas (in order to take noise under control, all the more – if the algorithm is iterative – while iterations proceed) and, contemporarily, to preserve as much as possible the spatial resolution in the reconstructed images. For this reason, the case with   x   0 has been subject of interest, since it corresponds to the so-called “total variation” function, which neither smooths nor enhances edges: this specific penalization provides an excellent edge preservation and smoothing of flat regions. Unfortunately, in regions with gradual variations (as often encountered in PET) it leads to unnatural staircasing effects. Many different penalty functions have been proposed [10]. The most used Bayesian regularization strategy is based on the Gaussian prior, which corresponds to the variational penalty function G  x   x 2 2 . Since

an

G  x   1 and considering that the function is used with argument x    , the prior smoothes all the

gradient intensities with equal relative strength, resulting in a strong smoothing effect. A modification of the Gaussian prior is the Huber prior [12] or, equivalently in the variational framework, the Gauss-Total Variation (GTV) prior [10]:

M

 x2 2 x    G x  x   GTV x     2  x   2 x   TV  x  x  

where the parameter  is a threshold to discriminate between the GR component – G x   1 for low x – for

Ac ce p

te

d

 x   0 for high x – for edge noise suppression in background regions and the TV component – TV preservation in signal regions. The generalized Gaussian or p-Gaussian prior (PR) [13] corresponds to the variational penalty function  P  x, p   x p p , with 1  p  2 . Since its second derivative is  P  x, p    p  1x p 2 , this regularization scheme has the capability of strongly penalizing low gradients, while reducing the smoothing effect on higher gradients. The proposed variational prior In this work a new prior was proposed as a generalization of  P  x,4 / 3 :

 3x 4 / 3 4 x   x     x  d d   ln x  d  c x  

with

c  d d   ln d        4 / 3 4



d  1 3 



to have  x  and  x  continuous in x   . Consequently, 3  x   x    x  d 1  d     xd

x  x 

and

1 3 27 x 2   x    d d      x  d 2 

x  x 

Page 4 of 20

ip t

The value p  4 3 was heuristically chosen by clinical evaluation of whole body PET images reconstructed with an algorithm using PSF and the generalized Gaussian prior, in which different values of the parameter p , spanning the range 1  p  2 , were used; the value p  4 3 gave the best compromise between preservation of spatial resolution, noise suppression and natural appearance of the clinical image. This value, which is also in agreement with indications from other studies (e.g. [19, 20]), was thus chosen as a reference point for the modified prior. The parameter , as described above, represents a threshold to discriminate between different regularization behaviours. Below the threshold ( x   ), since   x   0 the proposed prior smoothes with a strength independent of

 

the threshold but depending only on x. Moreover, in the interval 0  x   (   1 3 3 ), this smoothing

cr

effect is greater than the GR one (since   x   G  x  ) and, all the more, than the TV one (since

 x  ), while in the interval   x   the smoothing effect is intermediate between GR and TV.   x   TV

M

an

us

Above the threshold ( x   ), instead, the behaviour depends on the value of . If   1 , the prior smooths the image (   x   0 ), but less than the GR prior and with decreasing strength for larger x. Consequently, the result is a preservation of most edges, while globally maintaining “natural” transitions thanks to the very small smoothing effect. If   1 , the prior emphasises the edges (   x   0 ), contributing to a gain in spatial resolution. Consequently, the result is an enhancement of the contrast between regions of large x and regions with low x. In this case, however, it is important to note that the prior is not convex anymore. Both for   1 and   1 , for x   ,   x   0 , with a consequent asymptotic TV behaviour. In regions above the threshold, for every choice of  , the proposed prior smoothes less than the corresponding p-Gaussian prior  P  x,4 / 3 , while below the threshold the behaviour is identical.

d

Since the penalty function is used with argument x    , the strong smoothing is applied in background

te

regions (in which  is small), while the “natural” edge preservation (or edge enhancement) is used in

Ac ce p

signal regions (i.e. regions with high  ).

A summary of the prior characteristics is provided in table 1, while in figure 1 a graphical comparison of   x  and   x  for the Gauss-Total variation, the p-Gaussian and the proposed priors is presented.

 

0

min( ,  ) max( ,  )

  

  1

 1

 1

Smoothing stronger than GR and TV, equal to PR Smoothing lighter than GR, stronger than TV, equal to PR

Edges nearly preserved, with small smoothing effect lower than GR and PR

Edge preservation (TV)

Edge enhancement

 

Table 1. Summary of the new prior’s characteristics for different choices of    1 3 3  0.192450...

Page 5 of 20

ip t cr us an M d te Ac ce p Figure 1. Comparison of   x  and   x  for the Gauss-Total Variation, the generalized p-Gaussian and the

 

proposed priors for different choices of    1 3 3  0.192450... .

Page 6 of 20

Both the GTV and the proposed regularization schemes are based on two parameters:  (i.e. the contribution of the regularization to the resulting image) and  (i.e. the threshold to discriminate between background and signal regions). These two parameters have to be set properly to obtain the desired effects on the final image. In this work the tuning of the regularization parameters (andhas been performed by using an index of image quality (which can take into account both qualitative and quantitative information) to analyze reconstructed images obtained with different sets of parameters. The couple of parameters maximizing this index was considered the optimal one.

D  Cqnt C qlt

Cqnt  CR  100

S / B 1 R 1

cr

The quantitative factor was chosen to be the contrast recovery coefficient (CR)

ip t

The proposed index of image quality, here named detectability (D), was chosen equal to the product between two factors, Cqnt and Cqlt, representing respectively a quantitative and qualitative figure of merit:

S  B S  B

M

C qlt  ln

an

us

where R is the real signal-to-background ratio and S (B) is the mean content of the voxels inside a spherical volume of interest (VOI) drawn in the signal (background) region. The qualitative factor, instead, models the discrimination between signal and background regions (with mean values of activity concentration respectively of S and B) characterised by standard deviations of counts (chosen as a representation of “noise” perceived by the human operator) S and B, respectively. It assumes a logarithmic response of the human eye [21, 22]:

The complete expression for the detectability is therefore

D  100

S / B 1 S  B ln R 1  S  B

Validation of the D index

te

d

Larger values of D should correspond to better definition and visibility of the lesions and, thus, a higher probability of detection by a human operator.

Ac ce p

The index D was validated on experimental acquisitions of the NEMA IEC body phantom. It consists of a body (180 mm long), a cylindrical insert simulating lung (external diameter 51 mm) and an insert with six fillable spheres (inner diameters: 10 mm, 13 mm, 17 mm, 22 mm, 28 mm and 37 mm). Both the body and the spheres were filled with a 18F-FDG solution and two ratios between the activity concentrations in the spheres and in the background (S/B) were considered: 5 (high contrast, HC) and 3.3 (low contrast, LC). After centring the phantom in the scanner FOV, a CT scan for the attenuation correction was performed, followed by a 3D list mode emission scan. For each S/B ratio the raw data List-file was processed to generate two sets of data corresponding to 25 and 50 Mcounts. Each set of data was reconstructed with four different algorithms:  OSEM with PSF modelling (RwPSF).  OSEM with PSF modelling and a clinical post-filter (RwPSF-Filt). The post filter was a symmetric two-dimensional Gaussian with FWHM=4.28 mm in the transaxial plane and a three-point mean along the axial direction.  OSL-EM with PSF and Gaussian Total Variation Regularization (RwPSF-GTVR). Four sets of  and parameters were considered: =0.0025, =0.1; =0.025, =0.1; =0.0025, =1; =0.025, =1.  OSL-EM with PSF and the proposed Regularization (RwPSF-R). Four sets of  and parameters were considered: =0.0005, =0.2; =0.005, =0.2; =0.0005, =2; =0.005, =2. The choice of  and  for the two regularized reconstruction methods was done with the aim of exploring a large set of possibilities (low  low , high  low , low high  and high  high ) corresponding to images with different image quality. Common parameters for all the algorithms were: image matrix = 256x256, subsets = 18, reconstructed FOV = 60 cm, number of iterations = 1-50. Page 7 of 20

us

cr

ip t

The reconstructed images were analyzed by using 2D regions of interest (ROIs). For each sphere, a 2D ROI with the dimensions of the corresponding diameter was defined, while a big ROI was drawn in the background area. The choice of using two-dimensional instead of three-dimensional ROIs was motivated by the necessity of directly correlating the index D with the results of a visual inspection performed on 2D transaxial images. Two kinds of evaluations were performed: a) Comparison of D with the best image quality within the same set of data over the iterations The 40 sets (2 S/B x 2 count statistics x 10 algorithms) of 50 images (iterations 1:50) were independently analyzed by two expert operators. The operators were asked to choose the iteration or the range of iterations corresponding to the best lesion detectability, with particular reference to the smallest sphere. The evaluations of the two operators were compared with the measured D values. b) Comparison of D with the best image quality over different images set at the same iteration For each S/B ratio and count statistics, couples of images corresponding to the same number of iterations but obtained with different reconstruction algorithms and/or algorithm parameters were generated. The operator was asked to choose the best image inside each couple in terms of lesion detectability, again with particular reference to the smallest sphere. A total of 96 couples of images were evaluated by each operator. The evaluation of the operators was then compared with the computed D values.

an

Since both evaluations showed a good agreement between D and the operators’ analysis (see Results section), the D index was used to optimize the parameters  and of the regularized reconstruction algorithms.

M

ptimization of the parameters  and 

Ac ce p

te

d

The optimization of the  and  parameters to be used in the regularized algorithms was performed by reconstructing an acquisition of a NEMA IEC Body Phantom (acquired statistics: 52.4 Mcounts and S/B ratio 4.4:1 – see section “Validation of the D index” for further details about the phantom configuration) using different sets of parameters  ,   and choosing the set which maximized the detectability index D relative to the smallest sphere (diameter: 10 mm), i.e. representative of many oncological lesions. Larger lesions would suffer from less spatial resolution loss due to the regularization process, with consequent better results in terms of image quality and contrast. The optimization was performed using 10 iterations, considered as a good compromise between clinical requirements of reconstruction time and recovery of spatial resolution achieved by PSF modelling. In figure 2 (top) the detectability for the smallest sphere is plotted as a function of the parameters  and  for the proposed prior. The maximum detectability was obtained for   0.002 ,   0.3 With a similar strategy, it is possible to obtain the optimal parameters for the GTV prior. In figure 2 (bottom) the detectability for the smallest sphere using the GTV prior is plotted as a function of the parameters  and : the maximization was obtained for   0.015 ,   0.2

Page 8 of 20

ip t cr us an M Validation of the proposed prior

te

d

Figure 2. Detectability after 10 iterations for the smallest sphere (diameter: 10 mm) in the NEMA IEC Body phantom as a function of  and parameters for RwPSF-R (top) and for RwPSF-GTVR (bottom).

Ac ce p

The proposed prior has been validated quantitatively and qualitatively. In both cases, the validation consisted of a comparison among four different reconstructions: RwPSF, RwPSF-Filt, RwPSF-GTVR (with optimized and  parameters) and RwPSF-R (with optimized and  parameters) after 10 iterations. Phantom studies: Quantitative validation A quantitative validation was performed by using the NEMA IEC Body Phantom in two different configurations. The two configurations simulate representative oncologic whole body [18F]-2-fluoro-2deoxy-D-glucose (18F-FDG) studies. In the first experiment the phantom was used in its standard configuration with the 6 spheres, previously described in the Validation of the D index section. Phantom body and spheres were filled with a 18F-FDG solution in order to obtain a S/B ratio of 3.8. Two 3D emission scans (25 and 75 Mcounts) were acquired. In the second configuration the original set of spheres was substituted with a new one having diameters of: 16, 13, 10, 8, 6 and 5 mm. The body and the spheres were filled with a 18F-FDG solution to obtain a S/B ratio of 6.5. A 3D emission scan of 50 Mcounts was acquired. In both the experiments the cylindrical insert simulating the lung was removed from the phantom. Data were reconstructed with: image matrix = 256x256, subsets = 18, reconstructed FOV = 50 cm. For the quantitative analysis nine spherical VOIs were drawn, six centred on the spheres and three in the background region. The radii of the spheres’ VOIs were chosen according to the specifications of the manufacturer. Each background VOI had a diameter of 40 mm. Two Figures of Merit were computed: 1) the coefficient of variation (COV)

Page 9 of 20

COV 

STD 1   

1 2  b    N  1 bI

where , STD and N indicate respectively the average of counts inside the background VOI, the standard deviation and the number of voxels; 2) the hot contrast recovery (CR)

s

b  1 R 1

ip t

CR  100

cr

where R is the real signal-to-background ratio and S (B) is the mean content of the voxels inside the sphere (background) VOI.

Ac ce p

te

d

M

an

us

Clinical Data: Qualitative and Quantitative evaluation Two oncological PET/CT 18F-FDG studies were considered. The first oncological patient (height 169 cm, weight 80 kg) received an injection of 18F-FDG (407 MBq); the tracer uptake time was 60 min and the acquisition protocol consisted of a whole body CT scan (for anatomical localization and attenuation correction) and a 3D whole body emission scan (4 min for each bed containing the liver, 2 min per bed position otherwise, 7 bed positions, 9-slice overlap between adjacent beds). Data were reconstructed with RwPSF, RwPSF-Filt, RwPSF-GTVR (with optimized and  parameters), RwPSF-R (with optimized and  parameters). Common parameters to all algorithms were: image matrix = 256x256, subsets = 18, reconstructed FOV = 60 cm, number of iterations = 10. A qualitative evaluation of the reconstructed images was performed. The second oncological patient (height 183 cm, weight 73 kg) received an injection of 18F-FDG (351.5 MBq) and the acquisition protocol was the same adopted for the first patient. However, to allow a quantitative evaluation of clinical data, spheres of known activity and dimension were added to the patient data, following the method described in Schaefferkoetter et al. [23] In particular, three spheres (diameter 1.5 cm, 1.0 cm, and 0.8 cm) were acquired in air and in spatial positions corresponding to the liver and the spleen of the patient. The activity concentration in the sphere was set in order to generate a S/B ratio of 4.2 in the liver and 3.8 in the spleen. The sphere raw data were corrected for random coincidences, attenuation, decay, then attenuated with the patient attenuation map and lastly added to the patient raw data. Resulting sinograms were reconstructed with the same four algorithms and parameters used for the first patient. A quantitative evaluation of the results was performed by computing COV and CR for each of the three added “lesions”. Results

Validation of the D index In figure 3A a set of representative D curves, obtained for the 10 mm sphere, as a function of the iteration number, for different reconstruction algorithms are shown. As it can be seen, in some conditions a range of maximum values for D (Dmax) can be individuated, indicating the superiority in terms of lesion detectability of a specific iteration range, while in other conditions D has a plateau, indicating a wide range of iterations with similar characteristics.

Page 10 of 20

ip t cr us an

M

Figure 3A. NEMA IEC Body Phantom. Representative D curves, for the 10 mm diameter sphere, as a function of the iteration number for different reconstruction algorithms. In case of the RwPSF-R reconstruction  =0.0005,  =0.005 and = 0.2 were used.

Ac ce p

te

d

The validity of the D value as parameter for the characterization of the lesion detectability can be confirmed by a visual inspection of the corresponding images in figure 3B.

Figure 3B Representative images obtained at different iterations (from left to right: 1, 5, 10, 20, 30, 40, 50) using different types of reconstructions and acquisition statistics [from top to bottom: RwPSF-R (=0.0005, =0.2, 50 Mcounts), RwPSF-R (=0.005, =0.2, 50 Mcounts), RwPSF (50 Mcounts)].

Page 11 of 20

a) Comparison of D with the best image quality within the same set of data over the iterations In table 2, the results of the evaluation over the iterations are reported. The iterations corresponding to Dmax and those associated by the two operators to the best lesion detectability are shown. As it can be seen, images corresponding to Dmax always well overlap to those selected by the operators for the different types reconstructions as well as for the different levels of statistics.

Dmax 6

16 - 50

6-50

6-50

8-50

7

7-10

6-10

14-50

9-50

9-50

6 9 9

5-9 7-10 8-12

7-10 8-10 8-10

6-50

5-50

8 4-50

RwPSF-R =0.0005, =0.2 RwPSF-R =0.005, =0.2 RwPSF-GTVR =0.0025, =0.1 RwPSF-GTVR =0.025, =0.1 PSF PSF Filt RwPSF-R =0.0005, =2 RwPSF-R =0.005, =2 RwPSF-GTVR =0.0025, =1 RwPSF-GTVR =0.025, =1

5-8

6 - 50

5-50

6

6-9

6-8

8-50

8-50

6-50

5 6 5

4-8 6-9 5-8

5-7 6-9 6-9

7-50

5-50

5-50

5-50

6-10

6-10

6

6-9

6-8

4-50

4-50

4-50

3-50

6-50

LC S/B =3.3 : 1 Acq.Stat=50 Mcounts Op. 1 Op. 2

us

cr

5-8

M

LC S/B =3.3 : 1 Acq.Stat=50 Mcounts

ip t

7-12

HC S/B =5 : 1 Acq Stat=25 Mcounts Op.1 Op.2

an

8-12

Ac ce p

Type of Reconstruction

HC S/B =5 : 1 Acq.Stat=25 Mcounts

Dmax 9

d

RwPSF-R =0.0005, =0.2 RwPSF-R =0.005, =0.2 RwPSF-GTVR =0.0025, =0.1 RwPSF-GTVR =0.025, =0.1 PSF PSF Filt RwPSF-R =0.0005, =2 RwPSF-R =0.005, =2 RwPSF-GTVR =0.0025, =1 RwPSF-GTVR =0.025, =1

HC S/B =5:1 Acq.Stat=50 Mcounts Op.1 Op.2

HC S/B =5 : 1 Acq.Stat=50 Mcounts

te

Type of Reconstruction

LC S/B =3.3 : 1 Acq.Stat=25 Mcounts Dmax

LC S/B =3.3 : 1 Acq Stat=25 Mcounts Op. 1 Op. 2

Dmax 7

5-9

5-9

5

4-6

4-6

7- 50

6-50

6-50

6-50

6-50

6-50

7

5-9

5-6

5

4-7

4-6

10-50

10-50

8-50

8-50

7-50

8-50

6 7 9

5-8 7-10 7-10

5-7 6-9 6-10

4 7 5

3-7 4-8 4-6

3-6 6-9 4-6

7-50

5-50

5-50

8-50

6-50

5-50

8

6-10

6-10

5

5-8

4-6

4-50

3-50

6-50

7-50

5-50

4-50

Page 12 of 20

Table 2. Evaluation of lesion detectability over the 50 iterations for two acquisition statistics and for S/B ratio 5 (top table) and 3.3 (bottom table). Iterations corresponding to Dmax and those associated by the two operators to the best lesion detectability are shown.

Ac ce p

te

d

M

an

us

cr

ip t

b) Comparison of D with the best image quality over different images sets at the same iteration In figure 3C few representative couples of images with the corresponding D value, calculated for the smallest sphere of 10 mm and the operators judgement in terms of better lesion detectability are shown. Over the 96 couples of images an agreement of 93 % (operator 1) and 96 % (operator 2) was found between D and the operators’ choice.

Figure 3C Examples of image quality assessment on image couples: comparison between D and operator selection. Phantom studies: Quantitative validation of the proposed prior a) NEMA IQ: First configuration In Fig. 4 CR as a function of COV and D as a function of the iterations are shown for the different reconstructions. The graphs report the results obtained for two spheres with diameter 10 mm (A and B) and 37 mm (C and D) and for two acquisition statistics, 75 Mcounts (A and C) and 25 Mcounts (B and D). As it can be seen, for both spheres RwPSF always reaches the highest CR, but at the expense of the highest noise content. For the 10 mm sphere, the proposed prior results in a higher CR and in a lower noise when compared to RwPSF-Filt, while it shows a higher CR with a slightly higher noise when compared to RwPSFGTVR. For the 37 mm sphere, less affected by the degradation of the spatial resolution, the two regularized algorithms have higher CR than RwPSF-Filt, approaching that of RwPSF. Furthermore a significant lower noise is observed, in particular for RwPSF-GTVR. For what concerns the D values, for the 10 mm sphere RwPSF has the highest D at 75 Mcounts, while a slightly higher value is achieved by RwPSF-GTVR at 25 Mcounts. For the 37 mm sphere RwPSF-GTVR always reaches the highest D value. This can be explained by considering that, in presence of a slight degradation of the spatial resolution, the difference between the two priors in terms of contrast recovery

Page 13 of 20

us

cr

ip t

becomes relatively small; therefore RwPSF-GTVR is favoured by its smoothing effect, which limits COV and thus helps increasing D. The analysis performed on D is confirmed by the qualitative evaluation of the reconstructed images, shown in Fig. 5. On the 75 Mcounts images, the smallest sphere is clearly better detectable with the RwPSF-R reconstruction. RwPSF-GTVR follows, thanks to the absence of noise in the background which helps in revealing the object. In the RwPSF-Filt reconstruction the sphere is again well detectable, but the noise in the background acts as a confounding factor. In RwPSF the sphere is still visible, but some noise fluctuations are similar to the sphere activity. On the 25 Mcounts images the noise in the background is the prevailing factor. Therefore RwPSF-GTVR, which is able to reduce the amount of background noise, is the reconstruction which leads to the best visual detectability for the 10 mm sphere. In all the other reconstructions, particularly for RwPSF, the background noise becomes a confounding factor for the sphere detection. RwPSF-Filt and RwPSF-R have a similar appearance but RwPSF-R is preferable thanks to the slightly lower background noise. The qualitative evaluation of the biggest sphere revealed that RwPSF-GTVR always provides the best visual detectability. Furthermore, RwPSF-GTVR shows an appreciable uniformity within each sphere, both at high and low statistics. RwPSF-R has a comparable sphere uniformity at high statistics, but non uniformities are visible at low statistics. 30

an

60

25

50

30

20

RwPSF

M

CR

A

D

20

40

15

10

RwPSF-Filt RwPSF-R

RwPSF RwPSF-Filt

5

d

10

RwPSF-R RwPSF-GTVR

0 0,00

0,10

te

RwPSF-GTVR

0,20

0,30

0,40

0 0

0,50

2

4

Ac ce p

60

COV

D

CR

12

10

20

0,30

10

RwPSF-GTVR

15

0,20

8

RwPSF-R

30

0,10

12

RwPSF-Filt

25

20

0 0,00

10

RwPSF

40

10

8

30

50

B

6

Iteration

COV

RwPSF RwPSF-Filt

5

RwPSF-R RwPSF-GTVR 0

0,40

0,50

0

2

4

6

Iteration

Page 14 of 20

80

80

70 75

60 70

65

D

40 30

60

RwPSF

RwPSF

20

RwPSF-Filt 55

ip t

C

CR

50

RwPSF-Filt

RwPSF-R

RwPSF-R

10

RwPSF-GTVR

RwPSF-GTVR

0 0,10

0,20

0,30

0,40

0,50

0

COV

2

4

6

cr

50 0,00

8

10

12

8

10

12

Iteration

80

80

us

70

75

60

70

D

65

40

an

D

CR

50

30

60 RwPSF 55

RwPSF

20

RwPSF-Filt

M

RwPSF-R

RwPSF-Filt RwPSF-R

10

RwPSF-GTVR

RwPSF-GTVR 50 0,00

0,10

0,20

0,30

0,40

0,50

0

2

4

6

Iteration

d

COV

0

Ac ce p

te

Figure 4. NEMA IEC phantom, first configuration: CR as a function of COV (left column), D as a function of the iterations (right column). Row A) Sphere of 10 mm and 50 Mcounts scan. Row B) Sphere of 10 mm and 25 Mcounts scan. Row C) Sphere of 37 mm and 50 Mcounts scan. Row D) Sphere of 37 mm and 25 Mcounts scan.

Figure 5. NEMA IEC phantom, first configuration. From left to right: RwPSF, RwPSF-Filt, RwPSF-GTVR with =0.015, =0.2 and RwPSF-R with =0.002, =0.3; top row 75 Mcounts scan, bottom row 25 Mcounts scan. All the images are obtained after 10 iterations. [the images are shown using the same display parameters] b) Second configuration

Page 15 of 20

B

C

cr

A

ip t

In figure 6 and table 3 the images and the values of CR for the NEMA IEC Body Phantom with the second set of spheres are shown. Both images and quantitative parameters confirm the results obtained on the first setup configuration. It’s worth noticing that the two smallest spheres (5 and 6 mm) are not visible in any of the reconstructions here considered.

D

Lesion Diam. 13 mm CR 59.9 48.3 58.2 58.0

Lesion Diam. 10 mm CR 52.5 39.5 45.4 48.9

an

RwPSF RwPSF-Filt RwPSF-GTVR RwPSF-R

Lesion Diam. 16 mm CR 60.5 52.6 58.9 59.5

M

Reconstruction Algorithm

us

Figure 6. NEMA IEC phantom, second configuration (S/B 6.5:1). A (RwPSF – 10 iterations), B (RwPSFFilt – 10 iterations), C (Rw-PSF-GTVR – optimized parameters – 10 iterations) and D (RwPSF-R – optimized parameters – 10 iterations). Lesion Diam. 8 mm CR 48.3 31.3 30.3 38.0

d

Table 3. NEMA IEC phantom, second configuration (S/B = 6.5:1, 50 Mcounts). CR values calculated on the four largest spheres on images obtained using different reconstruction algorithms after 10 iterations.

Ac ce p

te

Clinical Data: Qualitative and Quantitative validation of the proposed prior In figures 7 and 8 some reconstructed images of the oncological patients, obtained with the different algorithms, are presented. In both patients RwPSF produces the best images in terms of spatial resolution and definition of the image, but the very high noise content could limit its clinical applicability. RwPSF-GTVR controls well the noise, but gives the image an unnatural appearance, also with some loss of spatial resolution. RwPSF-R results in lower visual noise than the clinical images (RwPSF-Filt), with a contemporary higher contrast and spatial resolution with respect to RwPSF-GTVR, even for smaller lesions. The quantitative results obtained for the patient in which three spheres were “added” in the liver and spleen regions well confirm, in terms of COV and CR (table 4), the qualitative evaluation.

Page 16 of 20

ip t cr

M

an

us

Figure 7. Qualitative validation on oncologic patient reconstructed using 10 iterations: top: Maximum Intensity Projection (MIP) coronal image; bottom (from left to right), RwPSF, RwPSF-Filt, RwPSF-GTVR (=0.015, =0.2) and RwPSF-R (=0.002, =0.3) [the images are shown using the same display parameters]

C

Ac ce p

B

te

d

A

D

E

Figure 8. Patient study. A) Original image of the patient (with no “lesions”) reconstructed with the clinical protocol. B, C, D and E: combined images showing the three added lesions. B (RwPSF), C (RwPSF-Filt), D (Rw-PSF-GTVR) and E (RwPSF-R) reconstuctions. Reconstruction algorithm

RwPSF RwPSF-Filt Rw-PSF-GTVR RwPSF-R

Liver

COV 0.26 0.14 0.05 0.08

Lesion 1 Diam.15 mm CR 40.1 34.7 35.8 36.5

Lesion 2 Diam.10 mm CR 41.7 32.8 29.1 35.3

Spleen COV 0.24 0.13 0.04 0.07

Lesion3 Diam. 8 mm CR 37,6 25.1 17.1 24.8

Table 4. Patient study: combined images. COV and CR values for the three spheres added to the patient data obtained using different reconstruction algorithms. Discussion Positron Emission Tomography has attracted increasing interest thanks to its capability of detecting metabolic information and performing in-vivo bio-imaging. Unfortunately, compared to the most common anatomical/morphological diagnostic techniques (CT or magnetic resonance), it is also characterized by poor spatial resolution, which reduces the quantitative accuracy and the image quality. To overcome this limitation, the physical effects leading to the spatial resolution loss may be accounted for by a “global” Point Spread Function (PSF), which can be introduced in the reconstruction algorithm to obtain images with higher quality and detail. Nowadays, the most used reconstruction algorithms in the clinical practice are iterative Page 17 of 20

Ac ce p

te

d

M

an

us

cr

ip t

ones, which converge to the final image by subsequent approximations: if the PSF information is introduced in such an algorithm, better spatial resolution is obtained for higher number of iterations. On the other hand, a common effect of iterative reconstruction techniques is the increase of noise as iterations proceed, due to the ill-posed nature of the reconstruction problem. Usually in clinical practice the quality of the images is privileged over their quantitative accuracy and the iterative algorithm is therefore stopped after few iterations and, consequently, far from convergence. The reduction of noise is usually obtained by convolving the final image with a filter kernel: this results in a worsened spatial resolution and underestimation of the activity, particularly in the smallest signal regions. The statistical reconstruction methods, albeit their increased computational requirements and time, allow implementing in-reconstruction denoising strategies with a reduced spatial resolution loss, since the image is repeatedly updated by direct comparison with the acquired data. Post-reconstruction methods, acting on the already reconstructed image, do not benefit from the continuous adjustments by the reconstruction algorithm itself: this limits the applicability of the methods, in particular if more refined image formation models are required. Among the different in-reconstruction solutions (e.g. see [24]), regularization techniques have been demonstrated to be useful for taking noise under control during the reconstruction and improving the benefits from the use of the PSF information by increasing the number of iterations used. Two important topics should be taken into account when dealing with regularization techniques:  since many regularization strategies exist (and each of them leads to different characteristics on the image), it is important to understand which effects are desired, in order to choose the most useful strategy for the scope under study;  since every regularization strategy depend on one or more parameters, it is useful to set some optimization criteria to obtain the best results. As far as the first issue is concerned, in PET it is usually desirable to suppress the noise (to obtain smooth, uniform background regions) while retaining as much spatial resolution as possible, all the more if the PSF information is taken into account. Two good candidates are the Huber (or Gauss-Total Variation) [10, 12] and the generalized p-Gaussian [13] priors. The former provides good preservation of spatial resolution thanks to the TV component for high gradients, but the Gaussian component for low gradients might be insufficient in controlling very noisy environments (as often encountered in PET, in particular when the number of iterations is increased to exploit the PSF action) unless the regularization strength is set to a very high value, obtaining an unnatural reconstructed image. The p-Gaussian prior provides a very strong smoothing on background regions, resulting in good noise control, while it smoothes much less in signal regions. In this work a new prior was proposed, based on a modification of the p-Gaussian prior, to maintain the smoothing effect for low gradients (i.e. in background regions) and to reduce the spatial resolution loss, while retaining “natural” transitions and appearance in the image. A 3D OSEM algorithm has been modified to include the proposed prior using a MAP OSL multiplicative approach. The second issue to be studied deals with the optimization criterion for the regularization parameters to be used. In this work a figure of merit, taking into account both the qualitative and the quantitative content, was proposed to evaluate the global “detectability” of a lesion. The validation of this index showed a very good correlation with the human response and, thus, justified its use to set the regularization parameters. The validation of the proposed prior was performed by comparing it with the results from OSEM with PSF information, from OSEM with PSF information and clinical post-filter and from OSEM with PSF information and GTV prior. The priors needed two regularization parameters: the regularization strength  (which controls the effect of the regularization strategy on the final image) and the signal-background threshold  (which distinguishes the different behaviours to be applied in the different regions of the image). Among them, the latter parameter is, in principle, more difficult and tricky to be set, since different regions in the same study might require different thresholds to maximize the general spatial resolution preservation. One conservative strategy may consist in setting the parameter for the worst conditions expected (e.g. for the smallest lesions and the lowest contrasts), assuming that larger lesions and/or higher contrasts lead to lower resolution loss introduced by the regularization. The regularization parameters were determined by maximizing the detectability index for each of the two priors considered in this work. This optimization was performed for a sphere with diameter 10 mm and 10 OSEM iterations. The maximization results hinted a different behaviour of the two priors: the GTV prior showed a narrow peak along the direction of the  parameter, suggesting quite different behaviours for small changes of this parameter and, consequently, indicating that setting this parameter could be a difficult

Page 18 of 20

an

us

cr

ip t

process. For the proposed prior, small changes of the  parameter around the maximum lead to very little differences in its behaviour, guaranteeing higher stability and “margins of security”. A qualitative and quantitative comparison of the performances of the different algorithms was performed on data acquired with a NEMA IEC Body Phantom and on data relative to two oncological patients. The qualitative evaluation showed that the regularized algorithms allow performing an effective noise control without any loss in spatial resolution. In fact, in all the considered experimental conditions, the two regularized algorithms allowed detecting the same smallest spheres as RwPSF (see NEMA phantom, second experiment). From a quantitative point of view, the RwPSF reconstruction provides the highest CR values but also the highest level of noise. The proposed prior provides instead the best compromise between CR and noise, since CR is higher than RwPSF-Filt and RwPSF-GTVR, while noise is lower than RwPSF-Filt and slightly higher than RwPSF-GTVR. It is important to underline that the priors were optimized considering 10 OSEM iterations as a limit for the image reconstruction in clinically acceptable times. If, in the future and/or for different clinical requirements, the number of iterations could be increased, different regularization parameters and consequent numerical results would be probably obtained and the comparison between the reconstruction algorithms might yield different conclusions. However, this goes beyond the scope of this work. The qualitative comparison on the images of the oncological patients confirmed the quantitative results. The algorithm with the proposed prior is capable of controlling the noise while retaining high enough spatial resolution and contrast, also in small lesions. Conclusions

M

The proposed prior represents a good compromise between noise control and preservation of spatial resolution and quantitative accuracy. The proposed strategy to set the regularization parameters appears to be useful in determining the optimal parameters for a specific clinical protocol. References

te

d

[1] Panin V, Kehren F, Michel C, Casey M. Fully 3-D PET reconstruction with system matrix derived from point source measurements. IEEE Trans Med Imag 2006; 25:907 21

Ac ce p

[2] Sureau F, Reader A, Comtat C, Leroy C, Ribeiro MJ, Buvat I, Trébossen R. Impact of image-space resolution modeling for studies with the High-Resolution Research Tomograph. J Nucl Med 2008; 49:1000 8 [3] Alessio A, Stearns C, Tong S, Ross S, Kohlmyer S, Ganin A, Kinahan P. Application and evaluation of a measured spatially variant system model of PET image reconstruction. IEEE Trans Med Imag 2010; 29:938 49 [4] Rapisarda E, Bettinardi V, Thielemans K, Gilardi MC. Image-based point spread function implementation in a fully 3D OSEM reconstruction algorithm for PET. Phys Med Biol 2010; 55:4131 51 [5] Vogel CR. Computational Methods for Inverse Problems. SIAM; 2002 [6] Thielemans K, Asma E, Ahn S, Manjeshwar RM, Deller T, Ross SG, Stearns CW, Ganin A. Impact of PSF Modeling on the Convergence Rate and Edge Behavior of EM Images in PET. Nucl Sci Symp Conf Rec IEEE 2010; 3267 72 [7] Qi J, Leahy R. Iterative reconstruction techniques in emission computed tomography. Phys Med Biol 2006; 51:R541 78 [8] Aubert G, Vese L. A variational method in image recovery J Numer Anal 1997; 34:1948 79 [9] Scherzer O, Grasmair M, Grossauer H, Haltmeier M, Lenzen F. Variational Methods in Imaging. Volume 167 of Applied Mathematical Sciences. Springer, 2009. [10] Keeling SL. Total Variation based convex filters for medical imaging. Appl Math Comput 2003; Page 19 of 20

139:101 19 [11] Dey N, Blanc-Féraud L, Zimmer C, Roux P, Kam Z, Olivo-Marin J, Zerubia J. 3D Microscopy Deconvolution using Richardson-Lucy Algorithm with Total Variation Regularization. 2004 INRIA Rapport de recherche n. 5272. [12] Huber PJ. Robust statistics. John Wiley and Sons; 1981

ip t

[13] Bouman C, Sauer K. A Generalized Gaussian Image Model for Edge-Preserving MAP Estimation. IEEE Trans Im Proc 1993; 2:296 310

cr

[14] Shepp L, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imag 1982; MI-2:113 22

us

[15] Green PJ. On Use of the EM for Penalized Likelihood Estimation. Journal of the Royal Statistical Society Series B (Methodological) 1990; 52:443 52

an

[16] Hudson H, Larkin R. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imag 1994; 13:601 9 [17] Mustafovic S, Thielemans K. Additive and Multiplicative Versions of the Maximum A Posteriori Algorithm With Median Root Prior. Nucl Sc Symp Conf Rec 2001 IEEE; Vol. 3:1783 5

M

[18] Bettinardi V, Presotto L, Rapisarda E, Picchio M, Gianolli L, Gilardi MC. Physical performance of the new hybrid PET∕CT Discovery-690. Med Phys 2011; 38:5394 411.

d

[19] Asma E, Manjeshwar R. Analysis of organ uniformity in low count density penalized likelihood PET images. Nucl Sci Symp Conf Rec 2007 IEEE; Vol. 6:4426 32

te

[20] Wilson JM, Ross SG, Deller T, Asma E, Manjeshwar R, Turkington TG. A Phantom Study of Regularized Image Reconstruction in PET. Nucl Sci Symp Conf Rec 2010 IEEE, pp. 3661 5

Ac ce p

[21] Weber EH. Der Tastsinn und das Gemeingefuhl. In “Handwörterbuch der Physiologie – vol. 3. Braunschweig”, 1850 [22] Fechner GT. Elemente der Psychophysik – vol. 2. Leipzig: Breitkopf und Haertel. 1860 [23] Schaefferkoetter J, Casey M, Townsend D, El Fakhri G. Clinical impact of time-of-flight and point response modeling in PET reconstructions: a lesion detection study. Phys Med Biol 2013; 58:1465 78 [24] Le Pogam A, Hanzouli H, Hatt M, Cheze Le Rest C, Visvikis D. Denoising of PET images by combining wavelets and curvelets for improved preservation of resolution and quantitation. Med Image Anal. 2013; 17:877 91

Page 20 of 20