Motion Dependent and Spatially Variant Resolution Modeling for PET Rigid Motion Correction

Recent advances in positron emission tomography (PET) have allowed to perform brain scans of freely moving animals by using rigid motion correction. One of the current challenges in these scans is that, due to the PET scanner spatially variant point spread function (SVPSF), motion corrected images have a motion dependent blurring since animals can move throughout the entire field of view (FOV). We developed a method to calculate the image-based resolution kernels of the motion dependent and spatially variant PSF (MD-SVPSF) to correct the loss of spatial resolution in motion corrected reconstructions. The resolution kernels are calculated for each voxel by sampling and averaging the SVPSF at all positions in the scanner FOV where the moving object was measured. In resolution phantom scans, the use of the MD-SVPSF resolution model improved the spatial resolution in motion corrected reconstructions and corrected the image deformation caused by the parallax effect consistently for all motion patterns, outperforming the use of a motion independent SVPSF or Gaussian kernels. Compared to motion correction in which the SVPSF is applied independently for every pose, our method performed similarly, but with more than two orders of magnitude faster computation time. Importantly, in scans of freely moving mice, brain regional quantification in motion-free and motion corrected images was better correlated when using the MD-SVPSF in comparison with motion independent SVPSF and a Gaussian kernel. The method developed here allows to obtain consistent spatial resolution and quantification in motion corrected images, independently of the motion pattern of the subject.

studied [1], [2]. For example, factors such as positron range, photon acollinearity and detector size, determine the spatial resolution properties of the PET scanner. Using this knowledge, several methods have been proposed to correct the loss of spatial resolution in PET scanners. Using iterative reconstruction algorithms, such as maximum-likelihood expectationmaximization (ML-EM), the loss of spatial resolution can be compensated for by introducing the resolution model in the system matrix [2], [3]. The scanner spatial resolution can be modeled by the system point spread function (PSF), which is the response of the imaging system to an impulse input. The PSF of the system can be estimated in several ways. For example, by experimental measurement of point sources distributed in the scanner field of view (FOV) [4], Monte Carlo simulations [5], and analytical modeling of the detection process [6]. A simplified model can be used assuming a spatially invariant Gaussian PSF with the same width of a point source measured at the center of the scanner FOV (CFOV). More realistic representations consider models which can capture the spatially variant and asymmetric shape of the true PSF [7], [8] usually observed in PET scanners. In addition, the resolution model can be implemented in the projection space [9] or in the image space [3]. However, there is minimal difference in image quality between implementation of both methods [4].
Recently, research has been carried out to perform brain PET scans of awake small animals to circumvent the confounding effect of anesthesia [10]. In addition, in scans of freely moving animals the behavior of the animal can be measured during the PET scan [11], [12]. These methods require to perform rigid motion correction on the PET data to obtain brain reconstructions unaffected by motion.
Despite of the fact that there is a great amount of research in PET motion correction [12]- [15], the vast majority of the spatial resolution correction methods have been developed to be implemented in motion-free PET image reconstruction. Given the advancement in the field to perform PET scans of freely moving animals, and its promise to change the paradigm of preclinical research [16], improvement of spatial resolution in motion corrected PET reconstructions is becoming more relevant.
In PET scans of freely moving animals, the image PSF of the reconstructed image is both spatially variant and motion dependent since the animal head can traverse the whole scanner FOV during the scan [11], [12], [14]. For example, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a scan of a freely moving animal that remains most of the scan time close to the scanner CFOV (where spatial resolution is usually better) will present a better spatial resolution than the scan of an animal remaining close to the edge of the FOV. This effect contributes to the observed loss of spatial resolution in the motion corrected images in comparison with motionfree reconstructions [13], leading to inaccurate quantification. Angelis et al. [17] have proposed to measure the motion dependent PSF in the image space by attaching a point source to the subject during the scan and fitting a mixture of Gaussian functions to the reconstruction of the point source. In this way the anisotropic nature of the PSF was captured, however it was assumed that the PSF is spatially invariant. Yao, et al. [18] proposed to model the spatial resolution for every line of response (LOR) for the motion compensation reconstruction algorithm MOLAR [19]. This approach performed similarly to the use of an isotropic Gaussian model in real data scans.
As pointed out by Chan et al. [20], two different approaches to perform motion correction in PET reconstruction using ML-EM can be considered. The first involves transformation of the reconstruction image from the reference (motion-free) pose to the original pose of measurement before forward projection of the image, followed by the inverse image transformation of the correction image. This approach does not require transformation of the system matrix, but the two image transformations need to be performed for every pose. The second approach instead involves transformation of the system matrix, and therefore the image is kept at the reference pose, requiring no image transformation for every pose. The second approach is particularly advantageous for brain rigid motion correction, where usually a great amount of object poses are measured, since only the LORs need to be transformed instead of the whole image for every pose.
Here we propose to analytically calculate the spatially variant and motion dependent PSF in the image space for rigid motion correction PET scans. The method is suited for the second motion correction approach mentioned above (system matrix transformation). We model the motion dependent and spatially variant PSF for every voxel in the motion corrected image by calculating the superposition of all the PSF's that the voxels traverse due to motion. Using the motion tracking information, the parametrized scanner PSF is sampled at the position each object voxel was originally detected at every time point. We performed phantom scans of two capillaries undergoing different motion patterns as well as scans of a moving resolution phantom for validation of the method. Finally, the method was applied to data from a freely moving mice experiment.

