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

@@ -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

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

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'],

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)

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: