diff --git a/skimage/util/noise.py b/skimage/util/noise.py index d0a3c5f6..9283f537 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. @@ -17,6 +17,8 @@ def random_noise(image, mode='gaussian', seed=None, **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. @@ -26,12 +28,20 @@ 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. 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 @@ -42,17 +52,51 @@ 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, 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 + 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) allowedtypes = { 'gaussian': 'gaussian_values', - 'poisson': '', + 'localvar': 'localvar_values', + 'poisson': 'poisson_values', 'salt': 'sp_values', 'pepper': 'sp_values', 's&p': 's&p_values', @@ -62,12 +106,15 @@ def random_noise(image, mode='gaussian', seed=None, **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']} + 's&p_values': ['amount', 'salt_vs_pepper'], + 'poisson_values': []} for key in kwargs: if key not in allowedkwargs[allowedtypes[mode]]: @@ -81,16 +128,32 @@ 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 == '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)) + 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.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) + + # 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) @@ -119,11 +182,15 @@ def random_noise(image, mode='gaussian', seed=None, **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, 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 87005d60..39477cb0 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 * 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(): @@ -72,10 +83,35 @@ 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 - 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 +120,78 @@ 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) + 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, np.clip(expected, 0., 1.)) + assert_allclose(cam_noisy2, expected) + + +def test_clip_poisson(): + seed = 42 + 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) + 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 = 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) + 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 = 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) + 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')