dipy logo

Site Navigation

NIPY Community

Previous topic

dipy.core.sphere_stats

Next topic

dipy.core.graph

dipy.core.geometry

Utility functions for algebra etc

dipy.core.geometry.cart2sphere(x, y, z)

Return angles for Cartesian 3D coordinates x, y, and z

See doc for sphere2cart for angle conventions and derivation of the formulae.

0 \le \theta \mathrm{(theta)} \le \pi and 0 \le \phi \mathrm{(phi)} \le 2 \pi

Parameters :

x : array-like

x coordinate in Cartesion space

y : array-like

y coordinate in Cartesian space

z : array-like

z coordinate

Returns :

r : array

radius

theta : array

inclination (polar) angle

phi : array

azimuth angle

dipy.core.geometry.cart_distance(pts1, pts2)

Cartesian distance between pts1 and pts2

If either of pts1 or pts2 is 2D, then we take the first dimension to index points, and the second indexes coordinate. More generally, we take the last dimension to be the coordinate dimension.

Parameters :

pts1 : (N,R) or (R,) array-like

where N is the number of points and R is the number of coordinates defining a point (R==3 for 3D)

pts2 : (N,R) or (R,) array-like

where N is the number of points and R is the number of coordinates defining a point (R==3 for 3D). It should be possible to broadcast pts1 against pts2

Returns :

d : (N,) or (0,) array

Cartesian distances between corresponding points in pts1 and pts2

See also

sphere_distance
distance between points on sphere surface

Examples

>>> cart_distance([0,0,0], [0,0,3])
3.0
dipy.core.geometry.circumradius(a, b, c)

a, b and c are 3-dimensional vectors which are the vertices of a triangle. The function returns the circumradius of the triangle, i.e the radius of the smallest circle that can contain the triangle. In the degenerate case when the 3 points are collinear it returns half the distance between the furthest apart points.

Parameters :

a, b, c : (3,) arraylike

the three vertices of the triangle

Returns :

circumradius : float

the desired circumradius

dipy.core.geometry.compose_matrix(scale=None, shear=None, angles=None, translate=None, perspective=None)

Return 4x4 transformation matrix from sequence of transformations.

Code modified from the work of Christoph Gohlke link provided here http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

This is the inverse of the decompose_matrix function.

Parameters :

scale : vector of 3 scaling factors

shear : list of shear factors for x-y, x-z, y-z axes

angles : list of Euler angles about static x, y, z axes

translate : translation vector along x, y, z axes

perspective : perspective partition of matrix

Returns :

matrix : 4x4 array

Examples

>>> import math
>>> import numpy as np
>>> import dipy.core.geometry as gm
>>> scale = np.random.random(3) - 0.5
>>> shear = np.random.random(3) - 0.5
>>> angles = (np.random.random(3) - 0.5) * (2*math.pi)
>>> trans = np.random.random(3) - 0.5
>>> persp = np.random.random(4) - 0.5
>>> M0 = gm.compose_matrix(scale, shear, angles, trans, persp)
dipy.core.geometry.decompose_matrix(matrix)

Return sequence of transformations from transformation matrix.

Code modified from the excellent work of Christoph Gohlke link provided here http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

Parameters :

matrix : array_like

Non-degenerative homogeneous transformation matrix

Returns :

tuple of: :

scale : vector of 3 scaling factors shear : list of shear factors for x-y, x-z, y-z axes angles : list of Euler angles about static x, y, z axes translate : translation vector along x, y, z axes perspective : perspective partition of matrix

Raise ValueError if matrix is of wrong type or degenerative. :

Examples

>>> import numpy as np
>>> T0=np.diag([2,1,1,1])
>>> scale, shear, angles, trans, persp = decompose_matrix(T0)
dipy.core.geometry.euler_matrix(ai, aj, ak, axes='sxyz')

Return homogeneous rotation matrix from Euler angles and axis sequence.

Code modified from the work of Christoph Gohlke link provided here http://www.lfd.uci.edu/~gohlke/code/transformations.py.html

Parameters :

ai, aj, ak : Euler’s roll, pitch and yaw angles

axes : One of 24 axis sequences as string or encoded tuple

Returns :

matrix: 4x4 numpy array :

Code modified from the work of Christoph Gohlke link provided here :

http://www.lfd.uci.edu/~gohlke/code/transformations.py.html :

Examples

