From de42ba831a125db2913545a364fd14989cb5db83 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Fri, 11 Oct 2013 01:53:15 -0500
Subject: [PATCH] FIX: Fix and improve Poisson random noise generator

The Poissson generator now works.

The improved Poisson generator now infers the bit depth of the image
after conversion to a floating point image, by analyzing the unique
values present and finding the next power of two. This value is then
used to scale the floating point image up, after which Poisson
noise is generated, and then image is then scaled back down.
---
 skimage/util/noise.py                   | 37 +++++++++++++++++++------
 skimage/util/tests/test_random_noise.py | 12 +++++++-
 2 files changed, 40 insertions(+), 9 deletions(-)

diff --git a/skimage/util/noise.py b/skimage/util/noise.py
index d0a3c5f6..5238a2e6 100644
--- a/skimage/util/noise.py
+++ b/skimage/util/noise.py
@@ -5,6 +5,25 @@ from .dtype import img_as_float
 __all__ = ['random_noise']
 
 
+def _next_pow2(n):
+    """
+    Returns next integer power of two.
+    """
+
+    next_pow = 0
+    if n == np.inf:
+        return np.inf
+
+    if n < 1:
+        raise ValueError("Unable to determine next power of two for %i" % (n))
+
+    while True:
+        if 2 ** next_pow >= n:
+            return next_pow
+        else:
+            next_pow += 1
+
+
 def random_noise(image, mode='gaussian', seed=None, **kwargs):
     """
     Function to add random noise of various types to a floating-point image.
@@ -52,7 +71,7 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs):
 
     allowedtypes = {
         'gaussian': 'gaussian_values',
-        'poisson': '',
+        'poisson': 'poisson_values',
         'salt': 'sp_values',
         'pepper': 'sp_values',
         's&p': 's&p_values',
@@ -67,7 +86,8 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs):
     allowedkwargs = {
         'gaussian_values': ['mean', 'var'],
         'sp_values': ['amount'],
-        's&p_values': ['amount', 'salt_vs_pepper']}
+        's&p_values': ['amount', 'salt_vs_pepper'],
+        'poisson_values': []}
 
     for key in kwargs:
         if key not in allowedkwargs[allowedtypes[mode]]:
@@ -84,13 +104,14 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs):
         out = np.clip(image + noise, 0., 1.)
 
     elif mode == 'poisson':
+        # Determine unique values present in image
+        vals = len(np.unique(image))
+
+        # Calculate the next lowest power of two
+        vals = 2 ** _next_pow2(vals)
+
         # Generating noise for each unique value in image.
-        out = np.zeros_like(image)
-        for val in np.unique(image):
-            # Generate mask for a unique value, replace w/values drawn from
-            # Poisson distribution about the unique value
-            mask = image == val
-            out[mask] = np.poisson(val, mask.sum())
+        out = np.random.poisson(image * vals) / float(vals)
 
     elif mode == 'salt':
         # Re-call function with mode='s&p' and p=1 (all salt noise)
diff --git a/skimage/util/tests/test_random_noise.py b/skimage/util/tests/test_random_noise.py
index 87005d60..fd05084a 100644
--- a/skimage/util/tests/test_random_noise.py
+++ b/skimage/util/tests/test_random_noise.py
@@ -75,7 +75,7 @@ def test_gaussian():
 def test_speckle():
     seed = 42
     data = np.zeros((128, 128)) + 0.1
-    np.random.seed(seed=42)
+    np.random.seed(seed=seed)
     noise = np.random.normal(0.1, 0.02 ** 0.5, (128, 128))
     expected = np.clip(data + data * noise, 0, 1)
 
@@ -84,6 +84,16 @@ def test_speckle():
     assert_allclose(expected, data_speckle)
 
 
+def test_poisson():
+    seed = 42
+    data = camera()  # 512x512 grayscale uint8
+    cam_noisy = random_noise(data, mode='poisson', seed=seed)
+
+    np.random.seed(seed=seed)
+    expected = np.random.poisson(img_as_float(data) * 256) / 256.
+    assert_allclose(cam_noisy, expected)
+
+
 def test_bad_mode():
     data = np.zeros((64, 64))
     assert_raises(KeyError, random_noise, data, 'perlin')