From f53a4e0764583a4f0f074efb466db4cf242b204c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Tue, 16 Jul 2013 12:40:52 +0200 Subject: [PATCH] unwrap: Add naive 1D unwrapper. The naive 1D unwrapper does not support masked arrays because the 1D unwrapping problem has an infite number of solutions when faced with missing data. Wrap around is not implemented because 1D phase unwrapping must start at a certain pixel, and there will always be a risk of a discontinuity there, wrap around or not. --- bento.info | 3 +++ skimage/exposure/_unwrap_1d.pyx | 22 ++++++++++++++++++ skimage/exposure/setup.py | 3 +++ skimage/exposure/tests/test_unwrap.py | 12 +++++++++- skimage/exposure/unwrap.py | 33 ++++++++++++++++++++------- 5 files changed, 64 insertions(+), 9 deletions(-) create mode 100644 skimage/exposure/_unwrap_1d.pyx diff --git a/bento.info b/bento.info index e8ae3b07..30367057 100644 --- a/bento.info +++ b/bento.info @@ -150,6 +150,9 @@ Library: Extension: skimage.exposure._unwrap_2d Sources: skimage/exposure/_unwrap_2d.pyx, skimage/exposure/unwrap_2d_ljmu.c + Extension: skimage.exposure._unwrap_1d + Sources: + skimage/exposure/_unwrap_1d.pyx Executable: skivi Module: skimage.scripts.skivi diff --git a/skimage/exposure/_unwrap_1d.pyx b/skimage/exposure/_unwrap_1d.pyx new file mode 100644 index 00000000..d39f9e57 --- /dev/null +++ b/skimage/exposure/_unwrap_1d.pyx @@ -0,0 +1,22 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +from libc.math cimport M_PI + + +def unwrap_1d(float[::1] image, float[::1] unwrapped_image): + '''Phase unwrapping using the naive approach.''' + cdef: + Py_ssize_t i + float difference + long periods = 0 + unwrapped_image[0] = image[0] + for i in range(1, image.shape[0]): + difference = image[i] - image[i - 1] + if difference > M_PI: + periods -= 1 + elif difference < -M_PI: + periods += 1 + unwrapped_image[i] = image[i] + 2 * M_PI * periods diff --git a/skimage/exposure/setup.py b/skimage/exposure/setup.py index d21f41f0..1e801088 100644 --- a/skimage/exposure/setup.py +++ b/skimage/exposure/setup.py @@ -13,9 +13,12 @@ def configuration(parent_package='', top_path=None): config = Configuration('exposure', parent_package, top_path) config.add_data_dir('tests') + cython(['_unwrap_1d.pyx'], working_path=base_path) cython(['_unwrap_2d.pyx'], working_path=base_path) cython(['_unwrap_3d.pyx'], working_path=base_path) + config.add_extension('_unwrap_1d', sources=['_unwrap_1d.c'], + include_dirs=[get_numpy_include_dirs()]) unwrap_sources_2d = ['_unwrap_2d.c', 'unwrap_2d_ljmu.c'] config.add_extension('_unwrap_2d', sources=unwrap_sources_2d, extra_compile_args=['-g'], diff --git a/skimage/exposure/tests/test_unwrap.py b/skimage/exposure/tests/test_unwrap.py index 1e1b0cd5..220bb177 100644 --- a/skimage/exposure/tests/test_unwrap.py +++ b/skimage/exposure/tests/test_unwrap.py @@ -2,7 +2,8 @@ from __future__ import print_function, division import numpy as np from numpy.testing import (run_module_suite, assert_array_almost_equal, - assert_almost_equal, assert_array_equal) + assert_almost_equal, assert_array_equal, + assert_raises) import warnings from skimage.exposure import unwrap_phase @@ -40,6 +41,15 @@ def check_unwrap(image, mask=None): assert_phase_almost_equal(image_unwrapped, image) +def test_unwrap_1d(): + image = np.linspace(0, 10 * np.pi, 100) + check_unwrap(image) + # Masked arrays are not allowed in 1D + assert_raises(ValueError, check_unwrap, image, True) + # wrap_around is not allowed in 1D + assert_raises(ValueError, unwrap_phase, image, True) + + def test_unwrap_2d(): x, y = np.ogrid[:8, :16] image = 2 * np.pi * (x * 0.2 + y * 0.1) diff --git a/skimage/exposure/unwrap.py b/skimage/exposure/unwrap.py index a38162bd..a2a8c42a 100644 --- a/skimage/exposure/unwrap.py +++ b/skimage/exposure/unwrap.py @@ -1,6 +1,7 @@ import numpy as np import warnings +from ._unwrap_1d import unwrap_1d from ._unwrap_2d import unwrap_2d from ._unwrap_3d import unwrap_3d from .._shared.six import string_types @@ -12,16 +13,18 @@ def unwrap_phase(image, wrap_around=False): Parameters ---------- - image : 2D or 3D ndarray of floats, optionally a masked array + image : 1D, 2D or 3D ndarray of floats, optionally a masked array The values should be in the range ``[-pi, pi)``. If a masked array is provided, the masked entries will not be changed, and their values will not be used to guide the unwrapping of neighboring, unmasked - values. + values. Masked 1D arrays are not allowed, and will raise a + ``ValueError``. wrap_around : bool or sequence of bool When an element of the sequence is ``True``, the unwrapping process will regard the edges along the corresponding axis of the image to be connected and use this connectivity to guide the phase unwrapping process. If only a single boolean is given, it will apply to all axes. + Wrap around is not supported for 1D arrays. Returns ------- @@ -29,6 +32,12 @@ def unwrap_phase(image, wrap_around=False): Unwrapped image of the same shape as the input. If the input ``image`` was a masked array, the mask will be preserved. + Raises + ------ + ValueError + If called with a masked 1D array or called with a 1D array and + ``wrap_around=True``. + Examples -------- >>> c0, c1 = np.ogrid[-1:1:128j, -1:1:128j] @@ -50,8 +59,8 @@ def unwrap_phase(image, wrap_around=False): C. Gorecki, & E. L. Novak (Eds.), Optical Metrology (2005) 32--40, International Society for Optics and Photonics. ''' - if image.ndim not in (2, 3): - raise ValueError('image must be 2 or 3 dimensional') + if image.ndim not in (1, 2, 3): + raise ValueError('image must be 1, 2 or 3 dimensional') if isinstance(wrap_around, bool): wrap_around = [wrap_around] * image.ndim elif (hasattr(wrap_around, '__getitem__') @@ -63,9 +72,15 @@ def unwrap_phase(image, wrap_around=False): else: raise ValueError('wrap_around must be a bool or a sequence with ' 'length equal to the dimensionality of image') - if image.ndim == 3 and 1 in image.shape: - warnings.warn('image is 3D and has a length 1 dimension; consider ' - 'using a 2D array to use the 2D unwrapping algorithm') + if image.ndim == 1: + if np.ma.isMaskedArray(image): + raise ValueError('1D masked images cannot be unwrapped') + if wrap_around[0]: + raise ValueError('wrap_around is not supported for 1D images') + if image.ndim in (2, 3) and 1 in image.shape: + warnings.warn('image has a length 1 dimension; consider using an ' + 'array of lower dimensionality to use a more efficient ' + 'algorithm') if np.ma.isMaskedArray(image): mask = np.require(image.mask, np.uint8, ['C']) @@ -74,7 +89,9 @@ def unwrap_phase(image, wrap_around=False): image_not_masked = np.asarray(image, dtype=np.float32, order='C') image_unwrapped = np.empty_like(image, dtype=np.float32, order='C') - if image.ndim == 2: + if image.ndim == 1: + unwrap_1d(image_not_masked, image_unwrapped) + elif image.ndim == 2: unwrap_2d(image_not_masked, mask, image_unwrapped, wrap_around) elif image.ndim == 3: