Sign determination methods for the respiratory signal in data-driven PET gating

Patient respiratory motion during PET image acquisition leads to blurring in the reconstructed images and may cause significant artifacts, resulting in decreased lesion detectability, inaccurate standard uptake value calculation and incorrect treatment planning in radiation therapy. To reduce these effects data can be regrouped into (nearly) ‘motion-free’ gates prior to reconstruction by selecting the events with respect to the breathing phase. This gating procedure therefore needs a respiratory signal: on current scanners it is obtained from an external device, whereas with data driven (DD) methods it can be directly obtained from the raw PET data. DD methods thus eliminate the use of external equipment, which is often expensive, needs prior setup and can cause patient discomfort, and they could also potentially provide increased fidelity to the internal movement. DD methods have been recently applied on PET data showing promising results. However, many methods provide signals whose direction with respect to the physical motion is uncertain (i.e. their sign is arbitrary), therefore a maximum in the signal could refer either to the end-inspiration or end-expiration phase, possibly causing inaccurate motion correction. In this work we propose two novel methods, CorrWeights and CorrSino, to detect the correct direction of the motion represented by the DD signal, that is obtained by applying principal component analysis (PCA) on the acquired data. They only require the PET raw data, and they rely on the assumption that one of the major causes of change in the acquired data related to the chest is respiratory motion in the axial direction, that generates a cranio-caudal motion of the internal organs. We also implemented two versions of a published registration-based method, that require image reconstruction. The methods were first applied on XCAT simulations, and later evaluated on cancer patient datasets monitored by the Varian Real-time Position ManagementTM (RPM) device, selecting the lower chest bed positions. For each patient different time intervals were evaluated ranging from 50 to 300 s in duration. The novel methods proved to be generally more accurate than the registration-based ones in detecting the correct sign of the respiratory signal, and their failure rates are lower than 3% when the DD signal is highly correlated with the RPM. They also have the advantage of faster computation time, avoiding reconstruction. Moreover, CorrWeights is not specifically related to PCA and considering its simple implementation, it could easily be applied together with any DD method in clinical practice.


