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')