diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 39a2f988..86463b69 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -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 diff --git a/bento.info b/bento.info index 84ddf4a8..30367057 100644 --- a/bento.info +++ b/bento.info @@ -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 diff --git a/doc/examples/plot_phase_unwrap.py b/doc/examples/plot_phase_unwrap.py new file mode 100644 index 00000000..39055995 --- /dev/null +++ b/doc/examples/plot_phase_unwrap.py @@ -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 +""" diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py index b873c339..c8f7b1be 100644 --- a/skimage/exposure/__init__.py +++ b/skimage/exposure/__init__.py @@ -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'] diff --git a/skimage/exposure/_unwrap_1d.pyx b/skimage/exposure/_unwrap_1d.pyx new file mode 100644 index 00000000..f23e5f40 --- /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(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 diff --git a/skimage/exposure/_unwrap_2d.pyx b/skimage/exposure/_unwrap_2d.pyx new file mode 100644 index 00000000..6e889729 --- /dev/null +++ b/skimage/exposure/_unwrap_2d.pyx @@ -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], + ) diff --git a/skimage/exposure/_unwrap_3d.pyx b/skimage/exposure/_unwrap_3d.pyx new file mode 100644 index 00000000..370d58be --- /dev/null +++ b/skimage/exposure/_unwrap_3d.pyx @@ -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], + ) diff --git a/skimage/exposure/setup.py b/skimage/exposure/setup.py new file mode 100644 index 00000000..b90b082c --- /dev/null +++ b/skimage/exposure/setup.py @@ -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()) + ) diff --git a/skimage/exposure/tests/test_unwrap.py b/skimage/exposure/tests/test_unwrap.py new file mode 100644 index 00000000..1b3ecca3 --- /dev/null +++ b/skimage/exposure/tests/test_unwrap.py @@ -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() diff --git a/skimage/exposure/unwrap.py b/skimage/exposure/unwrap.py new file mode 100644 index 00000000..86444aef --- /dev/null +++ b/skimage/exposure/unwrap.py @@ -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 diff --git a/skimage/exposure/unwrap_2d_ljmu.c b/skimage/exposure/unwrap_2d_ljmu.c new file mode 100644 index 00000000..6be7966a --- /dev/null +++ b/skimage/exposure/unwrap_2d_ljmu.c @@ -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 +#include +#include +#include + +#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); +} diff --git a/skimage/exposure/unwrap_3d_ljmu.c b/skimage/exposure/unwrap_3d_ljmu.c new file mode 100644 index 00000000..1874c051 --- /dev/null +++ b/skimage/exposure/unwrap_3d_ljmu.c @@ -0,0 +1,1054 @@ +// 3D phase unwrapping, modified for inclusion in scipy by Gregor Thalhammer +// Original file name: Hussein_3D_unwrapper_with_mask_and_wrap_around_option.c + +//This program was written by Hussein Abdul-Rahman and Munther Gdeisat to program the three-dimensional phase unwrapper +//entitled "Fast three-dimensional phase-unwrapping algorithm based on sorting by +//reliability following a noncontinuous path" +//by Hussein Abdul-Rahman, Munther A. Gdeisat, David R. Burton, and Michael J. Lalor, +//published in the Proceedings of SPIE - +//The International Society for Optical Engineering, Vol. 5856, No. 1, 2005, pp. 32-40 +//This program was written by Munther Gdeisat, Liverpool John Moores University, United Kingdom. +//Date 31st August 2007 +//The wrapped phase volume is assumed to be of floating point data type. The resultant unwrapped phase volume is also of floating point type. +//Read the data from the file frame by frame +//The mask is of byte data type. +//When the mask is 255 this means that the voxel is valid +//When the mask is 0 this means that the voxel is invalid (noisy or corrupted voxel) +//This program takes into consideration the image wrap around problem encountered in MRI imaging. + +#include +#include +#include +#include + +#define PI M_PI +#define TWOPI (2 * M_PI) + +#define NOMASK 0 +#define MASK 1 + +typedef struct +{ + double mod; + int x_connectivity; + int y_connectivity; + int z_connectivity; + int no_of_edges; +} params_t; + +//VOXELM information +struct VOXELM +{ + int increment; //No. of 2*pi to add to the voxel to unwrap it + int number_of_voxels_in_group;//No. of voxel in the voxel group + double value; //value of the voxel + double reliability; + unsigned char input_mask; //MASK voxel is masked. NOMASK voxel is not masked + unsigned char extended_mask; //MASK voxel is masked. NOMASK voxel is not masked + int group; //group No. + int new_group; + struct VOXELM *head; //pointer to the first voxel in the group in the linked list + struct VOXELM *last; //pointer to the last voxel in the group + struct VOXELM *next; //pointer to the next voxel in the group +}; + +typedef struct VOXELM VOXELM; + +//the EDGE is the line that connects two voxels. +//if we have S voxels, then we have S horizontal edges and S vertical edges +struct EDGE +{ + double reliab; //reliabilty of the edge and it depends on the two voxels + VOXELM *pointer_1; //pointer to the first voxel + VOXELM *pointer_2; //pointer to the second voxel + int increment; //No. of 2*pi to add to one of the + //voxels 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 voxels ---------------------------------- +//initiale voxels. See the explanation of the voxel class above. +//initially every voxel is assumed to belong to a group consisting of only itself +void initialiseVOXELs(double *WrappedVolume, unsigned char *input_mask, unsigned char *extended_mask, VOXELM *voxel, int volume_width, int volume_height, int volume_depth) +{ + VOXELM *voxel_pointer = voxel; + double *wrapped_volume_pointer = WrappedVolume; + unsigned char *input_mask_pointer = input_mask; + unsigned char *extended_mask_pointer = extended_mask; + int n, i, j; + + for (n=0; n < volume_depth; n++) + { + for (i=0; i < volume_height; i++) + { + for (j=0; j < volume_width; j++) + { + voxel_pointer->increment = 0; + voxel_pointer->number_of_voxels_in_group = 1; + voxel_pointer->value = *wrapped_volume_pointer; + voxel_pointer->reliability = 9999999 + rand(); + voxel_pointer->input_mask = *input_mask_pointer; + voxel_pointer->extended_mask = *extended_mask_pointer; + voxel_pointer->head = voxel_pointer; + voxel_pointer->last = voxel_pointer; + voxel_pointer->next = NULL; + voxel_pointer->new_group = 0; + voxel_pointer->group = -1; + voxel_pointer++; + wrapped_volume_pointer++; + input_mask_pointer++; + extended_mask_pointer++; + } + } + } +} +//-------------------end initialize voxels ----------- + +//gamma function in the paper +double wrap(double voxel_value) +{ + double wrapped_voxel_value; + if (voxel_value > PI) wrapped_voxel_value = voxel_value - TWOPI; + else if (voxel_value < -PI) wrapped_voxel_value = voxel_value + TWOPI; + else wrapped_voxel_value = voxel_value; + return wrapped_voxel_value; +} + +// voxelL_value is the left voxel, voxelR_value is the right voxel +int find_wrap(double voxelL_value, double voxelR_value) +{ + double difference; + int wrap_value; + difference = voxelL_value - voxelR_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 volume_width, int volume_height, int volume_depth, params_t *params) +{ + int n, i, j; + int vw = volume_width, vh = volume_height, vd = volume_depth; + int fs = volume_width * volume_height; //frame size + int frame_size = volume_width * volume_height; + int volume_size = volume_width * volume_height * volume_depth; //volume size + int vs = volume_size; + unsigned char *IMP = input_mask + frame_size + volume_width + 1; //input mask pointer + unsigned char *EMP = extended_mask + frame_size + volume_width + 1; //extended mask pointer + + //extend the mask for the volume except borders + for (n=1; n < volume_depth - 1; n++) + { + for (i=1; i < volume_height - 1; i++) + { + for (j=1; j < volume_width - 1; j++) + { + if( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) && + (*(IMP + vw) == NOMASK) && (*(IMP + vw - 1) == NOMASK) && (*(IMP + vw + 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP - vw - 1) == NOMASK) && (*(IMP - vw + 1) == NOMASK) && + (*(IMP + fs) == NOMASK) && (*(IMP + fs - 1) == NOMASK) && (*(IMP + fs + 1) == NOMASK) && + (*(IMP + fs - vw) == NOMASK) && (*(IMP + fs - vw - 1) == NOMASK) && (*(IMP + fs - vw + 1) == NOMASK) && + (*(IMP + fs + vw) == NOMASK) && (*(IMP + fs + vw - 1) == NOMASK) && (*(IMP + fs + vw + 1) == NOMASK) && + (*(IMP - fs) == NOMASK) && (*(IMP - fs - 1) == NOMASK) && (*(IMP - fs + 1) == NOMASK) && + (*(IMP - fs - vw) == NOMASK) && (*(IMP - fs - vw - 1) == NOMASK) && (*(IMP - fs - vw + 1) == NOMASK) && + (*(IMP - fs + vw) == NOMASK) && (*(IMP - fs + vw - 1) == NOMASK) && (*(IMP - fs + vw + 1) == NOMASK)) + { + *EMP = NOMASK; + } + ++EMP; + ++IMP; + } + EMP += 2; + IMP += 2; + } + EMP += 2 * volume_width; + IMP += 2 * volume_width; + } + + if (params->x_connectivity == 1) + { + //extend the mask to the front side of the phase volume + IMP = input_mask + frame_size + volume_width; //input mask pointer + EMP = extended_mask + frame_size + volume_width; //extended mask pointer + for (n=1; n < volume_depth - 1; n++) + { + for (i=1; i < volume_height - 1; i++) + { + if( (*IMP) == NOMASK && (*(IMP + vw - 1) == NOMASK) && (*(IMP + 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP + vw) == NOMASK) && + (*(IMP - fs) == NOMASK) && (*(IMP + fs) == NOMASK) && + (*(IMP - 1) == NOMASK) && (*(IMP + vw + 1) == NOMASK) && + (*(IMP - vw + 1) == NOMASK) && (*(IMP + 2 * vw - 1) == NOMASK) && + (*(IMP - fs - 1) == NOMASK) && (*(IMP + fs + vw + 1) == NOMASK) && + (*(IMP - fs - vw) == NOMASK) && (*(IMP + fs + vw) == NOMASK) && + (*(IMP - fs - vw + 1) == NOMASK) && (*(IMP + fs + 2 * vw - 1) == NOMASK) && + (*(IMP - fs + vw - 1) == NOMASK) && (*(IMP + fs + 1) == NOMASK) && + (*(IMP - fs + 1) == NOMASK) && (*(IMP + fs + vw - 1) == NOMASK) && + (*(IMP - fs + 2 * vw - 1) == NOMASK) && (*(IMP + fs - vw + 1) == NOMASK) && + (*(IMP - fs + vw) == NOMASK) && (*(IMP + fs - vw) == NOMASK) && + (*(IMP - fs + vw + 1) == NOMASK) && (*(IMP + fs - 1) == NOMASK) ) + { + *EMP = NOMASK; + } + EMP += vw; + IMP += vw; + } + EMP += 2 * vw; + IMP += 2 *vw; + } + + //extend the mask to the rear side of the phase volume + IMP = input_mask + frame_size + 2 * volume_width - 1; //input mask pointer + EMP = extended_mask + frame_size + 2 * volume_width - 1; //extended mask pointer + for (n=1; n < volume_depth - 1; n++) + { + for (i=1; i < volume_height - 1; i++) + { + if( (*IMP) == NOMASK && (*(IMP - vw + 1) == NOMASK) && (*(IMP - 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP + vw) == NOMASK) && + (*(IMP - fs) == NOMASK) && (*(IMP + fs) == NOMASK) && + (*(IMP - vw - 1) == NOMASK) && (*(IMP + 1) == NOMASK) && + (*(IMP + vw - 1) == NOMASK) && (*(IMP - 2 * vw + 1) == NOMASK) && + (*(IMP - fs - vw - 1) == NOMASK) && (*(IMP + fs + 1) == NOMASK) && + (*(IMP - fs - 2 * vw + 1) == NOMASK) && (*(IMP + fs + vw - 1) == NOMASK) && + (*(IMP - fs - 1) == NOMASK) && (*(IMP + fs - vw + 1) == NOMASK) && + (*(IMP - fs - vw + 1) == NOMASK) && (*(IMP + fs - 1) == NOMASK) && + (*(IMP - fs - vw) == NOMASK) && (*(IMP + fs + vw) == NOMASK) && + (*(IMP - fs + vw - 1) == NOMASK) && (*(IMP + fs - 2 * vw + 1) == NOMASK) && + (*(IMP - fs + vw) == NOMASK) && (*(IMP + fs - vw) == NOMASK) && + (*(IMP - fs + 1) == NOMASK) && (*(IMP + fs - vw - 1) == NOMASK) ) + { + *EMP = NOMASK; + } + EMP += vw; + IMP += vw; + } + EMP += 2 * vw; + IMP += 2 *vw; + } + } + + if (params->y_connectivity == 1) + { + //extend the mask to the left side of the phase volume + IMP = input_mask + frame_size + 1; + EMP = extended_mask + frame_size + 1; + for (n=1; n < volume_depth - 1; n++) + { + for (j=1; j < volume_width - 1; j++) + { + if( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) && + (*(IMP + fs - vw) == NOMASK) && (*(IMP + vw) == NOMASK) && + (*(IMP - fs) == NOMASK) && (*(IMP + fs) == NOMASK) && + (*(IMP + fs - vw - 1) == NOMASK) && (*(IMP + vw + 1) == NOMASK) && + (*(IMP + fs - vw + 1) == NOMASK) && (*(IMP + vw - 1) == NOMASK) && + (*(IMP - vw - 1) == NOMASK) && (*(IMP + fs + vw + 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP + fs + vw) == NOMASK) && + (*(IMP - vw + 1) == NOMASK) && (*(IMP + fs + vw - 1) == NOMASK) && + (*(IMP - fs - 1) == NOMASK) && (*(IMP + fs + 1) == NOMASK) && + (*(IMP - fs + 1) == NOMASK) && (*(IMP + fs - 1) == NOMASK) && + (*(IMP - fs + vw - 1) == NOMASK) && (*(IMP + 2 * fs - vw + 1) == NOMASK) && + (*(IMP - fs + vw) == NOMASK) && (*(IMP + 2 * fs - vw) == NOMASK) && + (*(IMP - fs + vw + 1) == NOMASK) && (*(IMP + 2 * fs - vw - 1) == NOMASK) ) + { + *EMP = NOMASK; + } + EMP++; + IMP++; + } + EMP += fs - vw + 2; + IMP += fs - vw + 2; + } + + //extend the mask to the right side of the phase volume + IMP = input_mask + 2 * frame_size - volume_width + 1; + EMP = extended_mask + 2 * frame_size - volume_width + 1; + for (n=1; n < volume_depth - 1; n++) + { + for (j=1; j < volume_width - 1; j++) + { + if( (*IMP) == NOMASK && (*(IMP + 1) == NOMASK) && (*(IMP - 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP - fs + vw) == NOMASK) && + (*(IMP - fs) == NOMASK) && (*(IMP + fs) == NOMASK) && + (*(IMP - vw - 1) == NOMASK) && (*(IMP - fs + vw + 1) == NOMASK) && + (*(IMP - vw + 1) == NOMASK) && (*(IMP - fs + vw - 1) == NOMASK) && + (*(IMP - fs - vw - 1) == NOMASK) && (*(IMP + vw + 1) == NOMASK) && + (*(IMP - fs - vw + 1) == NOMASK) && (*(IMP + vw - 1) == NOMASK) && + (*(IMP - fs - vw) == NOMASK) && (*(IMP + vw) == NOMASK) && + (*(IMP - fs - 1) == NOMASK) && (*(IMP + fs + 1) == NOMASK) && + (*(IMP - fs + 1) == NOMASK) && (*(IMP + fs - 1) == NOMASK) && + (*(IMP - 2 * fs + vw - 1) == NOMASK) && (*(IMP + fs - vw + 1) == NOMASK) && + (*(IMP - 2 * fs + vw) == NOMASK) && (*(IMP + fs - vw) == NOMASK) && + (*(IMP - 2 * fs + vw + 1) == NOMASK) && (*(IMP + fs - vw - 1) == NOMASK) ) + { + *EMP = NOMASK; + } + EMP++; + IMP++; + } + EMP += fs - vw + 2; + IMP += fs - vw + 2; + } + } + + if (params->z_connectivity == 1) + { + //extend the mask to the bottom side of the phase volume + IMP = input_mask + volume_width + 1; + EMP = extended_mask + volume_width + 1; + for (i=1; i < volume_height - 1; ++i) + { + for (j=1; j < volume_width - 1; ++j) + { + if( (*IMP) == NOMASK && (*(IMP - 1) == NOMASK) && (*(IMP + 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP + vw) == NOMASK) && + (*(IMP + fs) == NOMASK) && (*(IMP + vs - fs) == NOMASK) && + (*(IMP - vw - 1) == NOMASK) && (*(IMP + vw + 1) == NOMASK) && + (*(IMP - vw + 1) == NOMASK) && (*(IMP + vw - 1) == NOMASK) && + (*(IMP + vs - fs - vw - 1) == NOMASK) && (*(IMP + fs + vw + 1) == NOMASK) && + (*(IMP + vs - fs - vw) == NOMASK) && (*(IMP + fs + vw) == NOMASK) && + (*(IMP + vs - fs - vw + 1) == NOMASK) && (*(IMP + fs + vw - 1) == NOMASK) && + (*(IMP + vs - fs - 1) == NOMASK) && (*(IMP + fs + 1) == NOMASK) && + (*(IMP + vs - fs + 1) == NOMASK) && (*(IMP + fs - 1) == NOMASK) && + (*(IMP + vs - fs + vw - 1) == NOMASK) && (*(IMP + fs - vw + 1) == NOMASK) && + (*(IMP + vs - fs + vw) == NOMASK) && (*(IMP + fs - vw) == NOMASK) && + (*(IMP + vs - fs + vw + 1) == NOMASK) && (*(IMP + fs - vw - 1) == NOMASK) ) + { + *EMP = NOMASK; + } + EMP++; + IMP++; + } + EMP += 2; + IMP += 2; + } + + //extend the mask to the top side of the phase volume + IMP = input_mask + volume_size - frame_size + volume_width + 1; + EMP = extended_mask + volume_size - frame_size + volume_width + 1; + for (i=1; i < volume_height - 1; ++i) + { + for (j=1; j < volume_width - 1; ++j) + { + if( (*IMP) == NOMASK && (*(IMP + 1) == NOMASK) && (*(IMP - 1) == NOMASK) && + (*(IMP - vw) == NOMASK) && (*(IMP - fs + vw) == NOMASK) && + (*(IMP - fs) == NOMASK) && (*(IMP - vs + fs) == NOMASK) && + (*(IMP - vw - 1) == NOMASK) && (*(IMP + vw + 1) == NOMASK) && + (*(IMP - vw + 1) == NOMASK) && (*(IMP + vw - 1) == NOMASK) && + (*(IMP - fs - vw - 1) == NOMASK) && (*(IMP - vs + fs + vw + 1) == NOMASK) && + (*(IMP - fs - vw + 1) == NOMASK) && (*(IMP - vs + fs + vw - 1) == NOMASK) && + (*(IMP - fs - vw) == NOMASK) && (*(IMP - vs + fs + vw) == NOMASK) && + (*(IMP - fs - 1) == NOMASK) && (*(IMP - vs + fs + 1) == NOMASK) && + (*(IMP - fs + 1) == NOMASK) && (*(IMP - vs + fs - 1) == NOMASK) && + (*(IMP - fs + vw - 1) == NOMASK) && (*(IMP - vs + fs - vw + 1) == NOMASK) && + (*(IMP - fs + vw) == NOMASK) && (*(IMP - vs + fs - vw) == NOMASK) && + (*(IMP - fs + vw + 1) == NOMASK) && (*(IMP - vs + fs - vw - 1) == NOMASK) ) + { + *EMP = NOMASK; + } + EMP++; + IMP++; + } + EMP += 2; + IMP += 2; + } + } +} + +void calculate_reliability(double *wrappedVolume, VOXELM *voxel, int volume_width, int volume_height, int volume_depth, params_t *params) +{ + int frame_size = volume_width * volume_height; + int volume_size = volume_width * volume_height * volume_depth; + VOXELM *voxel_pointer; + double H, V, N, D1, D2, D3, D4, D5, D6, D7, D8, D9, D10; + double *WVP; + int n, i, j; + + WVP = wrappedVolume + frame_size + volume_width + 1; + voxel_pointer = voxel + frame_size + volume_width + 1; + for (n=1; n < volume_depth - 1; n++) + { + for (i=1; i < volume_height - 1; i++) + { + for (j=1; j < volume_width - 1; j++) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP - 1) - *WVP) - wrap(*WVP - *(WVP + 1)); + V = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + N = wrap(*(WVP - frame_size) - *WVP) - wrap(*WVP - *(WVP + frame_size)); + D1 = wrap(*(WVP - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + volume_width + 1)); + D2 = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + volume_width - 1)); + D3 = wrap(*(WVP - frame_size - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width + 1)); + D4 = wrap(*(WVP - frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width)); + D5 = wrap(*(WVP - frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width - 1)); + D6 = wrap(*(WVP - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 1)); + D7 = wrap(*(WVP - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 1)); + D8 = wrap(*(WVP - frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width + 1)); + D9 = wrap(*(WVP - frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width)); + D10 = wrap(*(WVP - frame_size + volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer++; + WVP++; + } + voxel_pointer += 2; + WVP += 2; + } + voxel_pointer += 2 * volume_width; + WVP += 2 * volume_width; + } + + if (params->x_connectivity == 1) + { + //calculating reliability for the front side of the phase volume...add volume_width + WVP = wrappedVolume + frame_size + volume_width; + voxel_pointer = voxel + frame_size + volume_width; + for (n=1; n < volume_depth - 1; ++n) + { + for (i=1; i < volume_height - 1; ++i) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + 1)); + V = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + N = wrap(*(WVP - frame_size) - *WVP) - wrap(*WVP - *(WVP + frame_size)); + D1 = wrap(*(WVP - 1) - *WVP) - wrap(*WVP - *(WVP + volume_width + 1)); + D2 = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + 2 * volume_width - 1)); + D3 = wrap(*(WVP - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width + 1)); + D4 = wrap(*(WVP - frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width)); + D5 = wrap(*(WVP - frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 2 * volume_width - 1)); + D6 = wrap(*(WVP - frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 1)); + D7 = wrap(*(WVP - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width - 1)); + D8 = wrap(*(WVP - frame_size + 2 * volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width + 1)); + D9 = wrap(*(WVP - frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width)); + D10 = wrap(*(WVP - frame_size + volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer += volume_width; + WVP += volume_width; + } + voxel_pointer += 2 * volume_width; + WVP += 2 * volume_width; + } + + //calculating reliability for the rear side of the phase volume..... subtract volume_width + WVP = wrappedVolume + frame_size + 2 * volume_width - 1; + voxel_pointer = voxel + frame_size + 2 * volume_width - 1; + for (n=1; n < volume_depth - 1; ++n) + { + for (i=1; i < volume_height - 1; ++i) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP - 1)); + V = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + N = wrap(*(WVP - frame_size) - *WVP) - wrap(*WVP - *(WVP + frame_size)); + D1 = wrap(*(WVP - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + 1)); + D2 = wrap(*(WVP + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP - 2 * volume_width + 1)); + D3 = wrap(*(WVP - frame_size - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 1)); + D4 = wrap(*(WVP - frame_size - 2 * volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width - 1)); + D5 = wrap(*(WVP - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width + 1)); + D6 = wrap(*(WVP - frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 1)); + D7 = wrap(*(WVP - frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width)); + D8 = wrap(*(WVP - frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 2 * volume_width + 1)); + D9 = wrap(*(WVP - frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width)); + D10 = wrap(*(WVP - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer += volume_width; + WVP += volume_width; + } + voxel_pointer += 2 * volume_width; + WVP += 2 * volume_width; + } + } + + if (params->y_connectivity == 1) + { + //calculating reliability for the left side of the phase volume...add frame_size + WVP = wrappedVolume + frame_size + 1; + voxel_pointer = voxel + frame_size + 1; + for (n=1; n < volume_depth - 1; ++n) + { + for (j=1; j < volume_width - 1; ++j) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP - 1) - *WVP) - wrap(*WVP - *(WVP + 1)); + V = wrap(*(WVP + frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + N = wrap(*(WVP - frame_size) - *WVP) - wrap(*WVP - *(WVP + frame_size)); + D1 = wrap(*(WVP + frame_size - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + volume_width + 1)); + D2 = wrap(*(WVP + frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + volume_width - 1)); + D3 = wrap(*(WVP - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width + 1)); + D4 = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width)); + D5 = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width - 1)); + D6 = wrap(*(WVP - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 1)); + D7 = wrap(*(WVP - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 1)); + D8 = wrap(*(WVP - frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + 2 * frame_size - volume_width + 1)); + D9 = wrap(*(WVP - frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP + 2 * frame_size - volume_width)); + D10 = wrap(*(WVP - frame_size + volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + 2 * frame_size - volume_width - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer++; + WVP++; + } + voxel_pointer += frame_size - volume_width + 2; + WVP += frame_size - volume_width + 2; + } + + //calculating reliability for the right side of the phase volume...subtract frame_size + WVP = wrappedVolume + 2 * frame_size - volume_width + 1; + voxel_pointer = voxel + 2 * frame_size - volume_width + 1; + for (n=1; n < volume_depth - 1; ++n) + { + for (j=1; j < volume_width - 1; ++j) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP + 1) - *WVP) - wrap(*WVP - *(WVP - 1)); + V = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP - frame_size + volume_width)); + N = wrap(*(WVP - frame_size) - *WVP) - wrap(*WVP - *(WVP + frame_size)); + D1 = wrap(*(WVP - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP - frame_size + volume_width + 1)); + D2 = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP - frame_size + volume_width - 1)); + D3 = wrap(*(WVP - frame_size - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + volume_width + 1) ); + D4 = wrap(*(WVP - frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + volume_width - 1)); + D5 = wrap(*(WVP - frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + D6 = wrap(*(WVP - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 1)); + D7 = wrap(*(WVP - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 1)); + D8 = wrap(*(WVP - 2 * frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width + 1)); + D9 = wrap(*(WVP - 2 * frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width)); + D10 = wrap(*(WVP - 2 * frame_size + volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer++; + WVP++; + } + voxel_pointer += frame_size - volume_width + 2; + WVP += frame_size - volume_width + 2; + } + } + + if (params->z_connectivity == 1) + { + //calculating reliability for the bottom side of the phase volume...add volume_size + WVP = wrappedVolume + volume_width + 1; + voxel_pointer = voxel + volume_width + 1; + for (i=1; i < volume_height - 1; ++i) + { + for (j=1; j < volume_width - 1; ++j) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP - 1) - *WVP) - wrap(*WVP - *(WVP + 1)); + V = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + N = wrap(*(WVP + frame_size) - *WVP) - wrap(*WVP - *(WVP + volume_size - frame_size)); + D1 = wrap(*(WVP - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + volume_width + 1)); + D2 = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + volume_width - 1)); + D3 = wrap(*(WVP + volume_size - frame_size - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width + 1)); + D4 = wrap(*(WVP + volume_size - frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width)); + D5 = wrap(*(WVP + volume_size - frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + volume_width - 1)); + D6 = wrap(*(WVP + volume_size - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size + 1)); + D7 = wrap(*(WVP + volume_size - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - 1)); + D8 = wrap(*(WVP + volume_size - frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width + 1)); + D9 = wrap(*(WVP + volume_size - frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width)); + D10 = wrap(*(WVP + volume_size - frame_size + volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + frame_size - volume_width - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer++; + WVP++; + } + voxel_pointer += 2; + WVP += 2; + } + + //calculating reliability for the top side of the phase volume...subtract volume_size + WVP = wrappedVolume + volume_size - frame_size + volume_width + 1; + voxel_pointer = voxel + volume_size - frame_size + volume_width + 1; + for (i=1; i < volume_height - 1; ++i) + { + for (j=1; j < volume_width - 1; ++j) + { + if (voxel_pointer->extended_mask == NOMASK) + { + H = wrap(*(WVP + 1) - *WVP) - wrap(*WVP - *(WVP - 1)); + V = wrap(*(WVP - volume_width) - *WVP) - wrap(*WVP - *(WVP + volume_width)); + N = wrap(*(WVP - frame_size) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size)); + D1 = wrap(*(WVP - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP + volume_width + 1)); + D2 = wrap(*(WVP - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP + volume_width - 1)); + D3 = wrap(*(WVP - frame_size - volume_width - 1) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size + volume_width + 1)); + D4 = wrap(*(WVP - frame_size - volume_width + 1) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size + volume_width - 1)); + D5 = wrap(*(WVP - frame_size - volume_width) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size + volume_width)); + D6 = wrap(*(WVP - frame_size - 1) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size + 1)); + D7 = wrap(*(WVP - frame_size + 1) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size - 1)); + D8 = wrap(*(WVP - frame_size + volume_width - 1) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size - volume_width + 1)); + D9 = wrap(*(WVP - frame_size + volume_width) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size - volume_width)); + D10 = wrap(*(WVP - frame_size + volume_width + 1) - *WVP) - wrap(*WVP - *(WVP - volume_size + frame_size - volume_width - 1)); + voxel_pointer->reliability = H*H + V*V + N*N + D1*D1 + D2*D2 + D3*D3 + D4*D4 + D5*D5 + D6*D6 + + D7*D7 + D8*D8 + D9*D9 + D10*D10; + } + voxel_pointer++; + WVP++; + } + voxel_pointer += 2; + WVP += 2; + } + } +} + +//calculate the reliability of the horizontal edges of the volume. it +//is calculated by adding the reliability of voxel and the relibility +//of its right neighbour. edge is calculated between a voxel and its +//next neighbour +void horizontalEDGEs(VOXELM *voxel, EDGE *edge, int volume_width, int volume_height, int volume_depth, params_t *params) +{ + int n, i, j; + EDGE *edge_pointer = edge; + VOXELM *voxel_pointer = voxel; + int no_of_edges = params->no_of_edges; + + for (n=0; n < volume_depth; n++) + { + for (i = 0; i < volume_height; i++) + { + for (j = 0; j < volume_width - 1; j++) + { + if (voxel_pointer->input_mask == NOMASK && (voxel_pointer + 1)->input_mask == NOMASK ) + { + edge_pointer->pointer_1 = voxel_pointer; + edge_pointer->pointer_2 = (voxel_pointer+1); + edge_pointer->reliab = voxel_pointer->reliability + (voxel_pointer + 1)->reliability; + edge_pointer->increment = find_wrap(voxel_pointer->value, (voxel_pointer + 1)->value); + edge_pointer++; + no_of_edges++; + } + voxel_pointer++; + } + voxel_pointer++; + } + } + if (params->x_connectivity == 1) + { + voxel_pointer = voxel + volume_width - 1; + for (n=0; n < volume_depth; n++) + { + for (i = 0; i < volume_height; i++) + { + if (voxel_pointer->input_mask == NOMASK && (voxel_pointer - volume_width + 1)->input_mask == NOMASK ) + { + edge_pointer->pointer_1 = voxel_pointer; + edge_pointer->pointer_2 = (voxel_pointer - volume_width + 1); + edge_pointer->reliab = voxel_pointer->reliability + (voxel_pointer - volume_width + 1)->reliability; + edge_pointer->increment = find_wrap(voxel_pointer->value, (voxel_pointer - volume_width + 1)->value); + edge_pointer++; + no_of_edges++; + } + voxel_pointer += volume_width; + } + } + } + params->no_of_edges = no_of_edges; +} + +void verticalEDGEs(VOXELM *voxel, EDGE *edge, int volume_width, int volume_height, int volume_depth, params_t *params) +{ + int n, i, j; + int no_of_edges = params->no_of_edges; + VOXELM *voxel_pointer = voxel; + EDGE *edge_pointer = edge + no_of_edges; + int frame_size = volume_width * volume_height; + int next_voxel = frame_size - volume_width; + + for (n=0; n < volume_depth; n++) + { + for (i=0; iinput_mask == NOMASK && (voxel_pointer + volume_width)->input_mask == NOMASK ) + { + edge_pointer->pointer_1 = voxel_pointer; + edge_pointer->pointer_2 = (voxel_pointer + volume_width); + edge_pointer->reliab = voxel_pointer->reliability + (voxel_pointer + volume_width)->reliability; + edge_pointer->increment = find_wrap(voxel_pointer->value, (voxel_pointer + volume_width)->value); + edge_pointer++; + no_of_edges++; + } + voxel_pointer++; + } + } + voxel_pointer += volume_width; + } + + if (params->y_connectivity == 1) + { + voxel_pointer = voxel + frame_size - volume_width; + for (n=0; n < volume_depth; n++) + { + for (i = 0; i < volume_width; i++) + { + if (voxel_pointer->input_mask == NOMASK && (voxel_pointer - next_voxel)->input_mask == NOMASK ) + { + edge_pointer->pointer_1 = voxel_pointer; + edge_pointer->pointer_2 = (voxel_pointer - next_voxel); + edge_pointer->reliab = voxel_pointer->reliability + (voxel_pointer - next_voxel)->reliability; + edge_pointer->increment = find_wrap(voxel_pointer->value, (voxel_pointer - next_voxel)->value); + edge_pointer++; + no_of_edges++; + } + voxel_pointer++; + } + voxel_pointer += next_voxel + 1; + } + } + params->no_of_edges = no_of_edges; +} + +void normalEDGEs(VOXELM *voxel, EDGE *edge, int volume_width, int volume_height, int volume_depth, params_t *params) +{ + int n, i, j; + int no_of_edges = params->no_of_edges; + int frame_size = volume_width * volume_height; + int volume_size = volume_width * volume_height * volume_depth; + VOXELM *voxel_pointer = voxel; + EDGE *edge_pointer = edge + no_of_edges; + int next_voxel = volume_size - frame_size; + + for (n=0; n < volume_depth - 1; n++) + { + for (i=0; iinput_mask == NOMASK && (voxel_pointer + frame_size)->input_mask == NOMASK ) + { + edge_pointer->pointer_1 = voxel_pointer; + edge_pointer->pointer_2 = (voxel_pointer + frame_size); + edge_pointer->reliab = voxel_pointer->reliability + (voxel_pointer + frame_size)->reliability; + edge_pointer->increment = find_wrap(voxel_pointer->value, (voxel_pointer + frame_size)->value); + edge_pointer++; + no_of_edges++; + } + voxel_pointer++; + } + } + } + + + if (params->z_connectivity == 1) + { + voxel_pointer = voxel + next_voxel; + for (i=0; i < volume_height; i++) + { + for (j = 0; j < volume_width; j++) + { + if (voxel_pointer->input_mask == NOMASK && (voxel_pointer - next_voxel)->input_mask == NOMASK ) + { + edge_pointer->pointer_1 = voxel_pointer; + edge_pointer->pointer_2 = (voxel_pointer - next_voxel); + edge_pointer->reliab = voxel_pointer->reliability + (voxel_pointer - next_voxel)->reliability; + edge_pointer->increment = find_wrap(voxel_pointer->value, (voxel_pointer - next_voxel)->value); + edge_pointer++; + no_of_edges++; + } + voxel_pointer++; + } + } + } + params->no_of_edges = no_of_edges; +} + +//gather the voxels of the volume into groups +void gatherVOXELs(EDGE *edge, params_t *params) +{ + int k; + VOXELM *VOXEL1; + VOXELM *VOXEL2; + VOXELM *group1; + VOXELM *group2; + EDGE *pointer_edge = edge; + int incremento; + + for (k = 0; k < params->no_of_edges; k++) + { + VOXEL1 = pointer_edge->pointer_1; + VOXEL2 = pointer_edge->pointer_2; + + //VOXELM 1 and VOXELM 2 belong to different groups + //initially each voxel is in a group by itself and one voxel can construct a group + //no else or else if to this if + if (VOXEL2->head != VOXEL1->head) + { + //VOXELM 2 is alone in its group + //merge this voxel with VOXELM 1 group and find the number of 2 pi to add + //to or subtract to unwrap it + if ((VOXEL2->next == NULL) && (VOXEL2->head == VOXEL2)) + { + VOXEL1->head->last->next = VOXEL2; + VOXEL1->head->last = VOXEL2; + (VOXEL1->head->number_of_voxels_in_group)++; + VOXEL2->head=VOXEL1->head; + VOXEL2->increment = VOXEL1->increment-pointer_edge->increment; + } + + //VOXELM 1 is alone in its group + //merge this voxel with VOXELM 2 group and find the number of 2 pi to add + //to or subtract to unwrap it + else if ((VOXEL1->next == NULL) && (VOXEL1->head == VOXEL1)) + { + VOXEL2->head->last->next = VOXEL1; + VOXEL2->head->last = VOXEL1; + (VOXEL2->head->number_of_voxels_in_group)++; + VOXEL1->head = VOXEL2->head; + VOXEL1->increment = VOXEL2->increment+pointer_edge->increment; + } + + //VOXELM 1 and VOXELM 2 both have groups + else + { + group1 = VOXEL1->head; + group2 = VOXEL2->head; + //if the no. of voxels in VOXELM 1 group is larger than the no. of voxels + //in VOXELM 2 group. Merge VOXELM 2 group to VOXELM 1 group + //and find the number of wraps between VOXELM 2 group and VOXELM 1 group + //to unwrap VOXELM 2 group with respect to VOXELM 1 group. + //the no. of wraps will be added to VOXELM 2 grop in the future + if (group1->number_of_voxels_in_group > group2->number_of_voxels_in_group) + { + //merge VOXELM 2 with VOXELM 1 group + group1->last->next = group2; + group1->last = group2->last; + group1->number_of_voxels_in_group = group1->number_of_voxels_in_group + group2->number_of_voxels_in_group; + incremento = VOXEL1->increment-pointer_edge->increment - VOXEL2->increment; + //merge the other voxels in VOXELM 2 group to VOXELM 1 group + while (group2 != NULL) + { + group2->head = group1; + group2->increment += incremento; + group2 = group2->next; + } + } + + //if the no. of voxels in VOXELM 2 group is larger than the no. of voxels + //in VOXELM 1 group. Merge VOXELM 1 group to VOXELM 2 group + //and find the number of wraps between VOXELM 2 group and VOXELM 1 group + //to unwrap VOXELM 1 group with respect to VOXELM 2 group. + //the no. of wraps will be added to VOXELM 1 grop in the future + else + { + //merge VOXELM 1 with VOXELM 2 group + group2->last->next = group1; + group2->last = group1->last; + group2->number_of_voxels_in_group = group2->number_of_voxels_in_group + group1->number_of_voxels_in_group; + incremento = VOXEL2->increment + pointer_edge->increment - VOXEL1->increment; + //merge the other voxels in VOXELM 2 group to VOXELM 1 group + while (group1 != NULL) + { + group1->head = group2; + group1->increment += incremento; + group1 = group1->next; + } // while + + } // else + } //else + } //if + pointer_edge++; + } +} + +//unwrap the volume +void unwrapVolume(VOXELM *voxel, int volume_width, int volume_height, int volume_depth) +{ + int i; + int volume_size = volume_width * volume_height * volume_depth; + VOXELM *voxel_pointer=voxel; + + for (i = 0; i < volume_size; i++) + { + voxel_pointer->value += TWOPI * (double)(voxel_pointer->increment); + voxel_pointer++; + } +} + +//set the masked voxels (mask = 0) to the minimum of the unwrapper phase +void maskVolume(VOXELM *voxel, unsigned char *input_mask, int volume_width, int volume_height, int volume_depth) +{ + int volume_width_plus_one = volume_width + 1; + int volume_height_plus_one = volume_height + 1; + int volume_width_minus_one = volume_width - 1; + int volume_height_minus_one = volume_height - 1; + + VOXELM *pointer_voxel = voxel; + unsigned char *IMP = input_mask; //input mask pointer + double min=99999999.; + int i, j; + int volume_size = volume_width * volume_height * volume_depth; + + //find the minimum of the unwrapped phase + for (i = 0; i < volume_size; i++) + { + if ((pointer_voxel->value < min) && (*IMP == NOMASK)) + min = pointer_voxel->value; + + pointer_voxel++; + IMP++; + } + + pointer_voxel = voxel; + IMP = input_mask; + + //set the masked voxels to minimum + for (i = 0; i < volume_size; i++) + { + if ((*IMP) == MASK) + { + pointer_voxel->value = min; + } + pointer_voxel++; + IMP++; + } +} + +//the input to this unwrapper is an array that contains the wrapped +//phase map. copy the volume on the buffer passed to this unwrapper +//to over-write the unwrapped phase map on the buffer of the wrapped +//phase map. +void returnVolume(VOXELM *voxel, double *unwrappedVolume, int volume_width, int volume_height, int volume_depth) +{ + int i; + int volume_size = volume_width * volume_height * volume_depth; + double *unwrappedVolume_pointer = unwrappedVolume; + VOXELM *voxel_pointer = voxel; + + for (i=0; i < volume_size; i++) + { + *unwrappedVolume_pointer = voxel_pointer->value; + voxel_pointer++; + unwrappedVolume_pointer++; + } +} + +//the main function of the unwrapper +void +unwrap3D(double* wrapped_volume, double* unwrapped_volume, unsigned char* input_mask, + int volume_width, int volume_height, int volume_depth, + int wrap_around_x, int wrap_around_y, int wrap_around_z) +{ + params_t params = {TWOPI, wrap_around_x, wrap_around_y, wrap_around_z, 0}; + unsigned char *extended_mask; + VOXELM *voxel; + EDGE *edge; + int volume_size = volume_height * volume_width * volume_depth; + int No_of_Edges_initially = 3 * volume_width * volume_height * volume_depth; + + extended_mask = (unsigned char *) calloc(volume_size, sizeof(unsigned char)); + voxel = (VOXELM *) calloc(volume_size, sizeof(VOXELM)); + edge = (EDGE *) calloc(No_of_Edges_initially, sizeof(EDGE));; + + extend_mask(input_mask, extended_mask, volume_width, volume_height, volume_depth, ¶ms); + initialiseVOXELs(wrapped_volume, input_mask, extended_mask, voxel, volume_width, volume_height, volume_depth); + calculate_reliability(wrapped_volume, voxel, volume_width, volume_height, volume_depth, ¶ms); + horizontalEDGEs(voxel, edge, volume_width, volume_height, volume_depth, ¶ms); + verticalEDGEs(voxel, edge, volume_width, volume_height, volume_depth, ¶ms); + normalEDGEs(voxel, edge, volume_width, volume_height, volume_depth, ¶ms); + + //sort the EDGEs depending on their reiability. The VOXELs with higher relibility (small value) first + quicker_sort(edge, edge + params.no_of_edges - 1); + + //gather VOXELs into groups + gatherVOXELs(edge, ¶ms); + + unwrapVolume(voxel, volume_width, volume_height, volume_depth); + maskVolume(voxel, input_mask, volume_width, volume_height, volume_depth); + + //copy the volume from VOXELM structure to the unwrapped phase array passed to this function + returnVolume(voxel, unwrapped_volume, volume_width, volume_height, volume_depth); + + free(edge); + free(voxel); + free(extended_mask); +} diff --git a/skimage/setup.py b/skimage/setup.py index 1082ba07..14f66897 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -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')