A. Scanner
All experiments were performed on a Siemens Inveon microPET scanner (Siemens Medical Solutions, Inc., Knoxville, USA). The scanner has 25600 detector crystals with a size of 1.5 × 1.5 × 10 mm arranged in 64 blocks of 20 × 20 elements. The FOV has an axial length of 127 mm and transaxial diameter of 99 mm. Therefore, the maximum radial distance from the center of the FOV is 49.5 mm. Images are reconstructed in a grid of 128 × 128 × 159 voxels with a size of 0.776 × 0.776 × 0.796 mm along the x, y and z directions respectively. The z direction coincides with the scanner axial direction. The spatial resolution is about 1.5 mm at the CFOV (full width at half maximum, FWHM) and larger than 3 mm at the edge of the FOV [21].

B. Motion Tracking
The motion of all moving phantom experiments and awake animal experiments, was tracked using the point source tracking (PST) method [22]. Briefly, sodium polyacrylate point sources were soaked in [ 18 F]FDG until they reached an activity greater than 222 kBq. Four point sources were pasted on the subject to track its motion. The list-mode PET data was processed in short time frames of 32 ms ( t = 32 ms) and the point sources were tracked in those frames in the image space. The six degrees of freedom pose of the subject T k (k = 1 . . . K ), represented as a homogenous rigid transformation, is determined for every frame. The PST has a tracking accuracy of 0.24 mm. Details and validation of the tracking algorithm can be found in [22].

C. Image Reconstruction
Motion-free scans are reconstructed using ordered subsets list-mode reconstruction (OSEM) without attenuation, scatter nor randoms correction using a line integral model: where λ q j is the number of annihilation events in voxel j ( j = 1 . . . J , voxels) at iteration q, g l j is the intersection path of detected LOR l with voxel j , s j is the sensitivity correction factor for voxel j , w is a weight factor for LOR ( = 1 . . . L, total number of LORs) to incorporate normalization (detector sensitivity) correction, and the detected events are subdivided in subsets S m (m = 1, . . . , M, subsets).
For the motion correction reconstruction, once the motion of the subject has been tracked over the entire scan, the average pose is calculated and is defined as the reference pose T re f , to which all LORs are transformed to. Motion correction reconstruction is performed using the list-mode reconstruction algorithm proposed in [23]: where g l j is the intersection path of LOR l with voxel j after transformation of the LOR l to the reference pose and s j is the motion compensated sensitivity correction factor for voxel j , calculated by interpolation of the voxel sensitivity factors from the sensisitivty image in (2) using all measures poses [23]. In the motion-free case, the scanner, represented as gray rectangular detectors around the pixelated image, it's coordinate system and the image coordinate systems (in red), are aligned. The shape of the scanner PSF varies across the scanner space, thus it is defined with respect to the scanner coordinate system. In this example, the pose of the object in the motion-free case will be defined as the reference pose T ref . b) In the moving scan case, the object moves anywhere inside the scanner space. Therefore, the PSF for every point in the object changes over time depending on its position in the scanner space. c) When motion compensation is performed, the object is transformed to the reference pose T ref using T d 0 . As a result, the scanner coordinate system (in blue) is transformed with the same transformation T d 0 . In the image coordinate system (in red), this transformation causes the PSF at every voxel, determined by the voxel position with respect to the scanner space, to change. d) The pose of the object is measured for every time point and e) motion compensation is performed for all time points. The spatially variant PSF for every voxel in the motion corrected image is therefore the weighted (by time duration) average of all transformed PSF's according to T d 0 .
D. Motion Dependent and Spatially Variant PSF 1) Motion Dependent Resolution Model: For a resolution model implemented in the image space, let us consider the activity distribution image of an object in 3D space λ(X) (X = [x yz] T ) and the spatially variant PSF of the PET scanner F(X). The blurred image λ (X) is then modelled as the linear transformation of the image using the matrix of resolution kernels: where the image λ(X) is a J × 1 vector and the matrix F(X) is a J × J matrix containing in every column the resolution kernel for each voxel, which can be different for every voxel, i.e. spatially variant. In the motion-free case the activity distribution λ(X) is considered to be constant over time (Fig. 1a), i.e. no tracer kinetics are considered. On the other hand, in scans where the object moves during the scan, the object activity in 3D space is time dependent (Fig. 1b, d).
The activity at every point in space is proportional to the object activity and the time the object remains in that point. The time dependent rigid motion of the object can be represented with a homogenous transformation T (t) = T t . This transformation defines the pose of the object at every time point. Therefore, the image of the moving object is obtained by integrating the activity of the moving object over time: where D is the total duration of the scan, and the matrix of resolution kernels is taken out of the integral since it does not depend on time. In PET rigid motion correction reconstruction, the object is tracked over time and the object poses T t are measured. These poses are used to compensate the motion of the object (Fig. 1c, e), moving back the object to a reference pose T re f at all time points. For this purpose, the differential transformation T d t that moves the object from its pose T t to a reference pose T re f is calculated: The moving object is then compensated for motion using this transformation: where λ MC (X) is the motion compensated object, stationary at T re f . After motion compensation, the activity of the object λ MC (X) is no longer time dependent, i.e. the activity distribution is constant at every point in space. However, after motion compensation, the object PSF becomes time dependent since the relative pose of the scanner with respect to the image coordinate system changes over time (Fig. 1c, e) according to the same transformation T d t . The measured blurred image of the objectλ MC (X) can be modelled by integrating over time the blurring of the moving object by the PSF, after both the object and the PSF have been transformed using T d t : Therefore, the motion dependent and spatially variant PSF (MD-SVPSF) is calculated by integrating the moving PSF over time: Finally, the motion compensated imageλ MC (X) is modelled as the product of the unblurred motion compensated object λ MC (X) with the MD-SVPSF kernel matrix: 2) Model of the Spatially Variant PSF: In this work we used the model proposed in [7] to represent the spatially variant PSF of the PET scanner. This model considers a radial internal, radial external, tangential and axial PSF width as a function of the radial distance. However, as the PSF width along the tangential and axial directions does not change much for the Inveon scanner [21], a fixed parameter was used for the PSF tangential and axial width. The model considers an asymmetric Gaussian with parameters: where ρ, α and τ are the radial, axial and tangential directions respectively, σ iρ (ρ) and σ eρ (ρ) are the internal and external radial width parameters, σ α and σ τ are the axial and tangential width parameters and μ ρ , μ α and μ τ are the mean radial, axial and tangential coordinates, respectively. These parameters define the shape and position of the scanner PSF: where H (•) is the Heaviside function and with the normalizing factor: Radial, axial and tangential coordinates in (12) are defined with respect to the radial, axial and tangential directions, which in turn are calculated from the Cartesian coordinates (x c , y c , z c ) of the PSF center. The radial direction has unit vector cos(θ c ) sin(θ c ) 0 , and the tangential direction −sin(θ c ) cos(θ c ) 0 , where θ c = atan(y c /x c ). The axial direction is aligned with the z axis. Fig. 2 shows the definition of the coordinate system and the relation between the Cartesian coordinates and the radial, axial and tangential coordinates.
The parameters in (11) where calculated from a non-linear least-square fit to 20 point sources distributed along the radial direction. The fit was performed in Matlab (The MathWorks, Inc., Massachusetts, USA) using the Levenberg-Marquardt algorithm. Parameters σ iρ and σ eρ where calculated as function of ρ by fitting a second order polynomial to the discrete values calculated on the 20 point sources, and σ α and σ τ were calculated as the average of all fitted values.

3) Calculation of the Motion Dependent and Spatially Variant
PSF: In the discrete case, calculation of (9) is performed by sampling the PSF for every pose k using (12), and summing all PSF's calculated at every pose. The PSF has to be sampled using the reconstructed image voxel size. To do so, instead of transforming the PSF using (6), as suggested in the derivation of the resolution model in (9), it is more convenient to transform the image space to the scanner space using the inverse transformation of (6). Then, the PSF can be sampled according to the image space sampling size at the corresponding scanner position. The calculation procedure is detailed below.
The MD-SVPSF is calculated for every voxel in the motion corrected image as follows. Initially, the voxels X j n in the neighborhood N j (n ∈ N j ) of every voxel j are considered to sample the PSF (Fig. 3a). For every pose k, the voxels in N j are transformed to their corresponding pose T k in the scanner FOV before motion compensation (Fig. 3b, d). Thus, the inverse of (6) is calculated: The voxel coordinates are then transformed as (Fig. 3b, d): where the central voxel has coordinates: and the central radial coordinate is: The coordinates of the central voxel, i.e. the center of the PSF, are used to calculate the radial and tangential coordinates of the voxels X k j n . Therefore, the radial and tangential directions are aligned with the x and y axes respectively, according to the angle θ k j with the x axis ( Fig. 3b, d): and then calculating the inverse rotation matrix about the z axis as (Fig. 3c, e): is used to rotate all voxels in N j : After performing this rotation, the x, y and z coordinates in X kr j n correspond to ρ, τ and α respectively. Finally, the PSF can be sampled at X kr j n using (12) with μ r = r k j , μ a = z k j , μ t = 0, and σ i (r k j ), σ e (r k j ), σ a and σ t defined by the fitted parameters, to obtain the kernel for voxel j at pose k, Q k j = f (ρ, τ, α). This procedure is performed for every pose and the average of all kernels at every pose is calculated, applying a weight factor to every kernel proportional to the time duration of the pose, i.e. t k : This kernel is normalized such that The MD-SVPSF is calculated for every voxel j in the motion compensated image. The voxels in N j are considered to sample the PSF. b) For every pose k, the voxels in N j are transformed to their original pose before motion compensation using T di k . c) Afterwards, voxels are rotated using R z (θ k j ). After rotation, voxels x, y and z coordinates correspond to radial, tangential and axial coordinates respectively. The PSF is sampled using these coordinates. d) Transformation to the original pose before motion compensation is performed for all poses, e) followed by rotation to sample the PSF. f) Finally, the sampled PSF's at every pose k are summed to calculate the MD-SVPSF.
where Q M D j is the motion dependent blurring kernel for voxel j in the motion corrected image. The algorithm to calculate the MD-SVPSF is summarized in algorithm 1.
Equation (10) can be written as: Therefore, the MD-SVPSF can be used in the motion correction reconstruction in the resolution modeling step as: where the product with the transpose of the blurring kernel matrix is calculated in the correction image as well as in the motion compensated sensitivity image. To reduce calculation time, MD-SVPSF kernels were calculated only for voxels in a predefined volume containing the object interest (phantom or mouse brain). Calculation of the MD-SVPSF kernels was performed on a NVIDIA GTX 1080 (NVIDIA Corporation) graphical processor unit (GPU) using CUDA [24]. The calculation was parallelized over the poses k which produced shorter calculation time than parallelizing over the voxels j .
The kernel size was set to 7×7×7 voxels as a compromise between calculation time and minimal truncation of the PSF in all directions. With this kernel size more than 4 ×σ (95% of the PSF) along the external radial, axial and tangential directions is considered in the entire FOV (see table 1). Only for the internal radial direction, at a radial distance > 40 mm, the kernel size is smaller than 4 × σ iρ . A kernel size of 9×9×9 voxels would consider more than 4 × σ iρ over the entire FOV, but the calculation time would increase about two times.

