Enrique Muñoz

Ane Etxebeste

Denis Dauvergne

Jean Michel Létang

David Sarrut

Voichita Maxim
email: voichita.maxim@creatis.insa-lyon.fr

Etienne Testa

Imaging of polychromatic sources through Compton spectral reconstruction

Objective: Study the performance of a spectral reconstruction method for Compton imaging of polychromatic sources and compare it to standard Compton reconstruction based on the selection of photopeak events. Approach: The proposed spectral and the standard photopeak reconstruction methods are used to reconstruct images from simulated sources emitting simultaneously photons of 140, 245, 364 and 511 keV. Data are simulated with perfect and realistic energy resolutions and including Doppler broadening. We compare photopeak and spectral reconstructed images both qualitatively and quantitatively by means of activity recovery coefficient and spatial resolution. Main results: The presented method allows improving the images of polychromatic sources with respect to standard reconstruction methods. The main reasons for this improvement are the increase of available statistics and the reduction of contamination from higher initial photon energies. The reconstructed images present lower noise, higher activity recovery coefficient and better spatial resolution. The improvements become more sensible as the energy resolution of the detectors decreases. Significance: Compton cameras have been studied for their capability of imaging polychromatic sources, thus allowing simultaneous imaging of multiple radiotracers. In such scenarios, Compton images are conventionally reconstructed for each emission energy independently, selecting only those measured events depositing a total energy within a fixed window around the known emission lines. We propose to employ a spectral image reconstruction method for polychromatic sources, which allows increasing the available statistics by using the information from events with partial energy deposition. The detector energy resolution influences the energy window used to select photopeak events and therefore the level of contamination by higher energies. The spectral method is expected to have a more important impact as the detector resolution worsens. In this paper we focus on energy ranges from nuclear medical imaging and we consider realistic energy resolutions.

Introduction

In nuclear medicine, single photon emission computed tomography (SPECT) is a well established technique for diagnostics imaging. This technique is commonly performed by means of collimated cameras, which detect γ-rays issued from a radiotracer previously injected in the patient's body and distributed in the organs or regions of interest. Each measurement position with a collimated camera allows the acquisition of a planar projection of the tracer distribution, so that acquisitions at different angles can be combined to reconstruct the full three-dimensional image. As an alternative to collimated cameras, Compton cameras have been proposed for this purpose. The main benefit of Compton cameras with respect to collimated cameras is the potential increase in detection efficiency derived from the absence of a physical collimator, which has been shown to be of at least one order of magnitude in the energy range of γ-rays emitted by the most commonly used radiotracers [START_REF] Han | Statistical performance evaluation and comparison of a compton medical imaging system and a collimated anger camera for higher energy photon imaging[END_REF][START_REF] Fontana | Compton camera study for high efficiency SPECT and benchmark with Anger system[END_REF]. Compton cameras can also be applied in a wider energy range, and are better suited for imaging of higher energy photons up to a few MeV, which present higher collimator transparency and would thus require a thicker septa.

The operation principle of Compton cameras is based on the detection in time coincidence of an initial Compton interaction and the subsequent photoelectric absorption of the scattered photon. In such events, the origin of the incident γ-ray is constricted to a conical surface with an aperture angle given by the Compton scattering formula

cos β = 1 - m e c 2 Ẽ1 E 0 E 0 -Ẽ1 , ( 1
)
where β is the scattering angle, m e c 2 = 511 keV, Ẽ1 is the energy transferred in the Compton interaction and E 0 is the incident γ-ray energy.

Recent studies have investigated the suitability of Compton cameras for polychromatic imaging, thereby allowing simultaneous imaging of more than one radiotracer [START_REF] Kishimoto | First demonstration of multi-color 3-d in vivo imaging using ultracompact compton camera[END_REF][START_REF] Sakai | In vivo simultaneous imaging with 99mTc and 18F using a Compton camera[END_REF][START_REF] Nakano | Imaging of 99mtc-dmsa and 18f-fdg in humans using a si/cdte compton camera[END_REF][START_REF] Uenomachi | Simultaneous in vivo imaging with pet and spect tracers using a compton-pet hybrid camera[END_REF]). In such scenario, γ-rays of different energies are emitted from the different tracer distributions, and images for each of them must be reconstructed independently. For this purpose, energy cuts are set in the measured spectra in order to select only those events whose total deposited energy corresponds to the known emission energies of the employed tracers. However, this approach implies that only events with full energy absorption in the second interaction can be properly reconstructed. This means that (a) events with partial energy deposition placed outside the chosen bounds in the spectra will not be used, thus diminishing the available statistics, and that (b) contamination events generated by a higher energy photon, with a measured energy falling within the bounds of a lower energy peak, will be reconstructed with a wrong initial energy, hence a wrong scattering angle, and will add noise to the image.

