David Morland
email: david.morland@reims.unicancer.fr

Sofiane Guendouzen

Edmond Rust

Dimitri Papathanassiou

Nicolas Passat

Fabrice Hubelé

Data-driven respiratory gating for ventilation/perfusion lung scan

Keywords:

BACKGROUND: Ventilation/perfusion lung scan is subject to blur due to respiratory motion whether with planar acquisition or single photon emission computed tomography (SPECT). We propose a data-driven gating method for extracting different respiratory phases from lung scan list-mode or dynamic data. METHODS: The algorithm derives a surrogate respiratory signal from an automatically detected diaphragmatic region of interest. The time activity curve generated is then filtered using a Savitzky-Golay filter. We tested this method on an oscillating phantom in order to evaluate motion blur decrease and on one lung SPECT. RESULTS: Our algorithm reduced motion blur on phantom acquisition: mean full width at half maximum 8.1 pixels on non-gated acquisition versus 5.3 pixels on gated acquisition and 4.1 pixels on reference image. Automated detection of the diaphragmatic region and time-activity curves generation were successful on patient acquisition. CONCLUSIONS: This algorithm is compatible with a clinical use considering its runtime. Further studies will be needed in order to validate this method.

I. INTRODUCTION

Ventilation/perfusion lung scan is subject to blur due to respiratory motion whether with planar acquisition or single photon emission computed tomography (SPECT). Diaphragm motion amplitude can indeed reach 10 cm during respiratory cycle [START_REF] Wang | Relationship between chest wall motion and diaphragmatic excursion in healthy adults in supine position[END_REF]. This phenomenon may have clinical consequences, particularly for lung bases. For instance, Suga et al. [START_REF] Suga | Clinical utility of coregistered respiratory-gated 99mTc-Technegas/MAA SPECT-CT images in the assessment of regional lung functional impairment in patients with lung cancer[END_REF] reported a lower sensitivity for detection of perfusion and ventilation defect when assessing regional functional impairment in patients with lung cancer by comparison with a gated approach. In order to overcome artifacts caused by respiratory motion, various methods have been proposed that aim to split the acquisition into time bins or gates during which motion is considered low [START_REF] Sanders | Fully automated data-driven respiratory signal extraction from SPECT images using Laplacian eigenmaps[END_REF]. Two different techniques have been described in the literature to subdivide the surrogate respiratory signal: phase gating (each respiratory cycle is split in a fixed number of gates) and amplitude gating (gates are only defined by the depth of the breath). Phase gating can result in large error in patients who have irregular respiration [START_REF] Chang | Implementation of an automated respiratory amplitude gating technique for PET/CT: Clinical evaluation[END_REF][START_REF] Chang | Design and performance of a respiratory amplitude gating device for PET/CT imaging[END_REF], amplitude gating can be more efficient in such cases [START_REF] Chang | Implementation of an automated respiratory amplitude gating technique for PET/CT: Clinical evaluation[END_REF][START_REF] Dawood | Respiratory gating in positron emission tomography: A quantitative comparison of different gating schemes[END_REF]. The acquisition of a signal describing the respiratory state over time is critical; it leads most of the time to the purchase of a dedicated sensor system. However such systems have their disadvantages besides the cost: large sensors systems prevent the detectors to come close to the patient and may lower spatial resolution; when placed directly on the patient, they may cause some discomfort. Data-driven approaches based on the variation of center of mass of target lesions have been previously described for Positron Emission Tomography (PET) [START_REF] Büther | List mode-driven cardiac and respiratory gating in PET[END_REF]. However detection of the center of mass on SPECT may be impaired by the sequential projection design of the acquisition [START_REF] Sanders | Fully automated data-driven respiratory signal extraction from SPECT images using Laplacian eigenmaps[END_REF]: some regions may not be visible from each projection angle due to attenuation. Count rate variations tracking methods have also been described on PET on phantom [START_REF] He | Respiratory motion gating based on list-mode data in 3D PET: A simulation study using the dynamic NCAT phantom[END_REF]. Methods based on diaphragm tracking were previously published on cone-beam CT [START_REF] Bögel | Respiratory motion compensation using diaphragm tracking for cone-beam C-Arm CT: A simulation and phantom study[END_REF][START_REF] Zijp | Extraction of the respiratory signal from sequential thorax Cone-Beam X-ray images.14th ICCR[END_REF]. Although SPECT projections have lower spatial resolution that might limit the detection of diaphragm, lung projections have high contrast and little background activity. In this paper we propose a gating algorithm based on the variation of the count rate of an automatically detected diaphragm region of interest.

II. METHODS

A. Input Data

We used a dynamic thorax phantom (Model 008A, Computerized Imaging Reference System, Inc.) with a spherical SPECT/CT insert of 8 ml (2.5 cm of diameter) filled with 2 MBq/ml of [99mTc] Pertechnetate. The phantom was positioned at the center of the field of view: the center of the insert was located 65 mm right from the center. The concentration was determined based on the mean lung capillary blood volume (89 ml) [START_REF] Rouatbi | Vieillissement du poumon profond[END_REF] and on the activity commonly used in patients (185 MBq). During respiratory motion, two portions are distinguishable in the lung: a diaphragmatic area in which the activity varies thoughout the respiratory cycle, and an extradiaphragmatic area in which activity is globally stable. Motion length of the phantom was set to 20 mm, lower than the diameter of the insert in order to create an area in which the radioactive sphere was always present at any given point, reproducing the extra-diaphragmatic portion of the lung. Three settings were tested using an infero-superior motion: 10, 20 and 25 cycles per minute, to encompass the range of human respiratory frequency. One patient referred for a suspicion of pulmonary embolism was retrospectively included and received 81 Krypton (81Kr) gas and an intravenous infusion of 187 MBq of [99mTc] Macroaggregated Albumine ([99mTc]MAA). This study was conducted in accordance to the principles set forth in the Helsinki Declaration.

B. Image Acquisitions

For phantom acquisitions, 30 projections (15 per detector, start angle: 0°, step angle: 12°) were acquired in dynamic mode (200 frames of 300 ms per projection) on a 128x128 matrix using a Symbia T2 system (Siemens Medical Solutions, USA). Patient acquisition was carried out using a General Electric Infinia SPECT/CT system (General Electric, Milwaukee, USA). 30 list-mode projections (15 per detector, start angle: 0°, step angle: 12°) of 50 seconds were acquired on a 128x128 matrix. Ventilation and perfusion events were extracted using appropriate energy windows (81Kr: 190 keV +/-10%, 99mTc: 140 keV +/-10%). This protocol was the one in use for clinical purpose at the time. Both phantom and patient SPECT were reconstructed using an OSEM algorithm (10 subsets, 2 iterations, post-filtering: Butterworth order 1, cut-off frequency 0.30 cycles/cm) on a Xeleris Workstation (General Electrics, Milwaukee, USA).

C. Respiratory Gating

1) Preprocessing

Two sets of images were reconstructed for each SPECT projection: a summed image consisting in the non-gated acquisition, and a second image series binned into 300 ms frames. No preprocessing was required for phantom acquisitions, which were already acquired with 300 ms frames. However for patient acquisition, list-mode data were temporally rebinned into 300 ms frames for each projection angle.

2) Time-activity curve calculation

A region of interest was automatically determined around the diaphragmatic region of the patient or the lower oscillating part of the phantom acquisition on the anterior summed projection of the SPECT. Each line of the image was summed to create a column vector that allowed us to define the upper bound of the region of interest: this upper bound corresponded to the last line whose value exceeds 75% of the maximum of the column vector. Lower, right and left bounds corresponded to the limits of the field of view. The same ROI was then used for each projection angle (Figure 1). For patient acquisition, [99mTc]MAA anterior projection was preferred to 81Kr anterior projection because of less scattering. This method remains valid as long as no extra pulmonary activity is observed in the field of view (e.g. injection site). A time-activity curve was then derived from the total number of counts inside the ROI on the 300ms frames. Geometric mean of both opposite detectors was used. The output signal was smoothed using a Savitzky-Golay filter (span: 7, order: 3). More precisely, the curve was filtered by fitting successive subsets of adjacent 7 data points with a 3-degree polynomial. This choice was motivated by the filter's tendency to preserve the magnitude of the peaks in a noisy signal as previously described [START_REF] Rouatbi | Vieillissement du poumon profond[END_REF][START_REF] Bromba | Application hints for Savitzky-Golay digital smoothing filters[END_REF]. Figure 1 3) Gating Amplitude-based gating was performed to obtain 5 different images per planar acquisition or SPECT projection angle. Minima and maxima of the time-activity curve were detected using the derivative of the curve. First quantile (Q1), third quantile (Q3) and interquantile interval (IQ) were calculated for minimal and maximal value respectively. Values outside the range of [Q1-IQ;Q3+IQ] were rejected. The algorithm then detected any two consecutive maximal (respectively minimal) values, and selected only the highest (respectively the lowest point) for each couple. The interval between maximal and minimal values was then split into 5 amplitude gates used to reconstruct the 5 different images per planar acquisition or projection. Each projection was composed of the same number of 300 ms frames to allow valid reconstruction. The 5 sets were imported in the SPECT/CT workstation and reconstructed.

