wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-29 20:02:22 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #772 from JDWarner/fix_poisson_noise

Browse Source 
FIX: Fix and improve Poisson random noise generator
This commit is contained in:
Johannes Schönberger
2013-10-13 10:48:03 -07:00
parent f3845e7e2f 1c5dc10f4d
commit 4b0cf8d27a
2 changed files with 192 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/util/noise.py
+81 -14
View File
@@ -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
skimage/util/tests/test_random_noise.py
+111 -3
View File
@@ -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')