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 -\pi\le\phi\mathrm{(phi)}\le\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,) array_like
|
---|---|
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 : (3,) array_like
shear : array_like
angles : array_like
translate : array_like
perspective : array_like
|
---|---|
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: | scale : (3,) ndarray
shear : (3,) ndarray
angles : (3,) ndarray
translate : (3,) ndarray
perspective : ndarray
|
Raises: | ValueError :
|
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 : ndarray (4, 4) 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
phi : array_like
|
---|---|
Returns: | y : (N,2) array
|
Least squares positive semi-definite tensor estimation
Parameters: | B : (3,3) array_like
|
---|---|
Returns: | npds : (3,3) array
|
References
[R2] | 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. |
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
(1, 3)
Rodrigues formula
Rotation matrix for rotation around axis r for angle theta.
The rotation matrix is given by the Rodrigues formula:
R = Id + sin(theta)*Sn + (1-cos(theta))*Sn^2
with:
0 -nz ny
Sn = nz 0 -nx
-ny nx 0
where n = r / ||r||
In case the angle ||r|| is very small, the above formula may lead to numerical instabilities. We instead use a Taylor expansion around theta=0:
R = I + sin(theta)/tetha Sr + (1-cos(theta))/teta2 Sr^2
leading to:
R = I + (1-theta2/6)*Sr + (1/2-theta2/24)*Sr^2
Parameters: | r : array_like shape (3,), axis theta : float, angle in degrees |
---|---|
Returns: | R : array, shape (3,3), rotation matrix |
Examples
>>> import numpy as np
>>> from dipy.core.geometry import rodrigues_axis_rotation
>>> v=np.array([0,0,1])
>>> u=np.array([1,0,0])
>>> R=rodrigues_axis_rotation(v,40)
>>> ur=np.dot(R,u)
>>> np.round(np.rad2deg(np.arccos(np.dot(ur,u))))
40.0
Convert spherical coordinates to latitude and longitude.
Returns: | lat, lon : ndarray
|
---|
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
rotation matrix from 2 unit vectors
u,v being unit 3d vectors return a 3x3 rotation matrix R than aligns u to v.
In general there are many rotations that will map u to v. If S is any rotation using v as an axis then R.S will also map u to v since (S.R)u = S(Ru) = Sv = v. The rotation R returned by vec2vec_rotmat leaves fixed the perpendicular to the plane spanned by u and v.
The transpose of R will align v to u.
Parameters: | u : array, shape(3,) v : array, shape(3,) |
---|---|
Returns: | R : array, shape(3,3) |
Examples
>>> import numpy as np
>>> from dipy.core.geometry import vec2vec_rotmat
>>> u=np.array([1,0,0])
>>> v=np.array([0,1,0])
>>> R=vec2vec_rotmat(u,v)
>>> np.dot(R,u)
array([ 0., 1., 0.])
>>> np.dot(R.T,v)
array([ 1., 0., 0.])
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 vector Euclidean (L2) norm
See unit vector and Euclidean norm
Parameters: | vec : array_like
axis : int
keepdims : bool
|
---|---|
Returns: | norm : array
|
Examples
>>> import numpy as np
>>> vec = [[8, 15, 0], [0, 36, 77]]
>>> vector_norm(vec)
array([ 17., 85.])
>>> vector_norm(vec, keepdims=True)
array([[ 17.],
[ 85.]])
>>> vector_norm(vec, axis=0)
array([ 8., 39., 77.])