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