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

Map coordinate fix #206

Closed
wants to merge 5 commits into from
Closed

Map coordinate fix #206

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
603eb3e
FIX: map_coordinate, fixed bug in Ticket #1378
jjhelmus May 8, 2012
55deb1e
FIX: map_coordinates, wrap mode, fixes Tickets #796 and #1520
jjhelmus May 8, 2012
3431e5a
new map_coordinate function and boundary unit tests
jjhelmus May 11, 2012
c26ea2c
stash
jjhelmus Jun 4, 2012
0ed6d46
Revert "stash"
jjhelmus Jun 4, 2012
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
99 changes: 94 additions & 5 deletions scipy/ndimage/src/ni_interpolation.c
View file
Open in desktop
Expand Up @@ -110,12 +110,101 @@ spline_coefficients(double x, int order, double *result)
}
}

#define EPSILON 1e-15
/* map a coordinate outside the borders, according to the requested
boundary condition: */
static double
map_coordinate(double in, npy_intp len, int mode)
{
if (in < 0) {
switch (mode) {
case NI_EXTEND_CONSTANT: // constant mode
if (in < 0 - EPSILON) {
in = -1;
break;
}
if (in > len - 1 + EPSILON) {
in = -1;
break;
}
case NI_EXTEND_NEAREST: // nearest mode
if (in < 0 - EPSILON) {
in = 0;
break;
}
if (in > len - 1 + EPSILON) {
in = len - 1;
break;
}
case NI_EXTEND_REFLECT: // reflect mode
if (in < 0 - EPSILON) {
if (len <= 1) {
in = 0;
} else {
npy_intp sz2 = 2 * len;
if (in < - sz2)
in = sz2 * (npy_intp)(-in / sz2) + in;
in = in < -len ? in + sz2 : -in - 1;
}
break;
}
if (in > len - 1 + EPSILON) {
if (len <= 1) {
in = 0;
} else {
npy_intp sz2 = 2 * len;
in -= sz2 * (npy_intp)(in / sz2);
if (in >= len)
in = sz2 - in - 1;
}
break;
}

case NI_EXTEND_MIRROR: // mirror mode
if (in < 0 - EPSILON) {
if (len <= 1) {
in = 0;
} else {
npy_intp sz2 = 2 * len - 2;
in = sz2 * (npy_intp)(-in / sz2) + in;
in = in <= 1 - len ? in + sz2 : -in;
}
break;
}
if (in > len - 1 + EPSILON) {
if (len <= 1) {
in = 0;
} else {
npy_intp sz2 = 2 * len - 2;
in -= sz2 * (npy_intp)(in / sz2);
if (in >= len)
in = sz2 - in;
}
break;
}
case NI_EXTEND_WRAP: // wrap mode
if (in < 0 - EPSILON) {
if (len <= 1) {
in = 0;
} else {
npy_intp sz = len;
in += sz * ((npy_intp)(-in / sz) + 1);
}
break;
}
if (in > len -1 + EPSILON) {
if (len <= 1) {
in = 0;
} else {
npy_intp sz = len;
in -= sz * (npy_intp)(in / sz);
}
break;
}

}
/*
// OLD CODE
if (in < -0.4999) {
switch (mode) {
case NI_EXTEND_MIRROR:
if (len <= 1) {
Expand All @@ -140,7 +229,7 @@ map_coordinate(double in, npy_intp len, int mode)
if (len <= 1) {
in = 0;
} else {
npy_intp sz = len - 1;
npy_intp sz = len;
// Integer division of -in/sz gives (-in mod sz)
// Note that 'in' is negative
in += sz * ((npy_intp)(-in / sz) + 1);
Expand All @@ -153,7 +242,7 @@ map_coordinate(double in, npy_intp len, int mode)
in = -1;
break;
}
} else if (in > len-1) {
} else if (in > len-0.5001) {
switch (mode) {
case NI_EXTEND_MIRROR:
if (len <= 1) {
Expand All @@ -179,7 +268,7 @@ map_coordinate(double in, npy_intp len, int mode)
if (len <= 1) {
in = 0;
} else {
npy_intp sz = len - 1;
npy_intp sz = len;
in -= sz * (npy_intp)(in / sz);
}
break;
Expand All @@ -191,7 +280,7 @@ map_coordinate(double in, npy_intp len, int mode)
break;
}
}

*/
return in;
}

Expand Down
135 changes: 101 additions & 34 deletions scipy/ndimage/tests/test_ndimage.py
View file
Open in desktop
Expand Up @@ -1306,40 +1306,81 @@ def test_extend10(self):
mode=mode, cval=0)
assert_array_equal(output, expected_value)

def test_boundaries(self):
"boundary modes"
def shift(x):
return (x[0] + 0.5,)

data = numpy.array([1,2,3,4.])
expected = {'constant': [1.5,2.5,3.5,-1,-1,-1,-1],
'wrap': [1.5,2.5,3.5,1.5,2.5,3.5,1.5],
'mirror' : [1.5,2.5,3.5,3.5,2.5,1.5,1.5],
'nearest' : [1.5,2.5,3.5,4,4,4,4]}

for mode in expected.keys():
assert_array_equal(expected[mode],
ndimage.geometric_transform(data,shift,
cval=-1,mode=mode,
output_shape=(7,),
order=1))

def test_boundaries2(self):
"boundary modes 2"
def shift(x):
return (x[0] - 0.9,)

data = numpy.array([1,2,3,4])
expected = {'constant': [-1,1,2,3],
'wrap': [3,1,2,3],
'mirror' : [2,1,2,3],
'nearest' : [1,1,2,3]}

