Merge pull request #1307 from ahojnnes/unwrap-rand

Remove random initialization from unwrap
This commit is contained in:
Stefan van der Walt
2014-12-16 20:08:48 +02:00
6 changed files with 59 additions and 41 deletions
+5 -3
View File
@@ -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)
+5 -3
View File
@@ -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)
+5 -4
View File
@@ -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])
+23 -17
View File
@@ -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)
+11 -7
View File
@@ -24,6 +24,7 @@
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#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, &params);
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,
&params);
horizontalEDGEs(pixel, edge, image_width, image_height, &params);
+10 -7
View File
@@ -28,6 +28,7 @@
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#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, &params);
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, &params);
horizontalEDGEs(voxel, edge, volume_width, volume_height, volume_depth,