To overcome these two problems, we propose to use a spectral image reconstruction method for polychromatic imaging of different radiotracers. The most important feature of the spectral method is that it does not assume that the second interaction of the measured event is a photoelectric absorption; instead, it calculates the probability for a photon of all the known emission energies to produce the measured outcome, and computes a possible cone of origin for all energies with non-zero probability. For each event, a row of the system matrix is constructed using four-dimensional (three spatial and one spectral) voxels, so that the reconstruction algorithm is able to find both the spatial distribution and the energy of emission of the measured sources. Spectral image reconstruction methods for Compton imaging devices have been previously developed as a means of reconstructing sources with an unknown emission spectrum or a continuous one [START_REF] Xu | Gamma-ray energy-imaging integrated spectral deconvolution[END_REF][START_REF] Gillam | A Compton imaging algorithm for on-line monitoring in hadron therapy[END_REF][START_REF] Muñoz | A spectral reconstruction algorithm for two-plane Compton cameras[END_REF]. The model employed in this work follows the one presented in [START_REF] Muñoz | A spectral reconstruction algorithm for two-plane Compton cameras[END_REF], with several modifications. In the first place, here we focus on the energy range of nuclear medical imaging. At these lower energies, Doppler broadening has a significant impact on the reconstructed images, so it was specifically included in the reconstruction process. In addition to Doppler broadening, realistic energy resolutions were simulated to model the performance of different detector materials. Another significant change is the reconstruction of only a few emission energies in the measured spectrum, as opposed to the reconstruction in a continuous energy range that was performed previously. This allows saving computing resources and speeding up computation time, and can be done in the context of multiple tracer imaging because the possible initial photon energies correspond to those emitted by the tracers, which are previously known. Finally, this work reports for the first time a direct comparison between images reconstructed with the spectral method and those obtained using a standard reconstruction method, selecting only fully absorbed events. Images are compared in terms of the activity recovery coefficient and spatial resolution using the Fourier ring correlation method.

The paper is organized as follows. Section 2 describes in detail the proposed spectral reconstruction method, as well as the simulations performed to carry out the study and the figures of merit used to quantify the improvement of the proposed method. Section 3 presents the results obtained with the proposed method compared to those obtained with standard reconstruction, including reconstructed images with different emission energies and varying statistics. Section 4 is dedicated to the discussion of the previously presented results, and section 5 is left for the conclusions and final remarks.

Materials and methods

Simulation of the imaging system

All simulations have been performed using the Compton camera module CC-Mod [START_REF] Etxebeste | Ccmod: a gate module for compton camera imaging simulation[END_REF] in GATE [START_REF] Sarrut | Advanced monte carlo simulations of emission tomography imaging systems with gate[END_REF]) version 9.0, employing the standard electromagnetic physics list. The simulated Compton camera corresponds to the geometry designed within the CLaRyS collaboration [START_REF] Fontana | Compton camera study for high efficiency SPECT and benchmark with Anger system[END_REF]). It consists of seven silicon scatterer layers with an area of 9.0×9.0 cm 2 and thickness of 0.2 cm and a 3 cm thick scintillator absorber with a surface of 28×21 cm 2 . There is a 1 cm gap between scatterer planes and the absorber is placed at 15 cm center-to-center from the last scatterer. In the original design, BGO was considered as the material for the scintillator absorber, although it could also be replaced by CeBr 3 to achieve a better energy resolution. In this work, both absorber materials have been considered.

The simulated source distribution is formed by six cylinders with a constant thickness of 4 mm and diameters of 10, 13, 17, 22, 28 and 37 mm, placed at a distance of 100 mm from the first scatterer (see Figure 1(b)). All cylinders emit γ-rays at four energies with the same probability and a constant activity density of 20.5 kBq/mL for each energy. The four simulated γ-ray energies chosen were 140, 245, 364 and 511 keV, which correspond respectively to photon energies emitted by 99m Tc, 111 In and 131 I and to annihilation photons found in the spectra of all β + emitters. No source material was included in the simulations, so the emitted photons cannot suffer from scatter or attenuation outside the detectors. In order to test the reconstruction method in the presence of background activity, an additional imaging test was performed including to the previous sources an homogeneous cylinder emitting at the same 4 energies and with an activity density of 25% of the sources. A final test in three-dimensional space was also carried out, for which the six cylindrical sources were replaced by spheres of the same radius. In this case, each sphere emitted in only one of the four energies, with lowest energies being emitted from the biggest spheres: 140 keV from the 37 mm sphere, 245 keV from the 28 mm one, 364 keV from those with 22 and 17 mm and 511 keV from the two smallest spheres. The activity densities and positions on the plane perpendicular to the detectors were the same as the initial cylinders.

The set of events selected for image reconstruction is formed only by those events with a forward Compton scatter in one of the scatterers and a subsequent interaction in the absorber. Events in which the photon is backscattered in the first interaction or those produced by the coincident detection of two different primary photons were not included. Due to the very asymmetric efficiencies of the simulated scatterer and absorber, the backscattered fraction of events is quite smaller than that of cameras using similar detectors as scatterer and absorber. The data were simulated with ideal energy resolution, including only Doppler broadening. In order to apply an energy resolution to the simulated data, the deposited energy values were convolved with a Gaussian distribution. The applied blurring corresponds to an energy resolution of 1.1 % full width half maximum (FWHM) at 1 MeV for the scatterers. For the absorber, two different resolution values were used: 7 % FWHM at 511 keV and 20 % FWHM at 662 keV. These resolution values correspond to the ones estimated for silicon double strip scatterers and an absorber of CeBr 3 [START_REF] Livingstone | Influence of sub-nanosecond time of flight resolution for online range verification in proton therapy using the line-cone reconstruction in compton imaging[END_REF] and BGO [START_REF] Fontana | Monitoring ion beam therapy with a compton camera: simulation studies of the clinical feasibility[END_REF] scintillator crystals. A minimum energy threshold of 10 keV was set for each detector plane. In all cases, the energy resolution was inversely proportional to the square root of the deposited energy. The spatial resolution of the detectors was not modeled in the simulations, and the interaction position of the measured events corresponds to that of the primary impact. Events employed for standard reconstruction were selected around the known photopeaks with an energy window of ±5 keV for the case of perfect resolution and ±σ of the modelled energy resolution, calculated for each emission energy in the total deposited energy spectrum, whereas all events without any energy restrictions were kept for the spectral method. One of the advantages of the spectral reconstruction method is that events with partial energy depositions can be used, thus increasing the available statistics and allowing for shorter measurement times. The performance of the two reconstruction methods with different measurement times has been studied by reconstructing images with different number of emitted primaries, which is equivalent to having different measurement times with a constant activity. In this study, images were reconstructed with a number of total emitted primary photons between 0.5 and 8×10 7 , evenly distributed amongst all four emission energies. This is equivalent to measurement times between 6.25 and 100 s with the simulated activity, which produced between 6.8×10 3 and 1.1×10 5 coincidences in the simulated system. Twenty independent realizations were reconstructed for each case.

Image reconstruction

Images have been reconstructed using the iterative list-mode Maximum Likelihood Expectation Maximization (MLEM) algorithm

λ ℓ+1 j = λ ℓ j s j i t ij k t ik λ ℓ k , (2)
where λ ℓ j is the image value at voxel j in iteration ℓ, t ij are the system matrix elements indexed over the i measured events and j image voxels and s j is the sensitivity matrix.

The system matrix elements t ij represent the probability that the measured event i is produced by a γ-ray emitted from voxel j. The calculation of this probability takes into account the physical interactions involved in the measurement process. The two compared reconstruction methods differ in the energy modelling of the photon causing the measurement, so their respective system matrices are computed differently. For the standard method, based on the reconstruction of events with full energy deposition, we use the model presented in equation ( 15) in [START_REF] Maxim | Probabilistic models and numerical calculation of system matrix and sensitivity in listmode MLEM 3D reconstruction of Compton camera images[END_REF]:

t ij ∝ | cos(θ --→ r 1 r j )| || --→ r 1 r j || 2 2 K( β|E 0 )p β| β (β| β, E 0 ) , (3)
where r j is the central position of voxel j, r 1 is the Compton interaction position, θ --→ r 1 r j is the angle formed by the vector --→ r 1 r j with the normal to the detector surface, K( β|E 0 ) is the Compton differential cross section from the Klein-Nishina equation and the factor p β| β (β| β, E 0 ) models the uncertainty of the aperture angle of the cone of response, being β the geometrical angle between voxel and cone axis at its apex and β the angle calculated using the measured energy and Compton scattering formula. This uncertainty is modelled as the linear combination of two Gaussian distributions, adjusted for each initial energy from simulated data to account for Doppler broadening [START_REF] Feng | Influence of doppler broadening model accuracy in compton camera list-mode mlem reconstruction[END_REF]). The linear combination of Gaussians was also investigated in [START_REF] Ida | Accurate modeling of event-byevent backprojection for a germanium semiconductor compton camera for system response evaluation in the lm-ml-em image reconstruction method[END_REF], as reflecting influences from scatter on electrons from the different orbits of the detectors atoms. Fitting functions as Lorentzian, Voigt or general parametric function [START_REF] Kim | Resolution recovery reconstruction for a compton camera[END_REF]) could also be used.

For the spectral reconstruction method, all measured events are reconstructed simultaneously with all the considered initial energies, as depicted in Figure 1(a). The expression in equation ( 3) is modified to include the factors relating the deposited energies with the different possible incident photon energies. After inclusion of these factors, as deduced in [START_REF] Muñoz | A spectral reconstruction algorithm for two-plane Compton cameras[END_REF], the calculation of the system matrix elements is given by the expression:

t ij,ε ∝ | cos(θ --→ r 1 r j )| || --→ r 1 r j || 2 2 e -µ 0 λ 0 K( β|E ε ) (E ε -Ẽ1 ) 2 × e -µ 1 λ 1 σ Ẽ2 ( Ẽ2 |E ε -Ẽ1 )p β| β (β| β, E ε ) , (4)
where the index ε runs over all considered initial energies E ε , Ẽk is the energy deposited in interaction k and µ k and λ k are the linear attenuation coefficient and distance inside the detectors traversed by the incident(0) and scattered(1) photons. It is assumed that Ẽk equals the energy transferred in the photon interaction, i.e., that there are no losses due to escaping electrons. The exponential factors account for the attenuation and must be included now because the µ depend on the photon energy. σ Ẽ2 ( Ẽ2 |E ε -Ẽ1 ) is the cross section of the interaction in the absorber, which has two possible values, according to the possible physical interactions in the absorber:

• Photoelectric: if |E ε -( Ẽ1 + Ẽ2 )| < ∆E/2 σ Ẽ2 ( Ẽ2 |E ε -Ẽ1 ) = σ P h (E ε -Ẽ1 ) , (5)
where σ P h (E) is the photoelectric cross section of the absorber material at energy E and ∆E is the width of the energy window around the emission energy for which an interaction in the absorber is treated as a photoelectric absorption.

• Compton: otherwise

σ Ẽ2 ( Ẽ2 |E ε -Ẽ1 ) = K( β2 |E ε -Ẽ1 )2πm e c 2 (E ε -Ẽ1 -Ẽ2 ) 2 ∆E , ( 6
)
where β2 is the scattering angle of a Compton interaction depositing Ẽ2 and with initial energy E 0 -Ẽ1 . In all cases, the sensitivity matrix employed is calculated as described in [START_REF] Muñoz | Study and comparison of different sensitivity models for a two-plane Compton camera[END_REF], including all the energy dependent terms and adapted to the simulated system geometry. The sensitivity is given by

s j,ε ∝ V j d 3 r j Sca d 3 r 1 e -µ 0 λ 0 |⃗ r 1 -⃗ r j | 2 Abs d 3 r 2 e -µ 1 λ 1 |⃗ r 2 -⃗ r 1 | 2 K(β|E ε )σ Ẽ2 ( Ẽ2 |E ε -Ẽ1 ), (7)
which is calculated by Monte Carlo integration by generating 10 4 random points in the scatterer and absorber detectors for each voxel position j and emission energy. Since the system consists of seven scatterer planes, the process is repeated for each of them and the final sensitivity matrix is computed as the sum of the seven sensitivity images.

Images have been reconstructed on a field of view (FoV) of 200×200×4 mm 3 , divided in 50×50×1 voxels. In the spectral method, the FoV contains 50×50×1×4 voxels, as it is extended with an additional dimension representing the possible emission energies. In all cases, images shown in this work correspond to the 20 th iteration of the MLEM algorithm. This number was empirically set after visual inspection of the results and ensures a good compromise between resolution and noise in all configurations.

Image analysis and quantification

Before image reconstruction, the goodness of the data sets obtained by selecting events around the photopeaks is quantified in terms of their precision and recall.Precision of a selected set is defined as the fraction of good events contained in it, i.e., the number of events triggered by a photon with initial energy corresponding to the photopeak energy over the total number of events in the set. Recall is defined as the number of good events selected in a photopeak set over the total number of events triggered by photons of that energy, i.e. the fraction of good events produced by that energy that are actually contained in the selected set around the photopeak. Mathematically, they can be written as:

Precision E i = G E i D E i ; Recall E i = G E i N E i , (8)
where G E i is the number of good events produced by photons emitted and detected at energy E i , D E i is the total number of events detected around E i and N E i is the number of events triggered by photons emitted with energy E i .

The quality of the reconstructed images is quantified by means of the Activity Recovery Coefficient (ARC) of the reconstructed sources and the images spatial resolution computed using the Fourier Ring Correlation (FRC) between reconstructed images and the ground truth (Figure 1(b)).

For each cylindrical source i, the ARC is defined as the ratio between the activity reconstructed in the simulated source position and its true activity:

ARC i = A i j V j V i A T , (9)
where A i and V i are the activity reconstructed in each cylindrical source i and its volume, and A T is the total reconstructed activity in the whole image, either inside or outside the true sources volumes. FRC is used to compute the correlation between a reconstructed image f (r) and the ground truth image g(r) of the simulated distribution at different spatial frequencies [START_REF] Banterle | Fourier ring correlation as a resolution criterion for super-resolution microscopy[END_REF]. It is calculated according to equation

F RC f,g (R) = i∈R f (r i ) • ĝ(r i ) * i∈R | f (r i )| 2 • i∈R |ĝ(r i )| 2 , ( 10
)
where f is the Fourier transform of f , f * denotes its conjugate and R is the ring in the Fourier domain with which the correlation is computed. The image resolution in the spatial domain can be calculated as the spatial equivalent of the intersection between the FRC and a threshold. As in [START_REF] Banterle | Fourier ring correlation as a resolution criterion for super-resolution microscopy[END_REF], the threshold was defined by the 2σ curve:

F 2σ (R) = 2 N p (R)/2 (11)
where N p is the number of pixels in ring R.

Results

Spectrum of simulated sources

The coincidence spectrum of the simulated sources is shown in Figure 2, as detected with realistic energy resolution values by a system based on silicon scatterers and a CeBr 3 (2(a)) or a BGO (2(b)) absorber. The different colors in the plotted spectra depict the contribution by the photons of the four emitted energies. Fully absorbed events are placed in the photopeak around the emission energy, whereas those with partial energy deposition form a continuum at lower energies. It can be appreciated from these spectra that events selected around the photopeaks are contaminated by the contribution of partial depositions from higher energies, especially in the case with worst energy resolution. The obtained precision and recall values are listed in Table 1 for the data sets around each energy with the three energy resolution models considered. Table 1. Precision and recall of events selected around photopeak for the different emission energies and energy resolution models. Precision is defined as the fraction of events produced by photons of the targeted energy contained in a selected set and recall as the number of good events in a photopeak set over the total number of events triggered by photons of that energy.

Energy

Reconstructed images at different energies

One realization of the images reconstructed for the four emission energies with the two reconstruction methods are shown in Figure 3. Figures 3(a), 3(b) and 3(c) show the images corresponding to a system with perfect energy resolution and with the considered energy resolutions of a system with silicon scatterer and CeBr 3 or BGO absorbers, respectively. In all cases, the six cylindrical sources can be identified, but the reconstructed images present different levels of noise. The noise component in the reconstructed images appears as sparse pixels with relatively high activity outside the true simulated sources. The relative activity reconstructed inside and outside the true emission volumes for each emission energy is plotted in Figure 4. Pixels placed inside and outside the true source volumes are selected based on the known geometry of the simulations, and those placed at the source boundaries are considered to be inside the source volume for this computation.

Reconstructed images with different number of primaries

Figure 5 shows one realization of the images reconstructed with different number of emitted primaries, summed over the four obtained images for the different emission energies. It can be appreciated that the increase of available statistics for the spectral method produces a more important enhancement in the images with fewer primaries, in terms of noise reduction, resolution enhancement and detectability of the smallest sources at the lowest energies. This is even more visible in the case of lower detector energy resolutions, as will be shown in sections 3.4 and 3.5 with the ARC and FRC metrics.

Activity recovery coefficient

The enhancement produced by the spectral method is also reflected in the ARC values calculated for the six cylindrical sources, plotted in Figures 6 and7. The ARC values were calculated from the images reconstructed for each emission energy (Figure 6) and from the combined image obtained by summing over the images of the four energies (Figure 7). The plotted values correspond to the mean ARC computed over 20 reconstructed images. In general, a higher ARC is obtained with the spectral method, which indicates a better concentration of the reconstructed activity in the true source volumes. As noted previously, the improvement is bigger for lower statistics and worse detector energy resolution.

The ratio of ARC obtained with the spectral and standard reconstruction methods is plotted in Figure 8 for the different source diameters. It can be seen that the gain is higher for lower statistics, with a tendency to stabilize at a constant value with increasing number of photons, and that the ratio is larger for smaller source diameters. The ratios obtained with highest statistics are listed in Table 2. This is the case that favors normal photopeak reconstruction the most, so that it can be interpreted as the expected minimum enhancement by the spectral method. In the case of perfect energy resolution, the ratio increases from 1.03 for the 37 mm cylinder to 1.11 for the 10 mm cylinder in the reconstructions with highest statistics. These values change to 1.03 and 1.21 using the considered energy resolutions of silicon scatterers and a CeBr 3 absorber, and to 1.23 and 1.39 considering a BGO absorber.

Fourier Ring Correlation

The FRC between the ground truth image and the reconstructed images was calculated according to equation ( 10), using the reconstructed images summed over the four emission energies. The mean results for the 20 realizations for each considered resolution and number of primaries are shown in Figure 9, where the 2σ curve (defined in equation ( 11)) is also plotted. The image resolution values are obtained by finding Table 2. Ratio of the ARC obtained with the spectral reconstruction over the one obtained with the standard method for the different cylinder sizes with 8 × 10 7 emitted primary photons.

Energy

Source diameter (mm) resolution 10 13 17 22 28 37 Perfect 1.11±0.01 1.08±0.01 1.072±0.005 1.048±0.004 1.032±0.003 1.029±0.001 Si+CeBr 3 1.21±0.02 1.16±0.01 1.13±0.01 1.086±0.004 1.053±0.004 1.028±0.002 Si+BGO 1.39±0.02 1.36±0.01 1.33±0.01 1.293±0.004 1.254±0.004 1.233±0.003 the first intersection position of the FRC and the 2σ curve, computed using linear interpolation between FRC points, and translating it into the spatial domain. The calculated resolution values are plotted in figure 10, showing a better image resolution when the spectral method is employed.

Sources in homogeneous background

Images obtained from the simulations including an homogeneous background cylinder with an activity density of 25% of the sources are shown in Figure 11. Shown images were reconstructed with data equivalent to 100 simulated seconds, which corresponds to the highest statistics in previously shown results. The two smallest sources are better reconstructed by the spectral method. More noise can be seen in the zero-activity background with the photopeak selection, which might be due to a lower signal-to-noise ratio in the data. Images obtained with the spectral reconstruction method present higher uniformity and smaller variance. Table 3 lists the variance over the mean, computed over a circular region of 32 mm diameter, in the center of the background cylinder and in the 37 mm source. The ghost activity outside the source is more visible on photopeak images and leads to a lower contrast.

Tests in 3D space

Three-dimensional images obtained from the simulations with spherical sources are shown in Figure 12. In this case, each sphere emitted only at one of the four considered energies, and results show that both methods are able to recover the activity distributions corresponding to each of them. In both reconstruction methods, the shape of the sources is well defined in the xy image plane, parallel to the detector surface, and presents an elongation towards the center of the detector in the xz plane due to the acquisition geometry and data truncation. The contours of the smallest sources are more precise for the spectral method.

Discussion

The proposed spectral reconstruction method aims at improving the performance of simultaneous Compton imaging of polychromatic sources by dealing with two sources of image degradation. The first of them is contamination in the photopeak from higher energy photons that are not completely absorbed. By generating a cone for each possible emission energy, the cone constructed with the true initial photon energy intersects its true emission position and the otherwise noisy event is correctly reconstructed with the spectral method. As seen from the precision values in Table 1, the ratio of such contamination events accepted in the standard reconstruction set depends on the overall system energy resolution, which dictates the width of the energy window of accepted events around the measured peaks. Consequently, the benefits of employing the spectral reconstruction method become more important when real resolution values are taken into account. This can also be appreciated from the obtained images, those reconstructed with the spectral method presenting in general a lower noise component.

The second feature that is improved by the spectral method is the available statistics for image reconstruction, since events with partial energy deposition are reconstructed. This is not important for low energy photons, which have a high probability of being completely absorbed in the detector. However, as the photon energy increases, the probability of being fully absorbed is reduced, resulting in a severe loss of statistics, as observed in the recall values in Table 1. In a potential clinical application, the increase of the available statistics could result in lower doses administered to the patient or reduced acquisition times.

Reconstructed images were quantified in terms of ARC and FRC. Regarding the ARC, higher values were obtained when the spectral method was employed. While this was true in all cases, the increase was more significant for worse detector energy resolutions and lower statistics. This can be seen in the ratio of ARC calculated for all cases tested, plotted in Figure 8. The ratio was consistently above one, which indicates a gain in ARC with the spectral method. Considering the results obtained with the FRC (Figure 9), it can be observed that the 2σ curve crosses the FRC at higher spatial frequencies when the spectral method is used, which indicates a higher correlation to the ground truth image. As expected, this intersection is also displaced to higher spatial frequencies as the number of primaries is increased and for detector models with better energy resolution, for which reconstructed images become more similar to the ground truth. In all cases, the spatial resolution computed from the FRC was better when the spectral method was employed. Figure 10 suggests that the spectral method could allow significant reduction of the radioactive dose to the patient, since roughly the same resolution can be obtained with the spectral method and half the dose required by the traditional method.