>>> import numpy
>>> R = euler_matrix(1, 2, 3, 'syxz')
>>> numpy.allclose(numpy.sum(R[0]), -1.34786452)
True
>>> R = euler_matrix(1, 2, 3, (0, 1, 0, 1))
>>> numpy.allclose(numpy.sum(R[0]), -0.383436184)
True
>>> ai, aj, ak = (4.0*math.pi) * (numpy.random.random(3) - 0.5)
>>> for axes in _AXES2TUPLE.keys():
...    R = euler_matrix(ai, aj, ak, axes)
>>> for axes in _TUPLE2AXES.keys():
...    R = euler_matrix(ai, aj, ak, axes)
dipy.core.geometry.lambert_equal_area_projection_cart(x, y, z)

Lambert Equal Area Projection from cartesian vector to plane

Return positions in (y_1,y_2) plane corresponding to the directions of the vectors with cartesian coordinates xyz under the Lambert Equal Area Projection mapping (see Mardia and Jupp (2000), Directional Statistics, p. 161).

The Lambert EAP maps the upper hemisphere to the planar disc of radius 1 and the lower hemisphere to the planar annulus between radii 1 and 2, The Lambert EAP maps the upper hemisphere to the planar disc of radius 1 and the lower hemisphere to the planar annulus between radii 1 and 2. and vice versa. See doc for sphere2cart for angle conventions

Parameters :

x : array-like

x coordinate in Cartesion space

y : array-like

y coordinate in Cartesian space

z : array-like

z coordinate

Returns :

y : (N,2) array

planar coordinates of points following mapping by Lambert’s EAP.

dipy.core.geometry.lambert_equal_area_projection_polar(theta, phi)

Lambert Equal Area Projection from polar sphere to plane

Return positions in (y1,y2) plane corresponding to the points with polar coordinates (theta, phi) on the unit sphere, under the Lambert Equal Area Projection mapping (see Mardia and Jupp (2000), Directional Statistics, p. 161).

See doc for sphere2cart for angle conventions

  • 0 \le \theta \le \pi and 0 \le \phi \le 2 \pi
  • |(y_1,y_2)| \le 2

The Lambert EAP maps the upper hemisphere to the planar disc of radius 1 and the lower hemisphere to the planar annulus between radii 1 and 2, and vice versa. Parameters ———- theta : array-like

theta spherical coordinates
phi : array-like
phi spherical coordinates
Returns :

y : (N,2) array

planar coordinates of points following mapping by Lambert’s EAP.

dipy.core.geometry.nearest_pos_semi_def(B)

Least squares positive semi-definite tensor estimation

Reference: Niethammer M, San Jose Estepar R, Bouix S, Shenton M, Westin CF. On diffusion tensor estimation. Conf Proc IEEE Eng Med Biol Soc. 2006;1:2622-5. PubMed PMID: 17946125; PubMed Central PMCID: PMC2791793.

Parameters :

B : (3,3) array-like

B matrix - symmetric. We do not check the symmetry.

Returns :

npds : (3,3) array

Estimated nearest positive semi-definite array to matrix B.

Examples

>>> B = np.diag([1, 1, -1])
>>> nearest_pos_semi_def(B)
array([[ 0.75,  0.  ,  0.  ],
       [ 0.  ,  0.75,  0.  ],
       [ 0.  ,  0.  ,  0.  ]])
dipy.core.geometry.normalized_vector(vec)

Return vector divided by its Euclidean (L2) norm

See unit vector and Euclidean norm

Parameters :

vec : array-like shape (3,)

Returns :

nvec : array shape (3,)

vector divided by L2 norm

Examples

>>> vec = [1, 2, 3]
>>> l2n = np.sqrt(np.dot(vec, vec))
>>> nvec = normalized_vector(vec)
>>> np.allclose(np.array(vec) / l2n, nvec)
True
>>> vec = np.array([[1, 2, 3]])
>>> vec.shape
(1, 3)
>>> normalized_vector(vec).shape
(3,)
dipy.core.geometry.sphere2cart(r, theta, phi)

Spherical to Cartesian coordinates

This is the standard physics convention where theta is the inclination (polar) angle, and phi is the azimuth angle.

