Utility functions for algebra etc
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
y : array-like
z : array-like
|
---|---|
Returns : | r : array
theta : array
phi : array
|
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
pts2 : (N,R) or (R,) array-like
|
---|---|
Returns : | d : (N,) or (0,) array
|
See also
Examples
>>> cart_distance([0,0,0], [0,0,3])
3.0
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
|
---|---|
Returns : | circumradius : float
|
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)
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
|
---|---|
Returns : | tuple of: :
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)
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)
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
y : array-like
z : array-like
|
---|---|
Returns : | y : (N,2) array
|
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
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
Returns : | y : (N,2) array
|
---|
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
|
---|---|
Returns : | npds : (3,3) array
|
Examples
>>> B = np.diag([1, 1, -1])
>>> nearest_pos_semi_def(B)
array([[ 0.75, 0. , 0. ],
[ 0. , 0.75, 0. ],
[ 0. , 0. , 0. ]])
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,)
|
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,)
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
theta : array-like
phi : array-like
|
---|---|
Returns : | x : array
y : array
z : array
|
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.
Distance across sphere surface between pts1 and pts2
Parameters : | pts1 : (N,R) or (R,) array-like
pts2 : (N,R) or (R,) array-like
radius : None or float, optional
check_radius : bool, optional
|
---|---|
Returns : | d : (N,) or (0,) array
|
See also
Examples
>>> print '%.4f' % sphere_distance([0,1],[1,0])
1.5708
>>> print '%.4f' % sphere_distance([0,3],[3,0])
4.7124
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
vecs1 : (N, R) or (R,) array-like
|
---|---|
Returns : | vcos : (N,) or (0,) array
|
Notes
The vector cosine will be the same as the correlation only if all the input vectors have zero mean.
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