Image-based modelling of residual blurring in motion corrected small animal PET imaging using motion dependent point spread functions

Optical motion tracking and motion compensation reconstruction algorithms enable the acquisition of quantitative measurements of brain function on conscious and freely moving rodents. However, motion corrected images often exhibit reduced resolution when compared with their stationary counterparts. This apparent loss of resolution can be attributed, among others, to jitter/noise in the measured motion estimates and brief periods of fast animal motion with insufficient motion sampling rate. In this paper we propose a novel methodology to experimentally characterise the residual blurring in the motion corrected images by measuring the motion-dependent point spread function (PSF) in image space using a point source rigidly attached on the moving object. We evaluated the proposed methodology using experimental phantom measurements acquired on the microPET Focus220 scanner. The motion dependent point spread function was extracted from the point source attached on the moving phantom, after motion correcting the images and modelling the point source in image space using an Expectation Maximisation algorithm as a weighted sum of two Gaussian distributions. Finally, the fitted blurring kernels were used within an iterative Lucy-Richardson algorithm to mitigate the deblurring in the motion corrected images. For motion typically encountered in an awake rat study, results showed that unprocessed motion corrected images suffer from lower resolution compared to a stationary acquisition. The shape of the measured blurring kernel, correlated well with the motion trajectory, while the width of the kernel was proportional to the speed/acceleration of the object. Post-processed images using the corresponding motion dependent blurring kernel appeared not only qualitatively, but also quantitatively (in terms of contrast) more similar to their stationary counterpart. We conclude that it is possible to experimentally measure the residual motion-dependent blurring kernel and use it within a post reconstruction deconvolution framework to improve resolution and quantification of motion corrected images.


Introduction
Dynamic small animal positron emission tomography (PET) imaging is a powerful research tool to gain unique insights into the physiology and molecular mechanisms underlying cognition and behaviour (Schulz and Vaska 2011). Conventionally, animals must be anaesthetized or restrained to avoid severe motion artefacts in the reconstructed images. However, the confounding effects of anaesthetic agents and stress induced by restraint are not always well understood and can alter or mask the physiological parameters of interest, such as brain metabolism (Spangler-Bickell et al 2016b), blood flow (Nakao et al 2001) and receptor binding (Alstrup and Smith 2013). Moreover, anaesthesia and physical restraint precludes the study of specific biological correlates of behaviour since the animals cannot (freely) move or respond to external stimuli (Rebec and Sun 2005).
For these reasons methods have been developed which allow brain function imaging while the animal is fully conscious and it is interacting with its environment. Such methods include the use of: a miniaturised PET detector ring which is surgically affixed to the animal's head , radioactive point sources attached to the animal's head, which are tracked and rigidly registered in image-space (Miranda et al 2017) and optical motion tracking of the animal's head using marker-based (Kyme et al 2011) or markerless  techniques. For those methods that involve any form of motion tracking, the extracted motion information is used to realign the measured lines of response (LORs) to where they would have been detected had the animal not moved at all (Rahmim et al 2004, Zhou et al 2013. Although the feasibility of reconstructing motion corrected brain images of a completely unrestrained rat using both marker-based and marker-less techniques has been demonstrated previously by our and other groups, the reconstructed images exhibit a noticeable reduction in contrast when compared to studies where the animal is completely stationary, i.e. under anaesthesia (Kyme et al 2012a). Nevertheless, exact knowledge and accurate modelling of animal motion within image reconstruction should, at least in theory, lead to improved resolution in the reconstructed images (Kennedy et al 2006). Several factors may be responsible for the resolution degradation observed with motion tracking and correction techniques, such as: (a) noise (jitter) in the measured motion parameters, (b) insufficient motion sampling (i.e. frame rate) relative to animal speed, (c) mis-synchronisation between the tracker and the emission data, (d) independent marker (or point) movement relative to the animal's brain and (e) systematic errors in trackerscanner cross calibration.
The impact of noise on the measured motion parameters can usually be minimised by temporally smoothing the estimated poses using a finite impulse response (Kyme et al 2011) or a polynomial (Kim et al 2015) pose filter. On the other hand, motion tracking techniques assume that animal motion within a recorded pose frame is negligible. This assumption (which directly depends on the animal's speed) can be strengthened by increasing the frequency of pose sampling and/or by incorporating an interpolation approach to correct the events between two recorded poses (Kim et al 2015, Spangler-Bickell et al 2016a. However, current optical tracking hardware may be limited in terms of sampling speed. For example, the Micron Tracker Sx60 has a maximum frame rate up to 48 Hz for a single tracker and up to 26 Hz when more than one trackers are running simultaneously on the same firewire bus. In addition, interpolation is based on the assumption that the observed motion follows the properties of the interpolant, which is not always easy to predict, particularly for awake and freely moving animal studies. Finally, synchronisation between motion and emission data streams can be refined by using unique, non-periodical wave-patterns to align the data, as well as experimentally measuring temporal fluctuations in the timings of the acquisition process (Spangler-Bickell et al 2016a). Despite all these improvements the contrast in motion corrected images appears to be worse compared to stationary ones, which indicates a residual (or systematic) blurring.
The nature of this residual blurring, caused by object motion (after correcting it using the measured motion parameters), can be broadly considered similar to the blurring caused by the system response function (Alessio et al 2006). Therefore, in the same way that one can characterise the geometric point spread function (PSF) either by experimental measurements (e.g Panin et al 2006) or analytic modelling (Qi et al 1998) and mitigate the blurring caused, for example, by parallax errors, via iterative deconvolution between the reconstructed images and the measured PSF (Rahmim et al 2013), we hypothesise that we can experimentally measure and characterise the motion dependent PSF and mitigate the residual blurring effects caused by imperfections in animal motion measurements. The motion dependent PSF has been previously analytically estimated using the acquired motion measurements (Miranda et al 2014). However, the estimated motion dependent PSFs may be vulnerable to systematic errors in the measured motion parameters, which may ultimately propagate to the deblurred images. Instead, in this work we propose to experimentally measure the study-specific, motiondependent PSF, in a similar way as the geometric PSF function; that is, by attaching a small point source (ideally with diameter smaller than the system resolution) either on the optical marker used for motion tracking or directly on the object being imaged. The advantage of the proposed image-based approach is that it inherently includes all factors that contribute to the apparent residual blurring and which cannot be easily modelled. Once estimated, the fitted motiondependent PSF can then be used within an iterative Lucy-Richardson deconvolution framework (Lucy 1974) to mitigate the residual blurring caused by the object motion.
Iterative deconvolution methods based on the Lucy-Richardson algorithm have been previously widely used in image processing to account for known (Khan et al 2013) or unknown (Fish et al 1995) imperfections of the imaging and acquisition process, as well as motion, that result in image blur. Likewise, within a PET context Lucy-Richardson deconvolution has been successfully applied on motion contaminated images to estimate the original non-corrupted image, mainly addressing intra-and inter-frame motion in dynamic brain imaging (Menke et al 1996, MohyudDin et al 2015. It has also been used to account for intraframe motion in whole-body imaging due to respiratory (Naqa et al 2006) and bulk patient motion (Karakatsanis and Tsoumpas 2017). Contrary to these publications, which aim to account for the overall motion blurring in the reconstructed images, in this work we use a post-reconstruction iterative Lucy-Richardson deconvolution algorithm to account for the residual motion blurring, after correcting the measured emission data, in an event-by-event basis, using motion tracking data.
Therefore, the overall aim of the present study was to investigate whether we can experimentally measure the motion-dependent PSF in image space and mitigate the residual blurring effects caused by imperfections of the motion tracking system. Since the residual blurring depends mostly on the observed motion rather than the object being imaged, we designed and performed two sets of controlled phantom experiments on the microPET Focus220 scanner, emulating simple, as well as realistic rat motion patterns. We demonstrate that the shape and width of the measured motion PSF is dependent on the observed motion pattern and the speed/acceleration of the object, respectively. Finally, we use the estimated motiondependent PSF within a post-reconstruction deconvolution framework to de-blur the motion corrected images and evaluate the improvements in attained resolution compared to the stationary counterpart.

Proposed methodology
Many studies have previously demonstrated the benefits of image-based modelling of those physical effects that degrade spatial resolution of the reconstructed images, such as positron range and parallax errors (Sureau et al 2008, Angelis et al 2015, Bickell et al 2016c. One way to characterise the overall spatially variant blurring caused by all resolution degrading effects is to experimentally measure the system response function in image space, using collimated or uncollimated radioactive point sources placed at multiple locations within the field of view (FOV) of the scanner (Rapisarda et al 2010, Kotasidis et al 2011, Rahmim et al 2013. The fitted and/or parameterised PSFs are then used during or after image reconstruction to improve contrast recovery in the reconstructed images (Rahmim et al 2013).
Along these lines, we propose a methodology (figure 1) to characterise the study dependent residual blurring in the motion corrected images using a small radioactive point source rigidly attached on the object being imaged. As such, the acquired list-mode emission data are reconstructed using an iterative maximum likelihood expectation maximisation (MLEM) motion compensation reconstruction algorithm (Rahmim et al 2008), without accounting for geometric or any other resolution degrading effects. After reconstruction the attached radioactive point source is isolated from the rest of the image and the study-dependent PSF is extracted by setting the background to zero (<1% of max value) and ensuring that the entire PSF comfortably fits within the isolated kernel. The point source is then fitted using an Expectation Maximisation (EM) algorithm to a Gaussian mixture model with two anisotropic components aiming to capture both the geometric and the residual motion blurring effects. In order to minimise the dependency on the initialisation parameters for the EM algorithm, 10 independent realisations were performed and the parameters that best fit the data (among the 10 realisations) were used to define the resolution kernel: where, j and l are the voxel linear indices for the image space and the blurring kernel respectively, m i and s i hold the mean and covariance for the radial, tangential and axial responses for the ith Gaussian distribution and w i is the mixing ratio of the two distributions. Once calculated, the fitted motion-dependent PSF is used within an iterative Lucy-Richardson algorithm (Lucy 1974) to de-blur the residual motion artefacts from the reconstructed images: where, are the reconstructed and deblurred images, respectively, with N number of voxels,  = Î{ } q q jl N N is the blurring kernel and k is the iteration number. The Lucy-Richardson deconvolution algorithm can be applied either on the final motion-corrected reconstructed image or in a nested approach within image reconstruction similar to Angelis et al (2013). In this study, the blurring kernel q jl was spatially invariant across all voxels in the reconstructed FOV, although a spatially variant kernel may be more appropriate (see discussion).

Phantom preparation
We carefully selected three molecular sieves (alkali metal alumino silicate) with almost spherical shape and diameter around 1 mm. The molecular sieves were soaked in a solution of 18 F and left to dry for a few minutes. Once dry, the molecular sieve with absorbed activity close to 200 kBq was selected and securely placed using a coating of super glue within a custom printed cubic enclosure (16 mm×16 mm×16 mm) (figure 2(a)), considered to be approximately equivalent to tissue density. The purpose of this printed enclosure was two-fold: to facilitate the secure mounting of the molecular sieves on the outer surface of the phantom away from the regions of interest and to provide a tissue equivalent material for positrons to annihilate. During an actual animal study the molecular sieve would be placed inside the plastic or perspex marker substrate and would be firmly glued on the animal's forehead, close to its brain (figure 2(c)). This marker enclosure provides not only adequate material for positrons to annihilate (thickness 4 mm), but it also it moves in synchrony with the animal's head and the marker.
In order to assess any improvements in contrast recovery for the deblurred motion corrected images using the proposed approach, we applied a range of robotically controlled motion patterns to a Micro-Deluxe hot-rod phantom (Data Spectrum Corporation). The phantom, which consists of six groups of hollow rods with internal diameters 4.2, 4.0, 3.2, 2.4, 1.6 and 1.2 mm, was filled with 65 MBq 18 F and scanned in the microPET Focus220 scanner (Preclinical Solutions, Siemens Healthcare Molecular Imaging, Knoxville, TN, USA) (Tai et al 2005). A composite marker consisting of two independent markers: a small 3-point marker (12 mm×8 mm) similar to those used for our awake animal studies and a large 4-point marker (36 mm×32 mm) was affixed to the phantom (figure 2(b)) to continuously track its pose. The larger marker was used to assess improvements in resolution due to the size of the marker, but results are not presented in this paper, since only the small marker relevant to small animal experiments was used.

Motion patterns
One of the aims of the present study was to demonstrate that the shape of the measured residual motion PSF is dependent on the trajectory of the motion and the width of the PSF is related to the speed and acceleration of the animal or object. For this purpose we designed two simple, robotically controlled, inplane motion patterns. The first pattern was a lateralonly movement, in which the phantom was allowed to Figure 1. Proposed image-based motion-dependent PSF deconvolution method: the emission data are first motion corrected and reconstructed into images. Then the motion-dependent PSF is extracted from the attached point source and is then fitted using a mixture of two anisotropic Gaussian distributions. Finally, the fitted PSF is used within an iterative Lucy-Richardson deconvolution step to deblur the motion-corrected reconstructed images. The hot-rod phantom used for all the experiments described in this paper (outer diameter 50 mm and length 96 mm). The molecular sieve (placed within the cubic enclosure) was securely attached on the phantom, close to the area of interest, while a large 4-point marker was rigidly attached on its front face. (c) An example of a typical small animal marker substrate (22 mm×27 mm) with a socket for the molecular sieve in the centre. The 3-point marker is stuck on the top preventing the rat from touching the radioactive point source. move radially from left to right only (figure 3(a)). The second was a star-like motion where the phantom was allowed to move in all directions within the same axial plane ( figure 3(b)). Both motion patterns involved only translations and they were specifically chosen to induce a distinct shape on the measured motiondependent PSF. For each of these motion patterns, 3 speed/acceleration settings were implemented (see table 1), which for this study were purposely selected to be faster than in routine practice. It is worth noting that in a typical awake rat study the average linear speed of the animal's head is around 15 mm/s, while animals may infrequently exhibit brief, rapid movements with maximum linear velocities up to 150-200 mm/s (Kyme et al 2012b). Therefore, including a stationary acquisition, which served as the gold standard, we performed 7 scans (6 with a motion pattern), each with a duration of 12 min.
Although these simple, in-plane motion patterns are convenient to demonstrate the dependencies of the measured PSF, they do not represent a typical animal motion. Therefore, the proposed methodology was also assessed in a more realistic motion scenario. Head movements of a previously tracked awake tube-bound rat were converted to robot coordinates and applied to the hot rod phantom using a high precision 6-dof robot. The simulated motion parameters included 3 translations and 3 rotations (figures 4(a) and (b)), with an average speed of 8 mm/s and maximum speeds about 60-80 mm/s across all three axes (figure 4(c)). To get an idea of the expected sampling error (and therefore an estimate of the blurring) due to sampling frequency for this particular animal study, the sampling error, in terms of average absolute distance travelled between two successive poses, for three sampling frequencies is shown in figure 4(d). These frequencies represent the sampling range of the motion tracking device used for this study, as well as the frequency that was actually used (31 Hz). In addition, to demonstrate that the behaviour of the animal can have an impact on the width of the measured motion dependent PSF, we increased the speed and acceleration of the robot to emulate a more agitated animal with an average head speed of 12 mm/s. Along with a completely stationary scan, three scans were performed, the duration of which was 20 min.
During data acquisition the phantom was continuously moved according to the predefined motion patterns using a high-precision 6-axis robot (Epson C3-A601S, Seiko Corp. Japan). The pose of the small marker attached to the phantom (which is routinely used for our awake animal experiments) was tracked using an optical tracking system (MicronTracker Sx60, ClaroNav, Toronto Canada) every 32 ms (i.e. 31 Hz) (Kyme et al 2011).

Image reconstruction and deconvolution
All scans were acquired in list mode and the emission data were reconstructed with an ordered subsets expectation maximisation (OSEM) list mode motion compensation reconstruction algorithm (Rahmim et al 2008). The acquired emission events were precorrected for the measured motion on an event-byevent basis (Zhou et al 2008) using the motion information obtained from the optical motion tracker. Corrections for normalisation and attenuation were included within the calculated time-averaged system matrix, in a factorised manner, similar to Angelis et al (2014). Corrections for scattered or random events were not applied, while reconstructed images were  corrected for isotope decay based on the starting time of each scan. All images were reconstructed in a fine image voxel grid 256×256×95, with voxel dimensions 0.47450×0.47450×0.796 mm 3 to ensure adequate sampling of the PSF. For each reconstructed image the point source that corresponds to the study-dependent PSF, was manually isolated from the rest of the imaging data and a threshold was applied to zero those pixels with intensity less than 1% of the maximum. As previously mentioned (see section 2.1), the extracted point source was modelled as a mixture of two Gaussian distributions, where the mean for each of the two Gaussian functions was a free parameter to account for potential asymmetry in the study-dependent PSF (either due to motion or geometry). The use of two Gaussian distributions, potentially with quite different kernel widths enables non-Gaussian distributions to be modelled. Finally, the fitted blurring kernel was used within a post-reconstruction iterative Lucy-Richardson deconvolution algorithm (Lucy 1974) to de-blur the motion corrected reconstructed images. For each study 4 iterations of the deconvolution algorithm were performed. Figure 5 shows the motion corrected reconstructed images for the three lateral motion schemes (top row, (b)-(d)) and the star-like motion pattern (bottom row, (e)-(g)), along with the stationary image as a reference (far left, a). For speed or acceleration typically encountered during an awake rat study (15 mm/s) images appear to be marginally degraded in terms of resolution (figures 5(b), (e)) when compared to their stationary counterpart ( figure 5(a)). However, residual motion blurring is clearly visible in the motion corrected images for moderate (figures 5(c), (f)) and fast motion speed (figures 5(d), (g)). It is interesting to note that although the speed/acceleration settings are very similar between the lateral and the star-like patterns, images corrupted with lateral motion appear visually less appealing. This is because resolution is degraded only along the direction of motion. In other words, the resolution loss is highly anisotropic for the lateral-only motion pattern, which makes the resolution degradation more obvious to the eye. On the other hand, the star-like motion pattern degrades image resolution by almost the same magnitude in all directions, which results in an image with reduced resolution, but no apparent artefacts. The isocontours of the fitted motion dependent PSFs, which were extracted from the molecular sieve attached on the moving phantom and encode the residual blurring after motion correction, are shown in figure 6. The isocontours for the lateral-only motion (top row, (b)-(d)) and the star-like motion (bottom row, (e)-(g)) clearly illustrate the dependencies of the shape of the measured PSF on the motion pattern and speed settings. As hypothesised, the shape of the PSF strongly depends on the nature of motion, with the lateral motion pattern leading to more elongated PSFs along the direction of motion (figures 6(b)-(d)). On the other hand, the star-like motion leads to more rounded PSFs, although there is still a strong diagonal component. In addition, the width of the distribution appears to become broader, for both motion patterns, as the speed of the phantom increases.

Results
The observation that for the lateral-only pattern the resolution is degraded only along the direction of motion is also supported by the line profiles through the fitted point source shown on the top row of figure 7. The horizontal profile (a) is along the x-axis, which is the same as the direction of motion for the lateral-only pattern, while the vertical profile (b) is  perpendicular to the direction of motion. As illustrated by the profiles, there is a clear deterioration of the FWHM only along the direction of motion, while there is almost no change in FWHM along the perpendicular direction. On the contrary, the degradation of the FWHM for the star-like motion, appears to be similar for both the horizontal and vertical directions, indicating a more isotropic loss of resolution.
The same observations are quantified on the plots shown in figure 8, where the radial and tangential FWHM of the first (figure 8(a)) and second ( figure 8(b)) Gaussian components are plotted against  the speed of the phantom. For the lateral-only case, which is represented by the solid lines, the FWHM of the first radial component increases almost linearly from 1.8 mm to almost 5 mm, while the first tangential component, which is perpendicular to the direction of motion stays below 2 mm. On the other hand, for the star-like motion, which is represented by the dashed lines, the FWHM for both radial and tangential components increases from 1.8 mm to 3.8 mm and 3.3 mm, respectively, since motion is multi-directional. The radial and tangential FWHM of the second Gaussian component for the lateral-only motion appear to be almost unaffected by the speed of motion ( figure 8(b)).
The estimated residual motion blurring kernels for each motion pattern and speed setting (figure 6) were used within an iterative deconvolution framework to deblur the motion corrected images post reconstruction. Figure 9 shows the unprocessed (first column) and de-blurred images (second column) for the stationary case (first row), the typical animal speed (second row) and moderate speed settings (third row) for the lateral-only motion pattern. As also shown in figure 5, the unprocessed motion corrected reconstructed images suffer from lower resolution due to the residual motion blurring. This loss of resolution appears to be marginal for the low speed setting (15 mm/s), while it is notably worse for the moderate speed setting (36 mm/s). In addition, due to the nature of the motion, some rods appear more elongated along the direction of motion (x-axis), especially for the moderate speed. However, the deblurred images appear to be qualitatively very similar, in terms of resolution and contrast, to their stationary counterpart. In particular, the 3 larger sets of rods (i.e. 4.2, 4.0, 3.2 mm) in the deblurred images appear more circular compared to the unprocessed images. As a note, the deblurring of the stationary acquisition is equivalent to conventional image-based resolution modelling using the experimentally measured spatially invariant blurring kernel, which leads (as expected) to improved contrast and potentially to unwanted Gibbs artefacts.
The line profiles through the 3.2 mm rods (third column) and 4.8 mm (fourth column) rods, shown in figure 9, quantify the differences in contrast before and after deblurring. The line profiles show a clear difference in contrast (peak-to-valley ratio) between the unprocessed stationary and motion corrected images (dashed curves), especially for the moderate speed (36 mm/s). After processing the reconstructed images for residual blurring line profiles between the stationary and motion corrected images (solid curves) appear very similar. Although line profiles between the stationary and motion corrected images indicate a better match compared to their unprocessed Figure 9. Lateral-only motion pattern. First column: transaxial slices of the motion corrected hot-rod phantom without postreconstruction processing, for no motion (top row), typical animal motion speed (middle row) and moderate speed (bottom row). Second column: the same transaxial slices, after de-blurring the motion corrected images using the corresponding estimated motion dependent PSF. Third and fourth columns: line profiles through the 3.2 mm and 4.8 mm rods, respectively, of the unprocessed (dotted lines) and the deblurred images (solid lines) for a given motion speed, compared against their stationary counterpart. counterparts, deblurred images may suffer from overshooting, which is common with deconvolution algorithms.
Similar images and line profiles are shown in figure 10 for the in-plane star-like motion. A clear deterioration in contrast can be seen for the unprocessed motion corrected reconstructed images (figures 10(b), (c)) when compared to the stationary image ( figure 10(a)). However, after post-processing the motion corrected images for residual motion blurring with the corresponding estimated blurring kernel, images appear not only sharper, but also qualitatively more similar to the stationary one. Line profiles shown on the right-hand side quantify the improvements of the proposed deblurring approach, with the deblurred motion corrected images being quantitatively very similar to their stationary counterparts in terms of contrast recovery. Figure 11 shows the motion corrected reconstructed images for the realistic rat motion, along with the isocontours of the estimated motion-dependent PSFs. In accordance with the results shown earlier, the motion corrected images appear to be marginally degraded, in terms of resolution, compared to the stationary acquisition, with resolution deterioration being worse for the fast moving case ( figure 11(c)). The isocontours of the fitted point source attached to the surface of the phantom clearly illustrate the dependency of the residual motion blurring on the speed of the object. Finally, the unprocessed and deblurred images for the realistic rat motion experiment are shown in figure 12, which demonstrate the improvements in contrast recovery offered by the post processing with the proposed deblurring methodology.

Discussion
Awake animal PET imaging is an important development towards the integration of functional neuroimaging and behavioural research. However, the statistical uncertainty introduced by factors such as noise/jitter and brief periods of insufficient sampling in existing motion tracking methods, can reduce the resolution of the motion corrected images. Such reduction in resolution may adversely affect quantification in small brain structures in rodents and therefore accounting for this residual motion blurring may be necessary. In this proof-of-principle study we presented a novel methodology to experimentally measure the residual blurring after motion correction and use the estimated motion-dependent blurring kernel within a post reconstruction Lucy-Richardson deconvolution algorithm to improve the contrast recovery in the motion corrected images. Figure 10. Star-like motion pattern. First column: transaxial slices of the motion corrected hot-rod phantom without postreconstruction processing, for no motion (top row), typical animal motion speed (middle row) and moderate speed (bottom row). Second column: the same transaxial slices, after de-blurring the motion corrected images using the corresponding estimated motion dependent PSF. Third and fourth columns: line profiles through the 3.2 mm and 4.8 mm rods, respectively, of the unprocessed (dotted lines) and the deblurred images (solid lines) for a given motion speed, compared against their stationary counterpart. Figure 11. Top row: transaxial slices through the motion corrected images of the hot rod phantom for the realistic rat motion experiment. Top row: reconstructed images without post-processing to account for motion blurring for the stationary scan (a), for typical animal motion, and (c) for a slightly agitated animal. Bottom row: isocontours of the fitted motion dependent PSFs superimposed on the raw reconstructed data from a single molecular sieve attached to the outer surface of the phantom. Figure 12. Realistic rat motion pattern. First column: transaxial slices of the motion corrected hot-rod phantom without postreconstruction processing, for no motion (top row), typical animal motion speed (middle row) and moderate speed (bottom row). Second column: the same transaxial slices, after de-blurring the motion corrected images using the corresponding estimated motion dependent PSF. Third and fourth columns: line profiles through the 3.2 mm and 4.8 mm rods, respectively, of the unprocessed (dotted lines) and the deblurred images (solid lines) for a given motion speed, compared against their stationary counterpart.
Contrary to a previously published method (Miranda et al 2014), which generates the motion dependent blurring kernels by estimating the displacement uncertainty in the measured motion parameters, the proposed method does not depend directly on the measured tracking data. Instead, it captures the blurring response as it is expressed in the final reconstructed image, in the same way as the geometric system response can be measured entirely in image space (Rahmim et al 2013). Although the proposed approach estimates an approximation of the motion dependent PSF (as opposed to an analytic derivation), there are numerous advantages to calculating the blurring kernel in image space. Firstly, the measured motion parameters are highly correlated and there may be a systematic bias (e.g. due to errors in trackerscanner cross calibration) that is challenging to identify and/or model in the estimated displacement uncertainty, which in turn may propagate to the deblurred image. Instead, the residual blurring in image space inherently includes all aspects of the motion tracking procedure that contribute to uncertainty in the realignment of the lines of response.
Secondly, although this paper focused on the motion component of the blurring kernel, the estimated motion dependent PSF also inherently includes a geometric component which accounts for resolution degrading effects (e.g. parallax error) as it would in a completely stationary scan. This observation is quite important, since the blurring in the final motion-corrected reconstructed images is straightforwardly treated as a single entity regardless of where this blurring originated from. However, since the blurring kernel is estimated using a single point source attached on the object being imaged, it means that the same spatiallyinvariant kernel will be applied to all voxels. Such an approach is certainly a crude approximation due to the spatially variant geometric response of most commercial PET scanners, as well as due to the rotational motion of the object, which amplifies motion away from the pivot point. Nevertheless, given that the size of a rodent's brain is rather small (i.e. approximately 25 mm for an adult rat), we assume that the geometric component should be effectively invariant (or marginally different) around the vicinity of the point source attached close to the rodent's brain. The same can be assumed for the voxel-wise motion parameters, where differences due to rotational motion are expected to be very small in the close neighbourhood of the rodent's brain. On the other hand, a major advantage of using a shift-invariant blurring kernel, particularly when the deconvolution is interleaved with tomographic iterations (Angelis et al 2015), is the ability to perform the convolution in the Fourier domain and speed up the overall reconstruction process.
Most importantly, the experimentally measured blurring kernel directly includes the time-varying geometric component due to the object moving to different locations within the FOV with substantially different resolution. This is particularly important for a freely moving animal scan, where the animal is free to roam within a small enclosure covering almost the entire radial extent of the scanner (Zhou et al 2013). The estimated motion dependent PSF is actually a time-weighted average of all individual PSFs along the trajectory of the animal's motion. This observation is valid, given the assumption that the blurring kernel, measured at the location of the point source, does not substantially change across the rest of the animal's brain. One way to overcome this approximation is to extract the motion-only component from the measured blurring kernel, by deconvolving the overall motion dependent kernel from a stationary blurring kernel. Nevertheless, this would require a second acquisition to obtain the stationary blurring kernel, as well as knowledge of the spatially variant nature of the PSF across all voxels in the FOV.
Deblurred images from experimental phantom data appear to be improved in terms of contrast recovery when compared to their unprocessed counterparts and also more similar to the stationary image. In addition, for a highly anisotropic residual motion blurring, such as that caused by the lateral-only pattern, the proposed deblurring methodology substantially improved the shape of the cylindrical rods, which appeared elongated along the direction of motion before processing (figure 9). However, since severe motion led to a rather broad blurring kernel, it did not provide any benefit in resolving the smaller rods in the phantom. Although these small rods challenge the intrinsic resolution of the scanner, lack of improvements in terms of resolution may also be due to the fact that the molecular sieve was placed relatively far from the region of the small rods (approximately 55 mm). As seen in figure 5(a) the molecular sieve was placed at a large radial distance, compared to the 1.6 mm rods, which were closer to the centre of the scanner, leading to an overall broader kernel. The use of a marginally overestimated motion-dependent blurring kernel leads to sub-optimal resolution recovery. However, in an actual rat experiment, ideally the point source will be placed very close to the brain (i.e. less than 20 mm from any part of the brain), leading to a more optimal resolution kernel and maximising the benefits of resolution modelling for all voxels. On the other hand, the close proximity of the point source to the object or animal requires careful consideration of the absorbed activity in the molecular sieve to avoid star artefacts or scatter contamination for substantially high and low activity in the sieve, respectively.
In this work we assumed that the motion dependent blurring kernel can be sufficiently characterised using a weighted sum of two Gaussian distributions. This assumption was based on the fact that the motion blurring kernel attempts to describe all resolution degrading effects, such as the geometric system response function and the residual blurring induced by the optical motion tracking. Previous studies (Sureau et al 2008, Kotasidis et al 2011 have demonstrated that a sum of two Gaussian distributions can adequately model the shape of the geometric PSF. Therefore, since the overall motion dependent PSF is a blurred version of the stationary geometric PSF (or a convolution between the stationary PSF and the motion uncertainty PSF) we assumed that the same two Gaussian model would also be sufficient to describe the motion dependent blurring kernel. According to the results presented in figure 8, the first Gaussian component (i.e. larger weighted contribution) was enough to describe the residual blurring caused by increasing speed of the phantom. On the contrary, the second Gaussian component (i.e. smaller weighted contribution) was almost independent of the phantom speed. The latter may indicate that a single Gaussian function may be enough to describe the motion-dependent blurring. However, in order to allow for non-Gaussian distributions to be modelled the proposed sum of Gaussian distributions with independent mean values may be more appropriate.
As already discussed, the image-based PSF is study-dependent and therefore a preliminary reconstruction is necessary to obtain the blurring kernel for the deconvolution step. Since Lucy-Richardson is a very efficient algorithm compared to image reconstruction, in order to avoid doubling the overall reconstruction time, we opted for a post-reconstruction deconvolution approach. However, the proposed methodology can be straightforwardly adapted to a nested deconvolution within tomographic image reconstruction similar to (Angelis et al 2015), potentially further improving image quality. In this case, however, in order to minimise reconstruction time and possibly the dependencies of the PSF on convergence, a preliminary reconstruction could be done with filtered back projection.
As with any other deconvolution approach the proposed method may be prone to overshooting, ringing artefacts and amplification of noise (Rahmim et al 2013), all of which may challenge the quantification of the deblurred images. The ringing artefacts are the result of the mismatch between the true and the estimated blurring kernel used within the Lucy-Richardson deconvolution algorithm and it becomes worse for increased iterations. Such overshoot, which was visible in the deblurred images, can be suppressed by using a slightly underestimated blurring kernel (Tong et al 2011). In this study, we opted for a relatively low number (4) of deconvolution iterations in order to achieve improvements in contrast, without amplifying overshooting artefacts and/or noise.
The proposed approach is easily adaptable to awake and freely moving animals. The molecular sieve can be placed within the printed marker substrate used to track the animal's pose (figure 2(c)), while it provides not only a material for positrons to annihilate, but also a simple way to rigidly attach the molecular sieve on the animal's head. In addition, this arrangement satisfies the assumption that the point source used to measure the motion-dependent PSF is very close to the region of interest to minimise errors due to the spatially variant nature of the resolution and the rotational motion of the animal. Although no animal scans were presented in this proof-of-principle study, we demonstrated the validity of the proposed method using realistic animal motion from a previous tube-bound rat study. Clearly, the performance of the proposed method depends on the measured motion, rather than the object being imaged. Therefore, the use of realistic animal motion to introduce motion effects on a resolution phantom, allowed for fair and controlled assessment of the proposed methodology.
As a final note, the computational time of the overall methodology depends among other things on factors, such as: (a) the specifications of the image processing workstation, (b) the deconvolution approach (FFT versus direct convolution), (c) the number of voxels in the reconstructed image and (d) the efficiency of the actual implementation (e.g. single versus multi-threaded). In addition, a computational bottleneck in the current methodology is the manual extraction of the point source from the reconstructed image, which adds a substantial overhead in the overall processing time. Despite this bottleneck, which can be fully automated in the future, the total post-reconstruction processing time is less than the actual image reconstruction.
In the future, we aim to further assess the proposed methodology in a nested deconvolution approach, within image reconstruction. Although such an approach will double the overall reconstruction time, we hypothesise that image quality will be further improved due to better handling the noise within the reconstruction algorithm. In addition, we currently investigate the option to include additional point sources around the object being imaged or the animal's head, in a similar way as in (Miranda et al 2017), to further improve the estimated blurring kernel for all the voxels in the reconstructed volume (due to rotational motion and geometric variability). Finally, we aim to apply the proposed deblurring approach to a small cohort of behaving rats to assess and quantify improvements in contrast recovery and resolution, for behaviour and motion patterns induced by different pharmacological interventions (e.g. amphetamine).

Conclusions
We conclude that it is possible to experimentally measure an estimate of the study-dependent residual motion blurring kernel using a small point source attached directly to the object being imaged. This blurring kernel can provide information about the uncertainty levels of motion tracking and can be used within a post reconstruction Lucy-Richardson deconvolution framework to improve resolution in motion corrected images. However, care has to be taken to avoid Gibbs artefacts and overshoot effects associated with deconvolution algorithms.