Introduction
For a diagnostic PET/CT exam a typical acquisition duration for the PET scan is 2 min per bed position, therefore making the patient respiratory motion unavoidable and the 'breath-hold' technique often applied to CT acquisitions (Fin et al 2008) unfeasible. The acquired data are therefore averaged over several breathing cycles and, if no correction is applied, respiratory motion will result in a blurring of the reconstructed images and in an overall degradation of the contrast (Nehmeh et al 2002b, Liu et al 2009. Moreover, in PET/CT respiratory motion may lead to mismatch between CT and PET datasets, possibly generating artifacts in the reconstructed PET images due to errors in the attenuation correction (Erdi et al 2004). There are therefore several reasons why accounting for respiratory motion when reconstructing PET data can improve PET image quality In order to reduce the artifacts induced by the patient breathing, PET data can be acquired in discrete bins in correlation with the respiratory cycle, therefore minimizing the motion within a single bin. As a result, multiple PET images are reconstructed, each related to a section of the respiratory cycle, and motion artefacts and image blurring are expected to decrease compared to the non-gated image. Studies have shown, both in phantom experiments and with lung cancer patient data, that the application of respiratory gating reduces the degrading effect of breathing motion in PET images, producing a reduction in the lesion volume up to 34% and an increase in SUV max up to 156.16% compared to the non-gated images (Nehmeh et al 2002a(Nehmeh et al , 2002b. In clinical practice respiratory gating techniques rely on external devices to measure the respiratory motion and several solutions have been implemented on commercial imaging systems (Visvikis et al 2006, Rahmim et al 2007, Pepin et al 2014. Common utilized devices are the Anzai AZ-733V system (Anzai Medical Corp., Tokyo, Japan), which is a pressure belt that measures the expansion of the patient's chest (Riedel et al 2006), and the Real-time Position Management (RPM) system (Varian Medical Systems, Palo Alto, California, USA), which tracks with a video camera the motion of an infrared marker placed on the patient's abdomen (Nehmeh et al 2002b, Otani et al 2010. Alternatively other methods have been proposed that measure the temperature (Boucher et al 2004) or flow of the air (Guivarc'h et al 2004) during patient respiration, via a thermistor and a spirometer respectively, then relating those changes to respiratory motion.
Considering the extra costs related to external tracking devices, the discomfort patients could experience because of them and also because of some evidence that the internal organ motion might differ from patient external surface motion (Ozhasoglu andMurphy 2002, Fayad et al 2011), data-driven (DD) post-acquisition approaches have been suggested, see Kesner et al (2014) for a review. In Nehmeh et al (2003) an 18 F-FDG point source is set on the patient's abdomen and its position used to track respiratory motion through the consecutive dynamic frames; in Visvikis et al (2003) time activity curves are obtained from ROIs selected on the series of dynamic images and a Fourier transform in time is performed in order to estimate the frequency of the motion, assuming periodicity; in Bundschuh et al (2007) a volume of interest (VOI) is manually defined in a summed image around the lesion, the center of mass (COM) of the activity distribution inside the VOI is determined and the z-coordinate of the COM as a function of time delivers the respiratory signal. These methods require image reconstruction whereas others only rely on the acquired raw data, and do not need any manual input (e.g. for drawing ROIs/VOIs) or the use of external radioactive sources: in Schleyer et al (2009) the variation of the counts within regions subject to respiratory motion is used to estimate the respiratory signal (Spectral Analysis Method, SAM); in Thielemans et al (2011) principal component analysis (PCA) is applied to the listmode PET data, and the respiratory signal is obtained as the principal component weight factor whose frequency spectrum has the highest peak in the respiratory band; in Wachinger et al (2012) an alternative dimensionality reduction technique is utilized, i.e. laplacian eigenmaps (LE) decomposition, in order to obtain the respiratory signal from the raw data; in He et al (2008) the non-uniformity of the geometric sensitivity of the scanner is exploited as a means of obtaining information on motion, since the count rate for a given organ will depend on the axial location of the organ within the scanner, the respiratory phase is determined from count rate changes in the listmode data (SENS method); in Kesner and Kuntner (2010) a respiratory trace is obtained from the rebinned sinogram by combining together information from the time activity curves extracted from all the voxels of the sinogram (sinogram region fluctuation method (SRF)).
However, many of these methods, and in particular SAM and PCA, provide respiratory signals that are only determined up to an arbitrary scale factor. As this scale factor can also be negative, this means that the direction of the represented motion is undefined and the maxima in the signal could refer either to the end-inspiration or end-expiration phase of the breathing cycle. The possibility of obtaining a signal that represents motion in the opposite direction compared to the truth could lead to inaccurate motion correction, inconsistencies between adjacent beds and an incorrect matching between PET and CT gates when performing attenuation correction, e.g. when using an end-expiration CT. In order to correctly identify the gates, the correspondence between the sign of the signal and the direction of motion has to be determined and fixed. This issue has been previously addressed in Schleyer et al (2011), where a method requiring the registration in the axial direction of the gated PET images (reconstructed without attenuation correction) is implemented.
In this work we propose two novel methods, CorrWeights and CorrSino, for the determination of the direction of motion of PCA respiratory signal, we implement a published method (Schleyer et al 2011) and compare the performance of the methods on simulation data and FDG lung cancer patients' data.

Extraction of respiratory signal
The data-driven (DD) method implemented in this work is PCA, as in Thielemans et al (2011), and in this section it will be described in detail as it is necessary for the later introduction of the proposed novel methods. To describe the method we will assume single slice rebinning (SSRB) (Daube-Witherspoon and Muehllehner 1987) has been performed on 3D PET data. The following coordinate system will be used (Bailey 2005): the radial position r is the radial displacement of the line of response (LOR) in the transaxial plane, the axial coordinate z, the azimuthal angle φ is the angle formed between the LOR (projected onto the transaxial plane) and the y-axis of the scanner, and the time index t. All coordinates are discretized. The dynamic sinogram data are therefore described as the vector φ d r z t , , , ( ) , and the empirical temporal mean on the total number of time frames T is given by: PCA is a 'dimensionality reduction' technique that attempts to find a linear transformation between the original data and a space of lower dimension (Pearson 1901), where the structure of the original data can be better observed. PCA represents the data as a combination of orthogonal basis vectors that are ordered with respect to descending variance, such that the first vector has the largest variance and thus represents the maximum observed variation. Given the data φ d r z t , , , ( ) , PCA provides the following linear expansion: ) are the principal components and w k (t) are the weight factors, that are the coordinates of the sinogram data described in the new basis, that has K vectors. The equality in (2b) derives from the orthonormality of the Principal Components, as they are basis vectors. Usually only the first three components contain motion related information so in our implementation K = 3.
PCA can be directly applied to PET dynamic sinograms as they are acquired by scanners but there are some issues related to this, e.g. the noise in the data and also the high spatial resolution, that provides more data than is actually needed by PCA to detect the biggest variations in time. In our implementation we therefore unlist the PET listmode data into dynamic sinograms (with time frame duration of 500 ms) with reduced spatial resolution. Furthermore, these are processed using the Freeman-Tukey variance stabilisation technique (Freeman and Tukey 1950) to obtain approximately normally distributed samples, as described in (Thielemans et al 2011): then spatial filtering is performed on φ d r z t , , , FT ( ) and finally PCA is applied. Note that this variance stabilization process effectively pre-whitens the sinogram data before the application of PCA.
The most respiratory-like signal is obtained by selecting the weight factor w k (t) whose main frequency most closely corresponds to respiration according to the following procedure (see figure 1): the power spectrum (where ŵ is the FFT of the signal at frequency f ) is computed for the three weight factors, the ratio R k between the peak in the respiratory frequency band ([0.1, 0.4] Hz) and the mean power above 0.4 Hz is evaluated, the w k (t) that provides the highest R k is selected (see figure 1(b)). The PC corresponding to the selected weight factor is defined as the respiratory PC (RPC) and contains information on the respiratory motion.
The sign of the respiratory signal thus obtained is not unequivocally determined: in (2a) the result of the product φ p r z w t , , therefore making the sign of the signal returned by PCA arbitrary, i.e. independent of noise and implementation details. Since the RPC weight factor (hereinafter referred to as RPC signal) needs to represent the physical respiratory motion, where maxima and minima should correspond to the end-expiration and end-inspiration phase respectively, the correct sign of w k (t) has to be determined prior to using it as a gating signal. Furthermore, consistency in the relationship between the signal and the direction of motion is required when combining adjacent bed positions.
As we are considering chest acquisitions, cardiac motion will also be present in the acquired data. Data-driven methods have proven to be capable of detecting cardiac motion (Büther et al 2009, Thielemans et al 2014, however the temporal resolution of the unlisted sinograms has to be high enough in order to be able to detect changes with cardiac frequencies, that are generally in the range of 1-2 Hz. In this work, the temporal resolution of 500 ms is designed to reduce the sensitivity of PCA to the cardiac contraction.

Sign determination methods
In order to determine the direction of the motion represented by the DD signal, Schleyer et al (2011) assume that respiration causes a cranio-caudal motion of the internal organs. The method analyses reconstructed images of the gated data. In this paper, we developed two new methods based on the same assumption that axial motion is the major cause of changes in the acquired data. However, these methods only require the PET sinograms and neither gating nor reconstruction are needed. All the following described methods take as input PET listmode data, that are then unlisted into dynamic sinograms φ d r z t , , , ( ) with a temporal resolution of 500 ms. The DD respiratory signal is obtained by applying PCA on sinograms that have undergone SSRB (see section 2.1), and radial and angular mashing (i.e. the summation of adjacent sinogram elements with respect to r and φ to reduce the dimension). (2011) performs registration via 1D rigid translation on gated images reconstructed without attenuation correction (to avoid issues with misalignment with the CT). The acquired data are gated with displacement gating based on the DD signal, normalized for detection efficiencies and reconstructed with FBP (with ramp filter). The registration process is performed on the sagittal and coronal Maximum Intensity Projection images (MIP) of the gates, in order to reduce noise. For each gate the displacement in axial direction (z coronal or z sagittal ) with respect to a reference gate is estimated. The final registration displacement is given as the average of the two.

Registration-based method. The method as proposed in Schleyer et al
We implemented the method with minor modifications. Each of the obtained gates is in turn taken as the reference for the registration and an average of the obtained displacements is obtained, rather than choosing as single reference the gate with most counts. We did this because our displacement gating was set to distribute the counts among the gates as evenly as possible to avoid large differences in noise between the gates. In addition, we have minimised the 1 norm of the sum of the absolute differences in pixel values, rather than the 2 norm, to increase robustness. This process is performed both with sagittal and coronal MIPs (obtaining z coronal,avg and z sagittal,avg ) and the direction of motion is defined by looking at: The weights w 1 , w 2 and w 3 related to the first three PCs, obtained from one patient dataset. (b) Power spectrum of the 3 weights. The peak within the respiratory frequency band is highlighted in red, the mean power above 0.4 Hz is displayed in lightblue and the value of their ratio for each weight is: R 1 = 15.0, R 2 = 2.9 and R 3 = 3.7. In this case, w 1 would be chosen as the respiratory signal.
if > sign 0 the sign of the PC respiratory signal is kept as it is, if < sign 0 then the respiratory signal is reversed. We also implemented a variation of this method by registering, as opposed to the MIPs, the images obtained by collapsing the 3D images in the sagittal or coronal plane (equivalent to summing together all the 2D sagittal or coronal planes). The two implementations are referred to as MIP and SUM.

Sinogram-based methods.
Under the assumption that motion occurs uniquely in the z-direction, we derived an approximation for the dynamic sinograms.
If we assume that motion occurs only in the z-direction, the change in the dynamic sinogram can be represented by d(z, t) (we omit the dependencies on r and φ in the following steps). Since the sensitivity of the scanner is constant in time and depends on z, in order to determine the time dependence of the sinograms, we will first consider how the sinograms normalised for the axial pattern in the sensitivity change over time: where s(z) is the sensitivity of the scanner (determined from the number of rings and axial acceptance of the scanner) and n(z, t) the normalized sinograms, that are going to be used in the following steps. The normalization of the sinograms is performed after SSRB. Assuming small translations δz t ( ) with respect to the normalized sinogram in t = 0, the normalized dynamic sinogram can be described as follows: and writing the Taylor expansion up to the first order we get: with ′ n z ( ) the gradient of the sinograms along the z-direction. Subtracting the temporal mean we get: where n z( ) and δz are the temporal mean of n(z, t) and δz t ( ). The left-hand side term of (8) can be obtained directly from the dynamic sinogram, whereas if ′ n z 0 ( ) were extracted from the initial temporal frame (t = 0) it would be too noisy. We therefore chose to approximate n 0 (z) with the temporal mean of the entire sinograms n z( ), consequently decreasing the noise. The gradient along z can be approximated using finite differences and is defined by z GradSino( ): As a result the 'component' ′ n z 0 ( ) in (8), that is multiplied by the function describing the axial motion, is approximated by GradSino: In order to reproduce the process performed on the sinograms that is ultimately utilised for PCA, we approximate the Freeman-Tukey transformation on the original dynamic sinogram d(z, t) and apply a Taylor expansion, given the assumption that the motion in the z-direction is small (see appendix for detailed description). The final result is that the dynamic sinogram that is given to PCA can be represented by the following (the temporal mean is subtracted): where: This description of the temporal changes of the sinograms from its mean resembles the result expected from PCA: it represents the dynamic sinogram with a product of a temporal one-dimensional signal with a component in the sinogram space. This similarity could be exploited to compare the assumed motion along the axial direction to the PCA output. The respiratory signal of interest is represented by the term δ δ − z t z ( ) , and can be described as the 'weight' of the component represented by WeightedGradSino. Figure 2(a) shows an example of the RPC obtained from the application of PCA on a time interval (i.e. temporal portion) of the acquisition of a patient, figure 2(b) shows the temporal mean of the sinograms and figure 2(c) shows the resulting WeightedGradSino (where it can be seen that 4 planes are set to zero at the top and at the bottom of the projection, in order to avoid issues with boundary conditions when applying (9)).
The above approximations motivate two new methods that only require the PET raw data (listmode data that are consequently unlisted into sinograms) to determine the sign of the RPC signal, and these are presented in the following sections.
CorrWeights method. The correspondence between the WeightedGradSino and the direction of motion is known as it derives from the application of the gradient in the z-direction, see (9), therefore its corresponding weight is related to the correct direction of the motion. An increase in the signal would correspond to motion towards the head ('up') while a decrease corresponds to motion towards the feet ('down'), thus maxima would correspond to the endexpiration phase and minima to the end-inspiration phase.
As a consequence of (11), the WeightedGradSino weight can be obtained making use of (2b). Therefore evaluating the inner-product between the WeightedGradSino and the dynamic sinogram on the left-hand side of (11): where the sums are over all sinogram elements. Examples of an RPC signal w and the WeightedGradSino weight u corresponding to the same time interval of PET data is shown in figure 3, where high similarity between the two signals can be observed. In practice, the signal u cannot usually be used as a good respiratory signal since (13) is based on various approximations (viz. that the changes in the dynamic sinogram are small and only in the axial direction, see (6) and (7)). Its importance lies in its known relationship with the direction of the axial motion, that makes it a comparison tool to determine whether w correctly represents the direction of the physical motion of the patient.
These arguments lead us to propose the following method to determine the direction of the DD respiratory signal: if the correlation between the two signals w and u is calculated, the sign of the correlation gives an indication about whether the DD respiratory signal and WeightedGradSino signal represent motion in the same direction. The Pearson correlation between the two signals is given by: If the correlation is negative it is an indicator that w represents a motion in the opposite direction compared to the motion related to u, thus the sign of w (i.e. the respiratory signal used for gating) has to be reversed. For the example shown in figure 3 (where only a 60 s interval of the total 300 s signals is displayed) the correlation is 0.98, while in other cases it can be lower. This method is not specifically related to PCA, therefore could be utilized to determine the correct direction of the motion of respiratory signals obtained with any other DD method.
CorrSino method. Considering the similarity between (2b) and (13), where the latter is expected to hold when the assumptions of small motion along the z-direction are met, it could be argued that WeightedGradSino should be similar to the RPC to a certain extent (see figure 2). WeightedGradSino could therefore be expected to be representative of the biggest changes that occur in the dynamic sinogram because of respiration, while additionally having a known relationship with the direction of axial motion. For this reason, to check if the sign of the RPC signal corresponds to axial motion in the desired direction, the RPC itself can be directly compared to the WeightedGradSino, evaluating the correlation between the two sinograms: The RPC and WeightedGradSino undergo a masking process prior to the evaluation of the correlation, in order to consider only non-zero bins. The absolute value of the correlation relates to the degree of agreement between the RPC and the WeightedGradSino sinogram, taking higher values in the case the gradient sinogram adequately represents the motion detected by the respiratory PC. If the correlation is negative the sign of the RPC signal has to be reversed, as in the CorrWeights method.
In figure 4 an example is shown of the appearance of the RPC when it has the wrong sign: the features and intensity patterns look similar in WeightedGradSino (4(b)) and in RPC (4(a)) and their correlation, from equation 15, is equal to 0.57, whereas 4(c) is the negative version of 4(a) (−RPC) and would be the result of the PCA related to the 'wrong' direction of motion. The correlation between 4(b) and 4(c) would result in −0.57 and the method would in this case reverse the corresponding respiratory signal w.
This method is specific to PCA, as it relies on the comparison of the RPC to the WeightedGradSino.

Phantom simulations.
As an initial test of the performance of the proposed methods, they were applied on simulated data produced with the XCAT phantom Segars et al (2010). In order to make the simulation realistic, the RPM respiratory trace of a normal-breathing patient was used as a template for the phantom motion. Diaphragm maximum displacement was set to 2 cm and anterior-posterior maximum displacement to 1.2 cm (default in the XCAT parameters). The duration was set to 50 s and the temporal sampling to 500 ms. Activity and attenuation images were produced. The activity images were projected through the use of STIR utilities (Thielemans et al 2012), simulating 3D projections in a GE Discovery STE PET/CT scanner (Teräs et al 2007), and their projections were attenuated using attenuation coefficients obtained from the attenuation images. Poisson noise was added to the projections, and no random nor scattered events were simulated (total number of counts per frame ∼ × 1.2 10 5 , which is similar to half the counts in clinical data, considering there are no random nor scatter events).

Patient data.
To test the performance of the four sign-determination methods we used FDG data of oncology patients. Patient datasets were acquired in 3D listmode on GE Discovery STE (21 datasets) and GE Discovery 690 (Bettinardi et  PET/CT scanners using activity levels for routine clinical protocols. The acquisition times ranged from 240 s to 720 s per bed position and only the lower chest bed positions were selected for this study. The acquisitions were monitored by the Varian Real-time Position Management TM (RPM) device, whose signal was used as the comparative standard for the evaluation process. For each patient the selected bed acquisition was also subdivided to independent shorter time intervals of 50, 100, 200 and 300 s in order to assess the performance of the methods on smaller amounts of data (i.e. for a 360 s acquisition, we selected 7 intervals of 50 s, 3 of 100 s, 1 of 200 s and 1 of 300 s).
2.3.3. Data processing. The 3D sinograms provided by the scanner (and those obtained from the XCAT simulation) consist of 553 2D sinograms of × 329 280 in φ r, for Discovery STE and 553 2D sinograms of × 381 288 in φ r, for Discovery 690, the axial field-of-view is for both scanners 157 mm and the listmode data were unlisted into sinograms with time frames of 500 ms. Before the application of PCA and consequently the sinogram-based sign determination methods, the sinograms underwent SSRB and radial and angular mashing: the number of summed angles was 40 for the Discovery STE and 32 for Discovery 690 and the number of summed radial bins was 4 in both cases, so that the final data format was × × 82 7 47 (radial positions × angles × transaxial planes) for Discovery STE datasets and × × 95 9 47 for Discovery 690 datasets.
For the registration-based method, the procedure of section 2.2.1 was used and images were reconstructed with voxel dimensions of × × 5.99 5.99 3.27 mm 3 for Discovery STE datasets and × × 5.33 5.33 3.27 mm 3 for Discovery 690 datasets.

Evaluation
The evaluation followed the following steps: (i) for the selected interval of the patient listmode file and for the projections of the simulation data, the PCA-method described in section 2.1 was used to produce the RPC and its respiratory signal, with a time frame duration of 500 ms; (ii) the four methods (CorrSino, CorrWeights, MIP and SUM ) were then applied on the selected intervals and the obtained RPC signal was corrected according to the method's response; (iii) the correlation was evaluated between each of the obtained respiratory signals and the corresponding RPM signal (downsampled to the same time resolution as the RPC signal). When the correlation is positive, the respiratory signal is correct, when the correlation is negative the respiratory signal is opposite to the RPM, therefore the method failed in correcting the sign.

Results
The application of the presented sign determination methods on XCAT data resulted in each case in the correct determination of the signal, and the correlation values obtained from the application of CorrWeights and CorrSino were 0.81 and 0.3 respectively (see (14) and (15)). For the patient data, examples of the appearance of the RPC and WeightedGradSino sinograms, and of the RPC signal and u(t) have been shown in figures 2 and 3. The application of the methods resulted in the failure rates reported in table 1, that includes all the acquired patients and all the different time duration intervals. Table 1 shows that methods CorrSino and CorrWeights have a more stable behavior with respect to the duration of the time intervals, compared to the registration-based methods MIP and SUM that show much higher failure rates in the 50 s and 100 s cases. The reduced amount of data that is utilized in such cases affects the quality of the reconstructed gated images, therefore the registration process is more likely to fail.
The DD output depends on the actual presence of respiratory motion and on the assumption that it represents one of the main causes of variation in the acquired data. The quality of the provided respiratory trace, considered as degree of agreement with the RPM signal, is therefore not predictable and can vary between different patient acquisitions and also between different time intervals of the same study. The same issues exist for the sign determination methods, as they too rely on the same assumptions as the DD method, and their outcome might vary with respect to different acquisitions and different respiratory signals.
In order to better understand the trend of the results with respect to the 'respiratorylikeness' of the PCA signal w(t), the failure rates have been analysed with respect to the correlation between w(t) and the RPM signal. The majority of the obtained signals w(t) (401 cases, equal to 93%) are highly similar to the RPM trace, showing correlations ranging from 0.75 to a maximum of 0.99, while in the remaining 7% of the intervals (30 cases) the quality of the PCA signal is sub-optimal, with 21 intervals with correlation included in To compare in more detail the performance and the reliability of CorrSino and CorrWeights, we examined the numerical results for each patient. In figure 5 the results obtained by the application of both methods are displayed per patient (for the Discovery STE datasets): each point in the plot corresponds to one evaluated interval and its absolute value corresponds to the absolute value of the methods' outcome (see (14) for CorrWeights and (15) for CorrSino). Each point is assigned a positive or negative sign, depending on whether the method succeeded or failed in detecting the correct direction of motion in that specific interval. For each patient, the displayed values correspond to the evaluated intervals ordered from the shortest to the longest in duration, and intervals of the same duration are ordered chronologically (i.e. [250,300], [300,350], [350,400],..). Figure 5 shows that the sinogram-based methods CorrWeights and CorrSino fail in the same cases (apart from two 50s intervals where CorrSino fails and CorrWeights succeeds). This could be expected considering that both methods rely on similar assumptions and therefore are potentially vulnerable to cases where the assumptions are not met, e.g. when the biggest changes in the acquired data are caused by bulk motion. Nevertheless it is clearly visible that CorrWeights values are considerably higher than those provided by CorrSino, therefore suggesting CorrWeights might be considered more reliable. For Patient 1, Patient 4 and Patient 5 both methods failed in several time intervals, possibly due to the patient moving during the acquisition. The plot also suggests that the results can vary considerably from patient to patient. This can be expected considering that the breathing patterns can be noticeably different between patients and also that interference due to other types of motion can be quite variable.

Discussion
When using the simulation data, all the presented sign-determination methods correctly detected the direction of axial motion, even in the presence of anterior-posterior motion. This suggests that the rationale for the two sinogram-based methods is reasonable. Moreover, the motion in the XCAT simulation data is non-rigid and relatively large (2 cm diaphragm movement and 1.2 cm chest expansion), therefore we have shown that the methods, that suppose motion is in axial direction, small and rigid, perform well even in this instance. Nevertheless, the structure and motion of the XCAT phantom is relatively simple and the evaluation of the methods on the patient data is therefore essential. When applying the methods on the available patients datasets, some failures occurred. CorrSino and CorrWeights proved to be generally more accurate than SUM and MIP in detecting the correct sign of the respiratory signal, providing an overall failure rate of 5.8% and 5.3% on the 431 evaluated intervals, compared to 18.0% and 12.8%. The former two methods also have the advantage of avoiding the reconstruction step making use of the raw PET data only. Table 1 also shows that the sinogram-based methods are more stable compared to the registration-based ones with respect to the amount of data that is used. In the majority of the cases in our cohort, which includes the lower chest bed position of FDG scans of oncology patients, the RPC signal was highly correlated with the external RPM signal (93% had correlation higher than 0.75).
The failure rates displayed in table 2 show that all sign-determination methods are more likely to fail as the correlation between RPC signal and RPM decreases, i.e. when the RPC signal is not adequately representing the external motion of the chest, although there is insufficient data (only 30 of the total 431 evaluated intervals) to draw significant conclusions about which method is best in these circumstances. For the cases with correlation higher than 0.75, the sinogram-based methods are clearly performing better compared to the registration-based methods, both failing in less than 3% of the cases.
The more detailed analysis carried out on CorrWeights and CorrSino shows that the correlation between the WeightedGradSino weight and the RPC signal is generally higher than the correlation between their corresponding sinograms (see figure 5). This indicates that CorrWeights could be a more reliable technique than CorrSino. Furthermore, the majority of the failures of CorrWeights arise when its result has an absolute value below 0.2, suggesting that a threshold could be defined to only use the method when it produces high values and is therefore highly likely to succeed. The determination of such a threshold will be the subject of future investigation.
In addition, it is worth noting that the CorrWeights method can be applied to fix the sign of the signal obtained from any DD method as it does not rely specifically on the use of PCA.
In general, failures in determining the sign of the DD respiratory motion with the presented techniques could be due to several reasons. Low contrast and noise in the data hinder the extraction of a good DD respiratory signal. Minimal motion along the z-direction contradicts the assumption at the core of the investigated methods. Also other types of patient movements during the acquisition might contaminate the RPC signal. Nevertheless figure 5 clearly shows that for the sinogram-based methods the strongest variability in the outcome is found between different patients and also that the duration of the interval of data taken into account does not have a noticeable impact.
In our implementation all of the data in the sinograms are used both for the extraction of the respiratory signal and for the determination of its sign. Selecting regions of interest (ROI) known to include respiratory motion might potentially improve both the quality of the signal and the performance of the methods. However, defining ROIs would require reconstruction of the data, as anatomical features are difficult to detect in sinogram space, and would therefore require additional processing time and external input from trained personnel. As the goal of our work is to provide a methodology which is exclusively dependent on the acquired PET raw data, we have not put this into practice.
The evaluation on patient data of the proposed methods has been performed via comparison with the RPM signal. This has been shown to potentially exhibit a time lag when compared to the internal organ motion (Ionascu et al 2007), while DD methods are expected to provide signals that are closely representative of the latter. However, as this time lag is expected to be small (below 0.2 s in Ionascu et al (2007)), we do not expect it to undermine the conclusions of this work, as it would only decrease the absolute value of the evaluated correlation between RPM and RPC signal but not its sign. Nevertheless a limitation of this study is that it relies on the RPM signal correctly representing the respiratory motion of the patient. A poor correlation between RPC signal and RPM could therefore be due to a deficient RPM rather than inaccurate RPC signal. However, the RPM signal is currently utilized in clinical practice and is usually considered to be a reliable method.
As for potential irregularity of the breathing in time, we expect the methods not to be affected by it, as long as the respiratory motion is similar throughout the acquisition. What could potentially change is the value of the ratio shown in figure 1 for an example of PCs: if the breathing is irregular the peak in the frequency band would be expected to decrease, and so would the resulting ratio value.
Future work will include a further analysis of the failures of CorrSino and CorrWeights and, in addition, the design of an informative metric for the goodness of the obtained signals, in order to provide confidence in the RPC signal in the absence of an external device (in our work RPC signals were compared to the RPM signal). The presented work is limited to the use of the PCA and the implemented methods were applied only to static studies in FDG oncology patients, therefore future developments and applications could include the extension to other DD algorithms and also to different types of acquisitions.

Conclusion
We presented two new methods, CorrWeights and CorrSino, for the determination of the direction of the data-driven respiratory signal obtained with PCA. The new methods provided lower failure rates compared to the previously published registration-based methods. When the extracted respiratory signal is highly correlated with the RPM they fail in less than 3% of the cases, compared to more than 10%. Among the two proposed methods, CorrWeights is not specifically related to PCA, therefore it could be utilized to determine the correct direction of the motion of respiratory signals obtained with any other DD method. Considering its simple implementation, it could easily be applied together with any DD method in clinical practice.