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

BF: Keep up with changes in scipy 0.16 #707

Merged
merged 7 commits into from Sep 14, 2015
Merged

BF: Keep up with changes in scipy 0.16 #707

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f97bdca
BF: Keep up with changes in scipy 0.16
arokem Sep 1, 2015
e9cc7fd
PEP8
arokem Sep 1, 2015
616d485
TST: Test more conditions.
arokem Sep 12, 2015
7a2302a
BF: Don't use Euclidean norm as distance function.
arokem Sep 14, 2015
e8617be
BF + PEP8: use the `warn` function. A bit of whitespace.
arokem Sep 14, 2015
71417bc
Test that a warning is raised when using euclidean_norm,
arokem Sep 14, 2015
d332cfb
PEP8
arokem Sep 14, 2015
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
130 changes: 67 additions & 63 deletions dipy/core/sphere.py
View file
Open in desktop
@@ -1,7 +1,5 @@
from __future__ import division, print_function, absolute_import

__all__ = ['Sphere', 'HemiSphere', 'faces_from_sphere_vertices', 'unique_edges']

import numpy as np
import warnings

Expand All @@ -11,6 +9,9 @@
from dipy.core.onetime import auto_attr
from dipy.reconst.recspeed import remove_similar_vertices

__all__ = ['Sphere', 'HemiSphere', 'faces_from_sphere_vertices',
'unique_edges']


def _all_specified(*args):
for a in args:
Expand Down Expand Up @@ -48,6 +49,7 @@ def faces_from_sphere_vertices(vertices):
else:
return faces


def unique_edges(faces, return_mapping=False):
"""Extract all unique edges from given triangular faces.

Expand Down Expand Up @@ -154,8 +156,8 @@ def __init__(self, x=None, y=None, z=None,

all_specified = _all_specified(x, y, z) + _all_specified(xyz) + \
_all_specified(theta, phi)
one_complete = _some_specified(x, y, z) + _some_specified(xyz) + \
_some_specified(theta, phi)
one_complete = (_some_specified(x, y, z) + _some_specified(xyz) +
_some_specified(theta, phi))

if not (all_specified == 1 and one_complete == 1):
raise ValueError("Sphere must be constructed using either "
Expand Down Expand Up @@ -347,11 +349,11 @@ def mirror(self):
vertices = np.vstack([self.vertices, -self.vertices])

edges = np.vstack([self.edges, n + self.edges])
_switch_vertex(edges[:,0], edges[:,1], vertices)
_switch_vertex(edges[:, 0], edges[:, 1], vertices)

faces = np.vstack([self.faces, n + self.faces])
_switch_vertex(faces[:,0], faces[:,1], vertices)
_switch_vertex(faces[:,0], faces[:,2], vertices)
_switch_vertex(faces[:, 0], faces[:, 1], vertices)
_switch_vertex(faces[:, 0], faces[:, 2], vertices)
return Sphere(xyz=vertices, edges=edges, faces=faces)

@auto_attr
Expand All @@ -370,7 +372,6 @@ def subdivide(self, n=1):
sphere = sphere.subdivide(n)
return HemiSphere.from_sphere(sphere)


def find_closest(self, xyz):
"""
Find the index of the vertex in the Sphere closest to the input vector,
Expand Down Expand Up @@ -429,11 +430,11 @@ def _get_forces(charges):
potential = 1. / r_mag

d = np.arange(len(charges))
force[d,d] = 0
force[d, d] = 0
force = force.sum(0)
force_r_comp = (charges*force).sum(-1)[:, None]
f_theta = force - force_r_comp*charges
potential[d,d] = 0
potential[d, d] = 0
potential = 2*potential.sum()
return f_theta, potential

Expand Down Expand Up @@ -498,8 +499,8 @@ def disperse_charges(hemi, iters, const=.2):


def interp_rbf(data, sphere_origin, sphere_target,
function='multiquadric', epsilon=None, smooth=0,
norm = "euclidean_norm"):
function='multiquadric', epsilon=None, smooth=0.1,
norm="angle"):
"""Interpolate data on the sphere, using radial basis functions.

Parameters
Expand All @@ -514,16 +515,18 @@ def interp_rbf(data, sphere_origin, sphere_target,
function : {'multiquadric', 'inverse', 'gaussian'}
Radial basis function.
epsilon : float
Radial basis function spread parameter.
Radial basis function spread parameter. Defaults to approximate average
distance between nodes.
a good start
smooth : float
values greater than zero increase the smoothness of the
approximation with 0 (the default) as pure interpolation.
approximation with 0 as pure interpolation. Default: 0.1
norm : str
A string indicating the function that returns the
"distance" between two points.
'angle' - The angle between two vectors
'euclidean_norm' - The Euclidean distance

Returns
-------
v : (M,) ndarray
Expand All @@ -535,22 +538,26 @@ def interp_rbf(data, sphere_origin, sphere_target,

"""
from scipy.interpolate import Rbf

def angle(x1, x2):
xx = np.arccos((x1 * x2).sum(axis=0))
xx[np.isnan(xx)] = 0
return xx

def euclidean_norm(x1, x2):
return np.sqrt(((x1 - x2)**2).sum(axis=0))
if norm is "angle":

if norm == "angle":
norm = angle
elif norm is "euclidean_norm":
elif norm == "euclidean_norm":
w_s = "The Eucldian norm used for interpolation is inaccurate "
w_s += "and will be deprecated in future versions. Please consider "
w_s += "using the 'angle' norm instead"
warnings.warn(w_s, DeprecationWarning)
norm = euclidean_norm
# Workaround for bug in SciPy that doesn't allow
# specification of epsilon None

