Comment on: ‘A Poisson resampling method for simulating reduced counts in nuclear medicine images’

In order to be able to calculate half-count images from already acquired data, White and Lawson published their method based on Poisson resampling. They verified their method experimentally by measurements with a Co-57 flood source. In this comment their results are reproduced and confirmed by a direct numerical simulation in Matlab. Not only Poisson resampling, but also two direct redrawing methods were investigated. Redrawing methods were based on a Poisson and a Gaussian distribution. Mean, standard deviation, skewness and excess kurtosis half-count/full-count ratios were determined for all methods, and compared to the theoretical values for a Poisson distribution. Statistical parameters showed the same behavior as in the original note and showed the superiority of the Poisson resampling method. Rounding off before saving of the half count image had a severe impact on counting statistics for counts below 100. Only Poisson resampling was not affected by this, while Gaussian redrawing was less affected by it than Poisson redrawing. Poisson resampling is the method of choice, when simulating half-count (or less) images from full-count images. It simulates correctly the statistical properties, also in the case of rounding off of the images.


Introduction
White and Lawson (2015) described a Poisson resampling scheme for simulating reduced counts in nuclear medicine images. The scheme was validated on images of a Co-57 flood source on a Skylight SPECT scanner (Philips Healthcare, Best, the Netherlands). It was compared with a scheme that was based on direct redrawing from a Poisson distribution, which showed severe problems regarding pixel values and noise and skewness below 10 counts per pixel on the original image. There are two mechanisms that can cause deviations in the process. The first one is the use of a statistical distribution in some way on the original data, which may result in statistical parameters that deviate from the theoretical correct values. The other mechanism is caused by the fact that the images are saved in DICOM (digital imaging and communication in medicine) format as integers, which may result in deviations of the statistical parameters due to the rounding off process. White & Lawson (2015) investigated mean m, standard deviation σ and skewness γ 1 for the sampled images. In this comment excess kurtosis γ 2 is included in this set of statistical parameters. For a Poisson distribution with parameter λ these are given by (Spiegel 1961, p 124): It is obvious that just dividing image by a factor of two is not sufficient, since that underestimates the standard deviation (noise). For the redrawing process two possibilities were considered. One, where samples were redrawn from a Poisson distribution, where parameter λ was chosen equal to the original amount of counts, and from a Gaussian distribution with mean and standard deviation equal to the original count data (assuming Poisson statistics). In this way, the resulting (relative) amount of noise in the half count images is correct. Dividing by 2 results in the proper mean and absolute standard deviation.
The process of saving as integers is simulated by rounding off the obtained values to the nearest integer, where 0.5 is rounded up.
The Poisson resampling scheme as described by White and Lawson (2015) results in integer values, so no rounding off is needed. Reversing the order of redrawing and rounding off does not lead to the correct amount of noise as shown by White and Lawson (2015).

Methods
A computer simulation was performed in Matlab R2012a (Mathworks, Natick, Massachusetts, USA), where large array was sampled from a Poisson distribution with parameter λ. A halfcount array was calculated by (1) redrawing from a Poisson distribution based on the full-count value, (2) redrawing from a Gauss distribution with standard deviation and mean based on the full-count value and assuming Poisson statistics, and (3) Poisson resampling as described in White and Lawson (2015). The crucial lines in the Matlab code for these calculations are shown in figure 1.
In order to verify the statistical relationships in equation (1) the half-count/full-count ratio of these statistical parameters is calculated from the simulations. Since excess kurtosis and skewness might be close to zero the half-count full-count ratios for the statistical parameters were determined by calculating the mean value for these separately over many repetitions and dividing the two mean values. Calculating the averaged ratio directly from the simulations may result in a numerically unstable solution, since the denominator can have a value close to zero, and a bias, since the ratio and the reciprocal of its inverse do not have the same expectation value (Moreno 1996, section 5 Relative statistical uncertainty was approximated as the square root of the summed relative variances (van Kempen 2000, equation (7) with zero covariance). The 5% double sided confidence interval was one digit or less wide for the last shown digit. For instance for a precision of 3 digits the maximal half confidence interval width was 0.000 5, equal to approximately 1.96 times the standard deviation. Simulations were run with and without rounding off. For skewness and excess kurtosis the ratio was also determined by using the theoretical full-count values of these parameters in the denominator in order to reduce the uncertainty in the calculated ratio.

Results
Without rounding off, the ratios did not depend on the number on counts in the range 0,1-100 (verified for 0.1, 1, 10, 100 and 1 000). All methods show the proper mean and standard deviation. However, skewness was not correct for the Poisson redrawing method. Excess kurtosis was only corrected for the Poisson resampling method. Results are shown in table 1 (calculated for a full-count pixel value of 1).
Results with rounding off for the four statistical parameters mean, standard deviation, skewness and excess kurtosis are shown in figure 2. The mean was correct for both Gaussian redrawing and Poisson resampling. Poisson redrawing resulted in a significant overestimation at counts below 100. At a full-count pixel value of 10 the overestimation was 5% in very good agreement with the results of White and Lawson (2015).
The graph of the standard deviation shows that only Poisson resampling gives the right result. Poisson redrawing overestimates the standard deviation below 10. At a full-count pixel value of one with 16%, which is in very good agreement with the results of White and Lawson (2015). Gaussian redrawing with rounding off performed better than Poisson redrawing, but still overestimated standard deviation by 5% at a full-count pixel value of one.

Comments
Phys. Med. Biol. 60 (2015) 5711 full-count pixel value of two, the ratio increases to the not rounded off value for skewness around 1,77 as also was shown in White and Lawson (2015). The overestimation in this simulation was a bit higher than in the original paper, 25% versus approximately 20%.
The behavior for the excess kurtosis was similar to the skewness. Above 100 the numerical simulations were no longer stable, probably because the excess kurtosis was close to zero.
Above 100 the results of White and Lawson (2015) are influenced by other sources of statistical fluctuation caused by the camera system. This is not the case for the simulation described in this comment and that is also why the values for the statistical values for pixel values above 100 are correct, if numerical simulations otherwise are stable.

Discussion and conclusions
With rounding off only Poisson resampling was correct for the four calculated statistical parameters. This is because the simulated half-count data only contains integer data and the statistical parameters are correct (see table 1). Gaussian redrawing performed reasonably with limited overestimation of the standard deviation, and values reasonably close to the correct values for skewness and excess kurtosis. Poisson redrawing failed at low counts (<10) for both mean and standard deviation, and failed for skewness and excess kurtosis for counts above approximately two.
It is shown that rounding off is the sole reason for the count dependent behavior and the main reason for large deviations from the correct values for the half-count/full-count ratios of the four investigated statistical parameters, but even without rounding off only Poisson resampling demonstrates the correct values for all four statistical parameters. For simulating half (or less)-count images Poisson resampling is the method of choice.