Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Make unwrap_phase fail for 2D input containing NaNs #7176

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Make unwrap_phase fail for 2D input containing NaNs #7176

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
20 changes: 13 additions & 7 deletions skimage/restoration/_unwrap_2d.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ cimport numpy as cnp


cdef extern from "unwrap_2d_ljmu.h":
void unwrap2D(
int unwrap2D(
cnp.float64_t *wrapped_image,
cnp.float64_t *UnwrappedImage,
unsigned char *input_mask,
Expand All @@ -27,15 +27,21 @@ def unwrap_2d(cnp.float64_t[:, ::1] image,
char use_seed
int wrap_around_x
int wrap_around_y
int error_code = 0

# convert from python types to C types so we can release the GIL
use_seed = seed is None
cseed = 0 if seed is None else seed
wrap_around_y, wrap_around_x = wrap_around
with nogil:
unwrap2D(&image[0, 0],
&unwrapped_image[0, 0],
&mask[0, 0],
image.shape[1], image.shape[0],
wrap_around_x, wrap_around_y,
use_seed, cseed)
error_code = unwrap2D(&image[0, 0],
&unwrapped_image[0, 0],
&mask[0, 0],
image.shape[1], image.shape[0],
wrap_around_x, wrap_around_y,
use_seed, cseed)
if error_code == 1:
raise RuntimeError(
"NaN encountered while sorting edges which may happen if `image` "
"contains unmasked NaNs."
)
15 changes: 15 additions & 0 deletions skimage/restoration/tests/test_unwrap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@
from skimage.restoration import unwrap_phase
import sys

import pytest

from skimage._shared import testing
from skimage._shared.testing import (assert_array_almost_equal_nulp,
assert_almost_equal, assert_array_equal,
Expand Down Expand Up @@ -215,3 +217,16 @@ def test_unwrap_3d_all_masked():
assert_(np.ma.isMaskedArray(unwrap))
assert_(np.sum(unwrap.mask) == 999) # all but one masked
assert_(unwrap[0, 0, 0] == 0)


def test_unwrap_2d_raise_on_nan():
"""Check that unwrap_phase fails for NaN input.

Regression test for
https://github.com/scikit-image/scikit-image/issues/3831.
"""
x = np.linspace(1, 100)
x[1] = np.nan
xx, yy = np.meshgrid(x, x)
with pytest.raises(RuntimeError, match="NaN encountered while sorting edges"):
unwrap_phase(xx)
3 changes: 3 additions & 0 deletions skimage/restoration/unwrap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,9 @@ def unwrap_phase(image, wrap_around=False, rng=None):
ValueError
If called with a masked 1D array or called with a 1D array and
``wrap_around=True``.
RuntimeError
If a NaN is encountered while sorting edges (abrupt phase steps). This
may happen if `image` contains unmasked NaNs.

Examples
--------
Expand Down
30 changes: 22 additions & 8 deletions skimage/restoration/unwrap_2d_ljmu.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
#define NOMASK 0
#define MASK 1

#define NAN_ERROR 1

typedef struct {
double mod;
Expand Down Expand Up @@ -138,15 +139,23 @@ EDGE *partition(EDGE *left, EDGE *right, double pivot) {
return left;
}

void quicker_sort(EDGE *left, EDGE *right) {
int quicker_sort(EDGE *left, EDGE *right) {
EDGE *p;
double pivot;
int error_code = 0;

if (find_pivot(left, right, &pivot) == yes) {
if (isnan(pivot)) {
return NAN_ERROR;
}
Comment on lines +148 to +150
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure if that is the best way to address the problem. I understand little of what's going on, but it's better than hanging with one core stuck at 100%. Performance wise I didn't see any significant impact for the following snippet taken from plot_phase_unwrap.py:

import numpy as np
from skimage import data, img_as_float, color, exposure
from skimage.restoration import unwrap_phase
image = color.rgb2gray(img_as_float(data.chelsea()))
image = exposure.rescale_intensity(image, out_range=(0, 4 * np.pi))
image_wrapped = np.angle(np.exp(1j * image))

%timeit unwrap_phase(image_wrapped)

Though I am a bit dubious because the paper makes it seem like there should be a huge amount of edges and I am not sure how often quicker_sort is called in the total scheme of things.

p = partition(left, right, pivot);
quicker_sort(left, p - 1);
quicker_sort(p, right);
error_code = quicker_sort(left, p - 1);
if (error_code != 0) {
return error_code;
}
error_code = quicker_sort(p, right);
}
return error_code;
}
//--------------end quicker_sort algorithm -----------------------------------

Expand Down Expand Up @@ -691,16 +700,17 @@ 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,
char use_seed, unsigned int seed) {
int unwrap2D(double *wrapped_image, double *UnwrappedImage,
unsigned char *input_mask, int image_width, int image_height,
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;
EDGE *edge;
int image_size = image_height * image_width;
int No_of_Edges_initially = 2 * image_width * image_height;
int error_code = 0;

extended_mask = (unsigned char *)calloc(image_size, sizeof(unsigned char));
pixel = (PIXELM *)calloc(image_size, sizeof(PIXELM));
Expand All @@ -717,7 +727,10 @@ void unwrap2D(double *wrapped_image, double *UnwrappedImage,
if (params.no_of_edges != 0) {
// sort the EDGEs depending on their reiability. The PIXELs with higher
// relibility (small value) first
quicker_sort(edge, edge + params.no_of_edges - 1);
error_code = quicker_sort(edge, edge + params.no_of_edges - 1);
if (error_code != 0) {
return error_code;
}
}
// gather PIXELs into groups
gatherPIXELs(edge, &params);
Expand All @@ -733,4 +746,5 @@ void unwrap2D(double *wrapped_image, double *UnwrappedImage,
free(edge);
free(pixel);
free(extended_mask);
return error_code;
}
2 changes: 1 addition & 1 deletion skimage/restoration/unwrap_2d_ljmu.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
void unwrap2D(
short unwrap2D(
double *wrapped_image,
double *UnwrappedImage,
unsigned char *input_mask,
Expand Down