Noise in smooth regions can be easily reduced by low-frequency filtering or by regularization. Homogeneity in the constant regions can be improved and variance reduced by low-frequency filtering or by regularization. In the experiments shown here total variation regularization could be used in this sense as in [START_REF] Feng | Influence of doppler broadening model accuracy in compton camera list-mode mlem reconstruction[END_REF]). However, these methods would have little impact on activity recovery and contrast and should not change qualitatively the obtained results. Resolution can be improved by deconvolving for the point spread function. Resolution enhancement aims to reinforce high frequencies in the image. They are easier to recover when the point spread function gets sharper, that is when the initial resolution of the images improves. The threedimensional point spread function, that can be roughly assimilated to the reconstructed image of the smaller source in Figure 12, seems sharper for the spectral method. Noise has a deleterious effect on deconvolution as it is mainly represented in high frequencies and tends to be reinforced. From Figure 11, it can be seen that the noise component in the zero-activity background is stronger for the photopeak reconstruction, especially when realistic energy resolutions are considered. We deduce that the spectral method provides initial estimates that could benefit more from advanced image processing.

The obtained results demonstrate the potential of the spectral method to ameliorate Compton imaging of polychromatic sources. However, there are still some limitations of the method that should be considered for its experimental application. There are known effects that degrade the image and were not taken into account in the present work. One of them is the spatial resolution of the detectors, which is usually modeled using Gaussian probabilities around the measured position [START_REF] Xu | Gamma-ray energy-imaging integrated spectral deconvolution[END_REF][START_REF] Muñoz | Study and comparison of different sensitivity models for a two-plane Compton camera[END_REF]) and can lead to additional blurring if it is not corrected. Backscattered events also represent an important fraction of the measured data, especially for low energy photons. Another source of signal degradation in Compton cameras is the fraction of backscattered events. In the simulated configuration, the ratio of backscattered events was between 0.42% and 0.72% for the four tested emission energies. Possible strategies to deal with these events include using event selection techniques to discard suspected backscattered events or to model all compatible events with an additional cone of response containing the backscatter trajectory. In the field of medical imaging, additional considerations arise, such as the need for three-dimensional imaging and photon attenuation and scatter in the body of patients. As shown in Figure 12, threedimensional Compton images reconstructed with a single camera present elongation artifacts in the direction of the camera due to data truncation, which can be solved by using several cameras simultaneously or by performing measurements in different angles [START_REF] Kishimoto | First demonstration of multi-color 3-d in vivo imaging using ultracompact compton camera[END_REF]. In clinical practice, a substantial component of the measured signal can be due to scattering inside the body, especially for low energy photons. Unlike events with partial energy deposition in the detectors, the scatteredphoton noise cannot be recovered by the spectral method as presented in this work. In the context of high scattering noise, it might be more adequate to set energy cuts to remove events with energy depositions outside the known emission peaks [START_REF] Sakai | Crosstalk reduction using a dual energy window scatter correction in compton imaging[END_REF]. In that approach, the spectral method would still be useful to recover unscattered events from higher energies with energy depositions falling in the limits of a lower peak. While all these effects are present in Compton imaging regardless of the employed reconstruction method, their impact on the proposed method should be evaluated in future works.

In this work only relatively low emission energies were simulated to study the proposed method with currently used radiotracers in SPECT and PET imaging. Nevertheless, there is a growing interest in hybrid imaging systems combining Compton reconstruction with other imaging modalities. These modalities aim at imaging sources emitting γ-rays of at least two different energies of interest, so the proposed method could also be advantageous. In addition, new radiotracers with higher energy photons are proposed in Compton-PET imaging; for instance, in [START_REF] Tashima | 3d compton image reconstruction method for whole gamma imaging[END_REF] , 909 keV photons issued from 89 Zr are imaged. In such cases, the spectral method is expected to produce a more substantial enhancement due to an effective use of events with partial energy deposition, which would highly benefit more compact systems.

Conclusion

The method presented allows performing Compton imaging of polychromatic sources and improving the images with respect to standard reconstruction methods. The main reasons for this improvement are the increase of available statistics and reduction of noise from incorrectly assigned initial photon energies by modelling the probability of events with partial energy depositions from all known incident energies in the system matrix. As a result, the reconstructed images present lower noise, higher activity recovery coefficient and better spatial resolution. The improvements become more sensible as the energy resolution of the detectors decreases and when the dose is reduced.

Figure 1 .

 1 Figure 1. (a) Representation of the set of cones of response generated in the spectral method for a single measurement with four different possible initial energies. Qualitatively, the initial energy increases from blue to red for fixed values of Ẽ1 and Ẽ2 . (b) Reference image of the simulated sources: six cylindrical sources with diameters of 10, 13, 17, 22, 28 and 37 mm. The color scale represents the simulated relative activity. Pixels with lower relative activity are only partially filled by the source volume.

Figure 2 .

 2 Figure 2. Spectrum of the simulated polychromatic sources with energy resolutions equivalent to Si+CeBr 3 (a) and Si+BGO (b). Note that histograms from lower energies are plotted stacked on top of those from higher energies.

Figure 3 .

 3 Figure 3. Reconstructed images at the 4 different emission energies with the spectral reconstruction method and with events selected around the photopeaks for the three energy resolution models: perfect resolution(a), Si+CeBr 3 (b) and Si+BGO(c). Images are normalized independently and shown in logarithmic color scale to highlight the noise component, with a minimum threshold of 1% of the image maximum.