D. Implementation and Run-Time

The algorithm was implemented inhouse using Matlab (The MathWorks, Inc., Natick, MA, USA). Processing was performed on a laptop with a 2.6 GHz Intel Core i7 processor.

E. Motion blur estimation

1) Blur estimation

SPECT volumes were summed to obtain one summed coronal slice. Each line of the image was then summed to create a column vector. A curve representing the variation of pixel value along the vertical axis was derived from the column vector. Full width at Half Maximum (FWHM) was calculated and compared to the FWHM of a reference image of the same duration acquired while motion was turned off. The blur index (BI) was calculated as:

BI = (FWHM -FWHMreference)/FWHMreference

A blur index of 0% indicated no blur.

2) Estimation of oscillation frequency

For each projection of SPECT acquisitions, oscillation frequency was derived from the Fourier transform of the time activity curve generated. Automatic detection of the most dominant frequency within a range of 5 to 30 cycles per minute was performed with a maximum detection and compared to the programmed input frequency.

III. RESULTS

A. Phantom tests

Data-derived frequency exactly corresponded to the input frequency whether at 10, 15 or 25 cycles per minute (Table 1). FWHM and BI are presented in Table 1. Non gated FWHM was 8.56, 8.58 and 7.21 pixels for 10, 20 and 25 cycles per minute acquisitions, respectively. Mean gated FWHM was reduced to 5.14, 5.90 and 5.02 for 10, 20 and 25 cycles per minute acquisitions, respectively. Reference FWHM was 4.1 pixels for all tested situations. An example of correspondence between input motion and calculated time-activity curve is shown Figure 2. Runtime was evaluated to less than 10 seconds. Smoothed time-activity curve and minima-maxima detection on the first projection of the SPECT is showed (Figure 3). A mean of 1 extremum per projection was removed by interquartile exclusion. Detected dominant frequency was 9 cycles/minute. SPECT acquisitions corresponding to the phase of maximal inspiration and maximal expiration (Gate 1 and 5) were reconstructed from the list-mode SPECT: perfusion coronal slices are presented in Figure 4. Our algorithm reduced motion blur on phantom acquisition. Blur index decreased between nongated and gated acquisition; FWHM was reduced by approximately one third. Data-derived respiratory frequency corresponded to the real motion frequency. Time-activity curves reliably described the input movement for the tested frequencies that encompassed the normal range of physiological adult respiratory frequency. This algorithm is furthermore compatible with a clinical use, considering its runtime. High lung signal to background ratio allows a good detection of diaphragmatic region, even on 300ms images. The threshold of 75% of the maximum line might however be less efficient when marked perfusion defects or hot spots are encountered.

Although no problem with outlier extrema exclusion is noted, the algorithm has not been tested on very irregular respiratory cycles.

