ENH: Add optional clip kwarg, docs, & handling of signed inputs.

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

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