Figure 4 .

 4 Figure 4. Normalized reconstructed activity at the different energies with both reconstruction methods. The lower bars with a darker color indicate the activity reconstructed in the voxels containing the true source volume. Error bars correspond to the standard deviation from 20 realizations.

Figure 5 .

 5 Figure 5. Images summed over the four energies with different levels of statistics, with the spectral reconstruction method and with events selected around the photopeaks, for the three considered energy resolutions: perfect resolution(a), Si+CeBr 3 (b) and Si+BGO(c). Titles on top of each column indicate the number of emitted photons. Images are normalized independently and shown in logarithmic color scale to highlight the noise component, with a minimum threshold of 1% of the image maximum.

Figure 6 .

 6 Figure6. Activity recovery coefficients calculated from images at each energy with statistics corresponding to 8 × 10 7 emitted primaries, using both reconstruction methods and different energy resolution models.

Figure 7 .

 7 Figure 7. Activity recovery coefficients calculated from the sum images with both reconstruction methods with different statistics and energy resolution models. Titles on top of each column indicate the number of emitted photons.

Figure 8 .

 8 Figure8. Ratio of ARC calculated with the spectral method and with events selected around the photopeaks.

Figure 9 .

 9 Figure9. FRC between ground truth and reconstructed images calculated with the spectral method and with events selected around the photopeaks. The continuous green plot depicts the 2σ curve used as threshold to calculate image resolution.

Figure 10 .

 10 Figure 10. Image resolutions calculated from the intersection of FRC and the 2σ curve. Points and error bars depict the mean resolutions and the standard error of the mean calculated from twenty different realizations with equivalent datasets.

Figure 11 .

 11 Figure 11. Reconstructed images of the sources in homogeneous background at the four different emission energies with the spectral reconstruction method and with events selected around the photopeaks for the three energy resolution models: perfect resolution(a), Si+CeBr 3 (b) and Si+BGO(c).

Figure 12 .

 12 Figure 12. Reconstructed three-dimensional images of spherical sources with perfect energy resolution, with views parallel to the xy (a) and xz (b) planes. Slices shown in (a) were extracted from images at a single energy, while images at all energies were summed to obtain the perpendicular slices shown in (b). The statistics used for reconstruction correspond to a simulation time of 50 s.

Table 3 .

 3 Variance over mean in background region (top) and in the 37 mm cylinder (bottom) with both reconstruction methods.

	Energy		Photopeak			Spectral
	resolution 140 keV 245 keV 364 keV 511 keV 140 keV 245 keV 364 keV 511 keV
	Perfect	0.18	0.17	0.22	0.26	0.16	0.11	0.12	0.17
	Si+CeBr	0.23	0.22	0.31	0.37	0.16	0.10	0.12	0.16
	Si+BGO	0.23	0.20	0.25	0.36	0.17	0.11	0.12	0.16
	Perfect	0.15	0.09	0.09	0.11	0.16	0.09	0.08	0.09
	Si+CeBr	0.18	0.12	0.12	0.13	0.17	0.09	0.08	0.09
	Si+BGO	0.18	0.13	0.13	0.13	0.17	0.11	0.09	0.10

Acknowledgements

This work was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007).


Une définition générale de l'influence entre processus stochastiques

Daniel Commenges between stochastic processes in a large class of semi-martingales that we call D . The definition is based on measurability condition on the characteristics of the decomposition of the involved semi-martingales. We also give a definition related to the same concept based on certain likelihood processes, using the Girsanov theorem. Under certain conditions, the 1 two definitions coincide on D . We present the graphical model associated to this notion of independence (or influence) in a set of stochastic processes. These results may be used in causal models. With these tools, influences of one component of a stochastic process on another can be described without ambiguity. From WCLI we can construct a concept of strong local conditional independence (SCLI). When WCLI does not hold, there is a direct influence while when SCLI does not hold there is direct or indirect influence.

We investigate whether WCLI and SCLI can be defined via conventional independence conditions and find that this is the case for the latter but not for the former.

Context

The issue of causality has attracted more and more interest from statisticians in recent years. One approach is directly based on dynamical models and has been developed, starting with [START_REF] Granger | Investigating causal relations by econometric models and crossspectral methods[END_REF] and [START_REF] Schweder | Composable Markov processes[END_REF], and more recently developed using the formalism of stochastic processes, by [START_REF] Aalen | Dynamic modelling and causality[END_REF], [START_REF] Florens | Noncausality in continuous time[END_REF] and [START_REF] Didelez | Graphical models for composable finite Markov processes[END_REF][START_REF] Didelez | Graphical models for marked point processes[END_REF].

Recently we have given more development to the dynamical models approach (Commenges and Gégout-Petit, 2009) using the basic idea of the Doob-Meyer decomposition proposed in [START_REF] Aalen | Dynamic modelling and causality[END_REF]. We have proposed a definition of weak local independence between processes (WCLI) for a certain class of special semi-martingales (called class D)

which involves the compensator of the Doob-Meyer decomposition of the studied semimartingale. Although it can be used in discrete time, this definition is especially adapted to continuous-time processes for which definitions based on conventional conditional independence may fail. The aim of this paper is to give an even more general definition of WCLI, and conversely of direct influence. What we call direct influence of one component