Extracted inspiration and expiration perfusion SPECT remains visually suitable for interpretation. Ventilation SPECT are, in our example, of lower quality due to radioaerosol leakage. Joint interpretation of both gated and ungated images might be useful to avoid lower signal to noise ratio due to shorter duration of each extracted SPECT. As [99mTc]MAA and 81Kr projections were acquired simultaneously, time-activity curves were only estimated on [99mTc]MAA and then applied on both [99mTc]MAA and 81Kr projections. Expiration Perfusion SPECT appears more contrasted than inspiration SPECT due to the concentration of the same amount of radiotracer in a smaller volume. On the other hand, lung expansion might allow us to detect small defects that could disappear due to partial volume effect on expiration study. This contrast pattern is inverted for ventilation study, probably because patients continuously breathe 81Kr gas during the acquisition: when lungs expand, more 81Kr is observed. Such inverted pattern may not be met using technetium aerosols, which require a single inhalation prior the acquisition. The question of which gate would be the most suitable for clinical interpretation remains unanswered. When considering the 1/8 data centered at peak inspiration on MAA/Technegas SPECT, Suga et al. [START_REF] Suga | Clinical utility of coregistered respiratory-gated 99mTc-Technegas/MAA SPECT-CT images in the assessment of regional lung functional impairment in patients with lung cancer[END_REF] reported a higher sensitivity for detection of perfusion and ventilation defect. Gated SPECT images detected 10.2% more ventilation defects (205 vs. 186) and 9% more perfusion defects (218 vs. 200) compared with ungated images. Other applications of this algorithm could include the therapy planning prior to hepatic 90Y radio embolization where diaphragm motion may impair the estimated dose to the lungs [START_REF] Garin | Personalized dosimetry with intensification using 90Y-loaded glass microsphere radioembolization induces prolonged overall survival in hepatocellular carcinoma patients with portal vein thrombosis[END_REF]. Although encouraging, the validation of our preliminary results on patients is still limited by the absence of ground truth and need to be tested on other subjects.

V. CONCLUSION

In this article, we described a fully automated data-driven procedure to perform an amplitude respiratory gating in ventilation/perfusion SPECT. Although promising, our preliminary results need to be tested on other subjects.

Table 1: Detected frequency, FWHM and Blur Index between gated and non-gated acquisitions. Figure 1: Example of automated diaphragmatic region of interest drawing. Upper limit corresponds to the last line whose sum value exceeds 75% of the maximum line (on the left : patient acquisition, on the right : phantom acquisition). Figure 2: Input motion (plain line) and calculated smoothed time-activity curve (dots) for 20 cycles per minute acquisition. X-axis corresponds to frame number (1 frame per 300 ms). Y-axis corresponds to amplitude. Figure 3: An example of maxima and minima analysis on the first projection of the lung SPECT. Top curve (A) corresponds to the smoothed time-activity curve (X-axis: frame number, Y-axis: amplitude). (B) corresponds to the maxima and minima detection (positive peak : maximum, minimal peak : minimum) before (top curve) and after (bottom curve) the removal of outlier points. Figure 4: coronal slices of reconstructed maximal expiration perfusion SPECT (A), maximal expiration ventilation SPECT (B), maximal inspiration perfusion SPECT (C) and maximal inspiration ventilation SPECT (D). Outline of lung fields drawn manually: light green outline for maximal inspiration and blue outline for maximal expiration.

Figure 2 B

 2 Figure 2 B. Patient AcquisitionRuntime was evaluated to 360 seconds on a laptop among which less than 10 seconds were used for the gating algorithm. The first 350 seconds corresponded to the import of the file and the generation of the perfusion and ventilation planar frames. Smoothed time-activity curve and minima-maxima detection on the first projection of the SPECT is showed (Figure3). A mean of 1 extremum per projection was removed by interquartile exclusion. Detected dominant frequency was 9 cycles/minute. SPECT acquisitions corresponding to the phase of maximal inspiration and maximal expiration (Gate 1 and 5) were reconstructed from the list-mode SPECT: perfusion coronal slices are presented in Figure4.

Figure 3 Figure 4 IV

 34 Figure 3

Table 1 : Detected frequency, FWHM and Blur Index between gated and non-gated acquisitions.

 1

	Input Frequency (cycles/min)	10	20	25
	Calculated frequency (cycles/min)*	10	20	25
	FWHM (pixels) / Blur index (%)
	-Non gated	8.56 (108.3%)	8.58 (108.8%)	7.21 (75.4%)
	-Gate 1	4.67 (13.6%)	6.04 (47.0%)	4.82 (33.1%)
	-Gate 2	4.74 (15.3%)	6.27 (52.6%)	5.07 (29.7%)
	-Gate 3	4.96 (20.7%)	6.37 (25.8%)	5.24 (27.5%)
	-Gate 4	6.13 (49.1%)	6.18 (50.4%)	5.31 (29.2%)
	-Gate 5	5.19 (26.3%)	4.61 (12.2%)	4.65 (13.1%)

* Mean frequency on all pairs of 15 projections.