for mode in expected.keys():
assert_array_equal(expected[mode],
ndimage.geometric_transform(data,shift,
cval=-1,mode=mode,
output_shape=(4,)))
def test_boundaries_constant(self):
"boundary constant"
input = numpy.array([0, 1, 2, 3.])
expected = {0.1 : [9., 0.9, 1.9, 2.9],
0.4 : [9., 0.6, 1.6, 2.6],
0.6 : [9., 0.4, 1.4, 2.4],
1.0 : [9., 0, 1, 2],
-0.1 : [0.1, 1.1, 2.1, 9.],
-0.4 : [0.4, 1.4, 2.4, 9.],
-0.6 : [0.6, 1.6, 2.6, 9.],
-1.0 : [1., 2., 3., 9.] }
for shift in expected.keys():
assert_almost_equal(expected[shift], ndimage.shift(input, shift,
mode='constant', order=1, cval=9))

def test_boundaries_nearest(self):
"boundary nearest"
input = numpy.array([0, 1, 2, 3.])
expected = {0.1 : [0., 0.9, 1.9, 2.9],
0.4 : [0., 0.6, 1.6, 2.6],
0.6 : [0., 0.4, 1.4, 2.4],
1.0 : [0., 0, 1, 2],
-0.1 : [0.1, 1.1, 2.1, 3.],
-0.4 : [0.4, 1.4, 2.4, 3.],
-0.6 : [0.6, 1.6, 2.6, 3.],
-1.0 : [1., 2., 3., 3.] }
for shift in expected.keys():
assert_almost_equal(expected[shift], ndimage.shift(input, shift,
mode='nearest', order=1))


def test_boundaries_reflect(self):
"boundary reflect"
input = numpy.array([0, 1, 2, 3.])
expected = {0.1 : [0., 0.9, 1.9, 2.9],
0.4 : [0., 0.6, 1.6, 2.6],
0.6 : [0., 0.4, 1.4, 2.4],
1.0 : [0., 0, 1, 2],
-0.1 : [0.1, 1.1, 2.1, 3.],
-0.4 : [0.4, 1.4, 2.4, 3.],
-0.6 : [0.6, 1.6, 2.6, 3.],
-1.0 : [1., 2., 3., 3.] }
for shift in expected.keys():
assert_almost_equal(expected[shift], ndimage.shift(input, shift,
mode='reflect', order=1))

def test_boundaries_mirror(self):
"boundary mirror"
input = numpy.array([0, 1, 2, 3.])
expected = {0.1 : [0.1, 0.9, 1.9, 2.9],
0.4 : [0.4, 0.6, 1.6, 2.6],
0.6 : [0.6, 0.4, 1.4, 2.4],
1.0 : [1., 0, 1, 2],
-0.1 : [0.1, 1.1, 2.1, 2.9],
-0.4 : [0.4, 1.4, 2.4, 2.6],
-0.6 : [0.6, 1.6, 2.6, 2.4],
-1.0 : [1., 2., 3., 2.] }
for shift in expected.keys():
assert_almost_equal(expected[shift], ndimage.shift(input, shift,
mode='mirror', order=1))

def test_boundaries_wrap(self):
"boundary wrap"
input = numpy.array([0, 1, 2, 3.])
expected = {#0.1 : [.3, 0.9, 1.9, 2.9],
#0.4 : [1.2, 0.6, 1.6, 2.6],
#0.6 : [1.8, 0.4, 1.4, 2.4],
1.0 : [3., 0, 1, 2],
#-0.1 : [0.1, 1.1, 2.1, 2.7],
#-0.4 : [0.4, 1.4, 2.4, 1.8],
#-0.6 : [0.6, 1.6, 2.6, 1.2],
-1.0 : [1., 2., 3., 0.] }
for shift in expected.keys():
assert_almost_equal(expected[shift], ndimage.shift(input, shift,
mode='wrap', order=1))

def test_fourier_gaussian_real01(self):
"gaussian fourier filter for real transforms 1"
Expand Down Expand Up @@ -2103,6 +2144,16 @@ def test_shift09(self):
[0, 4, 1, 3],
[0, 7, 6, 8]])

def test_shift10(self):
"Ticket #796"
data = numpy.arange(16).reshape(4, 4)
ref = np.array([[3, 0, 1, 2],
[7, 4, 5, 6],
[11, 8, 9, 10],
[15, 12, 13, 14]])
out = ndimage.shift(data, (0, 1), mode='wrap')
assert_array_almost_equal(out, ref)

def test_zoom1(self):
"zoom 1"
for order in range(0,6):
Expand Down Expand Up @@ -2260,6 +2311,22 @@ def test_rotate08(self):
out = ndimage.rotate(data, 90, axes=(0, 1),
reshape=False)
assert_array_almost_equal(out, expected)

def test_rotate09(self):
"Ticket #1378"
data = numpy.array([[[-1, -1, -1],
[-1, -1, -1],
[-1, -1, -1]],
[[ 0, 0, 0],
[ 0, 1, 0],
[ 0, 0, 0]],
[[ 0, 0, 0],
[ 0, 1, 0],
[ 0, 0, 0]]])
out = ndimage.rotate(data, 180, axes=(0, 1))
expected = ndimage.rotate(ndimage.rotate(data, 90, axes=(0, 1)),
90, axes=(0, 1))
assert_array_almost_equal(out, expected)

def test_watershed_ift01(self):
"watershed_ift 1"
Expand Down