From de42ba831a125db2913545a364fd14989cb5db83 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 11 Oct 2013 01:53:15 -0500 Subject: [PATCH 1/6] 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') From c2e4442eaff871e4ac205dcc35e4b13546aeb06a Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 11 Oct 2013 10:51:07 -0500 Subject: [PATCH 2/6] ENH: More concise next-power-of-2 calculation --- skimage/util/noise.py | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/skimage/util/noise.py b/skimage/util/noise.py index 5238a2e6..7d7c28d7 100644 --- a/skimage/util/noise.py +++ b/skimage/util/noise.py @@ -5,25 +5,6 @@ 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. @@ -108,7 +89,7 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): vals = len(np.unique(image)) # Calculate the next lowest power of two - vals = 2 ** _next_pow2(vals) + vals = 2 ** np.ceil(np.log2(vals)) # Generating noise for each unique value in image. out = np.random.poisson(image * vals) / float(vals) From f435afebd5453e679f1eb0325914cb3747215709 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 11 Oct 2013 16:39:22 -0500 Subject: [PATCH 3/6] ENH: Add optional `clip` kwarg, docs, & handling of signed inputs. --- skimage/util/noise.py | 62 ++++++++++++++++++++++--- skimage/util/tests/test_random_noise.py | 60 ++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 7 deletions(-) diff --git a/skimage/util/noise.py b/skimage/util/noise.py index 7d7c28d7..db470498 100644 --- a/skimage/util/noise.py +++ b/skimage/util/noise.py @@ -5,7 +5,7 @@ from .dtype import img_as_float __all__ = ['random_noise'] -def random_noise(image, mode='gaussian', seed=None, **kwargs): +def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): """ Function to add random noise of various types to a floating-point image. @@ -26,6 +26,11 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): seed : int If provided, this will set the random seed before generating noise, for valid pseudo-random comparisons. + clip : bool + If True (default), the output will be clipped after noise applied + for modes `'speckle'`, `'poisson'`, and `'gaussian'`. This is + needed to maintain the proper image data range. If False, clipping + is not applied, and the output may extend beyond the range [-1, 1]. mean : float Mean of random distribution. Used in 'gaussian' and 'speckle'. Default : 0. @@ -42,10 +47,42 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): Returns ------- out : ndarray - Output floating-point image data on range [0, 1]. + Output floating-point image data on range [0, 1] or [-1, 1] if the + input `image` was unsigned or signed, respectively. + + Notes + ----- + Speckle, Poisson, and Gaussian noise may generate noise outside the valid + image range. The default is to clip (not alias) these values, but they may + be preserved by setting `clip=False`. Note that in this case the output + may contain values outside the ranges [0, 1] or [-1, 1]. Use with care. + + Because of the prevalence of exclusively positive floating-point images in + intermediate calculations, it is not possible to intuit if an input is + signed based on dtype alone. Instead, negative values are explicity + searched for. Only if found does this function assume signed input. + Unexpected results only occur in rare, poorly exposes cases (e.g. if all + values are above 50 percent gray in a signed `image`). In this event, + manually scaling the input to the positive domain will solve the problem. + + The Poisson distribution is only defined for positive integers. To apply + this noise type, the number of unique values in the image is found and + the next round power of two is used to scale up the floating-point result, + after which it is scaled back down to the floating-point image range. + + To generate Poisson noise against a signed image, the signed image is + temporarily converted to an unsigned image in the floating point domain, + Poisson noise is generated, then it is returned to the original range. """ mode = mode.lower() + + # Detect if a signed image was input + if image.min() < 0: + low_clip = -1. + else: + low_clip = 0. + image = img_as_float(image) if seed is not None: np.random.seed(seed=seed) @@ -82,18 +119,25 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): if mode == 'gaussian': noise = np.random.normal(kwargs['mean'], kwargs['var'] ** 0.5, image.shape) - out = np.clip(image + noise, 0., 1.) + out = image + noise elif mode == 'poisson': - # Determine unique values present in image + # Determine unique values in image & calculate the next power of two vals = len(np.unique(image)) - - # Calculate the next lowest power of two vals = 2 ** np.ceil(np.log2(vals)) + # Ensure image is exclusively positive + if low_clip == -1.: + old_max = image.max() + image = (image + 1.) / (old_max + 1.) + # Generating noise for each unique value in image. out = np.random.poisson(image * vals) / float(vals) + # Return image to original range if input was signed + if low_clip == -1.: + out = out * (old_max + 1.) - 1. + elif mode == 'salt': # Re-call function with mode='s&p' and p=1 (all salt noise) out = random_noise(image, mode='s&p', seed=seed, @@ -126,6 +170,10 @@ def random_noise(image, mode='gaussian', seed=None, **kwargs): elif mode == 'speckle': noise = np.random.normal(kwargs['mean'], kwargs['var'] ** 0.5, image.shape) - out = np.clip(image + image * noise, 0., 1.) + out = image + image * noise + + # Clip back to original range, if necessary + if clip: + out = np.clip(out, low_clip, 1.0) return out diff --git a/skimage/util/tests/test_random_noise.py b/skimage/util/tests/test_random_noise.py index fd05084a..f5c9f9d0 100644 --- a/skimage/util/tests/test_random_noise.py +++ b/skimage/util/tests/test_random_noise.py @@ -94,6 +94,66 @@ def test_poisson(): assert_allclose(cam_noisy, expected) +def test_clip_poisson(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + data_signed = (data / 255.) * 2. - 1. # Same image, on range [-1, 1] + + # Signed and unsigned, clipped + cam_poisson = random_noise(data, mode='poisson', seed=seed, clip=True) + cam_poisson2 = random_noise(data_signed, mode='poisson', seed=seed, + clip=True) + assert (cam_poisson.max() == 1.) and (cam_poisson.min() == 0.) + assert (cam_poisson2.max() == 1.) and (cam_poisson2.min() == -1.) + + # Signed and unsigned, unclipped + cam_poisson = random_noise(data, mode='poisson', seed=seed, clip=False) + cam_poisson2 = random_noise(data_signed, mode='poisson', seed=seed, + clip=False) + assert (cam_poisson.max() > 1.15) and (cam_poisson.min() == 0.) + assert (cam_poisson2.max() > 1.3) and (cam_poisson2.min() == -1.) + + +def test_clip_gaussian(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + data_signed = (data / 255.) * 2. - 1. # Same image, on range [-1, 1] + + # Signed and unsigned, clipped + cam_gauss = random_noise(data, mode='gaussian', seed=seed, clip=True) + cam_gauss2 = random_noise(data_signed, mode='gaussian', seed=seed, + clip=True) + assert (cam_gauss.max() == 1.) and (cam_gauss.min() == 0.) + assert (cam_gauss2.max() == 1.) and (cam_gauss2.min() == -1.) + + # Signed and unsigned, unclipped + cam_gauss = random_noise(data, mode='gaussian', seed=seed, clip=False) + cam_gauss2 = random_noise(data_signed, mode='gaussian', seed=seed, + clip=False) + assert (cam_gauss.max() > 1.22) and (cam_gauss.min() < -0.36) + assert (cam_gauss2.max() > 1.219) and (cam_gauss2.min() < -1.337) + + +def test_clip_speckle(): + seed = 42 + data = camera() # 512x512 grayscale uint8 + data_signed = (data / 255.) * 2. - 1. # Same image, on range [-1, 1] + + # Signed and unsigned, clipped + cam_speckle = random_noise(data, mode='speckle', seed=seed, clip=True) + cam_speckle2 = random_noise(data_signed, mode='speckle', seed=seed, + clip=True) + assert (cam_speckle.max() == 1.) and (cam_speckle.min() == 0.) + assert (cam_speckle2.max() == 1.) and (cam_speckle2.min() == -1.) + + # Signed and unsigned, unclipped + cam_speckle = random_noise(data, mode='speckle', seed=seed, clip=False) + cam_speckle2 = random_noise(data_signed, mode='speckle', seed=seed, + clip=False) + assert (cam_speckle.max() > 1.219) and (cam_speckle.min() == 0.) + assert (cam_speckle2.max() > 1.219) and (cam_speckle2.min() < -1.306) + + def test_bad_mode(): data = np.zeros((64, 64)) assert_raises(KeyError, random_noise, data, 'perlin') From 9f7a2f4fbc17f1fb1987d7a7bfd58dbfe6b49464 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Fri, 11 Oct 2013 16:53:16 -0500 Subject: [PATCH 4/6] ENH: Tighten tests, all noise types now support signed I/O --- skimage/util/noise.py | 2 +- skimage/util/tests/test_random_noise.py | 19 ++++++++++++++++--- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/skimage/util/noise.py b/skimage/util/noise.py index db470498..aa8af300 100644 --- a/skimage/util/noise.py +++ b/skimage/util/noise.py @@ -165,7 +165,7 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): kwargs['amount'] * image.size * (1. - kwargs['salt_vs_pepper'])) coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape] - out[coords] = 0 + out[coords] = low_clip elif mode == 'speckle': noise = np.random.normal(kwargs['mean'], kwargs['var'] ** 0.5, diff --git a/skimage/util/tests/test_random_noise.py b/skimage/util/tests/test_random_noise.py index f5c9f9d0..b1dece36 100644 --- a/skimage/util/tests/test_random_noise.py +++ b/skimage/util/tests/test_random_noise.py @@ -23,12 +23,14 @@ def test_salt(): # Ensure approximately correct amount of noise was added proportion = float(saltmask.sum()) / (cam.shape[0] * cam.shape[1]) - assert 0.11 < proportion <= 0.18 + assert 0.11 < proportion <= 0.15 def test_pepper(): seed = 42 cam = img_as_float(camera()) + data_signed = (cam / 255.) * 2. - 1. # Same image, on range [-1, 1] + cam_noisy = random_noise(cam, seed=seed, mode='pepper', amount=0.15) peppermask = cam != cam_noisy @@ -37,7 +39,16 @@ def test_pepper(): # Ensure approximately correct amount of noise was added proportion = float(peppermask.sum()) / (cam.shape[0] * cam.shape[1]) - assert 0.11 < proportion <= 0.18 + assert 0.11 < proportion <= 0.15 + + # Check to make sure pepper gets added properly to signed images + orig_zeros = (data_signed == -1).sum() + cam_noisy_signed = random_noise(data_signed, seed=seed, mode='pepper', + amount=.15) + + proportion = (float((cam_noisy_signed == -1).sum() - orig_zeros) / + (cam.shape[0] * cam.shape[1])) + assert 0.11 < proportion <= 0.15 def test_salt_and_pepper(): @@ -88,10 +99,12 @@ def test_poisson(): seed = 42 data = camera() # 512x512 grayscale uint8 cam_noisy = random_noise(data, mode='poisson', seed=seed) + cam_noisy2 = random_noise(data, mode='poisson', seed=seed, clip=False) np.random.seed(seed=seed) expected = np.random.poisson(img_as_float(data) * 256) / 256. - assert_allclose(cam_noisy, expected) + assert_allclose(cam_noisy, np.clip(expected, 0., 1.)) + assert_allclose(cam_noisy2, expected) def test_clip_poisson(): From f20aa5e66191cb80b1e04725b85b4bf08da5cf75 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sat, 12 Oct 2013 09:32:46 -0500 Subject: [PATCH 5/6] FIX: Better handling of skimage.data.camera for Python3 compatibility --- skimage/util/tests/test_random_noise.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/skimage/util/tests/test_random_noise.py b/skimage/util/tests/test_random_noise.py index b1dece36..4d1752a0 100644 --- a/skimage/util/tests/test_random_noise.py +++ b/skimage/util/tests/test_random_noise.py @@ -29,7 +29,7 @@ def test_salt(): def test_pepper(): seed = 42 cam = img_as_float(camera()) - data_signed = (cam / 255.) * 2. - 1. # Same image, on range [-1, 1] + data_signed = cam * 2. - 1. # Same image, on range [-1, 1] cam_noisy = random_noise(cam, seed=seed, mode='pepper', amount=0.15) peppermask = cam != cam_noisy @@ -109,8 +109,8 @@ def test_poisson(): def test_clip_poisson(): seed = 42 - data = camera() # 512x512 grayscale uint8 - data_signed = (data / 255.) * 2. - 1. # Same image, on range [-1, 1] + data = camera() # 512x512 grayscale uint8 + data_signed = img_as_float(data) * 2. - 1. # Same image, on range [-1, 1] # Signed and unsigned, clipped cam_poisson = random_noise(data, mode='poisson', seed=seed, clip=True) @@ -129,8 +129,8 @@ def test_clip_poisson(): def test_clip_gaussian(): seed = 42 - data = camera() # 512x512 grayscale uint8 - data_signed = (data / 255.) * 2. - 1. # Same image, on range [-1, 1] + data = camera() # 512x512 grayscale uint8 + data_signed = img_as_float(data) * 2. - 1. # Same image, on range [-1, 1] # Signed and unsigned, clipped cam_gauss = random_noise(data, mode='gaussian', seed=seed, clip=True) @@ -149,8 +149,8 @@ def test_clip_gaussian(): def test_clip_speckle(): seed = 42 - data = camera() # 512x512 grayscale uint8 - data_signed = (data / 255.) * 2. - 1. # Same image, on range [-1, 1] + data = camera() # 512x512 grayscale uint8 + data_signed = img_as_float(data) * 2. - 1. # Same image, on range [-1, 1] # Signed and unsigned, clipped cam_speckle = random_noise(data, mode='speckle', seed=seed, clip=True) From 1c5dc10f4d8c8add889c79b6f2a6aed6808fa799 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sun, 13 Oct 2013 12:21:29 -0500 Subject: [PATCH 6/6] FEAT: Add 'localvar' mode to random_noise --- skimage/util/noise.py | 27 ++++++++++++++++++++----- skimage/util/tests/test_random_noise.py | 25 +++++++++++++++++++++++ 2 files changed, 47 insertions(+), 5 deletions(-) diff --git a/skimage/util/noise.py b/skimage/util/noise.py index aa8af300..9283f537 100644 --- a/skimage/util/noise.py +++ b/skimage/util/noise.py @@ -17,6 +17,8 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): One of the following strings, selecting the type of noise to add: 'gaussian' Gaussian-distributed additive noise. + 'localvar' Gaussian-distributed additive noise, with specified + local variance at each point of `image` 'poisson' Poisson-distributed noise generated from the data. 'salt' Replaces random pixels with 1. 'pepper' Replaces random pixels with 0. @@ -37,6 +39,9 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): var : float Variance of random distribution. Used in 'gaussian' and 'speckle'. Note: variance = (standard deviation) ** 2. Default : 0.01 + local_vars : ndarray + Array of positive floats, same shape as `image`, defining the local + variance at every image point. Used in 'localvar'. amount : float Proportion of image pixels to replace with noise on range [0, 1]. Used in 'salt', 'pepper', and 'salt & pepper'. Default : 0.05 @@ -52,10 +57,11 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): Notes ----- - Speckle, Poisson, and Gaussian noise may generate noise outside the valid - image range. The default is to clip (not alias) these values, but they may - be preserved by setting `clip=False`. Note that in this case the output - may contain values outside the ranges [0, 1] or [-1, 1]. Use with care. + Speckle, Poisson, Localvar, and Gaussian noise may generate noise outside + the valid image range. The default is to clip (not alias) these values, + but they may be preserved by setting `clip=False`. Note that in this case + the output may contain values outside the ranges [0, 1] or [-1, 1]. + Use this option with care. Because of the prevalence of exclusively positive floating-point images in intermediate calculations, it is not possible to intuit if an input is @@ -89,6 +95,7 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): allowedtypes = { 'gaussian': 'gaussian_values', + 'localvar': 'localvar_values', 'poisson': 'poisson_values', 'salt': 'sp_values', 'pepper': 'sp_values', @@ -99,10 +106,12 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): 'mean': 0., 'var': 0.01, 'amount': 0.05, - 'salt_vs_pepper': 0.5} + 'salt_vs_pepper': 0.5, + 'local_vars': np.zeros_like(image) + 0.01} allowedkwargs = { 'gaussian_values': ['mean', 'var'], + 'localvar_values': ['local_vars'], 'sp_values': ['amount'], 's&p_values': ['amount', 'salt_vs_pepper'], 'poisson_values': []} @@ -121,6 +130,14 @@ def random_noise(image, mode='gaussian', seed=None, clip=True, **kwargs): image.shape) out = image + noise + elif mode == 'localvar': + # Ensure local variance input is correct + if (kwargs['local_vars'] <= 0).any(): + raise ValueError('All values of `local_vars` must be > 0.') + + # Safe shortcut usage broadcasts kwargs['local_vars'] as a ufunc + out = image + np.random.normal(0, kwargs['local_vars'] ** 0.5) + elif mode == 'poisson': # Determine unique values in image & calculate the next power of two vals = len(np.unique(image)) diff --git a/skimage/util/tests/test_random_noise.py b/skimage/util/tests/test_random_noise.py index 4d1752a0..39477cb0 100644 --- a/skimage/util/tests/test_random_noise.py +++ b/skimage/util/tests/test_random_noise.py @@ -83,6 +83,31 @@ def test_gaussian(): assert 0.012 < data_gaussian.var() < 0.018 +def test_localvar(): + seed = 42 + data = np.zeros((128, 128)) + 0.5 + local_vars = np.zeros((128, 128)) + 0.001 + local_vars[:64, 64:] = 0.1 + local_vars[64:, :64] = 0.25 + local_vars[64:, 64:] = 0.45 + + data_gaussian = random_noise(data, mode='localvar', seed=seed, + local_vars=local_vars, clip=False) + assert 0. < data_gaussian[:64, :64].var() < 0.002 + assert 0.095 < data_gaussian[:64, 64:].var() < 0.105 + assert 0.245 < data_gaussian[64:, :64].var() < 0.255 + assert 0.445 < data_gaussian[64:, 64:].var() < 0.455 + + # Ensure local variance bounds checking works properly + bad_local_vars = np.zeros_like(data) + assert_raises(ValueError, random_noise, data, mode='localvar', seed=seed, + local_vars=bad_local_vars) + bad_local_vars += 0.1 + bad_local_vars[0, 0] = -1 + assert_raises(ValueError, random_noise, data, mode='localvar', seed=seed, + local_vars=bad_local_vars) + + def test_speckle(): seed = 42 data = np.zeros((128, 128)) + 0.1