mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-17 11:32:45 +08:00
Merge pull request #644 from josteinbf/phase-unwrap
Add phase unwrapping
This commit is contained in:
Stefan van der Walt
+14
-2
@@ -144,8 +144,8 @@
|
||||
Color separation (color deconvolution) for several stainings.
|
||||
|
||||
- Jostein Bø Fløystad
|
||||
Reconstruction circle mode for Radon transform
|
||||
Simultaneous Algebraic Reconstruction Technique for inverse Radon transform
|
||||
Tomography: radon/iradon improvements and SART implementation
|
||||
Phase unwrapping integration
|
||||
|
||||
- Matt Terry
|
||||
Color difference functions
|
||||
@@ -161,3 +161,15 @@
|
||||
|
||||
- Michael Hansen
|
||||
novice submodule
|
||||
|
||||
- Munther Gdeisat
|
||||
Phase unwrapping implementation
|
||||
|
||||
- Miguel Arevallilo Herraez
|
||||
Phase unwrapping implementation
|
||||
|
||||
- Hussein Abdul-Rahman
|
||||
Phase unwrapping implementation
|
||||
|
||||
- Gregor Thalhammer
|
||||
Phase unwrapping integration
|
||||
|
||||
@@ -144,6 +144,15 @@ Library:
|
||||
Extension: skimage.filter.rank.bilateral_cy
|
||||
Sources:
|
||||
skimage/filter/rank/bilateral_cy.pyx
|
||||
Extension: skimage.exposure._unwrap_3d
|
||||
Sources:
|
||||
skimage/exposure/_unwrap_3d.pyx, skimage/exposure/unwrap_3d_ljmu.c
|
||||
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
|
||||
|
||||
@@ -0,0 +1,116 @@
|
||||
"""
|
||||
================
|
||||
Phase Unwrapping
|
||||
================
|
||||
|
||||
Some signals can only be observed modulo 2*pi, and this can also apply to
|
||||
two- and three dimensional images. In these cases phase unwrapping is
|
||||
needed to recover the underlying, unwrapped signal. In this example we will
|
||||
demonstrate an algorithm [1]_ implemented in ``skimage`` at work for such a
|
||||
problem. One-, two- and three dimensional images can all be unwrapped using
|
||||
skimage. Here we will demonstrate phase unwrapping in the two dimensional case.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
from skimage import data, img_as_float, color, exposure
|
||||
from skimage.exposure import unwrap_phase
|
||||
|
||||
|
||||
# Load an image as a floating-point grayscale
|
||||
image = color.rgb2gray(img_as_float(data.chelsea()))
|
||||
# Scale the image to [0, 4*pi]
|
||||
image = exposure.rescale_intensity(image, out_range=(0, 4 * np.pi))
|
||||
# Create a phase-wrapped image in the interval [-pi, pi)
|
||||
image_wrapped = np.angle(np.exp(1j * image))
|
||||
# Perform phase unwrapping
|
||||
image_unwrapped = unwrap_phase(image_wrapped)
|
||||
|
||||
plt.figure()
|
||||
plt.subplot(221)
|
||||
plt.title('Original')
|
||||
plt.imshow(image, cmap='gray', vmin=0, vmax=4 * np.pi)
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(222)
|
||||
plt.title('Wrapped phase')
|
||||
plt.imshow(image_wrapped, cmap='gray', vmin=-np.pi, vmax=np.pi)
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(223)
|
||||
plt.title('After phase unwrapping')
|
||||
plt.imshow(image_unwrapped, cmap='gray')
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(224)
|
||||
plt.title('Unwrapped minus original')
|
||||
plt.imshow(image_unwrapped - image, cmap='gray')
|
||||
plt.colorbar()
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
The unwrapping procedure accepts masked arrays, and can also optionally
|
||||
assume cyclic boundaries to connect edges of an image. In the example below,
|
||||
we study a simple phase ramp which has been split in two by masking
|
||||
a row of the image.
|
||||
"""
|
||||
|
||||
# Create a simple ramp
|
||||
image = np.ones((100, 100)) * np.linspace(0, 8 * np.pi, 100).reshape((-1, 1))
|
||||
# Mask the image to split it in two horizontally
|
||||
mask = np.zeros_like(image, dtype=np.bool)
|
||||
mask[image.shape[0] // 2, :] = True
|
||||
|
||||
image_wrapped = np.ma.array(np.angle(np.exp(1j * image)), mask=mask)
|
||||
# Unwrap image without wrap around
|
||||
image_unwrapped_no_wrap_around = unwrap_phase(image_wrapped,
|
||||
wrap_around=(False, False))
|
||||
# Unwrap with wrap around enabled for the 0th dimension
|
||||
image_unwrapped_wrap_around = unwrap_phase(image_wrapped,
|
||||
wrap_around=(True, False))
|
||||
|
||||
plt.figure()
|
||||
plt.subplot(221)
|
||||
plt.title('Original')
|
||||
plt.imshow(np.ma.array(image, mask=mask), cmap='jet')
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(222)
|
||||
plt.title('Wrapped phase')
|
||||
plt.imshow(image_wrapped, cmap='jet', vmin=-np.pi, vmax=np.pi)
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(223)
|
||||
plt.title('Unwrapped without wrap_around')
|
||||
plt.imshow(image_unwrapped_no_wrap_around, cmap='jet')
|
||||
plt.colorbar()
|
||||
|
||||
plt.subplot(224)
|
||||
plt.title('Unwrapped with wrap_around')
|
||||
plt.imshow(image_unwrapped_wrap_around, cmap='jet')
|
||||
plt.colorbar()
|
||||
|
||||
plt.show()
|
||||
|
||||
"""
|
||||
.. image:: PLOT2RST.current_figure
|
||||
|
||||
In the figures above, the masked row can be seen as a white line across
|
||||
the image. The difference between the two unwrapped images in the bottom row
|
||||
is clear: Without unwrapping (lower left), the regions above and below the
|
||||
masked boundary do not interact at all, resulting in an offset between the
|
||||
two regions of an arbitrary integer times two pi. We could just as well have
|
||||
unwrapped the regions as two separate images. With wrap around enabled for the
|
||||
vertical direction (lower rigth), the situation changes: Unwrapping paths are
|
||||
now allowed to pass from the bottom to the top of the image and vice versa, in
|
||||
effect providing a way to determine the offset between the two regions.
|
||||
|
||||
References
|
||||
----------
|
||||
|
||||
.. [1] Miguel Arevallilo Herraez, David R. Burton, Michael J. Lalor,
|
||||
and Munther A. Gdeisat, "Fast two-dimensional phase-unwrapping
|
||||
algorithm based on sorting by reliability following a noncontinuous
|
||||
path", Journal Applied Optics, Vol. 41, No. 35, pp. 7437, 2002
|
||||
"""
|
||||
@@ -3,6 +3,7 @@ from .exposure import histogram, equalize, equalize_hist, \
|
||||
adjust_gamma, adjust_sigmoid, adjust_log
|
||||
|
||||
from ._adapthist import equalize_adapthist
|
||||
from .unwrap import unwrap_phase
|
||||
|
||||
__all__ = ['histogram',
|
||||
'equalize',
|
||||
@@ -12,4 +13,5 @@ __all__ = ['histogram',
|
||||
'cumulative_distribution',
|
||||
'adjust_gamma',
|
||||
'adjust_sigmoid',
|
||||
'adjust_log']
|
||||
'adjust_log',
|
||||
'unwrap_phase']
|
||||
|
||||
@@ -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(double[::1] image, double[::1] unwrapped_image):
|
||||
'''Phase unwrapping using the naive approach.'''
|
||||
cdef:
|
||||
Py_ssize_t i
|
||||
double 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
|
||||
@@ -0,0 +1,16 @@
|
||||
cdef extern void unwrap2D(double* wrapped_image,
|
||||
double* unwrapped_image,
|
||||
unsigned char* input_mask,
|
||||
int image_width, int image_height,
|
||||
int wrap_around_x, int wrap_around_y)
|
||||
|
||||
def unwrap_2d(double[:, ::1] image,
|
||||
unsigned char[:, ::1] mask,
|
||||
double[:, ::1] unwrapped_image,
|
||||
wrap_around):
|
||||
unwrap2D(&image[0, 0],
|
||||
&unwrapped_image[0, 0],
|
||||
&mask[0, 0],
|
||||
image.shape[1], image.shape[0],
|
||||
wrap_around[1], wrap_around[0],
|
||||
)
|
||||
@@ -0,0 +1,16 @@
|
||||
cdef extern void unwrap3D(double* wrapped_volume,
|
||||
double* unwrapped_volume,
|
||||
unsigned char* input_mask,
|
||||
int image_width, int image_height, int volume_depth,
|
||||
int wrap_around_x, int wrap_around_y, int wrap_around_z)
|
||||
|
||||
def unwrap_3d(double[:, :, ::1] image,
|
||||
unsigned char[:, :, ::1] mask,
|
||||
double[:, :, ::1] unwrapped_image,
|
||||
wrap_around):
|
||||
unwrap3D(&image[0, 0, 0],
|
||||
&unwrapped_image[0, 0, 0],
|
||||
&mask[0, 0, 0],
|
||||
image.shape[2], image.shape[1], image.shape[0], #TODO: check!!!
|
||||
wrap_around[2], wrap_around[1], wrap_around[0],
|
||||
)
|
||||
@@ -0,0 +1,40 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
import os
|
||||
|
||||
from skimage._build import cython
|
||||
|
||||
base_path = os.path.abspath(os.path.dirname(__file__))
|
||||
|
||||
|
||||
def configuration(parent_package='', top_path=None):
|
||||
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
|
||||
|
||||
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,
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
unwrap_sources_3d = ['_unwrap_3d.c', 'unwrap_3d_ljmu.c']
|
||||
config.add_extension('_unwrap_3d', sources=unwrap_sources_3d,
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy.distutils.core import setup
|
||||
setup(maintainer='scikit-image Developers',
|
||||
author='scikit-image Developers',
|
||||
maintainer_email='scikit-image@googlegroups.com',
|
||||
description='Exposure corrections',
|
||||
url='https://github.com/scikit-image/scikit-image',
|
||||
license='SciPy License (BSD Style)',
|
||||
**(configuration(top_path='').todict())
|
||||
)
|
||||
@@ -0,0 +1,141 @@
|
||||
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_raises)
|
||||
import warnings
|
||||
|
||||
from skimage.exposure import unwrap_phase
|
||||
|
||||
|
||||
def assert_phase_almost_equal(a, b, *args, **kwargs):
|
||||
'''An assert_almost_equal insensitive to phase shifts of n*2*pi.'''
|
||||
shift = 2 * np.pi * np.round((b.mean() - a.mean()) / (2 * np.pi))
|
||||
with warnings.catch_warnings():
|
||||
warnings.simplefilter("ignore")
|
||||
print('assert_phase_allclose, abs', np.max(np.abs(a - (b - shift))))
|
||||
print('assert_phase_allclose, rel',
|
||||
np.max(np.abs((a - (b - shift)) / a)))
|
||||
if np.ma.isMaskedArray(a):
|
||||
assert np.ma.isMaskedArray(b)
|
||||
assert_array_equal(a.mask, b.mask)
|
||||
au = np.asarray(a)
|
||||
bu = np.asarray(b)
|
||||
with warnings.catch_warnings():
|
||||
warnings.simplefilter("ignore")
|
||||
print('assert_phase_allclose, no mask, abs',
|
||||
np.max(np.abs(au - (bu - shift))))
|
||||
print('assert_phase_allclose, no mask, rel',
|
||||
np.max(np.abs((au - (bu - shift)) / au)))
|
||||
assert_array_almost_equal(a + shift, b, *args, **kwargs)
|
||||
|
||||
|
||||
def check_unwrap(image, mask=None):
|
||||
image_wrapped = np.angle(np.exp(1j * image))
|
||||
if not mask is None:
|
||||
print('Testing a masked image')
|
||||
image = np.ma.array(image, mask=mask)
|
||||
image_wrapped = np.ma.array(image_wrapped, mask=mask)
|
||||
image_unwrapped = unwrap_phase(image_wrapped)
|
||||
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)
|
||||
yield check_unwrap, image
|
||||
mask = np.zeros(image.shape, dtype=np.bool)
|
||||
mask[4:6, 4:8] = True
|
||||
yield check_unwrap, image, mask
|
||||
|
||||
|
||||
def test_unwrap_3d():
|
||||
x, y, z = np.ogrid[:8, :12, :16]
|
||||
image = 2 * np.pi * (x * 0.2 + y * 0.1 + z * 0.05)
|
||||
yield check_unwrap, image
|
||||
mask = np.zeros(image.shape, dtype=np.bool)
|
||||
mask[4:6, 4:6, 1:3] = True
|
||||
yield check_unwrap, image, mask
|
||||
|
||||
|
||||
def check_wrap_around(ndim, axis):
|
||||
# create a ramp, but with the last pixel along axis equalling the first
|
||||
elements = 100
|
||||
ramp = np.linspace(0, 12 * np.pi, elements)
|
||||
ramp[-1] = ramp[0]
|
||||
image = ramp.reshape(tuple([elements if n == axis else 1
|
||||
for n in range(ndim)]))
|
||||
image_wrapped = np.angle(np.exp(1j * image))
|
||||
|
||||
index_first = tuple([0] * ndim)
|
||||
index_last = tuple([-1 if n == axis else 0 for n in range(ndim)])
|
||||
# unwrap the image without wrap around
|
||||
with warnings.catch_warnings():
|
||||
# We do not want warnings about length 1 dimensions
|
||||
warnings.simplefilter("ignore")
|
||||
image_unwrap_no_wrap_around = unwrap_phase(image_wrapped)
|
||||
print('endpoints without wrap_around:',
|
||||
image_unwrap_no_wrap_around[index_first],
|
||||
image_unwrap_no_wrap_around[index_last])
|
||||
# without wrap around, the endpoints of the image should differ
|
||||
assert abs(image_unwrap_no_wrap_around[index_first]
|
||||
- image_unwrap_no_wrap_around[index_last]) > np.pi
|
||||
# unwrap the image with wrap around
|
||||
wrap_around = [n == axis for n in range(ndim)]
|
||||
with warnings.catch_warnings():
|
||||
# We do not want warnings about length 1 dimensions
|
||||
warnings.simplefilter("ignore")
|
||||
image_unwrap_wrap_around = unwrap_phase(image_wrapped, wrap_around)
|
||||
print('endpoints with wrap_around:',
|
||||
image_unwrap_wrap_around[index_first],
|
||||
image_unwrap_wrap_around[index_last])
|
||||
# with wrap around, the endpoints of the image should be equal
|
||||
assert_almost_equal(image_unwrap_wrap_around[index_first],
|
||||
image_unwrap_wrap_around[index_last])
|
||||
|
||||
|
||||
def test_wrap_around():
|
||||
for ndim in (2, 3):
|
||||
for axis in range(ndim):
|
||||
yield check_wrap_around, ndim, axis
|
||||
|
||||
|
||||
def test_mask():
|
||||
length = 100
|
||||
ramps = [np.linspace(0, 4 * np.pi, length),
|
||||
np.linspace(0, 8 * np.pi, length),
|
||||
np.linspace(0, 6 * np.pi, length)]
|
||||
image = np.vstack(ramps)
|
||||
mask_1d = np.ones((length,), dtype=np.bool)
|
||||
mask_1d[0] = mask_1d[-1] = False
|
||||
for i in range(len(ramps)):
|
||||
# mask all ramps but the i'th one
|
||||
mask = np.zeros(image.shape, dtype=np.bool)
|
||||
mask |= mask_1d.reshape(1, -1)
|
||||
mask[i, :] = False # unmask i'th ramp
|
||||
image_wrapped = np.ma.array(np.angle(np.exp(1j * image)), mask=mask)
|
||||
image_unwrapped = unwrap_phase(image_wrapped)
|
||||
image_unwrapped -= image_unwrapped[0, 0] # remove phase shift
|
||||
# The end of the unwrapped array should have value equal to the
|
||||
# endpoint of the unmasked ramp
|
||||
assert_array_almost_equal(image_unwrapped[:, -1], image[i, -1])
|
||||
|
||||
# Same tests, but forcing use of the 3D unwrapper by reshaping
|
||||
image_wrapped_3d = image_wrapped.reshape((1,) + image_wrapped.shape)
|
||||
image_unwrapped_3d = unwrap_phase(image_wrapped_3d)
|
||||
image_unwrapped_3d -= image_unwrapped_3d[0, 0, 0] # remove phase shift
|
||||
assert_array_almost_equal(image_unwrapped_3d[:, :, -1], image[i, -1])
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
run_module_suite()
|
||||
@@ -0,0 +1,104 @@
|
||||
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
|
||||
|
||||
|
||||
def unwrap_phase(image, wrap_around=False):
|
||||
'''From ``image``, wrapped to lie in the interval [-pi, pi), recover the
|
||||
original, unwrapped image.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
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. 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
|
||||
-------
|
||||
image_unwrapped : array_like, float32
|
||||
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]
|
||||
>>> image = 12 * np.pi * np.exp(-(c0**2 + c1**2))
|
||||
>>> image_wrapped = np.angle(np.exp(1j * image))
|
||||
>>> image_unwrapped = unwrap_phase(image_wrapped)
|
||||
>>> np.std(image_unwrapped - image) < 1e-6 # A constant offset is normal
|
||||
True
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Miguel Arevallilo Herraez, David R. Burton, Michael J. Lalor,
|
||||
and Munther A. Gdeisat, "Fast two-dimensional phase-unwrapping
|
||||
algorithm based on sorting by reliability following a noncontinuous
|
||||
path", Journal Applied Optics, Vol. 41, No. 35 (2002) 7437,
|
||||
.. [2] Abdul-Rahman, H., Gdeisat, M., Burton, D., & Lalor, M., "Fast
|
||||
three-dimensional phase-unwrapping algorithm based on sorting by
|
||||
reliability following a non-continuous path. In W. Osten,
|
||||
C. Gorecki, & E. L. Novak (Eds.), Optical Metrology (2005) 32--40,
|
||||
International Society for Optics and Photonics.
|
||||
'''
|
||||
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__')
|
||||
and not isinstance(wrap_around, string_types)):
|
||||
if len(wrap_around) != image.ndim:
|
||||
raise ValueError('Length of wrap_around must equal the '
|
||||
'dimensionality of image')
|
||||
wrap_around = [bool(wa) for wa in wrap_around]
|
||||
else:
|
||||
raise ValueError('wrap_around must be a bool or a sequence with '
|
||||
'length equal to the dimensionality of image')
|
||||
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'])
|
||||
else:
|
||||
mask = np.zeros_like(image, dtype=np.uint8, order='C')
|
||||
image_not_masked = np.asarray(image, dtype=np.float64, order='C')
|
||||
image_unwrapped = np.empty_like(image, dtype=np.float64, order='C')
|
||||
|
||||
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:
|
||||
unwrap_3d(image_not_masked, mask, image_unwrapped,
|
||||
wrap_around)
|
||||
|
||||
if np.ma.isMaskedArray(image):
|
||||
return np.ma.array(image_unwrapped, mask=mask)
|
||||
else:
|
||||
return image_unwrapped
|
||||
@@ -0,0 +1,725 @@
|
||||
// 2D phase unwrapping, modified for inclusion in scipy by Gregor Thalhammer
|
||||
// Original file name: Miguel_2D_unwrapper_with_mask_and_wrap_around_option.c
|
||||
|
||||
//This program was written by Munther Gdeisat and Miguel Arevallilo Herraez to program the two-dimensional unwrapper
|
||||
//entitled "Fast two-dimensional phase-unwrapping algorithm based on sorting by
|
||||
//reliability following a noncontinuous path"
|
||||
//by Miguel Arevallilo Herraez, David R. Burton, Michael J. Lalor, and Munther A. Gdeisat
|
||||
//published in the Journal Applied Optics, Vol. 41, No. 35, pp. 7437, 2002.
|
||||
//This program was written by Munther Gdeisat, Liverpool John Moores University, United Kingdom.
|
||||
//Date 26th August 2007
|
||||
//The wrapped phase map is assumed to be of floating point data type. The resultant unwrapped phase map is also of floating point type.
|
||||
//The mask is of byte data type.
|
||||
//When the mask is 255 this means that the pixel is valid
|
||||
//When the mask is 0 this means that the pixel is invalid (noisy or corrupted pixel)
|
||||
//This program takes into consideration the image wrap around problem encountered in MRI imaging.
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
|
||||
#define PI M_PI
|
||||
#define TWOPI (2 * M_PI)
|
||||
|
||||
//TODO: remove global variables
|
||||
//TODO: make thresholds independent
|
||||
|
||||
#define NOMASK 0
|
||||
#define MASK 1
|
||||
|
||||
typedef struct
|
||||
{
|
||||
double mod;
|
||||
int x_connectivity;
|
||||
int y_connectivity;
|
||||
int no_of_edges;
|
||||
} params_t;
|
||||
|
||||
//PIXELM information
|
||||
struct PIXELM
|
||||
{
|
||||
int increment; //No. of 2*pi to add to the pixel to unwrap it
|
||||
int number_of_pixels_in_group;//No. of pixel in the pixel group
|
||||
double value; //value of the pixel
|
||||
double reliability;
|
||||
unsigned char input_mask; //0 pixel is masked. NOMASK pixel is not masked
|
||||
unsigned char extended_mask; //0 pixel is masked. NOMASK pixel is not masked
|
||||
int group; //group No.
|
||||
int new_group;
|
||||
struct PIXELM *head; //pointer to the first pixel in the group in the linked list
|
||||
struct PIXELM *last; //pointer to the last pixel in the group
|
||||
struct PIXELM *next; //pointer to the next pixel in the group
|
||||
};
|
||||
|
||||
typedef struct PIXELM PIXELM;
|
||||
|
||||
//the EDGE is the line that connects two pixels.
|
||||
//if we have S pixels, then we have S horizontal edges and S vertical edges
|
||||
struct EDGE
|
||||
{
|
||||
double reliab; //reliabilty of the edge and it depends on the two pixels
|
||||
PIXELM *pointer_1; //pointer to the first pixel
|
||||
PIXELM *pointer_2; //pointer to the second pixel
|
||||
int increment; //No. of 2*pi to add to one of the pixels to
|
||||
//unwrap it with respect to the second
|
||||
};
|
||||
|
||||
typedef struct EDGE EDGE;
|
||||
|
||||
//---------------start quicker_sort algorithm --------------------------------
|
||||
#define swap(x,y) {EDGE t; t=x; x=y; y=t;}
|
||||
#define order(x,y) if (x.reliab > y.reliab) swap(x,y)
|
||||
#define o2(x,y) order(x,y)
|
||||
#define o3(x,y,z) o2(x,y); o2(x,z); o2(y,z)
|
||||
|
||||
typedef enum {yes, no} yes_no;
|
||||
|
||||
yes_no find_pivot(EDGE *left, EDGE *right, double *pivot_ptr)
|
||||
{
|
||||
EDGE a, b, c, *p;
|
||||
|
||||
a = *left;
|
||||
b = *(left + (right - left) /2 );
|
||||
c = *right;
|
||||
o3(a,b,c);
|
||||
|
||||
if (a.reliab < b.reliab)
|
||||
{
|
||||
*pivot_ptr = b.reliab;
|
||||
return yes;
|
||||
}
|
||||
|
||||
if (b.reliab < c.reliab)
|
||||
{
|
||||
*pivot_ptr = c.reliab;
|
||||
return yes;
|
||||
}
|
||||
|
||||
for (p = left + 1; p <= right; ++p)
|
||||
{
|
||||
if (p->reliab != left->reliab)
|
||||
{
|
||||
*pivot_ptr = (p->reliab < left->reliab) ? left->reliab : p->reliab;
|
||||
return yes;
|
||||
}
|
||||
return no;
|
||||
}
|
||||
}
|
||||
|
||||
EDGE *partition(EDGE *left, EDGE *right, double pivot)
|
||||
{
|
||||
while (left <= right)
|
||||
{
|
||||
while (left->reliab < pivot)
|
||||
++left;
|
||||
while (right->reliab >= pivot)
|
||||
--right;
|
||||
if (left < right)
|
||||
{
|
||||
swap (*left, *right);
|
||||
++left;
|
||||
--right;
|
||||
}
|
||||
}
|
||||
return left;
|
||||
}
|
||||
|
||||
void quicker_sort(EDGE *left, EDGE *right)
|
||||
{
|
||||
EDGE *p;
|
||||
double pivot;
|
||||
|
||||
if (find_pivot(left, right, &pivot) == yes)
|
||||
{
|
||||
p = partition(left, right, pivot);
|
||||
quicker_sort(left, p - 1);
|
||||
quicker_sort(p, right);
|
||||
}
|
||||
}
|
||||
//--------------end quicker_sort algorithm -----------------------------------
|
||||
|
||||
//--------------------start initialize pixels ----------------------------------
|
||||
//initialize pixels. See the explination of the pixel class above.
|
||||
//initially every pixel is assumed to belong to a group consisting of only itself
|
||||
void initialisePIXELs(double *wrapped_image, unsigned char *input_mask, unsigned char *extended_mask, PIXELM *pixel, int image_width, int image_height)
|
||||
{
|
||||
PIXELM *pixel_pointer = pixel;
|
||||
double *wrapped_image_pointer = wrapped_image;
|
||||
unsigned char *input_mask_pointer = input_mask;
|
||||
unsigned char *extended_mask_pointer = extended_mask;
|
||||
int i, j;
|
||||
|
||||
for (i=0; i < image_height; i++)
|
||||
{
|
||||
for (j=0; j < image_width; j++)
|
||||
{
|
||||
pixel_pointer->increment = 0;
|
||||
pixel_pointer->number_of_pixels_in_group = 1;
|
||||
pixel_pointer->value = *wrapped_image_pointer;
|
||||
pixel_pointer->reliability = 9999999. + rand();
|
||||
pixel_pointer->input_mask = *input_mask_pointer;
|
||||
pixel_pointer->extended_mask = *extended_mask_pointer;
|
||||
pixel_pointer->head = pixel_pointer;
|
||||
pixel_pointer->last = pixel_pointer;
|
||||
pixel_pointer->next = NULL;
|
||||
pixel_pointer->new_group = 0;
|
||||
pixel_pointer->group = -1;
|
||||
pixel_pointer++;
|
||||
wrapped_image_pointer++;
|
||||
input_mask_pointer++;
|
||||
extended_mask_pointer++;
|
||||
}
|
||||
}
|
||||
}
|
||||
//-------------------end initialize pixels -----------
|
||||
|
||||
//gamma function in the paper
|
||||
double wrap(double pixel_value)
|
||||
{
|
||||
double wrapped_pixel_value;
|
||||
if (pixel_value > PI) wrapped_pixel_value = pixel_value - TWOPI;
|
||||
else if (pixel_value < -PI) wrapped_pixel_value = pixel_value + TWOPI;
|
||||
else wrapped_pixel_value = pixel_value;
|
||||
return wrapped_pixel_value;
|
||||
}
|
||||
|
||||
// pixelL_value is the left pixel, pixelR_value is the right pixel
|
||||
int find_wrap(double pixelL_value, double pixelR_value)
|
||||
{
|
||||
double difference;
|
||||
int wrap_value;
|
||||
difference = pixelL_value - pixelR_value;
|
||||
|
||||
if (difference > PI) wrap_value = -1;
|
||||
else if (difference < -PI) wrap_value = 1;
|
||||
else wrap_value = 0;
|
||||
|
||||
return wrap_value;
|
||||
}
|
||||
|
||||
void extend_mask(unsigned char *input_mask, unsigned char *extended_mask,
|
||||
int image_width, int image_height,
|
||||
params_t *params)
|
||||
{
|
||||
int i,j;
|
||||
int image_width_plus_one = image_width + 1;
|
||||
int image_width_minus_one = image_width - 1;
|
||||
unsigned char *IMP = input_mask + image_width + 1; //input mask pointer
|
||||
unsigned char *EMP = extended_mask + image_width + 1; //extended mask pointer
|
||||
|
||||
//extend the mask for the image except borders
|
||||
for (i=1; i < image_height - 1; ++i)
|
||||
{
|
||||
for (j=1; j < image_width - 1; ++j)
|
||||
{
|
||||
if ( (*IMP) == NOMASK && (*(IMP + 1) == NOMASK) && (*(IMP - 1) == NOMASK) &&
|
||||
(*(IMP + image_width) == NOMASK) && (*(IMP - image_width) == NOMASK) &&
|
||||
(*(IMP - image_width_minus_one) == NOMASK) && (*(IMP - image_width_plus_one) == NOMASK) &&
|
||||
(*(IMP + image_width_minus_one) == NOMASK) && (*(IMP + image_width_plus_one) == NOMASK) )
|
||||
{
|
||||
*EMP = NOMASK;
|
||||
}
|
||||
++EMP;
|
||||
++IMP;
|
||||
}
|
||||
EMP += 2;
|
||||
IMP += 2;
|
||||
}
|
||||
|
||||
if (params->x_connectivity == 1)
|
||||
{
|
||||
//extend the mask for the right border of the image
|
||||
IMP = input_mask + 2 * image_width - 1;
|
||||
EMP = extended_mask + 2 * image_width -1;
|
||||
for (i=1; i < image_height - 1; ++ i)
|
||||
{
|
||||
if ( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) &&
|
||||
(*(IMP + image_width) == NOMASK) && (*(IMP - image_width) == NOMASK) &&
|
||||
(*(IMP - image_width - 1) == NOMASK) && (*(IMP - image_width + 1) == NOMASK) &&
|
||||
(*(IMP + image_width - 1) == NOMASK) && (*(IMP - 2 * image_width + 1) == NOMASK) )
|
||||
{
|
||||
*EMP = NOMASK;
|
||||
}
|
||||
EMP += image_width;
|
||||
IMP += image_width;
|
||||
}
|
||||
|
||||
//extend the mask for the left border of the image
|
||||
IMP = input_mask + image_width;
|
||||
EMP = extended_mask + image_width;
|
||||
for (i=1; i < image_height - 1; ++i)
|
||||
{
|
||||
if ( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) &&
|
||||
(*(IMP + image_width) == NOMASK) && (*(IMP - image_width) == NOMASK) &&
|
||||
(*(IMP - image_width + 1) == NOMASK) && (*(IMP + image_width + 1) == NOMASK) &&
|
||||
(*(IMP + image_width - 1) == NOMASK) && (*(IMP + 2 * image_width - 1) == NOMASK) )
|
||||
{
|
||||
*EMP = NOMASK;
|
||||
}
|
||||
EMP += image_width;
|
||||
IMP += image_width;
|
||||
}
|
||||
}
|
||||
|
||||
if (params->y_connectivity == 1)
|
||||
{
|
||||
//extend the mask for the top border of the image
|
||||
IMP = input_mask + 1;
|
||||
EMP = extended_mask + 1;
|
||||
for (i=1; i < image_width - 1; ++i)
|
||||
{
|
||||
if ( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) &&
|
||||
(*(IMP + image_width) == NOMASK) && (*(IMP + image_width * (image_height - 1)) == NOMASK) &&
|
||||
(*(IMP + image_width + 1) == NOMASK) && (*(IMP + image_width - 1) == NOMASK) &&
|
||||
(*(IMP + image_width * (image_height - 1) - 1) == NOMASK) && (*(IMP + image_width * (image_height - 1) + 1) == NOMASK) )
|
||||
{
|
||||
*EMP = NOMASK;
|
||||
}
|
||||
EMP++;
|
||||
IMP++;
|
||||
}
|
||||
|
||||
//extend the mask for the bottom border of the image
|
||||
IMP = input_mask + image_width * (image_height - 1) + 1;
|
||||
EMP = extended_mask + image_width * (image_height - 1) + 1;
|
||||
for (i=1; i < image_width - 1; ++i)
|
||||
{
|
||||
if ( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) &&
|
||||
(*(IMP - image_width) == NOMASK) && (*(IMP - image_width - 1) == NOMASK) && (*(IMP - image_width + 1) == NOMASK) &&
|
||||
(*(IMP - image_width * (image_height - 1) ) == NOMASK) &&
|
||||
(*(IMP - image_width * (image_height - 1) - 1) == NOMASK) &&
|
||||
(*(IMP - image_width * (image_height - 1) + 1) == NOMASK) )
|
||||
{
|
||||
*EMP = NOMASK;
|
||||
}
|
||||
EMP++;
|
||||
IMP++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void calculate_reliability(double *wrappedImage, PIXELM *pixel,
|
||||
int image_width, int image_height,
|
||||
params_t *params)
|
||||
{
|
||||
int image_width_plus_one = image_width + 1;
|
||||
int image_width_minus_one = image_width - 1;
|
||||
PIXELM *pixel_pointer = pixel + image_width_plus_one;
|
||||
double *WIP = wrappedImage + image_width_plus_one; //WIP is the wrapped image pointer
|
||||
double H, V, D1, D2;
|
||||
int i, j;
|
||||
|
||||
for (i = 1; i < image_height -1; ++i)
|
||||
{
|
||||
for (j = 1; j < image_width - 1; ++j)
|
||||
{
|
||||
if (pixel_pointer->extended_mask == NOMASK)
|
||||
{
|
||||
H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1));
|
||||
V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width));
|
||||
D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one));
|
||||
D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one));
|
||||
pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2;
|
||||
}
|
||||
pixel_pointer++;
|
||||
WIP++;
|
||||
}
|
||||
pixel_pointer += 2;
|
||||
WIP += 2;
|
||||
}
|
||||
|
||||
if (params->x_connectivity == 1)
|
||||
{
|
||||
//calculating the reliability for the left border of the image
|
||||
PIXELM *pixel_pointer = pixel + image_width;
|
||||
double *WIP = wrappedImage + image_width;
|
||||
|
||||
for (i = 1; i < image_height - 1; ++i)
|
||||
{
|
||||
if (pixel_pointer->extended_mask == NOMASK)
|
||||
{
|
||||
H = wrap(*(WIP + image_width - 1) - *WIP) - wrap(*WIP - *(WIP + 1));
|
||||
V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width));
|
||||
D1 = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one));
|
||||
D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP + 2* image_width - 1));
|
||||
pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2;
|
||||
}
|
||||
pixel_pointer += image_width;
|
||||
WIP += image_width;
|
||||
}
|
||||
|
||||
//calculating the reliability for the right border of the image
|
||||
pixel_pointer = pixel + 2 * image_width - 1;
|
||||
WIP = wrappedImage + 2 * image_width - 1;
|
||||
|
||||
for (i = 1; i < image_height - 1; ++i)
|
||||
{
|
||||
if (pixel_pointer->extended_mask == NOMASK)
|
||||
{
|
||||
H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP - image_width_minus_one));
|
||||
V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP + image_width));
|
||||
D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP + 1));
|
||||
D2 = wrap(*(WIP - 2 * image_width - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one));
|
||||
pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2;
|
||||
}
|
||||
pixel_pointer += image_width;
|
||||
WIP += image_width;
|
||||
}
|
||||
}
|
||||
|
||||
if (params->y_connectivity == 1)
|
||||
{
|
||||
//calculating the reliability for the top border of the image
|
||||
PIXELM *pixel_pointer = pixel + 1;
|
||||
double *WIP = wrappedImage + 1;
|
||||
|
||||
for (i = 1; i < image_width - 1; ++i)
|
||||
{
|
||||
if (pixel_pointer->extended_mask == NOMASK)
|
||||
{
|
||||
H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1));
|
||||
V = wrap(*(WIP + image_width*(image_height - 1)) - *WIP) - wrap(*WIP - *(WIP + image_width));
|
||||
D1 = wrap(*(WIP + image_width*(image_height - 1) - 1) - *WIP) - wrap(*WIP - *(WIP + image_width_plus_one));
|
||||
D2 = wrap(*(WIP + image_width*(image_height - 1) + 1) - *WIP) - wrap(*WIP - *(WIP + image_width_minus_one));
|
||||
pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2;
|
||||
}
|
||||
pixel_pointer++;
|
||||
WIP++;
|
||||
}
|
||||
|
||||
//calculating the reliability for the bottom border of the image
|
||||
pixel_pointer = pixel + (image_height - 1) * image_width + 1;
|
||||
WIP = wrappedImage + (image_height - 1) * image_width + 1;
|
||||
|
||||
for (i = 1; i < image_width - 1; ++i)
|
||||
{
|
||||
if (pixel_pointer->extended_mask == NOMASK)
|
||||
{
|
||||
H = wrap(*(WIP - 1) - *WIP) - wrap(*WIP - *(WIP + 1));
|
||||
V = wrap(*(WIP - image_width) - *WIP) - wrap(*WIP - *(WIP -(image_height - 1) * (image_width)));
|
||||
D1 = wrap(*(WIP - image_width_plus_one) - *WIP) - wrap(*WIP - *(WIP - (image_height - 1) * (image_width) + 1));
|
||||
D2 = wrap(*(WIP - image_width_minus_one) - *WIP) - wrap(*WIP - *(WIP - (image_height - 1) * (image_width) - 1));
|
||||
pixel_pointer->reliability = H*H + V*V + D1*D1 + D2*D2;
|
||||
}
|
||||
pixel_pointer++;
|
||||
WIP++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//calculate the reliability of the horizontal edges of the image
|
||||
//it is calculated by adding the reliability of pixel and the relibility of
|
||||
//its right-hand neighbour
|
||||
//edge is calculated between a pixel and its next neighbour
|
||||
void horizontalEDGEs(PIXELM *pixel, EDGE *edge,
|
||||
int image_width, int image_height,
|
||||
params_t *params)
|
||||
{
|
||||
int i, j;
|
||||
EDGE *edge_pointer = edge;
|
||||
PIXELM *pixel_pointer = pixel;
|
||||
int no_of_edges = params->no_of_edges;
|
||||
|
||||
for (i = 0; i < image_height; i++)
|
||||
{
|
||||
for (j = 0; j < image_width - 1; j++)
|
||||
{
|
||||
if (pixel_pointer->input_mask == NOMASK && (pixel_pointer + 1)->input_mask == NOMASK)
|
||||
{
|
||||
edge_pointer->pointer_1 = pixel_pointer;
|
||||
edge_pointer->pointer_2 = (pixel_pointer+1);
|
||||
edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer + 1)->reliability;
|
||||
edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer + 1)->value);
|
||||
edge_pointer++;
|
||||
no_of_edges++;
|
||||
}
|
||||
pixel_pointer++;
|
||||
}
|
||||
pixel_pointer++;
|
||||
}
|
||||
//construct edges at the right border of the image
|
||||
if (params->x_connectivity == 1)
|
||||
{
|
||||
pixel_pointer = pixel + image_width - 1;
|
||||
for (i = 0; i < image_height; i++)
|
||||
{
|
||||
if (pixel_pointer->input_mask == NOMASK && (pixel_pointer - image_width + 1)->input_mask == NOMASK)
|
||||
{
|
||||
edge_pointer->pointer_1 = pixel_pointer;
|
||||
edge_pointer->pointer_2 = (pixel_pointer - image_width + 1);
|
||||
edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer - image_width + 1)->reliability;
|
||||
edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer - image_width + 1)->value);
|
||||
edge_pointer++;
|
||||
no_of_edges++;
|
||||
}
|
||||
pixel_pointer+=image_width;
|
||||
}
|
||||
}
|
||||
params->no_of_edges = no_of_edges;
|
||||
}
|
||||
|
||||
//calculate the reliability of the vertical edges of the image
|
||||
//it is calculated by adding the reliability of pixel and the relibility of
|
||||
//its lower neighbour in the image.
|
||||
void verticalEDGEs(PIXELM *pixel, EDGE *edge,
|
||||
int image_width, int image_height,
|
||||
params_t *params)
|
||||
{
|
||||
int i, j;
|
||||
int no_of_edges = params->no_of_edges;
|
||||
PIXELM *pixel_pointer = pixel;
|
||||
EDGE *edge_pointer = edge + no_of_edges;
|
||||
|
||||
for (i=0; i < image_height - 1; i++)
|
||||
{
|
||||
for (j=0; j < image_width; j++)
|
||||
{
|
||||
if (pixel_pointer->input_mask == NOMASK && (pixel_pointer + image_width)->input_mask == NOMASK)
|
||||
{
|
||||
edge_pointer->pointer_1 = pixel_pointer;
|
||||
edge_pointer->pointer_2 = (pixel_pointer + image_width);
|
||||
edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer + image_width)->reliability;
|
||||
edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer + image_width)->value);
|
||||
edge_pointer++;
|
||||
no_of_edges++;
|
||||
}
|
||||
pixel_pointer++;
|
||||
} //j loop
|
||||
} // i loop
|
||||
|
||||
//construct edges that connect at the bottom border of the image
|
||||
if (params->y_connectivity == 1)
|
||||
{
|
||||
pixel_pointer = pixel + image_width *(image_height - 1);
|
||||
for (i = 0; i < image_width; i++)
|
||||
{
|
||||
if (pixel_pointer->input_mask == NOMASK && (pixel_pointer - image_width *(image_height - 1))->input_mask == NOMASK)
|
||||
{
|
||||
edge_pointer->pointer_1 = pixel_pointer;
|
||||
edge_pointer->pointer_2 = (pixel_pointer - image_width *(image_height - 1));
|
||||
edge_pointer->reliab = pixel_pointer->reliability + (pixel_pointer - image_width *(image_height - 1))->reliability;
|
||||
edge_pointer->increment = find_wrap(pixel_pointer->value, (pixel_pointer - image_width *(image_height - 1))->value);
|
||||
edge_pointer++;
|
||||
no_of_edges++;
|
||||
}
|
||||
pixel_pointer++;
|
||||
}
|
||||
}
|
||||
params->no_of_edges = no_of_edges;
|
||||
}
|
||||
|
||||
//gather the pixels of the image into groups
|
||||
void gatherPIXELs(EDGE *edge, params_t *params)
|
||||
{
|
||||
int k;
|
||||
PIXELM *PIXEL1;
|
||||
PIXELM *PIXEL2;
|
||||
PIXELM *group1;
|
||||
PIXELM *group2;
|
||||
EDGE *pointer_edge = edge;
|
||||
int incremento;
|
||||
|
||||
for (k = 0; k < params->no_of_edges; k++)
|
||||
{
|
||||
PIXEL1 = pointer_edge->pointer_1;
|
||||
PIXEL2 = pointer_edge->pointer_2;
|
||||
|
||||
//PIXELM 1 and PIXELM 2 belong to different groups
|
||||
//initially each pixel is a group by it self and one pixel can construct a group
|
||||
//no else or else if to this if
|
||||
if (PIXEL2->head != PIXEL1->head)
|
||||
{
|
||||
//PIXELM 2 is alone in its group
|
||||
//merge this pixel with PIXELM 1 group and find the number of 2 pi to add
|
||||
//to or subtract to unwrap it
|
||||
if ((PIXEL2->next == NULL) && (PIXEL2->head == PIXEL2))
|
||||
{
|
||||
PIXEL1->head->last->next = PIXEL2;
|
||||
PIXEL1->head->last = PIXEL2;
|
||||
(PIXEL1->head->number_of_pixels_in_group)++;
|
||||
PIXEL2->head=PIXEL1->head;
|
||||
PIXEL2->increment = PIXEL1->increment-pointer_edge->increment;
|
||||
}
|
||||
|
||||
//PIXELM 1 is alone in its group
|
||||
//merge this pixel with PIXELM 2 group and find the number of 2 pi to add
|
||||
//to or subtract to unwrap it
|
||||
else if ((PIXEL1->next == NULL) && (PIXEL1->head == PIXEL1))
|
||||
{
|
||||
PIXEL2->head->last->next = PIXEL1;
|
||||
PIXEL2->head->last = PIXEL1;
|
||||
(PIXEL2->head->number_of_pixels_in_group)++;
|
||||
PIXEL1->head = PIXEL2->head;
|
||||
PIXEL1->increment = PIXEL2->increment+pointer_edge->increment;
|
||||
}
|
||||
|
||||
//PIXELM 1 and PIXELM 2 both have groups
|
||||
else
|
||||
{
|
||||
group1 = PIXEL1->head;
|
||||
group2 = PIXEL2->head;
|
||||
//if the no. of pixels in PIXELM 1 group is larger than the
|
||||
//no. of pixels in PIXELM 2 group. Merge PIXELM 2 group to
|
||||
//PIXELM 1 group and find the number of wraps between PIXELM 2
|
||||
//group and PIXELM 1 group to unwrap PIXELM 2 group with respect
|
||||
//to PIXELM 1 group. the no. of wraps will be added to PIXELM 2
|
||||
//group in the future
|
||||
if (group1->number_of_pixels_in_group > group2->number_of_pixels_in_group)
|
||||
{
|
||||
//merge PIXELM 2 with PIXELM 1 group
|
||||
group1->last->next = group2;
|
||||
group1->last = group2->last;
|
||||
group1->number_of_pixels_in_group = group1->number_of_pixels_in_group + group2->number_of_pixels_in_group;
|
||||
incremento = PIXEL1->increment-pointer_edge->increment - PIXEL2->increment;
|
||||
//merge the other pixels in PIXELM 2 group to PIXELM 1 group
|
||||
while (group2 != NULL)
|
||||
{
|
||||
group2->head = group1;
|
||||
group2->increment += incremento;
|
||||
group2 = group2->next;
|
||||
}
|
||||
}
|
||||
|
||||
//if the no. of pixels in PIXELM 2 group is larger than the
|
||||
//no. of pixels in PIXELM 1 group. Merge PIXELM 1 group to
|
||||
//PIXELM 2 group and find the number of wraps between PIXELM 2
|
||||
//group and PIXELM 1 group to unwrap PIXELM 1 group with respect
|
||||
//to PIXELM 2 group. the no. of wraps will be added to PIXELM 1
|
||||
//group in the future
|
||||
else
|
||||
{
|
||||
//merge PIXELM 1 with PIXELM 2 group
|
||||
group2->last->next = group1;
|
||||
group2->last = group1->last;
|
||||
group2->number_of_pixels_in_group = group2->number_of_pixels_in_group + group1->number_of_pixels_in_group;
|
||||
incremento = PIXEL2->increment + pointer_edge->increment - PIXEL1->increment;
|
||||
//merge the other pixels in PIXELM 2 group to PIXELM 1 group
|
||||
while (group1 != NULL)
|
||||
{
|
||||
group1->head = group2;
|
||||
group1->increment += incremento;
|
||||
group1 = group1->next;
|
||||
} // while
|
||||
|
||||
} // else
|
||||
} //else
|
||||
} //if
|
||||
pointer_edge++;
|
||||
}
|
||||
}
|
||||
|
||||
//unwrap the image
|
||||
void unwrapImage(PIXELM *pixel, int image_width, int image_height)
|
||||
{
|
||||
int i;
|
||||
int image_size = image_width * image_height;
|
||||
PIXELM *pixel_pointer=pixel;
|
||||
|
||||
for (i = 0; i < image_size; i++)
|
||||
{
|
||||
pixel_pointer->value += TWOPI * (double)(pixel_pointer->increment);
|
||||
pixel_pointer++;
|
||||
}
|
||||
}
|
||||
|
||||
//set the masked pixels (mask = 0) to the minimum of the unwrapper phase
|
||||
void maskImage(PIXELM *pixel, unsigned char *input_mask, int image_width, int image_height)
|
||||
{
|
||||
int image_width_plus_one = image_width + 1;
|
||||
int image_height_plus_one = image_height + 1;
|
||||
int image_width_minus_one = image_width - 1;
|
||||
int image_height_minus_one = image_height - 1;
|
||||
|
||||
PIXELM *pointer_pixel = pixel;
|
||||
unsigned char *IMP = input_mask; //input mask pointer
|
||||
double min=99999999;
|
||||
int i;
|
||||
int image_size = image_width * image_height;
|
||||
|
||||
//find the minimum of the unwrapped phase
|
||||
for (i = 0; i < image_size; i++)
|
||||
{
|
||||
if ((pointer_pixel->value < min) && (*IMP == NOMASK))
|
||||
min = pointer_pixel->value;
|
||||
|
||||
pointer_pixel++;
|
||||
IMP++;
|
||||
}
|
||||
|
||||
pointer_pixel = pixel;
|
||||
IMP = input_mask;
|
||||
|
||||
//set the masked pixels to minimum
|
||||
for (i = 0; i < image_size; i++)
|
||||
{
|
||||
if ((*IMP) == MASK)
|
||||
{
|
||||
pointer_pixel->value = min;
|
||||
}
|
||||
pointer_pixel++;
|
||||
IMP++;
|
||||
}
|
||||
}
|
||||
|
||||
//the input to this unwrapper is an array that contains the wrapped
|
||||
//phase map. copy the image on the buffer passed to this unwrapper to
|
||||
//over-write the unwrapped phase map on the buffer of the wrapped
|
||||
//phase map.
|
||||
void returnImage(PIXELM *pixel, double *unwrapped_image, int image_width, int image_height)
|
||||
{
|
||||
int i;
|
||||
int image_size = image_width * image_height;
|
||||
double *unwrapped_image_pointer = unwrapped_image;
|
||||
PIXELM *pixel_pointer = pixel;
|
||||
|
||||
for (i=0; i < image_size; i++)
|
||||
{
|
||||
*unwrapped_image_pointer = pixel_pointer->value;
|
||||
pixel_pointer++;
|
||||
unwrapped_image_pointer++;
|
||||
}
|
||||
}
|
||||
|
||||
//the main function of the unwrapper
|
||||
void
|
||||
unwrap2D(double* wrapped_image, double* UnwrappedImage, unsigned char* input_mask,
|
||||
int image_width, int image_height,
|
||||
int wrap_around_x, int wrap_around_y)
|
||||
{
|
||||
params_t params = {TWOPI, wrap_around_x, wrap_around_y, 0};
|
||||
unsigned char *extended_mask;
|
||||
PIXELM *pixel;
|
||||
EDGE *edge;
|
||||
int image_size = image_height * image_width;
|
||||
int No_of_Edges_initially = 2 * image_width * image_height;
|
||||
|
||||
extended_mask = (unsigned char *) calloc(image_size, sizeof(unsigned char));
|
||||
pixel = (PIXELM *) calloc(image_size, sizeof(PIXELM));
|
||||
edge = (EDGE *) calloc(No_of_Edges_initially, sizeof(EDGE));
|
||||
|
||||
extend_mask(input_mask, extended_mask, image_width, image_height, ¶ms);
|
||||
initialisePIXELs(wrapped_image, input_mask, extended_mask, pixel, image_width, image_height);
|
||||
calculate_reliability(wrapped_image, pixel, image_width, image_height, ¶ms);
|
||||
horizontalEDGEs(pixel, edge, image_width, image_height, ¶ms);
|
||||
verticalEDGEs(pixel, edge, image_width, image_height, ¶ms);
|
||||
|
||||
//sort the EDGEs depending on their reiability. The PIXELs with higher
|
||||
//relibility (small value) first
|
||||
quicker_sort(edge, edge + params.no_of_edges - 1);
|
||||
|
||||
//gather PIXELs into groups
|
||||
gatherPIXELs(edge, ¶ms);
|
||||
|
||||
unwrapImage(pixel, image_width, image_height);
|
||||
maskImage(pixel, input_mask, image_width, image_height);
|
||||
|
||||
//copy the image from PIXELM structure to the unwrapped phase array
|
||||
//passed to this function
|
||||
//TODO: replace by (cython?) function to directly write into numpy array ?
|
||||
returnImage(pixel, UnwrappedImage, image_width, image_height);
|
||||
|
||||
free(edge);
|
||||
free(pixel);
|
||||
free(extended_mask);
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
@@ -10,6 +10,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_subpackage('color')
|
||||
config.add_subpackage('data')
|
||||
config.add_subpackage('draw')
|
||||
config.add_subpackage('exposure')
|
||||
config.add_subpackage('feature')
|
||||
config.add_subpackage('filter')
|
||||
config.add_subpackage('graph')
|
||||
|
||||
Reference in New Issue
Block a user