E. Validation Experiments 1) Moving Capillaries Phantom:
To assess the performance of the MD-SVPSF, a phantom with 2 capillaries was scanned while applying motion manually. Three different motion patterns were considered. Two thin glass wall capillaries (internal diameter 1.5 mm) were placed parallel to each other on a 2×10 cm foam platform with a spacing of 1 cm. The capillaries were filled with [ 18 F]FDG and scanned for 10 min for each of the three scans. Four point sources were fixed on the platform to track the phantom rigid motion.
Initially, the long axis of the capillaries was aligned with the scanner axial axis and the phantom was manually moved during the entire scan. For all motion patterns, motion consisted mainly of translation in the transaxial plane (different radial positions) and rotations around the axial axis. The first motion pattern considered translation and rotation confined to a small region midway between the CFOV and the edge of the transaxial FOV. The second motion pattern considered translation motion along the x axis over the full length. Finally, the third motion pattern consisted of translations along a circular trajectory at a radial distance close to the edge of the transaxial FOV. Heat maps showing the motion pattern for each scan are presented in Fig. 4 (a-c).
Motion corrected reconstructions were performed using (24) with 6 different resolution models. Three spatially invariant isotropic Gaussian kernels were used. The Gaussian kernels width were calculated from a Gaussian fit performed on Algorithm 1 Algorithm to Calculate the MD-SVPSF Kernels for Every Voxel in the Motion Corrected Imagee For  point sources located at a radial distance of 0 mm (CFOV), 24.5 mm (midway along the radial direction) and 49 mm (edge of the FOV). They had a FWHM equal to 1.60, 1.90 and 2.21 mm respectively. These kernels are referred to as Gs16, Gs19 and Gs22 respectively. In addition to spatially invariant kernels, the scanner's SVPSF calculated with (12), i.e. not motion dependent (motion independent, MI-SVPSF), and the MD-SVPSF were used. Finally, the reconstruction without resolution modeling (noRM) was performed. Reconstructions were performed with 16 subsets and 16 iterations. The motion corrected reconstructions were analyzed by first aligning the reconstructed capillaries with the axial axis. The average of 45 transaxial planes was calculated and the image was upsampled to a pixel size of 0.194 mm. The FWHM of the 2 capillaries was calculated in these images as follows. For each capillary the contour at 50% of the maximum intensity was calculated and the eigenvectors of the contour were determined. Then, the chords of the contour along the 2 eigenvectors were determined and the average of the 2 chords length was calculated as the capillary FWHM. In addition, we quantified the deformation caused by the parallax effect, as the deviation from the, ideally, circular shape of the capillaries in the transaxial plane. To do this, the ratio of the magnitude of the largest and smallest contour eigenvalues was calculated. This ratio is referred to as the eigenvalues ratio. A larger eigenvalues ratio reflects a more prominent deformation of the shape for round objects (eigenvalues ratio = 1). The analysis was performed individually for both capillaries, referred as "left" and "right" for the capillary placed towards the positive x axis and negative x axis, respectively. For all motion patterns, the left and right capillaries remained most of the time towards the positive and negative x axis respectively. The FWHM of both capillaries is reported as a function of the iteration number for every resolution model. The FWHM difference between both capillaries for each motion scan, as well as the eigenvalues ratio (average of both capillaries) are reported at iteration 16.

2) Comparison With a Pose-by-Pose SVPSF Reconstruction:
In order to compare the use of the MD-SVPSF resolution kernels with motion corrected reconstructions in which the SVPSF is performed for every motion pose, i.e. the SVPSF kernels are not averaged, the motion correction method proposed in [25] was implemented. For this approach, instead of transforming every LOR according to the subject motion (our approach), the image in the reference pose is transformed to the original pose in which the PET data was measured, and the LOR's are maintained in their original pose. The forward projection model for every LOR F l is then: where q j n is the SVPSF (motion independent) kernel, and M t (•) is a time dependent image transformation operator which transforms the image from the reference pose to the original pose of measurement at time t. Therefore, this transformation is performed for every pose k according to the time of measurement of LOR l. Moreover, the blurring operation using q j n has to be performed for every transformed image at every pose k. When the transposed operation is performed to calculate the correction image, the inverse of M t as well as the transpose blurring operation is performed for every pose k as well. In our implementation M t (•) consisted on a rigid transformation calculated with trilinear interpolation. We refer to the reconstructions using the forward model in (26) as PP-SVPSF (pose-by-pose SVPSF). The 3 moving capillaries phantom scans (detailed in the previous section) were reconstructed using PP-SVPSF, and were compared with reconstructions using the MD-SVPSF. Due to the long computation time using (26) only one minute of continuous motion data, corresponding to 1875 poses, was selected for each of the 3 moving capillaries phantom scans. The capillaries were analyzed in both reconstructions as described in the previous section.
3) Moving Resolution Phantom: A resolution phantom was scanned undergoing manual motion at two different positions; close to the CFOV and at an off-center position. The phantom was filled with [ 18 F]FDG and was scanned during 10 min at each position. The heat maps of the motion patterns are shown in Fig. 4 (d-e). The resolution phantom has six groups of rods, with diameters of 1.2, 1.6, 2.4, 3.2, 4 and 4.8 mm and spacing between rods is twice the diameter size. The two motion scans were corrected for motion and reconstructed using (24) with the same 6 resolution models detailed in the previous section. The number of iterations was increased to 40 (16 subsets) to allow convergence in the more complex phantom structures.
In addition to the qualitative assessment of the reconstructions, the deformation of the phantom rods was quantified using the eigenvalues ratio. After alignment of the rods long axis with the axial axis, the average of 20 transaxial slices was calculated and the image was upsampled to a pixel size of 0.194 mm. Then for every rod, in the groups of 2.4, 3.2, 4 and 4.8 mm, the eigenvalues ratio was calculated from the contour at 50% of the maximum rod intensity. The average eigenvalues ratio for every group of rods is reported for the centered and off-center motion corrected phantom using the 6 resolution models. 4) Freely Moving Mice Brain Scans: For this experiment data from a previous [ 18 F]FDG brain PET study in freely moving mice by our group [11] was used. The experiments followed the European Ethics Committee recommendations (Decree 86/609/CEE) and were approved by the Animal Experimental Ethical Committee of the University of Antwerp, Antwerp, Belgium (ECD 2016-89). Eight mice were scanned in awake state using the point source tracking in test-retest scans as well as in a challenge condition using the NMDA antagonist stimulant drug memantine. Awake animals were scanned for 20 minutes after an [ 18 F]FDG uptake period of 30 min. The animals were scanned in a cylindrical cage with a platform of 10×9 cm where they could move freely. For the memantine scans, the animals were injected with a memantine doses of 30 mg/kg (intraperitoneal) before [ 18 F]FDG injection. The dataset consisted of 24 scans (8 test, 8 retest and 8 memantine scans). In addition, the same protocol (test-retest and memantine) was performed in another 8 mice but using anesthesia during the PET scan to obtain motion-free reference images. The average injected activity was 19.0±0.44 and 17.7±1.99 MBq for the test-retest and memantine scans in awake animals, and 18.5±0.66 and 18.6±0.50 MBq for anesthetized animals respectively. Since in both awake and anesthesia scans, an awake [ 18 F]FDG uptake period was considered, brain uptake in both conditions should be similar. A detailed description of the study can be found in [11].
Awake animal scans were reconstructed with (24) using the motion dependent method (MD-SVPSF) and compared with motion corrected reconstructions using the motion independent method (MI-SVPSF) and considering resolution modeling using a spatially invariant Gaussian kernel with FWHM of 1.6 mm (Gs16). Anesthesia (motion-free) scans were reconstructed considering the same 1.6 mm FWHM Gaussian kernel. Reconstructions were calculated with 40 iterations.
Image processing was performed in PMOD 3.6 (PMOD Technologies Ltd.). The reconstructed images were upsampled to a pixel size of 0.2 mm and the brain was manually aligned with an [ 18 F]FDG mouse brain template. After manual alignment, a non-rigid registration (brain normalization) to the template was performed. The average of all brain reconstructions for each of the 2 conditions (test-retest and memantine) was calculated. Regional brain standard uptake value (SUV) quantification was performed in all brain regions, including hippocampus, thalamus and cortex, and the correlation between awake motion corrected (using Gs16, MI-SVPSF and MD-SVPSF) and anesthesia motion-free reconstructions quantification was calculated.

A. Parametrization of the Spatially Variant PSF
The contour of the point source image (experimental measurement) in the transaxial plane, as well as the model fit (12), are shown for the radial distances of 20.5 and 45.7 mm in Fig. 5. Values of σ iρ and σ eρ at the CFOV, midway between the CFOV and the edge of the FOV, and at the edge of the FOV are presented in table 1. The model replicates the shape of the point sources at both radial distances, specially the asymmetric elongation along the radial distance. However, at r = 45.7 mm, the shape of the point source slightly differs from a Gaussian shape, with a wider profile along the tangential direction (y axis in Fig. 5) at the radial external direction (positive x axis in Fig. 5).

B. Calculation of the MD-SVPSF Kernels
Calculation time of the MD-SVPSF on the GPU was on average 508 kernels (7×7×7) per second using 5500 poses. As an example, Fig. 6 shows the contour of the MI-SVPSF and MD-SVPSF kernel for the same voxel, in the 3 orthogonal planes. The MD-SVPSF was calculated for the motion corrected image of one of the awake mouse memantine challenge scans. The MD-SVPSF kernel is wider in the x − y and  x − z planes compared to the MI-SVPSF and becomes more elongated in the z − y plane. In all 3 planes the direction of the kernel's principal axes changes for the MD-SVPSF in comparison with the MI-SVPSF.

C. Moving Capillaries Phantom
The average speed of the phantom was 3.20, 3.12 and 3.33 cm/s for motion scans 1,2 and 3 respectively. Fig. 7 shows the FWHM as a function of the iteration number for both left and right capillaries in the 3 motion scans using the 6 different resolution models.
In all motion scans, the use of resolution modeling greatly improves the FWHM of the capillaries compared to the motion corrected reconstruction without resolution modeling (noRM). For motion scan 1, 2 and 3 the reconstructed FWHM is decreased when considering wider Gaussian kernels. However, there is some overcompensation for the 2 widest kernels (Gs19, Gs22), e.g. for Gs22 at 16 iterations there are 5 out of 6 FWHM values smaller than the capillary inner diameter (1.5 mm). For motion scan 1 using the MI-SVPSF and MD-SVPSF, the FWHM of both capillaries converges towards 1.5 mm. However, for motion scans 2 and 3 the FWHM using the MI-SVPSF does not converge to 1.5 mm, while for MD-SVPSF the FWHM converges to a value close to 1.5 mm. The reconstructed resolution is motion dependent as can be seen from the different reconstructed FWHM of the capillaries at iteration 16 for the different motion patterns for all resolution models. The smallest variation is found when using the MD-SVPSF where the maximal FWHM difference between the different motion patterns is 0.192 mm. In contrast, using MI-SVPSF the range is 0.363 mm. Among the Gaussian resolution models, Gs22 has the smallest variation at iteration 16, where the FWHM range is 0.239 mm.
The reconstructed resolution is also spatially variant as can be seen from the difference between the FWHM of the left and right capillaries for the same motion pattern. The difference between the left and right capillaries FWHM at the final iteration is shown in Fig. 8a. For all reconstructions using resolution modeling, difference is greater than noRM for all motion scans. However, as can be seen in Fig. 7, FWHM is much larger for noRM compared to the resolution models. When considering only the resolution models, it can be seen that for motion pattern 1 (red in Fig. 8) the differences between left and right FWHM is smallest for the MI-SVPSF and MD-SVPSF (<9%) and much larger when using the Gaussian kernels, (left-right differences ranging from 13% to 19%). For the other motion patterns all differences are similar among the different reconstructions and below 8% and 3% for motion patterns 2 and 3 respectively. Fig. 8b shows the results for the deformation quantification. Overall, the deformation of the transaxial circular shape of the capillaries is smallest when using the MD-SVPSF as seen in Fig. 8b. For motion scan 1 both MI-SVPSF and MD-SVPSF perform better than the other resolution models in terms of the eigenvalues ratio. However, for motion pattern 2 MI-SVPSF performs poorly, while MD-SVPSF still gives the smallest deformation. For motion scan 3 deformation is minimal for all resolution models with MD-SVPSF having the lowest eigenvalues ratio.   Table 2 shows the difference in capillaries FWHM between PP-SVPSF and MD-SVPSF at iteration 16. The minimum and maximum difference is 2.7 % and 7.1% respectively. Fig. 10a shows the FWHM difference between left and right capillaries using PP-SVPSF and MD-SVPSF for the 3 motion patterns. The minimum-maximum FWHM difference is 1.7-6.2% and 0.02-4.4% for PP-SVPSF and MD-SVPSF respectively. The capillaries eigenvalues ratio is similar between PP-SVPSF and MD-SVPSF as shown in Fig. 10b and table 3. For the 3 motion patterns and for both capillaries the difference in eigenvalues ratio between PP-SVPSF and  IV  TOTAL COMPUTATION TIME (SECONDS) FOR A MOTION CORRECTION  RECONSTRUCTION, CONSIDERING 16 SUBSETS, 16 ITERATIONS  AND 1875 POSES, FOR THE DIFFERENT RECONSTRUCTION  STEPS IN THE MD-SVPF Fig. 11 shows the motion corrected reconstructions of the resolution phantom using the different resolution models, at the center and off-center position. The average speed of the phantom was 1.85 and 1.23 cm/s for the centered and off-center phantom respectively. For all reconstructions, rods with a size of 1.6 mm and larger are recovered. For the 1.6 mm rods, Gs22 shows the most blurred image. In addition, edge overshoot artifacts are visible in the cluster of 1.2 mm rods using Gs22.

E. Moving Resolution Phantom
For rods with a size above 2.4 mm, performance is similar among all resolution models. For the off-center position, the 1.2 mm rods are slightly recovered using MI-SVPSF, MD-SVPSF, Gs16 and noRM, but not with Gs19 and Gs22. Deformation of the rods shape is visible in the off-center position, in the 2.4, 3.2 and 4 mm rods (visible as elongation) using the Gaussian kernels but not in reconstructions using the MI-SVPSF and MD-SVPSF.
The average eigenvalues ratio of the 2.4, 3.2, 4 and 4.8 mm rods is shown in Fig. 12 using the different resolution models,  for the center and off-center motion scans. For the centered phantom, rods eigenvalues ratios are lower than 1.12 in all cases, and the smallest and largest ranges are found using Gs19 (1.02-1.06) and MD-SVPSF (1.03-1.12) respectively. For the off-center phantom, the rod eigenvalues ratios substantially increase for all resolution models, except for MI-SVPSF and MD-SVPSF, for which only a modest increase or even a reduction in the 2.4 and 4.8 mm rods is observed. The smallest and largest difference in eigenvalues ratio between the centered and off-centered phantom rods is found using MD-SVPSF and Gs19 respectively, where the average eigenvalues ratio difference (considering all rods sizes) is 0.028 and 0.11.

F. Freely Moving Mice Brain Scans
The average mouse head speed during the test-retest and memantine challenge scans was 2.09 and 4.25 cm/s respectively. In test-retest scans animals moved sporadically, remaining for long time periods in a single position. In memantine challenge scans, due to the administration of the stimulant drug memantine, mice moved constantly during the entire scan, moving to regions close to the edge of the FOV. Fig. 13 shows the average motion corrected reconstructions of the mouse brain using the spatially invariant 1.6 mm Gaussian kernel (Gs16), the MI-SVPSF, and the MD-SVPSF, as well as the anesthesia motion-free reconstructions. Compared to the reconstructions using Gs16, reconstructions with MI-SVPSF and MD-SVPSF produce images with less spillover from hot regions of the brain to cold regions. Additionally, structures that were blurred in the motion corrected reconstructions using Gs16, such as the cortex in the test-retest scans and the hippocampus in the memantine scans (indicated with white arrows), were better recovered with the MD-SVPSF.
In both test-retest and memantine conditions, the Person's r correlation in regional brain quantification between awake motion corrected and anesthesia motion-free reconstructions increases when MI-SVPSF and MD-SVPSF is used compared to the use of Gs16. Diference between MI-DVPSF and MD-SVPSF is larger in memantine reconstructions compared to test-retest reconstructions. In the test-retest condition, r is 0.805, 0.830 and 0.836 (Fig. 14a, b, c), using Gs16, MI-SVPSF and MD-SVPSF respectively, while in the memantine condition r is 0.897, 0.938 and 0.964 (Fig. 14d, e, f).