# Workaround for bug in older versions of SciPy that don't allow
# specification of epsilon None:
if epsilon is not None:
kwargs = {'function': function,
'epsilon': epsilon,
Expand Down Expand Up @@ -608,13 +615,12 @@ def euler_characteristic_check(sphere, chi=2):


octahedron_vertices = np.array(
[[ 1.0 , 0.0, 0.0],
[-1.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0],
[ 0.0, -1.0, 0.0],
[ 0.0, 0.0, 1.0],
[ 0.0, 0.0, -1.0],
])
[[1.0, 0.0, 0.0],
[-1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, -1.0, 0.0],
[0.0, 0.0, 1.0],
[0.0, 0.0, -1.0], ])
octahedron_faces = np.array(
[[0, 4, 2],
[1, 5, 3],
Expand All @@ -623,47 +629,45 @@ def euler_characteristic_check(sphere, chi=2):
[1, 4, 3],
[0, 5, 2],
[0, 4, 3],
[1, 5, 2],
], dtype='uint16')
[1, 5, 2], ], dtype='uint16')

t = (1 + np.sqrt(5)) / 2
icosahedron_vertices = np.array(
[[ t, 1, 0], # 0
[ -t, 1, 0], # 1
[ t, -1, 0], # 2
[ -t, -1, 0], # 3
[ 1, 0, t], # 4
[ 1, 0, -t], # 5
[ -1, 0, t], # 6
[ -1, 0, -t], # 7
[ 0, t, 1], # 8
[ 0, -t, 1], # 9
[ 0, t, -1], # 10
[ 0, -t, -1], # 11
])
[[t, 1, 0], # 0
[-t, 1, 0], # 1
[t, -1, 0], # 2
[-t, -1, 0], # 3
[1, 0, t], # 4
[1, 0, -t], # 5
[-1, 0, t], # 6
[-1, 0, -t], # 7
[0, t, 1], # 8
[0, -t, 1], # 9
[0, t, -1], # 10
[0, -t, -1], ]) # 11

icosahedron_vertices /= vector_norm(icosahedron_vertices, keepdims=True)
icosahedron_faces = np.array(
[[ 8, 4, 0],
[ 2, 5, 0],
[ 2, 5, 11],
[ 9, 2, 11],
[ 2, 4, 0],
[ 9, 2, 4],
[[8, 4, 0],
[2, 5, 0],
[2, 5, 11],
[9, 2, 11],
[2, 4, 0],
[9, 2, 4],
[10, 8, 1],
[10, 8, 0],
[10, 5, 0],
[ 6, 3, 1],
[ 9, 6, 3],
[ 6, 8, 1],
[ 6, 8, 4],
[ 9, 6, 4],
[ 7, 10, 1],
[ 7, 10, 5],
[ 7, 3, 1],
[ 7, 3, 11],
[ 9, 3, 11],
[ 7, 5, 11],
], dtype='uint16')
[6, 3, 1],
[9, 6, 3],
[6, 8, 1],
[6, 8, 4],
[9, 6, 4],
[7, 10, 1],
[7, 10, 5],
[7, 3, 1],
[7, 3, 11],
[9, 3, 11],
[7, 5, 11], ], dtype='uint16')

unit_octahedron = Sphere(xyz=octahedron_vertices, faces=octahedron_faces)
unit_icosahedron = Sphere(xyz=icosahedron_vertices, faces=icosahedron_faces)
Expand Down
32 changes: 19 additions & 13 deletions dipy/core/tests/test_sphere.py
View file
Open in desktop
Expand Up @@ -358,21 +358,27 @@ def test_disperse_charges():
nt.assert_array_almost_equal(norms, 1)

def test_interp_rbf():
def data_func(s, a, b):
return a * np.cos(s.theta) + b * np.sin(s.phi)

from dipy.core.sphere import Sphere, interp_rbf
from dipy.core.subdivide_octahedron import create_unit_hemisphere
import numpy as np

s0 = create_unit_hemisphere(2)
s1 = create_unit_hemisphere(3)

data = np.cos(s0.theta) + np.sin(s0.phi)
expected = np.cos(s1.theta) + np.sin(s1.phi)
interp_data_en = interp_rbf(data, s0, s1, norm = "euclidean_norm")
interp_data_a = interp_rbf(data, s0, s1, norm = "angle")

nt.assert_(np.mean(np.abs(interp_data_en - expected)) < 0.1)
nt.assert_(np.mean(np.abs(interp_data_a - expected)) < 0.1)

s0 = create_unit_sphere(3)
s1 = create_unit_sphere(4)
for a, b in zip([1, 2, 0.5], [1, 0.5, 2]):
data = data_func(s0, a, b)
expected = data_func(s1, a, b)
interp_data_a = interp_rbf(data, s0, s1, norm="angle")
nt.assert_(np.mean(np.abs(interp_data_a - expected)) < 0.1)

# Test that using the euclidean norm raises a warning
# (following https://docs.python.org/2/library/warnings.html#testing-warnings)
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always")
interp_data_en = interp_rbf(data, s0, s1, norm ="euclidean_norm")
nt.assert_(len(w) == 1)
nt.assert_(issubclass(w[-1].category, DeprecationWarning))
nt.assert_("deprecated" in str(w[-1].message))

if __name__ == "__main__":
nt.run_module_suite()