X j on another component X k of a multivariate stochastic process X (noted X j -→ X X k )
is that X k is not WCLI of X j (we use WCLI both as the name of the condition and as an adjective, that is the "I" may mean "independence" or "independent" according to the P 0 -characteristics (B k 0 , C k 0 , ν k 0 ) of the semi-martingale X k are deterministic. We say that X k is LWCLI of X j in X on [0, t] if and only if the likelihood ratio process

Z P/P 0 t = L P/P 0 Ft is (F -jt )-measurable on [0, t]. We have denoted F -jt = H ∨ X -jt and X -jt = ∨ l =j X -lt . 2. X k is LWCLI of X j in X on [r, s] if and only if the process Z P/P 0 t Z P/P 0 r is (F -jt )-predictable
for all t ∈ [r, s] for all the probabilities P 0 as above.

These definitions allow us to define a strong local conditional independence (SCLI) and to use graphical model. It describes a mathematical framework which is well suited for formalizing causality. However we prefer to speak of influence rather than causal influence. A causal interpretation needs an epistemological act to link the mathematical model to a physical reality.

  1,2 , & Anne Gégout-Petit2,3

	(1) INSERM, U 897, Equipe de biostatistique
	(2) Université Victor Segalen Bordeaux 2
	(3) IMB UMR 5251, INRIA Bordeaux Sud-Ouest, Projet CQFD
	présentons le modèle graphique associé à cette notion d'indépendance (ou d'influence)
	dans un ensemble de processus stochastiques. Ces résultats peuvent être utilisés dans
	des modèles de causalité. Nous pensons que cette définition permet de décrire l'influence
	d'une composante d'un processus stochastique sur une autre sans ambiguïté. De WCLI
	on peut construire un concept d'indépendance conditionnelle locale forte (SCLI). Lorsque

Résumé: Nous présentons la notion de indépendance locale faible conditionnelle (WCLI) entre processus stochastiques dans une large classe de semi-martingales, que nous appelons D . La définition est fondée sur la mesurabilité des caractéristiques de la décomposition de la semi-martingale cible. Nous donnons également une définition liée à la notion de vraisemblance sur la base de certains processus, en utilisant le théorème de Girsanov. Sous certaines conditions, les deux définitions coïncident sur D . Nous WCLI n'est pas vérifié, il y a une influence directe ; lorsque SCLI n'est pas vérifié, il y a une influence directe ou indirecte. Nous étudions le lien entre WCLI et les conditions classiques d'indépendance. Abstract: We present the notion of weak local conditional independence (WCLI)

context). This concept of influence is a good starting point for defining causal influence.

In the perspective of extending WCLI to a larger class of processes, we see two ways.

The first one is to stay in the class of semi-martingales and try to be more general about the conditions. In particular we could use the triplet of the characteristics of a semi-martingales. For an exact definition of the characteristics of a semi-martingale, see [START_REF] Jacod | Limit Theorems for Stochastic Processes[END_REF]. Roughly speaking, the characteristics of a semi-martingale are represented by the triplet (B, C, ν) where ν is the compensator of the jump part of the semi-martingale, B the finite variation part not included in ν, and C is the angle bracket process of the continuous martingale. The second way is to work with the likelihood of the process which is also tightly linked with the characteristics of the semi-martingale. In this paper we explore these two ways, extending the WCLI definition to a very large class of processes that we call D , and showing that another definition of WCLI is possible by the use of likelihood processes. Another issue that we explore is the link between WCLI and analogous definitions based on conventional conditional independence; this angle of attack is closer to [START_REF] Granger | Investigating causal relations by econometric models and crossspectral methods[END_REF] proposal for time series.

Definitions and Results

In order to give the general definition of WCLI, using the characteristics of the semimartingales, we present the notation Consider a filtered space (Ω, F, (F t ), P ) and a multivariate stochastic process X = (X t ) t≥0 ; X t takes values in m , and the whole process X takes values in D( m ), the Skorohod space of all cadlag functions: + → m . We suppose that all the filtrations satisfy "the usual conditions". We have X = (X j , j = 1, . . . , m) where X j = (X jt ) t≥0 . We denote by X t the history of X up to time t, that is X t is the σ-field σ(X u , 0 ≤ u ≤ t), and by (X t ) = (X t ) t≥0 the families of these histories, that is the filtration generated by X.

Similarly we shall denote by X jt and (X jt ) the histories and the filtration associated to X j .

Let F t = H ∨ X t ; H may contain information known at t = 0, in addition to the initial value of X. We shall consider the class of special semi-martingales in the filtration (F t ).

We denote by (B, C, ν) the characteristics of the semi-martingale X under probability P , by M j the martingale part of X j , and by M c j the continuous part of this martingale. We denote by (B k , C k , ν k ) the characteristics of the semi-martingale X k under probability P .

We shall assume two conditions on X: A1 M j and M k are square integrable orthogonal martingales, for all j = k.

Under assumption A1, the jumping parts of the martingales M j and M k are orthogonal. Moreover, the characteristic C of X (the angle bracket of the continuous part of the martingale) is a diagonal matrix. Indeed by definition of orthogonality of semi-

We call D the class of all special semi-martingales satisfying A1 and A2'. In fact, A1 and A2' coud be merged into a single compact assuption: the characteristic C of X is a deterministic diagonal matrix. D is stable by change of absolutely continuous probability (C does not change with the probability). D is a very large class of processes: it includes random measures, marked point processes, diffusions and diffusions with jumps. Let X = (X j , j = 1, . . . , m) be in the class D .

1. Suppose the existence of a probability P 0 such that (i) P P 0 , (ii) the characteristics of the semi-martingales X i 's with i = k are the same under P and P 0 and (iii) the