Imagine a sphere with center (0,0,0). Orient it with the z axis running south->north, the y axis running west-east and the x axis from posterior to anterior. theta (the inclination angle) is the angle to rotate from the z-axis (the zenith) around the y-axis, towards the x axis. Thus the rotation is counter-clockwise from the point of view of positive y. phi (azimuth) gives the angle of rotation around the z-axis towards the y axis. The rotation is counter-clockwise from the point of view of positive z.

Equivalently, given a point P on the sphere, with coordinates x, y, z, theta is the angle between P and the z-axis, and phi is the angle between the projection of P onto the XY plane, and the X axis.

Geographical nomenclature designates theta as ‘co-latitude’, and phi as ‘longitude’

Parameters :

r : array-like

radius

theta : array-like

inclination or polar angle

phi : array-like

azimuth angle

Returns :

x : array

x coordinate(s) in Cartesion space

y : array

y coordinate(s) in Cartesian space

z : array

z coordinate

Notes

See these pages:

for excellent discussion of the many different conventions possible. Here we use the physics conventions, used in the wikipedia page.

Derivations of the formulae are simple. Consider a vector x, y, z of length r (norm of x, y, z). The inclination angle (theta) can be found from: cos(theta) == z / r -> z == r * cos(theta). This gives the hypotenuse of the projection onto the XY plane, which we will call Q. Q == r*sin(theta). Now x / Q == cos(phi) -> x == r * sin(theta) * cos(phi) and so on.

We have deliberately named this function sphere2cart rather than sph2cart to distinguish it from the Matlab function of that name, because the Matlab function uses an unusual convention for the angles that we did not want to replicate. The Matlab function is trivial to implement with the formulae given in the Matlab help.

dipy.core.geometry.sphere_distance(pts1, pts2, radius=None, check_radius=True)

Distance across sphere surface between pts1 and pts2

Parameters :

pts1 : (N,R) or (R,) array-like

where N is the number of points and R is the number of coordinates defining a point (R==3 for 3D)

pts2 : (N,R) or (R,) array-like

where N is the number of points and R is the number of coordinates defining a point (R==3 for 3D). It should be possible to broadcast pts1 against pts2

radius : None or float, optional

Radius of sphere. Default is to work out radius from mean of the length of each point vector

check_radius : bool, optional

If True, check if the points are on the sphere surface - i.e check if the vector lengths in pts1 and pts2 are close to radius. Default is True.

Returns :

d : (N,) or (0,) array

Distances between corresponding points in pts1 and pts2 across the spherical surface, i.e. the great circle distance

See also

cart_distance
cartesian distance between points
vector_cosine
cosine of angle between vectors

Examples

>>> print '%.4f' % sphere_distance([0,1],[1,0])
1.5708
>>> print '%.4f' % sphere_distance([0,3],[3,0])
4.7124
dipy.core.geometry.vector_cosine(vecs1, vecs2)

Cosine of angle between two (sets of) vectors

The cosine of the angle between two vectors v1 and v2 is given by the inner product of v1 and v2 divided by the product of the vector lengths:

v_cos = np.inner(v1, v2) / (np.sqrt(np.sum(v1**2)) *
                            np.sqrt(np.sum(v2**2)))
Parameters :

vecs1 : (N, R) or (R,) array-like

N vectors (as rows) or single vector. Vectors have R elements.

vecs1 : (N, R) or (R,) array-like

N vectors (as rows) or single vector. Vectors have R elements. It should be possible to broadcast vecs1 against vecs2

Returns :

vcos : (N,) or (0,) array

Vector cosines. To get the angles you will need np.arccos

Notes

The vector cosine will be the same as the correlation only if all the input vectors have zero mean.

dipy.core.geometry.vector_norm(data, axis=None, out=None)

Return length, i.e. euclidean norm, of ndarray along axis.

Examples

>>> import numpy
>>> v = numpy.random.random(3)
>>> n = vector_norm(v)
>>> numpy.allclose(n, numpy.linalg.norm(v))
True
>>> v = numpy.random.rand(6, 5, 3)
>>> n = vector_norm(v, axis=-1)
>>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=2)))
True
>>> n = vector_norm(v, axis=1)
>>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1)))
True
>>> v = numpy.random.rand(5, 4, 3)
>>> n = numpy.empty((5, 3), dtype=numpy.float64)
>>> vector_norm(v, axis=1, out=n)
>>> numpy.allclose(n, numpy.sqrt(numpy.sum(v*v, axis=1)))
True
>>> vector_norm([])
0.0
>>> vector_norm([1.0])
1.0