Home
Search
Collections
Journals
About
Contact us
My IOPscience
A Poisson resampling method for simulating reduced counts in nuclear medicine images
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Med. Biol. 60 N167 (http://iopscience.iop.org/0031-9155/60/9/N167) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 157.89.65.129 This content was downloaded on 23/04/2015 at 12:30
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) N167–N176
Physics in Medicine & Biology doi:10.1088/0031-9155/60/9/N167
Note
A Poisson resampling method for simulating reduced counts in nuclear medicine images Duncan White1,3,4 and Richard S Lawson2,3 1
Barnsley Hospital NHS Foundation Trust, Barnsley, UK Central Manchester University Hospitals NHS Foundation Trust, Manchester, UK 3 The Institute of Physics and Engineering in Medicine, Nuclear Medicine Software Quality Group 2
E-mail: [email protected] Received 10 October 2014, revised 6 March 2015 Accepted for publication 11 March 2015 Published 16 April 2015 Abstract
Nuclear medicine computers now commonly offer resolution recovery and other software techniques which have been developed to improve image quality for images with low counts. These techniques potentially mean that these images can give equivalent clinical information to a full-count image. Reducing the number of counts in nuclear medicine images has the benefits of either allowing reduced activity to be administered or reducing acquisition times. However, because acquisition and processing parameters vary, each user should ideally evaluate the use of images with reduced counts within their own department, and this is best done by simulating reduced-count images from the original data. Reducing the counts in an image by division and rounding off to the nearest integer value, even if additional Poisson noise is added, is inadequate because it gives incorrect counting statistics. This technical note describes how, by applying Poisson resampling to the original raw data, simulated reduced-count images can be obtained while maintaining appropriate counting statistics. The authors have developed manufacturer independent software that can retrospectively generate simulated data with reduced counts from any acquired nuclear medicine image. Keywords: Poisson sampling, reduced count image, software, computer methods (Some figures may appear in colour only in the online journal) 4 Author to whom any correspondence should be addressed. Department of Nuclear Medicine, Barnsley Hospital, Gawber Road, Barnsley, S75 2EP, UK
0031-9155/15/09N167+10$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
N167
Note
Phys. Med. Biol. 60 (2015) N167
1. Introduction Currently there is considerable interest in using resolution recovery and other software techniques to improve image quality for low-count images such that the resulting image provides equivalent clinical information to a full-count image (Ebdon-Jackson 2011, DePuey 2012). Reducing the number of counts in nuclear medicine images has the benefits of either allowing reduced activity to be administered or reducing acquisition times. Proper evaluation of whether clinically equivalent images can be acquired with counts reduced below the normal number requires a comparison between normally acquired images and the reduced-count versions. In many situations re-acquiring patient images with reduced counts is not a viable option, because of the additional time it would take, and because repeat studies might differ in other ways such as in the patient position or changes in radiopharmaceutical distribution. In addition, ethics approval and patient consent would need to be obtained in order to repeat a patient acquisition. Acquiring static images as a series of dynamic images which can be summed afterwards is also a possibility but it would be difficult to utilise in whole body acquisitions. Either of these methods also means that any investigation must be prospective, rather than being applied retrospectively to an existing archive of image studies. A method of creating simulated reduced-count images from the original raw acquired data would overcome these problems and, if used for service evaluation purposes, would not generally require research ethics approval. Radioactive decay is statistical in nature, so given a sufficiently large population of radioactive atoms, it can be determined that an average number will decay per unit time, but the exact number that disintegrate in a given time fluctuates around this average value. Although other physical processes, such as attenuation of photons by the collimator and the detector efficiency, determine the overall sensitivity of a gamma camera image, the underlying statistical nature of the decay process means that the noise for counts in nuclear medicine images obeys a Poisson distribution (Parker et al 1978). The Poisson distribution is characterised by only one parameter, the mean, m. The Poisson distribution has the special characteristic that the standard deviation σ and skewness of the distribution γ1 are directly related to the mean: 1 σ(1) = m ; γ1 = . m
Simulation of a reduced-count image by simple mathematical division of the counts in each pixel of the raw data, rounded off to maintain integer values, will not retain these relationships and will therefore produce simulated images with incorrect noise relationships. Adding back the appropriate amount of Poisson noise into a scaled down image needs to be done cautiously, because the original raw data will already have some noise, but the amount to be added can be estimated from the scaled down counts. However, the correct noise relationships in images as counts are reduced can be more elegantly maintained by Poisson sampling (Chatfield 1970). A discussion about this is included in the Discussion section. In Poisson sampling, each element i in a population is included in the sample independently of all other units with probability πi. This selection can be carried out computationally using real time sampling. Real time sampling is a sequential method where each member of a population passes the sampler once, and as it does so the sampler determines whether that member should be included in the sample. For general Poisson sampling each element of the population may have a different probability of being included in the sample, but for the purposes of simulating a reduced-count image, we want each count in the image to have the same probability of being selected for N168
Note
Phys. Med. Biol. 60 (2015) N167
the sample. In this case Bernoulli sampling is used. Bernoulli sampling is a special case of Poisson sampling where each element of the population sampled is subjected to the same independent Bernoulli (or Binomial) trial to determine whether the element becomes part of the sample during the drawing of a single sample, in other words probability πi is the same for all elements. 2. Method Each count in an acquired image is considered to be an independent variable, or element in the population of total counts. To simulate a reduced-count image one first defines the probability πi of selection in the new image, for example 0.5 for a half-count image, or 0.75 for a three quarter-count image. Each count, c, in the original image is then allocated an independent uniform random number U1 … UN on [0,1] where N is the number of counts in the pixel. On a pixel by pixel basis, the random number associated with each count is subjected to the Bernoulli trial to see if it is less than or equal to the probability of selection (Uc ≤ πi). If true the count is selected for inclusion in the new image, otherwise it is discarded (Ghosh and Vogt 2002). See figure 1 for a flow chart of this process. One of the authors (DW) has written a program using IDL to sample in this manner. The basic code works on any 2D integer matrix and is called as a subroutine from the interface program that extracts the images from a DICOM (Digital Imaging and Communications in Medicine) file and returns a clone of the original file with the resampled images replacing the original. It adjusts the total acquired counts and image maximum tags in the DICOM header but leaves all other parameters, including any private fields, untouched so that it can still be processed normally. Any type of raw acquired nuclear medicine image can be sampled, including dynamic, SPECT and gated. It will produce an output image of the same type that simulates what would have been obtained if the administered activity had been reduced by a given factor. If the imaged activity distribution remains constant this would be equivalent to using the same administered activity but with reduced acquisition time. The program cannot be used with processed data, such as reconstructed SPECT slices or filtered images because Poisson statistics will no longer apply. The program uses IDL with the add-on DICOM package, so it can only be run on a computer that has a licence to use the DICOM package. The wrapper program allows sampling in either single image or batch mode. In batch mode all nuclear medicine images contained in the given directory or subdirectories will be processed, while non-nuclear medicine files such as CT attenuation images or screenshots will be copied but not processed. Running on a typical PC the basic code can generate a simulated image from a 100 million-count original with a 512 × 512 pixel matrix in approximately 6.5 s and will normally take much less than 1 s to complete a typical clinical image. 3. Validation To validate our sampling program, we obtained sets of uniform flood images on a dual headed Philips Skylight gamma camera, working comfortably within clinical parameters for uniformity of response. An advantage of the Philips Skylight is that it can acquire ‘concurrent imaging’ meaning that up to 15 different acquisitions can be made simultaneously with different stop conditions. Three sets of flood images were acquired using a Co-57 flood source and a matrix size of 512 × 512 pixels. The high-count set contained 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90 and 100 million total counts per image. The medium-count set contained N169
Note
Phys. Med. Biol. 60 (2015) N167
Decide probability of selection, πi
Start For each pixel in original image, p
For each count in the pixel, c Generate uniform random number on [0,1], Upc Next p Next c
Upc ≤ πi
No
Discard count
Yes Retain count in sampled image
End Sampled image created
Figure 1. Flow chart showing the process for obtaining a sampled image using a predetermined Bernoulli selection probability.
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0 and 10.0 million total counts. The low-count set contained 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100 and 200 thousand total counts. N170
Note
Phys. Med. Biol. 60 (2015) N167
Using the 5 and 10 million-count images gives a direct comparison of a pair of full-count and true half-count images that might have been obtained in the same time using half the activity. Each set of acquisitions yielded 9 pairs of full-count and corresponding half-count images, giving 27 pairs in total. Average counts per pixel in the full-count images ranged from 700 down to 0.2. Using a uniform source allowed us to quickly obtain several thousand replicate counts (one from each pixel) in a short time. Each pair of images contained 119 thousand pixels within the central field of view which can be treated as independent measures of the same count and hence should follow a Poisson distribution. Any small non-uniformities in the source or the gamma camera would not matter at the lowest pixel counts because Poisson noise dominates. At the highest pixel counts Poisson noise would be 4%, which is greater than the camera non-uniformity but camera effects might be expected to become apparent. From the original raw data of each pair, half-count images were generated by Poisson sampling. The mean, standard deviation and skewness of the distribution of all pixels were calculated within a standardised region of interest containing the useful field of view and excluding pixels near the edge of the field of view. Results for head 1 and head 2 were combined so a total of 238 thousand pixels were included in each measurement. It follows from equation (1) that if the pixel counts vary only because of noise then the standard deviation of counts in each half-count image should be reduced by a factor of 1/ 2 and the skewness increased by a factor of 2 compared with the full-count image. Therefore the mean, standard deviation and skewness of the counts in each acquired and simulated half-count image were compared with the values expected due to random noise only. Results for the validation studies are shown in figure 2. It can be seen from figure 2(a) that the mean counts per pixel for both true and Poisson resampled half counts differ from the expected value by less than 0.2%. None of the differences between true and simulated mean counts were significant when compared with 95% confidence limits for the standard error of the difference. Figure 2(b) shows that as the mean count per pixel increases, the standard deviation of counts in both the true and Poisson resampled half-count images becomes less than the expected value. This is because, in addition to the random noise, the image actually contains some variability due to real camera non-uniformity which becomes increasingly apparent as the counts increase. However it is clear that the simulated and true half-count images both show the same deviation, demonstrating that the simulation has worked correctly. The difference between the true and simulated images is always less than 0.5 percentage points. Fourteen out of the 27 simulated standard deviations did not differ significantly from the true values by an F-ratio test. Figure 2(c) shows that the skewness of the count distribution also deviates from that expected from random noise only at high counts, because of real camera non-uniformity, but there is still reasonable agreement between the real and simulated images. It is attractive to think that simpler and more natural methods of simulating reduced count data, such as direct Poisson redrawing, could obtain images with the correct statistics, but unfortunately these fail, especially at low count levels. During our validation we considered two alternative options for comparison. Alternative A is an attempt to double the noise by sampling from a Poisson distribution with a mean equal to the pixel value, and then dividing the result by 2 and rounding off to the nearest integer (remainders of 0.5 are rounded up). Alternative B first divides each pixel value by 2 (without rounding) and then adds noise by sampling from a Poisson distribution with mean equal to the halved pixel value, which guarantees an integer result but will add too much noise. Without rounding, both alternatives will result in a correct mean value. Alternative A will also give a correct magnitude for the noise but alternative B will theoretically result in noise that is too high by a factor of 1.225 N171
Note
Phys. Med. Biol. 60 (2015) N167
(a) Percentage difference from expected in mean pixel count of half-count images 10%
Percentage difference from expected
9% 8% 7% True half counts
6%
Poisson resampled half counts
5%
Alternative A
4%
Alternative B
3% 2% 1% 0% -1% -2% 0.1
1.0
10.0
100.0
1000.0
Counts per pixel in original image
(b) Percentage difference from expected in standard deviation of half-count images 30%
Percentage difference from expected
25% 20% True half counts
15%
Poisson resampled half counts Alternative A
10%
Alternative B
5% 0% -5% -10% 0.1
1.0
10.0
100.0
1000.0
Counts per pixel in original image
(c) Percentage difference from expected in skewness of half-count images 50%
Percentage difference from expected
40% True half counts Poisson resampled half counts
30%
Alternative A
20%
Alternative B
10% 0% -10% -20% -30% 0.1
1.0
10.0
100.0
1000.0
Counts per pixel in original image
Figure 2. Charts showing the percentage difference from expected values due to
random noise in a flood image for (a) mean, (b) standard deviation and (c) skewness for true half-count images and simulated images generated by Poisson resampling and two alternative methods which are described in the text. Error bars in (a) represent the standard error of the mean difference. N172
Note
Phys. Med. Biol. 60 (2015) N167
Figure 3. Clinical images demonstrating the equivalence of reduced acquisition time and Poisson resampled counts.
(see the Discussion section for the derivation of this value). The results for these two alternative methods are also shown in figure 2 for comparison with true half-count and Poisson resampled data. Figure 2(a) demonstrates that, solely due to rounding errors, alternative A only gives an accurate estimate of mean counts at very high count levels (0.5% deviation at a full count of 100), whereas alternative B, which does not involve a rounding step, gives a good estimate of the mean for all count levels. Figure 2(b) shows that alternative A reproduces the correct standard deviation at high counts but significantly overestimates the standard deviation at levels of less than 10 counts per pixel, again solely due to rounding off. Alternative B overestimates the standard deviation at all levels and, until the camera non-uniformity becomes significant, the overestimation is close to the theoretical value of 22.5%. Figure 2(c) shows that neither of the alternative methods are as good as Poisson resampling at reproducing the skewness, although alternative A comes close at very low counts. To demonstrate that our Poisson resampling method maintains the texture and distribution of noise in a clinical setting when compared to true acquisitions we acquired a static bone scan image on a Philips Skylight gamma camera, using the concurrent imaging facility to acquire images simultaneously for 200 s and 50 s. We then simulated the 50 s acquisition by Poisson resampling the original with a cut-off of 25%. The results are shown in figure 3. Although the true 50 s image and that simulated by Poisson resampling are not identical, the images have a very similar count distribution and the texture of noise is preserved for the high and low-count regions of the image. 4. Discussion Simulating reduced counts in an image by simple division of counts in each pixel will obviously not give the correct result because the standard deviation will be too low, as the noise distribution is based on the noise from an image with the higher number of counts. However, adding additional Poisson noise either before division by two (alternative A) or after division (alternative B) still does not result in an image with the correct statistical properties. As can be seen in figure 2(a) alternative A only gives an accurate estimate of mean counts at high count levels. This occurs because when any pixel is divided by two, pixels with even counts can be correctly halved, but all odd numbers will have the fractional part rounded up, so the mean is overestimated by approximately 0.25, and thus an error of 5% at a count level of 10. When the average counts per pixel are less than 1 the rounding effect is even worse, because most pixels will contain only 0 or 1 count and neither of these values are N173
Note
Phys. Med. Biol. 60 (2015) N167
changed after division by two with rounding and so counts in these pixels are not halved at all. Alternative B does not involve a rounding step and so it gives a good estimate of the mean for all count levels. The effect of rounding off could be reduced by applying a Bernoulli process to the rounding but that would make the method more complicated than the straightforward Bernoulli procedure described in figure 1 and would probably introduce other undesirable effects on the counting statistics. Alternative A accurately reproduces the correct standard deviation for images with more than about 10 counts per pixel but fails for lower counts (figure 2(b)). This can be explained by observing that if the mean count in the original image is m it will have a variance due to Poisson noise of m and if we add some additional Poisson noise this will increase to 2m. When we divide all pixels by 2 the variance reduces by a factor of 4 and so the variance becomes m/2 which is correct. This works in practice if m is large due to the relatively small effect of rounding errors. However, the method does not work properly for low counts per pixel where the effects of rounding off are significant. Note also that it wouldn’t work for division by any factor other than 2. For alternative B we start by dividing the original data by 2, which simply halves the mean and reduces the variance to m/4. We then add more Poisson noise to give an additional variance of m/2 which makes the total variance 3m/4 instead of the expected m/2, yielding a standard deviation which is too high by a factor of 1.5 = 1.225. This accounts for why alternative B overestimates the correct standard deviation even at high counts (figure 2(b)). The skewness of a Poisson distribution with mean m is 1/√ m so the skewness of the halfcount distribution should be increased by a factor of √2. However figure 2(c) shows that both alternative methods produce a skewness that is greater than expected. Dividing by two will not change the skewness, so this overestimate is presumably due to the effect of adding additional Poisson noise. Ideally the noise would be based on the true signal in each pixel but that is unknown; all we have is an already noisy distribution of pixel values. So large values (in the right hand tail of the distribution) will get more noise added than smaller values (in the left hand tail of the distribution) and therefore the skewness becomes increased. However alternative A includes a rounding step and at very low levels of only a few counts per pixel the rounding errors associated with pixels containing counts of 0 and 1 have the opposite effect of enhancing the left hand tail of the distribution. This appears to mitigate the overestimate and results in an approximately correct skewness when there are on average fewer than 2 counts per pixel in the original image. In contrast, Poisson sampling accurately simulates the statistical nature of a real acquisition and so maintains the correct mean, standard deviation and skewness, over the full range of counts, even when overall counts are so low that most pixels have no counts in them. Using this method we have shown that the true and simulated half-counts for multiple replications of a single pixel follow each other very closely over the range from 0.2 to 700 counts per pixel. Therefore the results can be applied to any arbitrary activity distribution with pixel counts in this range, and there is no reason to suppose that it will not continue to be applicable beyond this range. It is possible to produce true half-count images in other ways. Repeating the acquisition with half the acquisition time is not generally desirable because of the additional patient time required and the possibility of other differences, such as patient movement, between the two studies. It is possible to acquire a static image as a dynamic sequence and subsequently add some frames to simulate different acquisition times. Whole body studies, however cannot generally be acquired as a dynamic sequence so it would be necessary to acquire several separate studies and add them, which could lead to additional differences due to patient movement between studies. Dynamic SPECT acquisitions can be used to acquire several successive N174
Note
Phys. Med. Biol. 60 (2015) N167
rotations which can be summed afterwards, but this also introduces the possibility of patient movement between rotations. It is better to acquire a SPECT study in gated mode using an ECG simulator as a trigger, because this allows the data at each angle to be acquired into 8 equivalent bins which can be added subsequently. This method is preferable to acquiring more frames per rotation and rejecting or summing alternate frames as this introduces variations in angular resolution. Some Philips gamma cameras include a concurrent imaging facility that allows full-time and half-time images to be acquired simultaneously which is ideal for this purpose. Data could also be acquired in list mode and rebinned afterwards, but this is an option that is no longer available on most gamma cameras (although it is often available on PET scanners). These methods are very convenient but they must be set up in advance so they cannot be applied retrospectively to existing data. Some gamma camera manufacturers offer software similar to that described here that can perform Poisson resampling for particular applications, but it is usually specific to one manufacturer. The software that we have developed has the advantage that it is manufacturer independent and can be applied retrospectively to data acquired on any gamma camera. The simulated images represent new acquired data which may be processed in the same way as the original or in a different way, depending on the application. For example different processing may be required when testing resolution recovery SPECT reconstruction software with half counts when compared with standard reconstruction software on the original data (Lawson et al 2014). The methods and programs described here have been used successfully for a multi-centre, multi-vendor audit to evaluate the use of resolution recovery software from several manufactures when applied to myocardial perfusion data with half the normal counts acquired under a variety of clinical protocols in a range of departments (Lawson et al 2014). In total 16 centres, using cameras from three manufactures, and processing computers from four manufacturers, processed and compared a total of 769 patient studies with good results from Poisson resampled half-count images in all cases. Simulating half-count images from existing data is an easy way to investigate the possibility of either using a lower administered activity for current studies or reducing acquisition time. In the former case acquisition parameters would be the same so the simulation should be exact, but in the latter case the shorter acquisition time would also give the added benefit of less chance of patient movement. At a time when there are increasing pressures to improve patient throughput and also minimise the administered activity of 99mTc, we have demonstrated that Poisson resampling is an effective method for simulating images with reduced counts retrospectively from existing acquired image sets in order to optimise clinical acquisition parameters. Conflicts of interest and source of funding The Nuclear Medicine Software Quality Group is a sub-committee of the Institute of Physics and Engineering in Medicine Nuclear Medicine Special Interest Group. The group’s web site is hosted by Link Medical. No conflicts were declared for the authors. References Chatfield C 1970 Statistics for Technology (Hammondsworth: Penguin Books) DePuey E G 2012 Advances in SPECT camera software and hardware: currently available and new on the horizon J. Nucl. Cardiol. 19 551–81 N175
Note
Phys. Med. Biol. 60 (2015) N167
Ebdon-Jackson S 2011 The use of resolution recovery software in nuclear medicine from an ALARA perspective, European ALARA Network Newsletter 29 (www.eu-alara.net/index.php/newslettersmainmenu-37/65-newsletter-29/269-ebdon.html) Ghosh D and Vogt A 2002 Sampling methods related to Bernoulli and Poisson sampling Proc. of the American Statistical Association, Survey Research Methods Section pp 3569–70 (www.amstat.org/ sections/srms/proceedings/y2002/files/JSM2002-001080.pdf) Lawson R S, White D, Nijran K, Cade S C, Hall D O, Kenny B, Knight A, Livieratos L, Murray A and Towey D 2014 An audit of half-count myocardial perfusion imaging using resolution recovery software Nucl. Med. Commun. 35 511–21 Parker P, Smith P and Taylor D 1978 Basic Science of Nuclear Medicine (Edinburgh: Churchill Livingstone)
N176
Search
Collections
Journals
About
Contact us
My IOPscience
A Poisson resampling method for simulating reduced counts in nuclear medicine images
This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Med. Biol. 60 N167 (http://iopscience.iop.org/0031-9155/60/9/N167) View the table of contents for this issue, or go to the journal homepage for more
Download details: IP Address: 157.89.65.129 This content was downloaded on 23/04/2015 at 12:30
Please note that terms and conditions apply.
Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) N167–N176
Physics in Medicine & Biology doi:10.1088/0031-9155/60/9/N167
Note
A Poisson resampling method for simulating reduced counts in nuclear medicine images Duncan White1,3,4 and Richard S Lawson2,3 1
Barnsley Hospital NHS Foundation Trust, Barnsley, UK Central Manchester University Hospitals NHS Foundation Trust, Manchester, UK 3 The Institute of Physics and Engineering in Medicine, Nuclear Medicine Software Quality Group 2
E-mail: [email protected] Received 10 October 2014, revised 6 March 2015 Accepted for publication 11 March 2015 Published 16 April 2015 Abstract
Nuclear medicine computers now commonly offer resolution recovery and other software techniques which have been developed to improve image quality for images with low counts. These techniques potentially mean that these images can give equivalent clinical information to a full-count image. Reducing the number of counts in nuclear medicine images has the benefits of either allowing reduced activity to be administered or reducing acquisition times. However, because acquisition and processing parameters vary, each user should ideally evaluate the use of images with reduced counts within their own department, and this is best done by simulating reduced-count images from the original data. Reducing the counts in an image by division and rounding off to the nearest integer value, even if additional Poisson noise is added, is inadequate because it gives incorrect counting statistics. This technical note describes how, by applying Poisson resampling to the original raw data, simulated reduced-count images can be obtained while maintaining appropriate counting statistics. The authors have developed manufacturer independent software that can retrospectively generate simulated data with reduced counts from any acquired nuclear medicine image. Keywords: Poisson sampling, reduced count image, software, computer methods (Some figures may appear in colour only in the online journal) 4 Author to whom any correspondence should be addressed. Department of Nuclear Medicine, Barnsley Hospital, Gawber Road, Barnsley, S75 2EP, UK
0031-9155/15/09N167+10$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
N167
Note
Phys. Med. Biol. 60 (2015) N167
1. Introduction Currently there is considerable interest in using resolution recovery and other software techniques to improve image quality for low-count images such that the resulting image provides equivalent clinical information to a full-count image (Ebdon-Jackson 2011, DePuey 2012). Reducing the number of counts in nuclear medicine images has the benefits of either allowing reduced activity to be administered or reducing acquisition times. Proper evaluation of whether clinically equivalent images can be acquired with counts reduced below the normal number requires a comparison between normally acquired images and the reduced-count versions. In many situations re-acquiring patient images with reduced counts is not a viable option, because of the additional time it would take, and because repeat studies might differ in other ways such as in the patient position or changes in radiopharmaceutical distribution. In addition, ethics approval and patient consent would need to be obtained in order to repeat a patient acquisition. Acquiring static images as a series of dynamic images which can be summed afterwards is also a possibility but it would be difficult to utilise in whole body acquisitions. Either of these methods also means that any investigation must be prospective, rather than being applied retrospectively to an existing archive of image studies. A method of creating simulated reduced-count images from the original raw acquired data would overcome these problems and, if used for service evaluation purposes, would not generally require research ethics approval. Radioactive decay is statistical in nature, so given a sufficiently large population of radioactive atoms, it can be determined that an average number will decay per unit time, but the exact number that disintegrate in a given time fluctuates around this average value. Although other physical processes, such as attenuation of photons by the collimator and the detector efficiency, determine the overall sensitivity of a gamma camera image, the underlying statistical nature of the decay process means that the noise for counts in nuclear medicine images obeys a Poisson distribution (Parker et al 1978). The Poisson distribution is characterised by only one parameter, the mean, m. The Poisson distribution has the special characteristic that the standard deviation σ and skewness of the distribution γ1 are directly related to the mean: 1 σ(1) = m ; γ1 = . m
Simulation of a reduced-count image by simple mathematical division of the counts in each pixel of the raw data, rounded off to maintain integer values, will not retain these relationships and will therefore produce simulated images with incorrect noise relationships. Adding back the appropriate amount of Poisson noise into a scaled down image needs to be done cautiously, because the original raw data will already have some noise, but the amount to be added can be estimated from the scaled down counts. However, the correct noise relationships in images as counts are reduced can be more elegantly maintained by Poisson sampling (Chatfield 1970). A discussion about this is included in the Discussion section. In Poisson sampling, each element i in a population is included in the sample independently of all other units with probability πi. This selection can be carried out computationally using real time sampling. Real time sampling is a sequential method where each member of a population passes the sampler once, and as it does so the sampler determines whether that member should be included in the sample. For general Poisson sampling each element of the population may have a different probability of being included in the sample, but for the purposes of simulating a reduced-count image, we want each count in the image to have the same probability of being selected for N168
Note
Phys. Med. Biol. 60 (2015) N167
the sample. In this case Bernoulli sampling is used. Bernoulli sampling is a special case of Poisson sampling where each element of the population sampled is subjected to the same independent Bernoulli (or Binomial) trial to determine whether the element becomes part of the sample during the drawing of a single sample, in other words probability πi is the same for all elements. 2. Method Each count in an acquired image is considered to be an independent variable, or element in the population of total counts. To simulate a reduced-count image one first defines the probability πi of selection in the new image, for example 0.5 for a half-count image, or 0.75 for a three quarter-count image. Each count, c, in the original image is then allocated an independent uniform random number U1 … UN on [0,1] where N is the number of counts in the pixel. On a pixel by pixel basis, the random number associated with each count is subjected to the Bernoulli trial to see if it is less than or equal to the probability of selection (Uc ≤ πi). If true the count is selected for inclusion in the new image, otherwise it is discarded (Ghosh and Vogt 2002). See figure 1 for a flow chart of this process. One of the authors (DW) has written a program using IDL to sample in this manner. The basic code works on any 2D integer matrix and is called as a subroutine from the interface program that extracts the images from a DICOM (Digital Imaging and Communications in Medicine) file and returns a clone of the original file with the resampled images replacing the original. It adjusts the total acquired counts and image maximum tags in the DICOM header but leaves all other parameters, including any private fields, untouched so that it can still be processed normally. Any type of raw acquired nuclear medicine image can be sampled, including dynamic, SPECT and gated. It will produce an output image of the same type that simulates what would have been obtained if the administered activity had been reduced by a given factor. If the imaged activity distribution remains constant this would be equivalent to using the same administered activity but with reduced acquisition time. The program cannot be used with processed data, such as reconstructed SPECT slices or filtered images because Poisson statistics will no longer apply. The program uses IDL with the add-on DICOM package, so it can only be run on a computer that has a licence to use the DICOM package. The wrapper program allows sampling in either single image or batch mode. In batch mode all nuclear medicine images contained in the given directory or subdirectories will be processed, while non-nuclear medicine files such as CT attenuation images or screenshots will be copied but not processed. Running on a typical PC the basic code can generate a simulated image from a 100 million-count original with a 512 × 512 pixel matrix in approximately 6.5 s and will normally take much less than 1 s to complete a typical clinical image. 3. Validation To validate our sampling program, we obtained sets of uniform flood images on a dual headed Philips Skylight gamma camera, working comfortably within clinical parameters for uniformity of response. An advantage of the Philips Skylight is that it can acquire ‘concurrent imaging’ meaning that up to 15 different acquisitions can be made simultaneously with different stop conditions. Three sets of flood images were acquired using a Co-57 flood source and a matrix size of 512 × 512 pixels. The high-count set contained 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90 and 100 million total counts per image. The medium-count set contained N169
Note
Phys. Med. Biol. 60 (2015) N167
Decide probability of selection, πi
Start For each pixel in original image, p
For each count in the pixel, c Generate uniform random number on [0,1], Upc Next p Next c
Upc ≤ πi
No
Discard count
Yes Retain count in sampled image
End Sampled image created
Figure 1. Flow chart showing the process for obtaining a sampled image using a predetermined Bernoulli selection probability.
0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0 and 10.0 million total counts. The low-count set contained 10, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80, 90, 100 and 200 thousand total counts. N170
Note
Phys. Med. Biol. 60 (2015) N167
Using the 5 and 10 million-count images gives a direct comparison of a pair of full-count and true half-count images that might have been obtained in the same time using half the activity. Each set of acquisitions yielded 9 pairs of full-count and corresponding half-count images, giving 27 pairs in total. Average counts per pixel in the full-count images ranged from 700 down to 0.2. Using a uniform source allowed us to quickly obtain several thousand replicate counts (one from each pixel) in a short time. Each pair of images contained 119 thousand pixels within the central field of view which can be treated as independent measures of the same count and hence should follow a Poisson distribution. Any small non-uniformities in the source or the gamma camera would not matter at the lowest pixel counts because Poisson noise dominates. At the highest pixel counts Poisson noise would be 4%, which is greater than the camera non-uniformity but camera effects might be expected to become apparent. From the original raw data of each pair, half-count images were generated by Poisson sampling. The mean, standard deviation and skewness of the distribution of all pixels were calculated within a standardised region of interest containing the useful field of view and excluding pixels near the edge of the field of view. Results for head 1 and head 2 were combined so a total of 238 thousand pixels were included in each measurement. It follows from equation (1) that if the pixel counts vary only because of noise then the standard deviation of counts in each half-count image should be reduced by a factor of 1/ 2 and the skewness increased by a factor of 2 compared with the full-count image. Therefore the mean, standard deviation and skewness of the counts in each acquired and simulated half-count image were compared with the values expected due to random noise only. Results for the validation studies are shown in figure 2. It can be seen from figure 2(a) that the mean counts per pixel for both true and Poisson resampled half counts differ from the expected value by less than 0.2%. None of the differences between true and simulated mean counts were significant when compared with 95% confidence limits for the standard error of the difference. Figure 2(b) shows that as the mean count per pixel increases, the standard deviation of counts in both the true and Poisson resampled half-count images becomes less than the expected value. This is because, in addition to the random noise, the image actually contains some variability due to real camera non-uniformity which becomes increasingly apparent as the counts increase. However it is clear that the simulated and true half-count images both show the same deviation, demonstrating that the simulation has worked correctly. The difference between the true and simulated images is always less than 0.5 percentage points. Fourteen out of the 27 simulated standard deviations did not differ significantly from the true values by an F-ratio test. Figure 2(c) shows that the skewness of the count distribution also deviates from that expected from random noise only at high counts, because of real camera non-uniformity, but there is still reasonable agreement between the real and simulated images. It is attractive to think that simpler and more natural methods of simulating reduced count data, such as direct Poisson redrawing, could obtain images with the correct statistics, but unfortunately these fail, especially at low count levels. During our validation we considered two alternative options for comparison. Alternative A is an attempt to double the noise by sampling from a Poisson distribution with a mean equal to the pixel value, and then dividing the result by 2 and rounding off to the nearest integer (remainders of 0.5 are rounded up). Alternative B first divides each pixel value by 2 (without rounding) and then adds noise by sampling from a Poisson distribution with mean equal to the halved pixel value, which guarantees an integer result but will add too much noise. Without rounding, both alternatives will result in a correct mean value. Alternative A will also give a correct magnitude for the noise but alternative B will theoretically result in noise that is too high by a factor of 1.225 N171
Note
Phys. Med. Biol. 60 (2015) N167
(a) Percentage difference from expected in mean pixel count of half-count images 10%
Percentage difference from expected
9% 8% 7% True half counts
6%
Poisson resampled half counts
5%
Alternative A
4%
Alternative B
3% 2% 1% 0% -1% -2% 0.1
1.0
10.0
100.0
1000.0
Counts per pixel in original image
(b) Percentage difference from expected in standard deviation of half-count images 30%
Percentage difference from expected
25% 20% True half counts
15%
Poisson resampled half counts Alternative A
10%
Alternative B
5% 0% -5% -10% 0.1
1.0
10.0
100.0
1000.0
Counts per pixel in original image
(c) Percentage difference from expected in skewness of half-count images 50%
Percentage difference from expected
40% True half counts Poisson resampled half counts
30%
Alternative A
20%
Alternative B
10% 0% -10% -20% -30% 0.1
1.0
10.0
100.0
1000.0
Counts per pixel in original image
Figure 2. Charts showing the percentage difference from expected values due to
random noise in a flood image for (a) mean, (b) standard deviation and (c) skewness for true half-count images and simulated images generated by Poisson resampling and two alternative methods which are described in the text. Error bars in (a) represent the standard error of the mean difference. N172
Note
Phys. Med. Biol. 60 (2015) N167
Figure 3. Clinical images demonstrating the equivalence of reduced acquisition time and Poisson resampled counts.
(see the Discussion section for the derivation of this value). The results for these two alternative methods are also shown in figure 2 for comparison with true half-count and Poisson resampled data. Figure 2(a) demonstrates that, solely due to rounding errors, alternative A only gives an accurate estimate of mean counts at very high count levels (0.5% deviation at a full count of 100), whereas alternative B, which does not involve a rounding step, gives a good estimate of the mean for all count levels. Figure 2(b) shows that alternative A reproduces the correct standard deviation at high counts but significantly overestimates the standard deviation at levels of less than 10 counts per pixel, again solely due to rounding off. Alternative B overestimates the standard deviation at all levels and, until the camera non-uniformity becomes significant, the overestimation is close to the theoretical value of 22.5%. Figure 2(c) shows that neither of the alternative methods are as good as Poisson resampling at reproducing the skewness, although alternative A comes close at very low counts. To demonstrate that our Poisson resampling method maintains the texture and distribution of noise in a clinical setting when compared to true acquisitions we acquired a static bone scan image on a Philips Skylight gamma camera, using the concurrent imaging facility to acquire images simultaneously for 200 s and 50 s. We then simulated the 50 s acquisition by Poisson resampling the original with a cut-off of 25%. The results are shown in figure 3. Although the true 50 s image and that simulated by Poisson resampling are not identical, the images have a very similar count distribution and the texture of noise is preserved for the high and low-count regions of the image. 4. Discussion Simulating reduced counts in an image by simple division of counts in each pixel will obviously not give the correct result because the standard deviation will be too low, as the noise distribution is based on the noise from an image with the higher number of counts. However, adding additional Poisson noise either before division by two (alternative A) or after division (alternative B) still does not result in an image with the correct statistical properties. As can be seen in figure 2(a) alternative A only gives an accurate estimate of mean counts at high count levels. This occurs because when any pixel is divided by two, pixels with even counts can be correctly halved, but all odd numbers will have the fractional part rounded up, so the mean is overestimated by approximately 0.25, and thus an error of 5% at a count level of 10. When the average counts per pixel are less than 1 the rounding effect is even worse, because most pixels will contain only 0 or 1 count and neither of these values are N173
Note
Phys. Med. Biol. 60 (2015) N167
changed after division by two with rounding and so counts in these pixels are not halved at all. Alternative B does not involve a rounding step and so it gives a good estimate of the mean for all count levels. The effect of rounding off could be reduced by applying a Bernoulli process to the rounding but that would make the method more complicated than the straightforward Bernoulli procedure described in figure 1 and would probably introduce other undesirable effects on the counting statistics. Alternative A accurately reproduces the correct standard deviation for images with more than about 10 counts per pixel but fails for lower counts (figure 2(b)). This can be explained by observing that if the mean count in the original image is m it will have a variance due to Poisson noise of m and if we add some additional Poisson noise this will increase to 2m. When we divide all pixels by 2 the variance reduces by a factor of 4 and so the variance becomes m/2 which is correct. This works in practice if m is large due to the relatively small effect of rounding errors. However, the method does not work properly for low counts per pixel where the effects of rounding off are significant. Note also that it wouldn’t work for division by any factor other than 2. For alternative B we start by dividing the original data by 2, which simply halves the mean and reduces the variance to m/4. We then add more Poisson noise to give an additional variance of m/2 which makes the total variance 3m/4 instead of the expected m/2, yielding a standard deviation which is too high by a factor of 1.5 = 1.225. This accounts for why alternative B overestimates the correct standard deviation even at high counts (figure 2(b)). The skewness of a Poisson distribution with mean m is 1/√ m so the skewness of the halfcount distribution should be increased by a factor of √2. However figure 2(c) shows that both alternative methods produce a skewness that is greater than expected. Dividing by two will not change the skewness, so this overestimate is presumably due to the effect of adding additional Poisson noise. Ideally the noise would be based on the true signal in each pixel but that is unknown; all we have is an already noisy distribution of pixel values. So large values (in the right hand tail of the distribution) will get more noise added than smaller values (in the left hand tail of the distribution) and therefore the skewness becomes increased. However alternative A includes a rounding step and at very low levels of only a few counts per pixel the rounding errors associated with pixels containing counts of 0 and 1 have the opposite effect of enhancing the left hand tail of the distribution. This appears to mitigate the overestimate and results in an approximately correct skewness when there are on average fewer than 2 counts per pixel in the original image. In contrast, Poisson sampling accurately simulates the statistical nature of a real acquisition and so maintains the correct mean, standard deviation and skewness, over the full range of counts, even when overall counts are so low that most pixels have no counts in them. Using this method we have shown that the true and simulated half-counts for multiple replications of a single pixel follow each other very closely over the range from 0.2 to 700 counts per pixel. Therefore the results can be applied to any arbitrary activity distribution with pixel counts in this range, and there is no reason to suppose that it will not continue to be applicable beyond this range. It is possible to produce true half-count images in other ways. Repeating the acquisition with half the acquisition time is not generally desirable because of the additional patient time required and the possibility of other differences, such as patient movement, between the two studies. It is possible to acquire a static image as a dynamic sequence and subsequently add some frames to simulate different acquisition times. Whole body studies, however cannot generally be acquired as a dynamic sequence so it would be necessary to acquire several separate studies and add them, which could lead to additional differences due to patient movement between studies. Dynamic SPECT acquisitions can be used to acquire several successive N174
Note
Phys. Med. Biol. 60 (2015) N167
rotations which can be summed afterwards, but this also introduces the possibility of patient movement between rotations. It is better to acquire a SPECT study in gated mode using an ECG simulator as a trigger, because this allows the data at each angle to be acquired into 8 equivalent bins which can be added subsequently. This method is preferable to acquiring more frames per rotation and rejecting or summing alternate frames as this introduces variations in angular resolution. Some Philips gamma cameras include a concurrent imaging facility that allows full-time and half-time images to be acquired simultaneously which is ideal for this purpose. Data could also be acquired in list mode and rebinned afterwards, but this is an option that is no longer available on most gamma cameras (although it is often available on PET scanners). These methods are very convenient but they must be set up in advance so they cannot be applied retrospectively to existing data. Some gamma camera manufacturers offer software similar to that described here that can perform Poisson resampling for particular applications, but it is usually specific to one manufacturer. The software that we have developed has the advantage that it is manufacturer independent and can be applied retrospectively to data acquired on any gamma camera. The simulated images represent new acquired data which may be processed in the same way as the original or in a different way, depending on the application. For example different processing may be required when testing resolution recovery SPECT reconstruction software with half counts when compared with standard reconstruction software on the original data (Lawson et al 2014). The methods and programs described here have been used successfully for a multi-centre, multi-vendor audit to evaluate the use of resolution recovery software from several manufactures when applied to myocardial perfusion data with half the normal counts acquired under a variety of clinical protocols in a range of departments (Lawson et al 2014). In total 16 centres, using cameras from three manufactures, and processing computers from four manufacturers, processed and compared a total of 769 patient studies with good results from Poisson resampled half-count images in all cases. Simulating half-count images from existing data is an easy way to investigate the possibility of either using a lower administered activity for current studies or reducing acquisition time. In the former case acquisition parameters would be the same so the simulation should be exact, but in the latter case the shorter acquisition time would also give the added benefit of less chance of patient movement. At a time when there are increasing pressures to improve patient throughput and also minimise the administered activity of 99mTc, we have demonstrated that Poisson resampling is an effective method for simulating images with reduced counts retrospectively from existing acquired image sets in order to optimise clinical acquisition parameters. Conflicts of interest and source of funding The Nuclear Medicine Software Quality Group is a sub-committee of the Institute of Physics and Engineering in Medicine Nuclear Medicine Special Interest Group. The group’s web site is hosted by Link Medical. No conflicts were declared for the authors. References Chatfield C 1970 Statistics for Technology (Hammondsworth: Penguin Books) DePuey E G 2012 Advances in SPECT camera software and hardware: currently available and new on the horizon J. Nucl. Cardiol. 19 551–81 N175
Note
Phys. Med. Biol. 60 (2015) N167
Ebdon-Jackson S 2011 The use of resolution recovery software in nuclear medicine from an ALARA perspective, European ALARA Network Newsletter 29 (www.eu-alara.net/index.php/newslettersmainmenu-37/65-newsletter-29/269-ebdon.html) Ghosh D and Vogt A 2002 Sampling methods related to Bernoulli and Poisson sampling Proc. of the American Statistical Association, Survey Research Methods Section pp 3569–70 (www.amstat.org/ sections/srms/proceedings/y2002/files/JSM2002-001080.pdf) Lawson R S, White D, Nijran K, Cade S C, Hall D O, Kenny B, Knight A, Livieratos L, Murray A and Towey D 2014 An audit of half-count myocardial perfusion imaging using resolution recovery software Nucl. Med. Commun. 35 511–21 Parker P, Smith P and Taylor D 1978 Basic Science of Nuclear Medicine (Edinburgh: Churchill Livingstone)
N176
