From 1c5dc10f4d8c8add889c79b6f2a6aed6808fa799 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 13 Oct 2013 12:21:29 -0500
Subject: [PATCH] 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