IV. DISCUSSION
A method to calculate the spatially variant and motion dependent PSF in motion corrected reconstructions has been developed. The resolution kernels of the MD-SVPSF were calculated by transforming the scanner PSF according to the same transformation used in the rigid motion correction procedure to compensate the object motion. In this work we used the model proposed in [7] to represent the spatially variant PSF of the PET scanner. This model assumes an asymmetric and anisotropic Gaussian PSF. For the Inveon scanner this parameterization was a good approximation as shown in the comparison of the experimental measurement of a point source and the evaluation of the model at different radial distances. However, the developed method does not depend on the Gaussian model and it can be used in conjunction with other PSF models.
The robustness of the method in the presence of different motion patterns was evaluated using a moving phantom with 2 capillaries. The importance of the spatially variant correction was also tested by placing the 2 capillaries at a distance of 1 cm from each other. The method was compared with 3 spatially invariant Gaussian kernels with different widths, the spatially variant but motion independent kernel, and with the motion corrected reconstruction without resolution modeling. The first motion pattern confined the phantom to a Fig. 13. Average images of the awake mouse brain motion corrected (MC) reconstructions for the test-retest (top row) and memantine (bottom row) conditions, using a spatially invariant 1.6 mm (FWHM) Gaussian kernel, MI-SVPSF and MD-SVPSF. In addition, the anesthesia motion-free brain reconstructions using the 1.6 mm Gaussian kernel are shown for reference. All reconstructions are performed with 40 iterations. Image units are relative standard uptake value (SUV r ), normalized to the whole brain activity for visualization purposes.
small region in the FOV, similar to a motion-free case. In this case, the MI-SVPSF and MD-SVPSF gave similar image spatial resolution due to the small change of the PSF due to motion. For both cases the FWHM of the capillaries (1.53, 1.54 mm FWHM for MI and MD respectively) was close the capillaries diameter (1.5 mm) at 16 iterations. Using the Gaussian kernels, spatial resolution was improved proportional to the Gaussian width, showing undercompensation (1.74 mm) and overcompensation (1.42 mm) using the narrower and wider kernels respectively. In addition, the difference between left and right capillaries FWHM was the smallest using the MI-SVPSF and MD-SVPSF (lower than 9%), compared to the Gaussian kernels (between 13%-19%). The deformation of the capillaries circular shape, quantified with the eigenvalues ratio, was also the lowest using the spatially variant kernels.
The second motion pattern considered more extensive motion than pattern 1, spanning several radial positions, similar to the motion of an animal that moves on a horizontal platform. Contrary to motion pattern 1, MI-SVPF performed poorly in pattern 2 (1.69 mm FWHM) since, although the spatially variant nature of the PSF is considered, the motion dependence is not taken in consideration. Using the MD-SVPSF the motion dependence of the PSF is considered, resulting in an image spatial resolution similar to that of the motion pattern 1 (1.47 mm). Using the Gaussian kernels, again there is improvement of the spatial resolution in proportion to the kernel width (1.58, 1.42 and 1.32 mm FWHM from narrowest to widest Gaussian width) but with more severe overcompensation. For motion pattern 2 the difference between both capillaries FWHM is less prominent than in pattern 1. This can be explained by the fact that in motion pattern 2 both capillaries span approximately the same radial positions in the FOV, therefore undergoing similar loss of spatial resolution. The correction for shape deformation observed in motion pattern 1 (low eigenvalues ratio) using MI-SVPSF was not observed in pattern 2, in which even MI-SVPF produces the highest deformation among all resolution kernels. On the other hand, considering the motion dependence of the PSF, using MD-SVPSF produces the smallest deformation in the capillaries shape among all resolution models.
In motion pattern 3, the phantom was moved in a circular pattern in the transaxial plane to allow positioning the capillaries in similar radial positions, where the PSF width is also similar, but with different orientation. In this case MI-SVPSF gave the worst performance among all resolution models (1.84 mm FWHM). The MD-SVPSF performed better (1.61 mm FWHM) than the Gaussian kernel with the narrowest width (1.75 mm FWHM) but not as good as the Gaussian kernel with the widest width (1.48 mm). For this motion pattern the difference between capillaries FWHM was minimal (lower than 3%) in all cases, explained by the circular pattern motion which allows spanning of PSF's with different orientations for both capillaries. This effect also causes compensation of the capillaries shape deformation in all directions, producing minimal deformation in the capillaries shape (as quantified with the eigenvalues ratio) using all resolution models.
For all motion patterns, the use of resolution modeling greatly improved the spatial resolution of motion corrected images. In reconstructions without resolution modeling, the capillaries FWHM was larger than 2.3 mm in all cases.
In addition, deformation of the capillaries observed in reconstructions without resolution modeling was reduced using the MD-SVPSF for all motion patterns, but not for all cases using the Gaussian kernels.
In comparison with a motion correction reconstruction calculating the SVPSF resolution modeling for every pose (PP-SVPSF), the MD-SVPSF (motion dependent average SVPSF) performed similarly. The average difference in the capillaries FWHM between PP-SVPSF and MD-SVPSF was 4.6 %, while the difference in eigenvalues ratio was minimal (1.1% on average). However, motion correction reconstruction using PP-SVPSF is mostly practical to perform when few motion frames are present, such as in respiratory or cardiac motion [25]. For scans with non-repetitive motion, as in our case, the calculation time using PP-SVPSF becomes prohibitive for practical purposes due to the multiple image interpolation and blurring operations that need to be performed for every subject pose. For example, calculation time of the capillaries phantom scan 1 using PP-SVPSF, with 16 subsets and 16 iterations considering 1875 poses (1 min data), was 64 hours using parallelized code in a 10 core 3.3 GHz processor.
These experiments also show that using a spatially variant resolution model in motion corrected reconstructions is not recommendable if motion dependence is not considered in the resolution model. This is due to the fact that the reference pose in which the motion corrected reconstruction is performed, i.e. the pose that determines in which area in the image space the object will be reconstructed, can be chosen arbitrarily. In the spatially variant case, the area in the image space in which the motion corrected object is placed, can have a very different PSF than the area in which the object was actually measured. In our experiments, even considering the average pose as the reference pose, which in most cases reflect the most frequent pose of the object, using the motion independent PSF did not perform well (e.g. motion pattern 2 and 3) due to the difference between the average pose and the poses of the object during the scan.
In the moving resolution phantom experiment, the phantom was scanned at a center and off-center position to visualize the effects of the spatially variant PSF and to quantify the deformation of the phantom rods in the motion corrected reconstructions. In both positions, recovery of the rods of 1.6 mm diameter and larger was possible using all resolution models. However, using the widest Gaussian kernel, the cluster of 1.6 mm rods are visually more blurred compared to the rest of the resolution models. This can be due to slower convergence of PET reconstruction when applying resolution modeling [26], although the number of iterations was increased from 16 (capillaries experiment) to 40. In addition, edge artifacts are more pronounced using the widest Gaussian kernel compared to the other resolution models, as can be seen in the cluster of 1.2 mm rods.
The loss of resolution and deformation of the circular cross sections of the rods in the off-center position compared to the center position is visible using all resolution models. The eigenvalues ratio of the rods was higher (larger deformation) in the off-centered position compared to the centered position using all resolution models, except when using the MD-SVPSF which even produced a lower eigenvalues ratio in some rods in the off-center position. Moreover, the eigenvalues ratio difference between off-centered and centered positions is the smallest using the MD-SVPSF, which translates into more consistent shape of the object independently of its position in the FOV.
Finally, the use of the MD-SVPSF was demonstrated in scans of freely moving mice and were compared to motion correction using resolution modeling with a 1.6 mm FWHM Gaussian, and with the MI-SVPSF. The similarity with motion-free reconstruction using a 1.6 mm FWHM Gaussian was quantified. Since anesthetized (motion-free) mice were scanned at a position close to the CFOV, a 1.6 mm Gaussian kernel represent a good approximation to the PSF at the CFOV. Images using MD-SVPSF present more well-defined brain structures and less spill-over activity from hot regions to cold regions. In the [ 18 F]FDG brain test-retest scans, the cortex uptake is more uniform using the MI-SVPSF and MD-SVPSF compared to the use of the Gaussian kernel. In the memantine challenge scans the brain uptake distribution changed compared to that in the test-retest scans. This can be observed in the motion-free as well as in the motion corrected reconstructions. In memantine images, the increase uptake in hippocampus observed in the motion-free images was better recovered using the MD-SVPSF compared to the use of MI-SVPSF and the Gaussian kernel. In the memantine condition the difference between MI-SVPSF and MD-SVPSF was larger than in the test-retest condition. This can be explained by the higher motion speed and motion range of the animals in the memantine condition compared to the test-retest condition. When animals remain for long periods of time in the same position (test-retest), the MI-SVPSF is a good approximation, but when they move to different locations of the FOV (memantine), the MD-SVPSF is needed to model the scanner PSF with more accuracy. Notably, brain regional quantification was improved using the MD-SVPSF, in both test-retest and memantine challenge conditions. Pearson's r correlation between regional SUV quantification in motion corrected reconstructions and motion-free reconstructions increased from 0.805 (Gaussian kernel) to 0.836 (MD-SVPSF) in the test-retest condition and from 0.897 to 0.964 in the memantine condition.
In PET motion corrected reconstructions, resolution is not only determined by the motion dependent loss of resolution due to the spatially variant PSF. Other factors such as finite motion tracking sampling as well as tracking errors can affect the resolution in motion corrected images [13], [27]. After characterization of these effects, the blurring caused by these errors could additionally be introduced in the resolution model to further improve spatial resolution in motion corrected reconstructions.
The method presented here was developed and validated for brain scans of freely moving animals. However, it can also be useful for human brain rigid motion correction scans, especially when the patient motion is severe and accurate brain quantification in small regions is required. Moreover, although calculation of the MD-SVPSF was derived for rigid motion correction, it could also be calculated for non-rigid motion correction. Provided that one knows the non-rigid motion of every voxel in the reference image, e.g. using displacement fields, one can average the SVPSF over all motion frames at the original position of measurement to calculate the MD-SVPSF kernels.

V. CONCLUSIONS
A method to correct the loss of spatial resolution in PET rigid motion corrected reconstructions due to the spatially variant PSF was developed. Moving objects in PET scans can traverse areas in the FOV with different PSF's width and shape, thus being affected by a motion dependent PSF. Compared to the use of a spatially invariant Gaussian, the use of a motion dependent and spatially variant resolution model produced images with consistent spatial resolution as well as reduced deformation due to the parallax effect in phantoms undergoing different motion patterns. In in vivo experiments, mouse brain structures are better defined and less spill-over from hot regions to cold regions is observed using our method. As a conclusion, spatial resolution and quantification consistency in motion corrected reconstructions, independently of the subject motion characteristics, is improved using the proposed motion dependent spatially variant PSF.