Adding Poisson Noise to an Image

QUESTION: What is the best way to add random Poisson noise to an image in IDL? I would like to add five percent random noise to my image.

ANSWER: Yes, adding random Poisson noise to an image is a little harder than just adding random uniform noise to an image, since you have to compute the Poisson random deviate for each value in the image. (The fact that the Poisson keyword to RandomU only accepts a scalar Poisson mean is what makes this step more difficult.)

Wayne Landsman, curator of the NASA IDL Astronomy Library, recommends using the POIDEV routine from that library, which makes the process exceedingly easy.

   IDL> image = image + 0.5 * PoiDev(image)

However, you can also write your own function in IDL to do this, as shown in the code below. The function is named AddPoissonNoise, and it is adapted from a program offered by Wayne Landsman and modified by David Platten in a 12 April 2013 IDL newsgroup article on this subject.

   FUNCTION AddPoissonNoise, image, PERCENTNOISE=percentNoise, SEED=seed
       SetDefaultValue, image, cgDemoData(20)
       SetDefaultValue, percentNoise, 5.0
       h = cgHistogram(image, MIN=0, REVERSE_INDICES=ri)
       noiseImage = image
       FOR j=0,N_Elements(h)-1 DO BEGIN
           currentMean = j
           indices = ReverseIndices(ri, j, COUNT=count)
           IF count GT 0 THEN BEGIN
               noiseImage[indices] = RandomU(seed, count, POISSON=currentMean)
               newMean = SqRt(currentMean) / (percentNoise/100.0)
               
               ; Adjust the standard deviation to the required level.
               noiseImage[indices] += newMean - currentMean
               
               ; Now adjust it back to starting the mean, but with noise added.
               noiseImage[indices] *= currentMean / newMean
           ENDIF
       ENDFOR
       RETURN, noiseImage
   END

You can test the function like this, in which 15% noise is added to the image.

   image = cgDemoData(20)
   seed = -3L ; Only so that the example here matches your result!
   Print, "Mean and standard deviation of original image is: " $
          + StrTrim(Mean(image),2) + ", " + StrTrim(StdDev(image),2)
   noisyImage = AddPoissonNoise(image, PercentNoise=15, Seed=seed)
   Print, "Mean and standard deviation of noisy image is: " $
          + StrTrim(Mean(noisyImage ),2) + ", " + StrTrim(StdDev(noisyImage),2)
   !P.Multi=[0,2,1]
   cgDisplay, 700, 375, /Free, Title='Poisson Noise Added to Image'
   cgImage, image, MultiMargin=2
   cgImage, noisyImage, MultiMargin=2
   !P.Multi = 0

Executing these commands produces the following output, and the image plot shown in the figure below.

   Mean and standard deviation of original image is: 115.418, 21.1413
   Mean and standard deviation of noisy image is: 113.930, 27.4261

Original image on the left and image with 15% noise added on the right.

 

Version of IDL used to prepare this article: IDL 8.2.2.