From 3b0651c202f2a13387ddacd4a73bc2f0e0710977 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 16 Dec 2014 16:36:30 +0100 Subject: [PATCH 1/2] Expose seed parameter, fix doc string --- skimage/restoration/_unwrap_2d.pyx | 8 +++--- skimage/restoration/_unwrap_3d.pyx | 8 +++--- skimage/restoration/unwrap.py | 40 ++++++++++++++++------------ skimage/restoration/unwrap_2d_ljmu.c | 18 ++++++++----- skimage/restoration/unwrap_3d_ljmu.c | 17 +++++++----- 5 files changed, 54 insertions(+), 37 deletions(-) diff --git a/skimage/restoration/_unwrap_2d.pyx b/skimage/restoration/_unwrap_2d.pyx index 6e486439..040baa03 100644 --- a/skimage/restoration/_unwrap_2d.pyx +++ b/skimage/restoration/_unwrap_2d.pyx @@ -3,15 +3,17 @@ cdef extern from *: double* unwrapped_image, unsigned char* input_mask, int image_width, int image_height, - int wrap_around_x, int wrap_around_y) + int wrap_around_x, int wrap_around_y, + char use_seed, unsigned int seed) def unwrap_2d(double[:, ::1] image, unsigned char[:, ::1] mask, double[:, ::1] unwrapped_image, - wrap_around): + wrap_around, + seed): unwrap2D(&image[0, 0], &unwrapped_image[0, 0], &mask[0, 0], image.shape[1], image.shape[0], wrap_around[1], wrap_around[0], - ) + seed is None, 0 if seed is None else seed) diff --git a/skimage/restoration/_unwrap_3d.pyx b/skimage/restoration/_unwrap_3d.pyx index 3a1869a8..bbae61a3 100644 --- a/skimage/restoration/_unwrap_3d.pyx +++ b/skimage/restoration/_unwrap_3d.pyx @@ -3,15 +3,17 @@ cdef extern from *: 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) + int wrap_around_x, int wrap_around_y, int wrap_around_z, + char use_seed, unsigned int seed) def unwrap_3d(double[:, :, ::1] image, unsigned char[:, :, ::1] mask, double[:, :, ::1] unwrapped_image, - wrap_around): + wrap_around, + seed): 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], - ) + seed is None, 0 if seed is None else seed) diff --git a/skimage/restoration/unwrap.py b/skimage/restoration/unwrap.py index f0025775..c6810229 100644 --- a/skimage/restoration/unwrap.py +++ b/skimage/restoration/unwrap.py @@ -7,29 +7,34 @@ from ._unwrap_2d import unwrap_2d from ._unwrap_3d import unwrap_3d -def unwrap_phase(image, wrap_around=False): - '''From ``image``, wrapped to lie in the interval [-pi, pi), recover the +def unwrap_phase(image, wrap_around=False, seed=None): + '''Recover the original from a wrapped phase image. + + From an 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 + 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 + `ValueError`. + wrap_around : bool or sequence of bool, optional + 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. + seed : int, optional + Unwrapping 2D or 3D images uses random initialization. This sets the + seed of the PRNG to achieve deterministic behavior. Returns ------- - image_unwrapped : array_like, float32 - Unwrapped image of the same shape as the input. If the input ``image`` + image_unwrapped : array_like, double + Unwrapped image of the same shape as the input. If the input `image` was a masked array, the mask will be preserved. Raises @@ -60,25 +65,25 @@ def unwrap_phase(image, wrap_around=False): International Society for Optics and Photonics. ''' if image.ndim not in (1, 2, 3): - raise ValueError('image must be 1, 2 or 3 dimensional') + 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 ' + 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 ' + 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') + 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 ' + warnings.warn('Image has a length 1 dimension. Consider using an ' 'array of lower dimensionality to use a more efficient ' 'algorithm') @@ -87,17 +92,18 @@ def unwrap_phase(image, wrap_around=False): image = image.data 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') + + image_not_masked = np.asarray(image, dtype=np.double, order='C') + image_unwrapped = np.empty_like(image, dtype=np.double, 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) + wrap_around, seed) elif image.ndim == 3: unwrap_3d(image_not_masked, mask, image_unwrapped, - wrap_around) + wrap_around, seed) if np.ma.isMaskedArray(image): return np.ma.array(image_unwrapped, mask=mask) diff --git a/skimage/restoration/unwrap_2d_ljmu.c b/skimage/restoration/unwrap_2d_ljmu.c index c2c47525..605c7993 100644 --- a/skimage/restoration/unwrap_2d_ljmu.c +++ b/skimage/restoration/unwrap_2d_ljmu.c @@ -24,6 +24,7 @@ #include #include #include +#include #ifndef M_PI #define M_PI 3.1415926535897932384626433832795 @@ -154,22 +155,24 @@ void quicker_sort(EDGE *left, EDGE *right) { // itself void initialisePIXELs(double *wrapped_image, unsigned char *input_mask, unsigned char *extended_mask, PIXELM *pixel, - int image_width, int image_height) { + int image_width, int image_height, + char use_seed, unsigned int seed) { 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; - // Make the initialization deterministic - srand(0); + if (use_seed) { + srand(seed); + } 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->reliability = rand(); pixel_pointer->input_mask = *input_mask_pointer; pixel_pointer->extended_mask = *extended_mask_pointer; pixel_pointer->head = pixel_pointer; @@ -647,7 +650,7 @@ void maskImage(PIXELM *pixel, unsigned char *input_mask, int image_width, PIXELM *pointer_pixel = pixel; unsigned char *IMP = input_mask; // input mask pointer - double min = 99999999; + double min = DBL_MAX; int i; int image_size = image_width * image_height; @@ -694,7 +697,8 @@ void returnImage(PIXELM *pixel, double *unwrapped_image, int image_width, // 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) { + int wrap_around_x, int wrap_around_y, + char use_seed, unsigned int seed) { params_t params = {TWOPI, wrap_around_x, wrap_around_y, 0}; unsigned char *extended_mask; PIXELM *pixel; @@ -708,7 +712,7 @@ void unwrap2D(double *wrapped_image, double *UnwrappedImage, extend_mask(input_mask, extended_mask, image_width, image_height, ¶ms); initialisePIXELs(wrapped_image, input_mask, extended_mask, pixel, image_width, - image_height); + image_height, use_seed, seed); calculate_reliability(wrapped_image, pixel, image_width, image_height, ¶ms); horizontalEDGEs(pixel, edge, image_width, image_height, ¶ms); diff --git a/skimage/restoration/unwrap_3d_ljmu.c b/skimage/restoration/unwrap_3d_ljmu.c index 2f498f66..be47f926 100644 --- a/skimage/restoration/unwrap_3d_ljmu.c +++ b/skimage/restoration/unwrap_3d_ljmu.c @@ -28,6 +28,7 @@ #include #include #include +#include #ifndef M_PI #define M_PI 3.1415926535897932384626433832795 @@ -159,15 +160,17 @@ void quicker_sort(EDGE *left, EDGE *right) { // itself void initialiseVOXELs(double *WrappedVolume, unsigned char *input_mask, unsigned char *extended_mask, VOXELM *voxel, - int volume_width, int volume_height, int volume_depth) { + int volume_width, int volume_height, int volume_depth, + char use_seed, unsigned int seed) { 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; - // Make the initialization deterministic - srand(0); + if (use_seed) { + srand(seed); + } for (n = 0; n < volume_depth; n++) { for (i = 0; i < volume_height; i++) { @@ -175,7 +178,7 @@ void initialiseVOXELs(double *WrappedVolume, unsigned char *input_mask, 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->reliability = rand(); voxel_pointer->input_mask = *input_mask_pointer; voxel_pointer->extended_mask = *extended_mask_pointer; voxel_pointer->head = voxel_pointer; @@ -1060,7 +1063,7 @@ void maskVolume(VOXELM *voxel, unsigned char *input_mask, int volume_width, VOXELM *pointer_voxel = voxel; unsigned char *IMP = input_mask; // input mask pointer - double min = 99999999.; + double min = DBL_MAX; int i, j; int volume_size = volume_width * volume_height * volume_depth; @@ -1108,7 +1111,7 @@ void returnVolume(VOXELM *voxel, double *unwrappedVolume, int volume_width, 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) { + int wrap_around_z, char use_seed, unsigned int seed) { params_t params = {TWOPI, wrap_around_x, wrap_around_y, wrap_around_z, 0}; unsigned char *extended_mask; VOXELM *voxel; @@ -1124,7 +1127,7 @@ void unwrap3D(double *wrapped_volume, double *unwrapped_volume, 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); + volume_width, volume_height, volume_depth, use_seed, seed); calculate_reliability(wrapped_volume, voxel, volume_width, volume_height, volume_depth, ¶ms); horizontalEDGEs(voxel, edge, volume_width, volume_height, volume_depth, From 03be934e6d4f19d09fe698336fa9722abdf36a61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Tue, 16 Dec 2014 16:38:01 +0100 Subject: [PATCH 2/2] Set seed for test cases --- skimage/restoration/tests/test_unwrap.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/skimage/restoration/tests/test_unwrap.py b/skimage/restoration/tests/test_unwrap.py index 8c2a6487..39723894 100644 --- a/skimage/restoration/tests/test_unwrap.py +++ b/skimage/restoration/tests/test_unwrap.py @@ -37,7 +37,7 @@ def check_unwrap(image, mask=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) + image_unwrapped = unwrap_phase(image_wrapped, seed=0) assert_phase_almost_equal(image_unwrapped, image) @@ -47,7 +47,7 @@ def test_unwrap_1d(): # 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) + assert_raises(ValueError, unwrap_phase, image, True, seed=0) def test_unwrap_2d(): @@ -83,7 +83,7 @@ def check_wrap_around(ndim, axis): 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) + image_unwrap_no_wrap_around = unwrap_phase(image_wrapped, seed=0) print('endpoints without wrap_around:', image_unwrap_no_wrap_around[index_first], image_unwrap_no_wrap_around[index_last]) @@ -95,7 +95,8 @@ def check_wrap_around(ndim, axis): 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) + image_unwrap_wrap_around = unwrap_phase(image_wrapped, wrap_around, + seed=0) print('endpoints with wrap_around:', image_unwrap_wrap_around[index_first], image_unwrap_wrap_around[index_last])