Adding poisson noise to an image

Question

I have some images that I need to add incremental amounts of Poisson noise to in order to more thoroughly analyze them. I know you can do this in MATLAB, but how do you go about doing it in Python? Searches have yielded nothing so far.

Undersell

Answer 1

If numpy/scipy are available to you, then the following should help. I recommend that you cast the arrays to float for intermediate computations then cast back to uint8 for output/display purposes. As poisson noise is all >=0, you will need to decide how you want to handle overflow of your arrays as you cast back to uint8. You could scale or truncate depending on what your goals were.

filename = 'myimage.png'
imagea = (scipy.misc.imread(filename)).astype(float)

poissonNoise = numpy.random.poisson(imagea).astype(float)

noisyImage = imagea + poissonNoise

#here care must be taken to re cast the result to uint8 if needed or scale to 0-1 etc...

Proportional

It's confusing because the answer has been corrected but is still downvoted... Stackoverflow definitely needs a rewamp. – Nought

I did try to address the original concern for the downvote @Nought as it was something I wanted to keep up to date and address any concerns. My use of the poisson function was initially off but has been fixed. Glad you noticed. – Proportional

Of course, thanks for correcting it! I'm only upset at stackoverflow, their system just makes it more confusing to people trying to understand. – Nought

Answer 2

The answer of Helder is correct. I just want to add the fact that Poisson noise is not additive and you can not add it as Gaussian noise.

Depend on what you want to achieve, here is some suggestions:

• Simulate a low-light noisy image (if PEAK = 1, it will be really noisy)

  import numpy as np
  image = read_image("YOUR_IMAGE")  # need a rescale to be more realistic
  noisy = np.random.poisson(image / 255.0 * PEAK) / PEAK * 255  # noisy image
  

• Add a noise layer on top of the clean image

  import numpy as np
  image = read_image("YOUR_IMAGE") 
  noisemap = create_noisemap() 
  noisy = image + np.random.poisson(noisemap)  
  

Then you can crop the result to 0 - 255 if you like (I use PIL so I use 255 instead of 1).

Sexcentenary

Poisson distribution generates integers. If the image is a np.float array with a range [0,1], the correct code should be noisy = np.random.poisson(image * 255.0 * PEAK) / PEAK / 255 – Mispronounce

What is create_noisemap in this case? – Dysarthria

what is create_noisemap...? – Contiguous

This works for me: noisemap = np.ones((image_size_x, image_size_y)) * mean – Diehl

Answer 3

Actually the answer of Paul doesn't make sense.

Poisson noise is signal dependent! And using those commands, provided by him, the noise later added to the image is not signal dependent.

To make it signal dependent you should pass the image to the NumPy's poisson function:

filename = 'myimage.png'
img = (scipy.misc.imread(filename)).astype(float)

noise_mask = numpy.random.poisson(img)

noisy_img = img + noise_mask

Haerr

You r answer makes sense. Though how do you handle the saturation this causes? – Ribbonwood

Well, in this case you have to handle the saturation by yourself, setting to 0 the negative values and setting to the highest value of the number of bits of your image, those values that became highest than that. – Haerr

Thanks. I handled it by saturating the top 5% pixels – Ribbonwood

Can't there be intrinsic and extrinsic sources of poisson distributed noise? – Dena

I don't think Poisson Noise is additive like Gaussian noise. Please refer @Sexcentenary 's answer for reference. – Hakon

@Hakon you didnt get the idea. Poisson noise is not additive, yes we both know. But you need to generate a random mask (created by a Poisson process) where the parameter of poisson distribution is the pixel value. Later you need to add that mask, with the variation above the mean (pixel value) to the real image... – Haerr

For all you python noobs (such as myself), you can use noisyImage = numpy.clip(tempImage, 0.0, 255.0) to clamp – Otey

Answer 4

You could use skimage.util.random_noise:

from skimage.util import random_noise

noisy = random_noise(img, mode="poisson")

Dysarthria

Answer 5

From the item 1.4.4 - "Gaussian Approximation of the Poisson Distribution" of Chapter 1 of this book:

For large mean values, the Poisson distribution is well approximated by a Gaussian distribution with mean and variance equal to the mean of the Poisson random variable:

P(μ) ≈ N (μ,μ)

Then, we can generate Poisson noise from a normal distribution N (0,1), scale its standard deviation by the square root of μ and add it to the image which is the μ value:

# Image size
M, N = 1000, 1000

# Generate synthetic image
image = np.tile(np.arange(0,N,dtype='float64'),(M,1)) * 20

# -- sqrt(mu) * normal(0,1) --
poisson_noise = np.sqrt(image) * np.random.normal(0, 1, image.shape)

# Add the noise to the mu values
noisy_image = image + poisson_noise

plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.title('Image')
plt.imshow(image,'gray')
plt.subplot(2,2,2)
plt.title('Noisy image noise')
plt.imshow(noisy_image,'gray')  
plt.subplot(2,2,3)
plt.title('Image profile')
plt.plot(image[0,:])
plt.subplot(2,2,4)
plt.title('Noisy image profile')
plt.plot(noisy_image[0,:])

print("Synthetic image mean: {}".format(image[:,1].mean()))
print("Synthetic image variance: {}".format(image[:,1].var()))    
print("Noisy image mean: {}".format(noisy_image[:,1].mean()))
print("Noisy image variance: {}".format(noisy_image[:,1].var()))

As Poisson noise is signal-dependent, as we increase the underlying signal the noise variance also increases, as we can see in this row profiles:

Output for statistics in a single column:

Synthetic image mean: 20.0

Synthetic image variance: 0.0

Noisy image mean: 19.931120555821597

Noisy image variance: 19.39456713877459

Further references: [1][2]

Underage