wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-05 17:55:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'emmanuelle/tv_denoise'

Browse Source
This commit is contained in:
Stefan van der Walt
2011-05-31 12:48:18 +02:00
parent 430821a910 a0f7eb51cd
commit 2d4fc6e480
5 changed files with 378 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files
scikits/image/filter/__init__.py
+1
View File
@@ -2,3 +2,4 @@ from lpi_filter import *
from ctmf import median_filter
from canny import canny
from edges import sobel, hsobel, vsobel, hprewitt, vprewitt, prewitt
from tv_denoise import tv_denoise
scikits/image/filter/tests/test_tv_denoise.py
+58
View File
@@ -0,0 +1,58 @@
import numpy as np
from numpy.testing import run_module_suite
import scikits.image.filter as F
class TestTvDenoise():
def test_tv_denoise_2d(self):
"""
Apply the TV denoising algorithm on the lena image provided
by scipy
"""
import scipy
# lena image
lena = scipy.lena().astype(np.float)
# add noise to lena
lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
# denoise
denoised_lena = F.tv_denoise(lena, weight=60.0)
# which dtype?
assert denoised_lena.dtype in [np.float, np.float32, np.float64]
from scipy import ndimage
grad = ndimage.morphological_gradient(lena, size=((3,3)))
grad_denoised = ndimage.morphological_gradient(denoised_lena, size=((3,3)))
# test if the total variation has decreased
assert np.sqrt((grad_denoised**2).sum()) < np.sqrt((grad**2).sum())
denoised_lena_int = F.tv_denoise(lena.astype(np.int32), \
weight=60.0, keep_type=True)
assert denoised_lena_int.dtype is np.dtype('int32')
def test_tv_denoise_3d(self):
"""
Apply the TV denoising algorithm on a 3D image representing
a sphere.
"""
x, y, z = np.ogrid[0:40, 0:40, 0:40]
mask = (x -22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
mask = 100 * mask.astype(np.float)
mask += 60
mask += 20*np.random.randn(*mask.shape)
mask[mask < 0] = 0
mask[mask > 255] = 255
res = F.tv_denoise(mask.astype(np.uint8), weight=100, keep_type=True)
assert res.std() < mask.std()
assert res.dtype is np.dtype('uint8')
res = F.tv_denoise(mask.astype(np.uint8), weight=100)
assert res.std() < mask.std()
assert res.dtype is not np.dtype('uint8')
# test wrong number of dimensions
a = np.random.random((8, 8, 8, 8))
try:
res = F.tv_denoise(a)
except ValueError:
pass
if __name__ == "__main__":
run_module_suite()
scikits/image/filter/tv_denoise.py
+273
View File
@@ -0,0 +1,273 @@
import numpy as np
def _tv_denoise_3d(im, eps=2.e-4, weight=100, keep_type=False, n_iter_max=200):
"""
Perform total-variation denoising on 3-D arrays
Parameters
----------
im: ndarray
3-D input data to be denoised
eps: float, optional
relative difference of the value of the cost function that determines
the stop criterion. The algorithm stops when
(E_(n-1) - E_n) < eps * E_0
weight: float, optional
denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
keep_type: bool, optional (False)
whether the output has the same dtype as the input array.
keep_type is False by default, and the dtype of the output
is np.float
n_iter_max: int, optional
maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
denoised array
Notes
-----
Rudin, Osher and Fatemi algorithm
Examples
---------
First build synthetic noisy data
>>> x, y, z = np.ogrid[0:40, 0:40, 0:40]
>>> mask = (x -22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
>>> mask = mask.astype(np.float)
>>> mask += 0.2*np.random.randn(*mask.shape)
>>> res = tv_denoise_3d(mask, weight=100)
"""
im_type = im.dtype
if im_type is not np.float:
im = im.astype(np.float)
px = np.zeros_like(im)
py = np.zeros_like(im)
pz = np.zeros_like(im)
gx = np.zeros_like(im)
gy = np.zeros_like(im)
gz = np.zeros_like(im)
d = np.zeros_like(im)
i = 0
while i < n_iter_max:
d = - px - py - pz
d[1:] += px[:-1]
d[:, 1:] += py[:, :-1]
d[:, :, 1:] += pz[:, :, :-1]
out = im + d
E = (d**2).sum()
gx[:-1] = np.diff(out, axis=0)
gy[:, :-1] = np.diff(out, axis=1)
gz[:, :, :-1] = np.diff(out, axis=2)
norm = np.sqrt(gx**2 + gy**2 + gz**2)
E += weight * norm.sum()
norm *= 0.5 / weight
norm += 1.
px -= 1./6.*gx
px /= norm
py -= 1./6.*gy
py /= norm
pz -= 1/6.*gz
pz /= norm
E /= float(im.size)
print E
if i == 0:
E_init = E
E_previous = E
else:
if np.abs(E_previous - E) < eps * E_init:
print E_previous, E
break
else:
E_previous = E
i += 1
if keep_type:
return out.astype(im_type)
else:
return out
def _tv_denoise_2d(im, weight=50, eps=2.e-4, keep_type=False, n_iter_max=200):
"""
Perform total-variation denoising
Parameters
----------
im: ndarray
input data to be denoised
eps: float, optional
relative difference of the value of the cost function that determines
the stop criterion. The algorithm stops when
(E_(n-1) - E_n) < eps * E_0
weight: float, optional
denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
keep_type: bool, optional (False)
whether the output has the same dtype as the input array.
keep_type is False by default, and the dtype of the output
is np.float
n_iter_max: int, optional
maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
denoised array
Notes
-----
The principle of total variation denoising is explained in
http://en.wikipedia.org/wiki/Total_variation_denoising
This code is an implementation of the algorithm of Rudin, Fatemi and Osher
that was proposed by Chambolle in [1]_.
References
----------
.. [1] A. Chambolle, An algorithm for total variation minimization and
applications, Journal of Mathematical Imaging and Vision,
Springer, 2004, 20, 89-97.
Examples
---------
>>> import scipy
>>> lena = scipy.lena()
>>> import scipy
>>> lena = scipy.lena().astype(np.float)
>>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
>>> denoised_lena = tv_denoise(lena, weight=60.0)
"""
im_type = im.dtype
if im_type is not np.float:
im = im.astype(np.float)
px = np.zeros_like(im)
py = np.zeros_like(im)
gx = np.zeros_like(im)
gy = np.zeros_like(im)
d = np.zeros_like(im)
i = 0
while i < n_iter_max:
d = -px -py
d[1:] += px[:-1]
d[:, 1:] += py[:, :-1]
out = im + d
E = (d**2).sum()
gx[:-1] = np.diff(out, axis=0)
gy[:, :-1] = np.diff(out, axis=1)
norm = np.sqrt(gx**2 + gy**2)
E += weight * norm.sum()
norm *= 0.5 / weight
norm += 1
px -= 0.25*gx
px /= norm
py -= 0.25*gy
py /= norm
E /= float(im.size)
print E
if i == 0:
E_init = E
E_previous = E
else:
if np.abs(E_previous - E) < eps * E_init:
break
else:
E_previous = E
i += 1
print i
if keep_type:
return out.astype(im_type)
else:
return out
def tv_denoise(im, eps=2.e-4, weight=50, keep_type=False, n_iter_max=200):
"""
Perform total-variation denoising on 2-d and 3-d images
Parameters
----------
im: ndarray (2d or 3d) of ints, uints or floats
input data to be denoised. `im` can be of any numeric type,
but it is cast into an ndarray of floats for the computation
of the denoised image.
eps: float, optional
relative difference of the value of the cost function that
determines the stop criterion. The algorithm stops when
(E_(n-1) - E_n) < eps * E_0
weight: float, optional
denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
keep_type: bool, optional (False)
whether the output has the same dtype as the input array.
keep_type is False by default, and the dtype of the output
is np.float
n_iter_max: int, optional
maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
denoised array
Notes
-----
The principle of total variation denoising is explained in
http://en.wikipedia.org/wiki/Total_variation_denoising
The principle of total variation denoising is to minimize the
total variation of the image, which can be roughly described as
the integral of the norm of the image gradient. Total variation
denoising tends to produce "cartoon-like" images, that is,
piecewise-constant images.
This code is an implementation of the algorithm of Rudin, Fatemi and Osher
that was proposed by Chambolle in [1]_.
References
----------
.. [1] A. Chambolle, An algorithm for total variation minimization and
applications, Journal of Mathematical Imaging and Vision,
Springer, 2004, 20, 89-97.
Examples
---------
>>> import scipy
>>> # 2D example using lena
>>> lena = scipy.lena()
>>> import scipy
>>> lena = scipy.lena().astype(np.float)
>>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
>>> denoised_lena = tv_denoise(lena, weight=60)
>>> # 3D example on synthetic data
>>> x, y, z = np.ogrid[0:40, 0:40, 0:40]
>>> mask = (x -22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
>>> mask = mask.astype(np.float)
>>> mask += 0.2*np.random.randn(*mask.shape)
>>> res = tv_denoise_3d(mask, weight=100)
"""
if im.ndim == 2:
return _tv_denoise_2d(im, eps, weight, keep_type, n_iter_max)
elif im.ndim == 3:
return _tv_denoise_3d(im, eps, weight, keep_type, n_iter_max)
else:
raise ValueError('only 2-d and 3-d images may be denoised with this function')
scikits/image/graph/heap.pyx
+21 -19
View File
@@ -127,20 +127,10 @@ cdef class BinaryHeap:
#
# To calculate the capacity at a certain level:
# 2**l
def __init__(self, int initial_capcity=128):
"""__init__(initial_capacity=128)
Class constructor.
Takes an optional parameter 'initial_capacity' so that
if the required heap capacity is known or can be estimated in advance,
there will need to be fewer resize operations on the heap."""
def __cinit__(self, int initial_capacity=128, *args, **kws):
# calc levels from the default capacity
cdef int levels = 0
while 2**levels < initial_capcity:
while 2**levels < initial_capacity:
levels += 1
# set levels
self.min_levels = self.levels = levels
@@ -150,10 +140,19 @@ cdef class BinaryHeap:
# allocate arrays
cdef int number = 2**self.levels
cdef VALUE_T *values
values = self._values = <VALUE_T *>malloc( 2*number * sizeof(VALUE_T))
self._values = <VALUE_T *>malloc( 2*number * sizeof(VALUE_T))
self._references = <REFERENCE_T *>malloc(number * sizeof(REFERENCE_T))
def __init__(self, int initial_capacity=128):
"""__init__(initial_capacity=128)
Class constructor.
Takes an optional parameter 'initial_capacity' so that
if the required heap capacity is known or can be estimated in advance,
there will need to be fewer resize operations on the heap."""
if self._values is NULL or self._references is NULL:
raise MemoryError()
self.reset()
def reset(self):
@@ -512,19 +511,22 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
is lower than the current value in the heap. This is again useful for
pathfinding algorithms.
"""
def __cinit__(self, int initial_capacity=128, max_reference=None):
if max_reference is None:
max_reference = initial_capacity - 1
self.max_reference = max_reference
self._crossref = <REFERENCE_T *>malloc((max_reference+1) * \
sizeof(REFERENCE_T))
def __init__(self, int initial_capacity=128, max_reference=None):
"""__init__(initial_capacity=128, max_reference=None)
Class constructor.
"""
if max_reference is None:
max_reference = initial_capacity - 1
self.max_reference = max_reference
self._crossref = <REFERENCE_T *>malloc((max_reference+1) * \
sizeof(REFERENCE_T))
# below will call self.reset
BinaryHeap.__init__(self, initial_capacity)
def __dealloc__(self):
if self._crossref is not NULL:
free(self._crossref)
scikits/image/io/_plugins/freeimage_plugin.py
+25 -29
View File
@@ -16,7 +16,7 @@ if 'HOME' in os.environ:
lib_dirs.append(os.path.join(os.environ['HOME'], 'lib'))
API = {
'FreeImage_Load': (ctypes.c_voidp,
'FreeImage_Load': (ctypes.c_void_p,
[ctypes.c_int, ctypes.c_char_p, ctypes.c_int]),
'FreeImage_GetWidth': (ctypes.c_uint,
[ctypes.c_void_p]),
@@ -263,10 +263,12 @@ def read_multipage(filename, flags=0):
if not multibitmap:
raise ValueError('Could not open %s as multi-page image.' % filename)
try:
multibitmap = ctypes.c_void_p(multibitmap)
pages = _FI.FreeImage_GetPageCount(multibitmap)
arrays = []
for i in range(pages):
bitmap = _FI.FreeImage_LockPage(multibitmap, i)
bitmap = ctypes.c_void_p(bitmap)
try:
arrays.append(_array_from_bitmap(bitmap))
finally:
@@ -284,36 +286,28 @@ def _read_bitmap(filename, flags):
bitmap = _FI.FreeImage_Load(ftype, filename, flags)
if not bitmap:
raise ValueError('Could not load file %s' % filename)
return bitmap
return ctypes.c_void_p(bitmap)
def _wrap_bitmap_bits_in_array(bitmap, shape, dtype):
"""Return an ndarray view on the data in a FreeImage bitmap. Only
valid for as long as the bitmap is loaded (if single page) / locked
in memory (if multipage).
"""Return an ndarray view on the data in a FreeImage bitmap. Only
valid for as long as the bitmap is loaded (if single page) / locked
in memory (if multipage).
"""
pitch = _FI.FreeImage_GetPitch(bitmap)
height = shape[-1]
itemsize = dtype.itemsize
"""
pitch = _FI.FreeImage_GetPitch(bitmap)
height = shape[-1]
byte_size = height * pitch
itemsize = dtype.itemsize
if len(shape) == 3:
strides = (itemsize, shape[0]*itemsize, pitch)
else:
strides = (itemsize, pitch)
bits = _FI.FreeImage_GetBits(bitmap)
class DummyArray:
__array_interface__ = {
'data': (bits, False),
'strides': strides,
'typestr': dtype.str,
'shape': tuple(shape),
'version' : 3,
}
# Still segfaulting on 64-bit machine because of illegal memory access
return numpy.array(DummyArray(), copy=False)
if len(shape) == 3:
strides = (itemsize, shape[0]*itemsize, pitch)
else:
strides = (itemsize, pitch)
bits = _FI.FreeImage_GetBits(bitmap)
array = numpy.ndarray(shape, dtype=dtype,
buffer=(ctypes.c_char*byte_size).from_address(bits),
strides=strides)
return array
def _array_from_bitmap(bitmap):
"""Convert a FreeImage bitmap pointer to a numpy array
@@ -399,8 +393,9 @@ def write_multipage(arrays, filename, flags=0):
raise ValueError('Could not open %s for writing multi-page image.' %
filename)
try:
multibitmap = ctypes.c_void_p(multibitmap)
for array in arrays:
bitmap = _array_to_bitmap(array)
bitmap, fi_type = _array_to_bitmap(array)
_FI.FreeImage_AppendPage(multibitmap, bitmap)
finally:
_FI.FreeImage_CloseMultiBitmap(multibitmap, flags)
@@ -436,6 +431,7 @@ def _array_to_bitmap(array):
try:
def n(arr): # normalise to freeimage's in-memory format
return arr.T[:,::-1]
bitmap = ctypes.c_void_p(bitmap)
wrapped_array = _wrap_bitmap_bits_in_array(bitmap, w_shape, dtype)
# swizzle the color components and flip the scanlines to go to
# FreeImage's BGR[A] and upside-down internal memory format