223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'vectorize_rotation' into development

This commit is contained in:
Karo Sedighiani 2020-04-14 19:34:11 +02:00
parent d3eba41519 3d10266fbc
commit 03aa6f9c80
6 changed files with 617 additions and 222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit aa0b9fe992ce0bfc172989ab006893ddb8939c1e Subproject commit 232a094c715bcbbd1c6652c4dc4a4a50d402b82f

34
python/damask/_Lambert.py
View File

@ -51,19 +51,18 @@ def cube_to_ball(cube):
https://doi.org/10.1088/0965-0393/22/7/075013 https://doi.org/10.1088/0965-0393/22/7/075013
""" """
if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5: cube_ = np.clip(cube,None,np.pi**(2./3.) * 0.5) if np.isclose(np.abs(np.max(cube)),np.pi**(2./3.) * 0.5,atol=1e-6) else cube
raise ValueError
# transform to the sphere grid via the curved square, and intercept the zero point # transform to the sphere grid via the curved square, and intercept the zero point
if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(cube_,0.0,rtol=0.0,atol=1.0e-16):
ball = np.zeros(3) ball = np.zeros(3)
else: else:
# get pyramide and scale by grid parameter ratio # get pyramide and scale by grid parameter ratio
p = _get_order(cube) p = _get_order(cube_)
XYZ = cube[p] * sc XYZ = cube_[p[0]] * sc
# intercept all the points along the z-axis # intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300): if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16):
ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else: else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1]
@ -82,7 +81,7 @@ def cube_to_ball(cube):
ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ])
# reverse the coordinates back to the regular order according to the original pyramid number # reverse the coordinates back to the regular order according to the original pyramid number
ball = ball[p] ball = ball[p[1]]
return ball return ball
@ -102,15 +101,14 @@ def ball_to_cube(ball):
https://doi.org/10.1088/0965-0393/22/7/075013 https://doi.org/10.1088/0965-0393/22/7/075013
""" """
rs = np.linalg.norm(ball) ball_ = ball/np.linalg.norm(ball)*R1 if np.isclose(np.linalg.norm(ball),R1,atol=1e-6) else ball
if rs > R1: rs = np.linalg.norm(ball_)
raise ValueError
if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(ball_,0.0,rtol=0.0,atol=1.0e-16):
cube = np.zeros(3) cube = np.zeros(3)
else: else:
p = _get_order(ball) p = _get_order(ball_)
xyz3 = ball[p] xyz3 = ball_[p[0]]
# inverse M_3 # inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
@ -118,7 +116,7 @@ def ball_to_cube(ball):
# inverse M_2 # inverse M_2
qxy = np.sum(xyz2**2) qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300): if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16):
Tinv = np.zeros(2) Tinv = np.zeros(2)
else: else:
q2 = qxy + np.max(np.abs(xyz2))**2 q2 = qxy + np.max(np.abs(xyz2))**2
@ -132,7 +130,7 @@ def ball_to_cube(ball):
# inverse M_1 # inverse M_1
cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc
# reverse the coordinates back to the regular order according to the original pyramid number # reverse the coordinates back to the regular order according to the original pyramid number
cube = cube[p] cube = cube[p[1]]
return cube return cube
@ -157,10 +155,10 @@ def _get_order(xyz):
""" """
if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \ if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
(abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]): (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
return [0,1,2] return [[0,1,2],[0,1,2]]
elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \ elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \
(abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]): (abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]):
return [1,2,0] return [[1,2,0],[2,0,1]]
elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \ elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \
(abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]): (abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]):
return [2,0,1] return [[2,0,1],[1,2,0]]

362
python/damask/_rotation.py
View File

@ -2,7 +2,7 @@ import numpy as np
from ._Lambert import ball_to_cube, cube_to_ball from ._Lambert import ball_to_cube, cube_to_ball
P = -1 _P = -1
def iszero(a): def iszero(a):
return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
@ -89,7 +89,7 @@ class Rotation:
other_q = other.quaternion[0] other_q = other.quaternion[0]
other_p = other.quaternion[1:] other_p = other.quaternion[1:]
R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p), R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))) self_q*other_p + other_q*self_p + _P * np.cross(self_p,other_p)))
return R.standardize() return R.standardize()
elif isinstance(other, (tuple,np.ndarray)): elif isinstance(other, (tuple,np.ndarray)):
if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid
@ -97,7 +97,7 @@ class Rotation:
B = 2.0 * ( self.quaternion[1]*other[0] B = 2.0 * ( self.quaternion[1]*other[0]
+ self.quaternion[2]*other[1] + self.quaternion[2]*other[1]
+ self.quaternion[3]*other[2]) + self.quaternion[3]*other[2])
C = 2.0 * P*self.quaternion[0] C = 2.0 * _P*self.quaternion[0]
return np.array([ return np.array([
A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]), A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
@ -262,9 +262,9 @@ class Rotation:
if acceptHomomorph: if acceptHomomorph:
qu *= -1. qu *= -1.
else: else:
raise ValueError('Quaternion has negative first component.\n{}'.format(qu[0])) raise ValueError('Quaternion has negative first component: {}.'.format(qu[0]))
if not np.isclose(np.linalg.norm(qu), 1.0): if not np.isclose(np.linalg.norm(qu), 1.0):
raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu)) raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu))
return Rotation(qu) return Rotation(qu)
@ -276,7 +276,7 @@ class Rotation:
else np.array(eulers,dtype=float) else np.array(eulers,dtype=float)
eu = np.radians(eu) if degrees else eu eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi: if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu)) raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu))
return Rotation(Rotation.eu2qu(eu)) return Rotation(Rotation.eu2qu(eu))
@ -292,9 +292,9 @@ class Rotation:
if degrees: ax[ 3] = np.radians(ax[3]) if degrees: ax[ 3] = np.radians(ax[3])
if normalise: ax[0:3] /= np.linalg.norm(ax[0:3]) if normalise: ax[0:3] /= np.linalg.norm(ax[0:3])
if ax[3] < 0.0 or ax[3] > np.pi: if ax[3] < 0.0 or ax[3] > np.pi:
raise ValueError('Axis angle rotation angle outside of [0..π].\n{}'.format(ax[3])) raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3]))
if not np.isclose(np.linalg.norm(ax[0:3]), 1.0): if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[0:3])) raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3]))
return Rotation(Rotation.ax2qu(ax)) return Rotation(Rotation.ax2qu(ax))
@ -312,11 +312,11 @@ class Rotation:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition (U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.dot(U,Vh) om = np.dot(U,Vh)
if not np.isclose(np.linalg.det(om),1.0): if not np.isclose(np.linalg.det(om),1.0):
raise ValueError('matrix is not a proper rotation.\n{}'.format(om)) raise ValueError('matrix is not a proper rotation: {}.'.format(om))
if not np.isclose(np.dot(om[0],om[1]), 0.0) \ if not np.isclose(np.dot(om[0],om[1]), 0.0) \
or not np.isclose(np.dot(om[1],om[2]), 0.0) \ or not np.isclose(np.dot(om[1],om[2]), 0.0) \
or not np.isclose(np.dot(om[2],om[0]), 0.0): or not np.isclose(np.dot(om[2],om[0]), 0.0):
raise ValueError('matrix is not orthogonal.\n{}'.format(om)) raise ValueError('matrix is not orthogonal: {}.'.format(om))
return Rotation(Rotation.om2qu(om)) return Rotation(Rotation.om2qu(om))
@ -336,9 +336,9 @@ class Rotation:
if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1 if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1
if normalise: ro[0:3] /= np.linalg.norm(ro[0:3]) if normalise: ro[0:3] /= np.linalg.norm(ro[0:3])
if not np.isclose(np.linalg.norm(ro[0:3]), 1.0): if not np.isclose(np.linalg.norm(ro[0:3]), 1.0):
raise ValueError('Rodrigues rotation axis is not of unit length.\n{} {} {}'.format(*ro[0:3])) raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3]))
if ro[3] < 0.0: if ro[3] < 0.0:
raise ValueError('Rodrigues rotation angle not positive.\n{}'.format(ro[3])) raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3]))
return Rotation(Rotation.ro2qu(ro)) return Rotation(Rotation.ro2qu(ro))
@ -350,6 +350,9 @@ class Rotation:
else np.array(homochoric,dtype=float) else np.array(homochoric,dtype=float)
if P > 0: ho *= -1 # convert from P=1 to P=-1 if P > 0: ho *= -1 # convert from P=1 to P=-1
if np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9:
raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho))
return Rotation(Rotation.ho2qu(ho)) return Rotation(Rotation.ho2qu(ho))
@staticmethod @staticmethod
@ -358,6 +361,10 @@ class Rotation:
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \ cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
else np.array(cubochoric,dtype=float) else np.array(cubochoric,dtype=float)
if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9:
raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu))
ho = Rotation.cu2ho(cu) ho = Rotation.cu2ho(cu)
if P > 0: ho *= -1 # convert from P=1 to P=-1 if P > 0: ho *= -1 # convert from P=1 to P=-1
@ -383,7 +390,7 @@ class Rotation:
""" """
if not all(isinstance(item, Rotation) for item in rotations): if not all(isinstance(item, Rotation) for item in rotations):
raise TypeError("Only instances of Rotation can be averaged.") raise TypeError('Only instances of Rotation can be averaged.')
N = len(rotations) N = len(rotations)
if not weights: if not weights:
@ -441,34 +448,69 @@ class Rotation:
#---------- Quaternion ---------- #---------- Quaternion ----------
@staticmethod @staticmethod
def qu2om(qu): def qu2om(qu):
if len(qu.shape) == 1:
"""Quaternion to rotation matrix.""" """Quaternion to rotation matrix."""
qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2)
om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2)
om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) om[0,1] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3])
om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) om[1,0] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3])
om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) om[1,2] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1])
om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1])
om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2])
om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2])
return om if P > 0.0 else om.T else:
qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2)
om = np.block([qq + 2.0*qu[...,1:2]**2,
2.0*(qu[...,2:3]*qu[...,1:2]+qu[...,0:1]*qu[...,3:4]),
2.0*(qu[...,3:4]*qu[...,1:2]-qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,1:2]*qu[...,2:3]-qu[...,0:1]*qu[...,3:4]),
qq + 2.0*qu[...,2:3]**2,
2.0*(qu[...,3:4]*qu[...,2:3]+qu[...,0:1]*qu[...,1:2]),
2.0*(qu[...,1:2]*qu[...,3:4]+qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,2:3]*qu[...,3:4]-qu[...,0:1]*qu[...,1:2]),
qq + 2.0*qu[...,3:4]**2,
]).reshape(qu.shape[:-1]+(3,3))
return om if _P < 0.0 else np.swapaxes(om,(-1,-2))
@staticmethod @staticmethod
def qu2eu(qu): def qu2eu(qu):
"""Quaternion to Bunge-Euler angles.""" """Quaternion to Bunge-Euler angles."""
if len(qu.shape) == 1:
q03 = qu[0]**2+qu[3]**2 q03 = qu[0]**2+qu[3]**2
q12 = qu[1]**2+qu[2]**2 q12 = qu[1]**2+qu[2]**2
chi = np.sqrt(q03*q12) chi = np.sqrt(q03*q12)
if np.abs(q12) < 1.e-8:
if iszero(chi): eu = np.array([np.arctan2(-_P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0])
eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \ elif np.abs(q03) < 1.e-8:
np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0])
else: else:
eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), eu = np.array([np.arctan2((-_P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]-qu[2]*qu[3])*chi ),
np.arctan2( 2.0*chi, q03-q12 ), np.arctan2( 2.0*chi, q03-q12 ),
np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) np.arctan2(( _P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]+qu[2]*qu[3])*chi )])
else:
q02 = qu[...,0:1]*qu[...,2:3]
q13 = qu[...,1:2]*qu[...,3:4]
q01 = qu[...,0:1]*qu[...,1:2]
q23 = qu[...,2:3]*qu[...,3:4]
q03_s = qu[...,0:1]**2+qu[...,3:4]**2
q12_s = qu[...,1:2]**2+qu[...,2:3]**2
chi = np.sqrt(q03_s*q12_s)
# reduce Euler angles to definition range, i.e a lower limit of 0.0 eu = np.where(np.abs(q12_s) < 1.0e-8,
np.block([np.arctan2(-_P*2.0*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2),
np.zeros(qu.shape[:-1]+(2,))]),
np.where(np.abs(q03_s) < 1.0e-8,
np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
np.broadcast_to(np.pi,qu.shape[:-1]+(1,)),
np.zeros(qu.shape[:-1]+(1,))]),
np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi),
np.arctan2( 2.0*chi, q03_s-q12_s ),
np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)])
)
)
# reduce Euler angles to definition range
eu[np.abs(eu)<1.e-6] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu return eu
@ -479,38 +521,65 @@ class Rotation:
Modified version of the original formulation, should be numerically more stable Modified version of the original formulation, should be numerically more stable
""" """
if iszero(qu[1]**2+qu[2]**2+qu[3]**2): # set axis to [001] if the angle is 0/360 if len(qu.shape) == 1:
ax = [ 0.0, 0.0, 1.0, 0.0 ] if np.abs(np.sum(qu[1:4]**2)) < 1.e-6: # set axis to [001] if the angle is 0/360
elif not iszero(qu[0]): ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
elif qu[0] > 1.e-6:
s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2)
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ])
else: else:
ax = [ qu[1], qu[2], qu[3], np.pi] ax = ax = np.array([ qu[1], qu[2], qu[3], np.pi])
return np.array(ax) else:
with np.errstate(invalid='ignore',divide='ignore'):
s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2)
omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-6,qu.shape),
np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]),
np.block([qu[...,1:4]*s,omega]))
ax[np.sum(np.abs(qu[...,1:4])**2,axis=-1) < 1.0e-6,] = [0.0, 0.0, 1.0, 0.0]
return ax
@staticmethod @staticmethod
def qu2ro(qu): def qu2ro(qu):
"""Quaternion to Rodrigues-Frank vector.""" """Quaternion to Rodrigues-Frank vector."""
if len(qu.shape) == 1:
if iszero(qu[0]): if iszero(qu[0]):
ro = [qu[1], qu[2], qu[3], np.inf] ro = np.array([qu[1], qu[2], qu[3], np.inf])
else: else:
s = np.linalg.norm([qu[1],qu[2],qu[3]]) s = np.linalg.norm(qu[1:4])
ro = [0.0,0.0,P,0.0] if iszero(s) else \ ro = np.array([0.0,0.0,_P,0.0] if iszero(s) else \
[ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))])
return np.array(ro) else:
with np.errstate(invalid='ignore',divide='ignore'):
s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True)
ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape),
np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]),
np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s,
np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0)))
])
)
ro[np.abs(s).squeeze(-1) < 1.0e-12] = [0.0,0.0,_P,0.0]
return ro
@staticmethod @staticmethod
def qu2ho(qu): def qu2ho(qu):
"""Quaternion to homochoric vector.""" """Quaternion to homochoric vector."""
if len(qu.shape) == 1:
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
if np.abs(omega) < 1.0e-12:
if iszero(omega): ho = np.zeros(3)
ho = np.array([ 0.0, 0.0, 0.0 ])
else: else:
ho = np.array([qu[1], qu[2], qu[3]]) ho = np.array([qu[1], qu[2], qu[3]])
f = 0.75 * ( omega - np.sin(omega) ) f = 0.75 * ( omega - np.sin(omega) )
ho = ho/np.linalg.norm(ho) * f**(1./3.) ho = ho/np.linalg.norm(ho) * f**(1./3.)
else:
with np.errstate(invalid='ignore'):
omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
ho = np.where(np.abs(omega) < 1.0e-12,
np.zeros(3),
qu[...,1:4]/np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) \
* (0.75*(omega - np.sin(omega)))**(1./3.))
return ho return ho
@staticmethod @staticmethod
@ -532,37 +601,71 @@ class Rotation:
@staticmethod @staticmethod
def om2eu(om): def om2eu(om):
"""Rotation matrix to Bunge-Euler angles.""" """Rotation matrix to Bunge-Euler angles."""
if abs(om[2,2]) < 1.0: if len(om.shape) == 2:
if not np.isclose(np.abs(om[2,2]),1.0,1.e-4):
zeta = 1.0/np.sqrt(1.0-om[2,2]**2) zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
np.arccos(om[2,2]), np.arccos(om[2,2]),
np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) np.arctan2(om[0,2]*zeta, om[1,2]*zeta)])
else: else:
eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation
else:
# reduce Euler angles to definition range, i.e a lower limit of 0.0 with np.errstate(invalid='ignore',divide='ignore'):
zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2)
eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4),
np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]),
np.pi*0.5*(1-om[...,2,2:3]),
np.zeros(om.shape[:-2]+(1,)),
]),
np.block([np.arctan2(om[...,2,0:1]*zeta,-om[...,2,1:2]*zeta),
np.arccos(om[...,2,2:3]),
np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta)
])
)
eu[np.abs(eu)<1.e-6] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu return eu
@staticmethod @staticmethod
def om2ax(om): def om2ax(om):
"""Rotation matrix to axis angle pair.""" """Rotation matrix to axis angle pair."""
if len(om.shape) == 2:
ax=np.empty(4) ax=np.empty(4)
# first get the rotation angle # first get the rotation angle
t = 0.5*(om.trace() -1.0) t = 0.5*(om.trace() -1.0)
ax[3] = np.arccos(np.clip(t,-1.0,1.0)) ax[3] = np.arccos(np.clip(t,-1.0,1.0))
if np.abs(ax[3])<1.e-6:
if iszero(ax[3]): ax = np.array([ 0.0, 0.0, 1.0, 0.0])
ax = [ 0.0, 0.0, 1.0, 0.0]
else: else:
w,vr = np.linalg.eig(om) w,vr = np.linalg.eig(om)
# next, find the eigenvalue (1,0j) # next, find the eigenvalue (1,0j)
i = np.where(np.isclose(w,1.0+0.0j))[0][0] i = np.where(np.isclose(w,1.0+0.0j))[0][0]
ax[0:3] = np.real(vr[0:3,i]) ax[0:3] = np.real(vr[0:3,i])
diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) diagDelta[np.abs(diagDelta)<1.e-6] = 1.0
return np.array(ax) ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta))
else:
diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
om[...,2,0:1]-om[...,0,2:3],
om[...,0,1:2]-om[...,1,0:1]
])
diag_delta[np.abs(diag_delta)<1.e-6] = 1.0
t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
w,vr = np.linalg.eig(om)
# mask duplicated real eigenvalues
w[np.isclose(w[...,0],1.0+0.0j),1:] = 0.
w[np.isclose(w[...,1],1.0+0.0j),2:] = 0.
vr = np.swapaxes(vr,-1,-2)
ax = np.where(np.abs(diag_delta)<0,
np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \
*np.sign(diag_delta))
ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))])
ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0]
return ax
@staticmethod @staticmethod
def om2ro(om): def om2ro(om):
@ -584,32 +687,57 @@ class Rotation:
@staticmethod @staticmethod
def eu2qu(eu): def eu2qu(eu):
"""Bunge-Euler angles to quaternion.""" """Bunge-Euler angles to quaternion."""
if len(eu.shape) == 1:
ee = 0.5*eu ee = 0.5*eu
cPhi = np.cos(ee[1]) cPhi = np.cos(ee[1])
sPhi = np.sin(ee[1]) sPhi = np.sin(ee[1])
qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), qu = np.array([ cPhi*np.cos(ee[0]+ee[2]),
-P*sPhi*np.cos(ee[0]-ee[2]), -_P*sPhi*np.cos(ee[0]-ee[2]),
-P*sPhi*np.sin(ee[0]-ee[2]), -_P*sPhi*np.sin(ee[0]-ee[2]),
-P*cPhi*np.sin(ee[0]+ee[2]) ]) -_P*cPhi*np.sin(ee[0]+ee[2]) ])
if qu[0] < 0.0: qu*=-1 if qu[0] < 0.0: qu*=-1
else:
ee = 0.5*eu
cPhi = np.cos(ee[...,1:2])
sPhi = np.sin(ee[...,1:2])
qu = np.block([ cPhi*np.cos(ee[...,0:1]+ee[...,2:3]),
-_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]),
-_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]),
-_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])])
qu[qu[...,0]<0.0]*=-1
return qu return qu
@staticmethod @staticmethod
def eu2om(eu): def eu2om(eu):
"""Bunge-Euler angles to rotation matrix.""" """Bunge-Euler angles to rotation matrix."""
if len(eu.shape) == 1:
c = np.cos(eu) c = np.cos(eu)
s = np.sin(eu) s = np.sin(eu)
om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]],
[-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]],
[+s[0]*s[1], -c[0]*s[1], +c[1] ]]) [+s[0]*s[1], -c[0]*s[1], +c[1] ]])
else:
om[np.where(iszero(om))] = 0.0 c = np.cos(eu)
s = np.sin(eu)
om = np.block([+c[...,0:1]*c[...,2:3]-s[...,0:1]*s[...,2:3]*c[...,1:2],
+s[...,0:1]*c[...,2:3]+c[...,0:1]*s[...,2:3]*c[...,1:2],
+s[...,2:3]*s[...,1:2],
-c[...,0:1]*s[...,2:3]-s[...,0:1]*c[...,2:3]*c[...,1:2],
-s[...,0:1]*s[...,2:3]+c[...,0:1]*c[...,2:3]*c[...,1:2],
+c[...,2:3]*s[...,1:2],
+s[...,0:1]*s[...,1:2],
-c[...,0:1]*s[...,1:2],
+c[...,1:2]
]).reshape(eu.shape[:-1]+(3,3))
om[np.abs(om)<1.e-12] = 0.0
return om return om
@staticmethod @staticmethod
def eu2ax(eu): def eu2ax(eu):
"""Bunge-Euler angles to axis angle pair.""" """Bunge-Euler angles to axis angle pair."""
if len(eu.shape) == 1:
t = np.tan(eu[1]*0.5) t = np.tan(eu[1]*0.5)
sigma = 0.5*(eu[0]+eu[2]) sigma = 0.5*(eu[0]+eu[2])
delta = 0.5*(eu[0]-eu[2]) delta = 0.5*(eu[0]-eu[2])
@ -617,24 +745,45 @@ class Rotation:
alpha = np.pi if iszero(np.cos(sigma)) else \ alpha = np.pi if iszero(np.cos(sigma)) else \
2.0*np.arctan(tau/np.cos(sigma)) 2.0*np.arctan(tau/np.cos(sigma))
if iszero(alpha): if np.abs(alpha)<1.e-6:
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else: else:
ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front ax = -_P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front
ax = np.append(ax,alpha) ax = np.append(ax,alpha)
if alpha < 0.0: ax *= -1.0 # ensure alpha is positive if alpha < 0.0: ax *= -1.0 # ensure alpha is positive
else:
t = np.tan(eu[...,1:2]*0.5)
sigma = 0.5*(eu[...,0:1]+eu[...,2:3])
delta = 0.5*(eu[...,0:1]-eu[...,2:3])
tau = np.linalg.norm(np.block([t,np.sin(sigma)]),axis=-1,keepdims=True)
alpha = np.where(np.abs(np.cos(sigma))<1.e-12,np.pi,2.0*np.arctan(tau/np.cos(sigma)))
with np.errstate(invalid='ignore',divide='ignore'):
ax = np.where(np.broadcast_to(np.abs(alpha)<1.0e-12,eu.shape[:-1]+(4,)),
[0.0,0.0,1.0,0.0],
np.block([-_P/tau*t*np.cos(delta),
-_P/tau*t*np.sin(delta),
-_P/tau* np.sin(sigma),
alpha
]))
ax[(alpha<0.0).squeeze()] *=-1
return ax return ax
@staticmethod @staticmethod
def eu2ro(eu): def eu2ro(eu):
"""Bunge-Euler angles to Rodrigues-Frank vector.""" """Bunge-Euler angles to Rodrigues-Frank vector."""
if len(eu.shape) == 1:
ro = Rotation.eu2ax(eu) # convert to axis angle pair representation ro = Rotation.eu2ax(eu) # convert to axis angle pair representation
if ro[3] >= np.pi: # Differs from original implementation. check convention 5 if ro[3] >= np.pi: # Differs from original implementation. check convention 5
ro[3] = np.inf ro[3] = np.inf
elif iszero(ro[3]): elif iszero(ro[3]):
ro = np.array([ 0.0, 0.0, P, 0.0 ]) ro = np.array([ 0.0, 0.0, _P, 0.0 ])
else: else:
ro[3] = np.tan(ro[3]*0.5) ro[3] = np.tan(ro[3]*0.5)
else:
ax = Rotation.eu2ax(eu)
ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)])
ro[ax[...,3]>=np.pi,3] = np.inf
ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, _P, 0.0 ]
return ro return ro
@staticmethod @staticmethod
@ -652,17 +801,24 @@ class Rotation:
@staticmethod @staticmethod
def ax2qu(ax): def ax2qu(ax):
"""Axis angle pair to quaternion.""" """Axis angle pair to quaternion."""
if iszero(ax[3]): if len(ax.shape) == 1:
if np.abs(ax[3])<1.e-6:
qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) qu = np.array([ 1.0, 0.0, 0.0, 0.0 ])
else: else:
c = np.cos(ax[3]*0.5) c = np.cos(ax[3]*0.5)
s = np.sin(ax[3]*0.5) s = np.sin(ax[3]*0.5)
qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ])
return qu return qu
else:
c = np.cos(ax[...,3:4]*.5)
s = np.sin(ax[...,3:4]*.5)
qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s]))
return qu
@staticmethod @staticmethod
def ax2om(ax): def ax2om(ax):
"""Axis angle pair to rotation matrix.""" """Axis angle pair to rotation matrix."""
if len(ax.shape) == 1:
c = np.cos(ax[3]) c = np.cos(ax[3])
s = np.sin(ax[3]) s = np.sin(ax[3])
omc = 1.0-c omc = 1.0-c
@ -672,7 +828,20 @@ class Rotation:
q = omc*ax[idx[0]] * ax[idx[1]] q = omc*ax[idx[0]] * ax[idx[1]]
om[idx[0],idx[1]] = q + s*ax[idx[2]] om[idx[0],idx[1]] = q + s*ax[idx[2]]
om[idx[1],idx[0]] = q - s*ax[idx[2]] om[idx[1],idx[0]] = q - s*ax[idx[2]]
return om if P < 0.0 else om.T else:
c = np.cos(ax[...,3:4])
s = np.sin(ax[...,3:4])
omc = 1. -c
om = np.block([c+omc*ax[...,0:1]**2,
omc*ax[...,0:1]*ax[...,1:2] + s*ax[...,2:3],
omc*ax[...,0:1]*ax[...,2:3] - s*ax[...,1:2],
omc*ax[...,0:1]*ax[...,1:2] - s*ax[...,2:3],
c+omc*ax[...,1:2]**2,
omc*ax[...,1:2]*ax[...,2:3] + s*ax[...,0:1],
omc*ax[...,0:1]*ax[...,2:3] + s*ax[...,1:2],
omc*ax[...,1:2]*ax[...,2:3] - s*ax[...,0:1],
c+omc*ax[...,2:3]**2]).reshape(ax.shape[:-1]+(3,3))
return om if _P < 0.0 else np.swapaxes(om,(-1,-2))
@staticmethod @staticmethod
def ax2eu(ax): def ax2eu(ax):
@ -682,21 +851,35 @@ class Rotation:
@staticmethod @staticmethod
def ax2ro(ax): def ax2ro(ax):
"""Axis angle pair to Rodrigues-Frank vector.""" """Axis angle pair to Rodrigues-Frank vector."""
if iszero(ax[3]): if len(ax.shape) == 1:
ro = [ 0.0, 0.0, P, 0.0 ] if np.abs(ax[3])<1.e-6:
ro = [ 0.0, 0.0, _P, 0.0 ]
else: else:
ro = [ax[0], ax[1], ax[2]] ro = [ax[0], ax[1], ax[2]]
# 180 degree case # 180 degree case
ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \
[np.tan(ax[3]*0.5)] [np.tan(ax[3]*0.5)]
return np.array(ro) return np.array(ro)
else:
ro = np.block([ax[...,:3],
np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0),
np.inf,
np.tan(ax[...,3:4]*0.5))
])
ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0]
return ro
@staticmethod @staticmethod
def ax2ho(ax): def ax2ho(ax):
"""Axis angle pair to homochoric vector.""" """Axis angle pair to homochoric vector."""
if len(ax.shape) == 1:
f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
ho = ax[0:3] * f ho = ax[0:3] * f
return ho return ho
else:
f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0)
ho = ax[...,:3] * f
return ho
@staticmethod @staticmethod
def ax2cu(ax): def ax2cu(ax):
@ -723,27 +906,38 @@ class Rotation:
@staticmethod @staticmethod
def ro2ax(ro): def ro2ax(ro):
"""Rodrigues-Frank vector to axis angle pair.""" """Rodrigues-Frank vector to axis angle pair."""
ta = ro[3] if len(ro.shape) == 1:
if np.abs(ro[3]) < 1.e-6:
if iszero(ta): ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
ax = [ 0.0, 0.0, 1.0, 0.0 ] elif not np.isfinite(ro[3]):
elif not np.isfinite(ta): ax = np.array([ ro[0], ro[1], ro[2], np.pi ])
ax = [ ro[0], ro[1], ro[2], np.pi ]
else: else:
angle = 2.0*np.arctan(ta) angle = 2.0*np.arctan(ro[3])
ta = 1.0/np.linalg.norm(ro[0:3]) ta = np.linalg.norm(ro[0:3])
ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ])
return np.array(ax) else:
with np.errstate(invalid='ignore',divide='ignore'):
ax = np.where(np.isfinite(ro[...,3:4]),
np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]),
np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)]))
ax[np.abs(ro[...,3]) < 1.e-6] = np.array([ 0.0, 0.0, 1.0, 0.0 ])
return ax
@staticmethod @staticmethod
def ro2ho(ro): def ro2ho(ro):
"""Rodrigues-Frank vector to homochoric vector.""" """Rodrigues-Frank vector to homochoric vector."""
if iszero(np.sum(ro[0:3]**2.0)): if len(ro.shape) == 1:
ho = [ 0.0, 0.0, 0.0 ] if np.sum(ro[0:3]**2.0) < 1.e-6:
ho = np.zeros(3)
else: else:
f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi
ho = ro[0:3] * (0.75*f)**(1.0/3.0) ho = ro[0:3] * (0.75*f)**(1.0/3.0)
return np.array(ho) else:
f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi)
ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape),
np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0))
return ho
@staticmethod @staticmethod
def ro2cu(ro): def ro2cu(ro):
@ -778,6 +972,7 @@ class Rotation:
+0.0001703481934140054, -0.00012062065004116828, +0.0001703481934140054, -0.00012062065004116828,
+0.000059719705868660826, -0.00001980756723965647, +0.000059719705868660826, -0.00001980756723965647,
+0.000003953714684212874, -0.00000036555001439719544]) +0.000003953714684212874, -0.00000036555001439719544])
if len(ho.shape) == 1:
# normalize h and store the magnitude # normalize h and store the magnitude
hmag_squared = np.sum(ho**2.) hmag_squared = np.sum(ho**2.)
if iszero(hmag_squared): if iszero(hmag_squared):
@ -791,6 +986,17 @@ class Rotation:
hm *= hmag_squared hm *= hmag_squared
s += tfit[i] * hm s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
else:
hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True)
hm = hmag_squared.copy()
s = tfit[0] + tfit[1] * hmag_squared
for i in range(2,16):
hm *= hmag_squared
s += tfit[i] * hm
with np.errstate(invalid='ignore'):
ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-6,ho.shape[:-1]+(4,)),
[ 0.0, 0.0, 1.0, 0.0 ],
np.block([ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))]))
return ax return ax
@staticmethod @staticmethod
@ -801,7 +1007,10 @@ class Rotation:
@staticmethod @staticmethod
def ho2cu(ho): def ho2cu(ho):
"""Homochoric vector to cubochoric vector.""" """Homochoric vector to cubochoric vector."""
if len(ho.shape) == 1:
return ball_to_cube(ho) return ball_to_cube(ho)
else:
raise NotImplementedError
#---------- Cubochoric ---------- #---------- Cubochoric ----------
@ -833,4 +1042,7 @@ class Rotation:
@staticmethod @staticmethod
def cu2ho(cu): def cu2ho(cu):
"""Cubochoric vector to homochoric vector.""" """Cubochoric vector to homochoric vector."""
if len(cu.shape) == 1:
return cube_to_ball(cu) return cube_to_ball(cu)
else:
raise NotImplementedError

220
python/tests/test_Rotation.py
View File

@ -5,12 +5,82 @@ import numpy as np
from damask import Rotation from damask import Rotation
n = 1000 n = 1100
atol=1.e-4
scatter=1.e-2
@pytest.fixture @pytest.fixture
def default(): def default():
"""A set of n random rotations.""" """A set of n random rotations."""
return [Rotation.fromRandom() for r in range(n)] specials = np.array(
[np.array([ 1.0, 0.0, 0.0, 0.0]),
#-----------------------------------------------
np.array([0.0, 1.0, 0.0, 0.0]),
np.array([0.0, 0.0, 1.0, 0.0]),
np.array([0.0, 0.0, 0.0, 1.0]),
np.array([0.0,-1.0, 0.0, 0.0]),
np.array([0.0, 0.0,-1.0, 0.0]),
np.array([0.0, 0.0, 0.0,-1.0]),
#-----------------------------------------------
np.array([1.0, 1.0, 0.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0, 1.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0, 0.0, 1.0])/np.sqrt(2.),
np.array([0.0, 1.0, 1.0, 0.0])/np.sqrt(2.),
np.array([0.0, 1.0, 0.0, 1.0])/np.sqrt(2.),
np.array([0.0, 0.0, 1.0, 1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([1.0,-1.0, 0.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0,-1.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 1.0,-1.0, 0.0])/np.sqrt(2.),
np.array([0.0, 1.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 0.0, 1.0,-1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([0.0, 1.0,-1.0, 0.0])/np.sqrt(2.),
np.array([0.0, 1.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 0.0, 1.0,-1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([0.0,-1.0,-1.0, 0.0])/np.sqrt(2.),
np.array([0.0,-1.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 0.0,-1.0,-1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([1.0, 1.0, 1.0, 0.0])/np.sqrt(3.),
np.array([1.0, 1.0, 0.0, 1.0])/np.sqrt(3.),
np.array([1.0, 0.0, 1.0, 1.0])/np.sqrt(3.),
np.array([1.0,-1.0, 1.0, 0.0])/np.sqrt(3.),
np.array([1.0,-1.0, 0.0, 1.0])/np.sqrt(3.),
np.array([1.0, 0.0,-1.0, 1.0])/np.sqrt(3.),
np.array([1.0, 1.0,-1.0, 0.0])/np.sqrt(3.),
np.array([1.0, 1.0, 0.0,-1.0])/np.sqrt(3.),
np.array([1.0, 0.0, 1.0,-1.0])/np.sqrt(3.),
np.array([1.0,-1.0,-1.0, 0.0])/np.sqrt(3.),
np.array([1.0,-1.0, 0.0,-1.0])/np.sqrt(3.),
np.array([1.0, 0.0,-1.0,-1.0])/np.sqrt(3.),
#-----------------------------------------------
np.array([0.0, 1.0, 1.0, 1.0])/np.sqrt(3.),
np.array([0.0, 1.0,-1.0, 1.0])/np.sqrt(3.),
np.array([0.0, 1.0, 1.0,-1.0])/np.sqrt(3.),
np.array([0.0,-1.0, 1.0, 1.0])/np.sqrt(3.),
np.array([0.0,-1.0,-1.0, 1.0])/np.sqrt(3.),
np.array([0.0,-1.0, 1.0,-1.0])/np.sqrt(3.),
np.array([0.0,-1.0,-1.0,-1.0])/np.sqrt(3.),
#-----------------------------------------------
np.array([1.0, 1.0, 1.0, 1.0])/2.,
np.array([1.0,-1.0, 1.0, 1.0])/2.,
np.array([1.0, 1.0,-1.0, 1.0])/2.,
np.array([1.0, 1.0, 1.0,-1.0])/2.,
np.array([1.0,-1.0,-1.0, 1.0])/2.,
np.array([1.0,-1.0, 1.0,-1.0])/2.,
np.array([1.0, 1.0,-1.0,-1.0])/2.,
np.array([1.0,-1.0,-1.0,-1.0])/2.,
])
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1
return [Rotation.fromQuaternion(s) for s in specials] + \
[Rotation.fromQuaternion(s) for s in specials_scatter] + \
[Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))]
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
@ -22,35 +92,151 @@ class TestRotation:
def test_Eulers(self,default): def test_Eulers(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asQuaternion(), m = rot.asQuaternion()
Rotation.fromEulers(rot.asEulers()).asQuaternion()) o = Rotation.fromEulers(rot.asEulers()).asQuaternion()
ok = np.allclose(m,o,atol=atol)
if np.isclose(rot.asQuaternion()[0],0.0,atol=atol):
ok = ok or np.allclose(m*-1.,o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.isclose(np.linalg.norm(o),1.0)
def test_AxisAngle(self,default): def test_AxisAngle(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asEulers(), m = rot.asEulers()
Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers()) o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers()
u = np.array([np.pi*2,np.pi,np.pi*2])
ok = np.allclose(m,o,atol=atol)
ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol):
sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]])
ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol)
print(m,o,rot.asQuaternion())
assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
def test_Matrix(self,default): def test_Matrix(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asAxisAngle(), m = rot.asAxisAngle()
Rotation.fromMatrix(rot.asMatrix()).asAxisAngle()) o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle()
ok = np.allclose(m,o,atol=atol)
if np.isclose(m[3],np.pi,atol=atol):
ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9
def test_Rodriques(self,default): def test_Rodrigues(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asMatrix(), m = rot.asMatrix()
Rotation.fromRodrigues(rot.asRodrigues()).asMatrix()) o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix()
ok = np.allclose(m,o,atol=atol)
print(m,o)
assert ok and np.isclose(np.linalg.det(o),1.0)
def test_Homochoric(self,default): def test_Homochoric(self,default):
cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in default: for rot in default:
assert np.allclose(rot.asRodrigues(), m = rot.asRodrigues()
Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=1.e-4) o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0)
def test_Cubochoric(self,default): def test_Cubochoric(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asHomochoric(), m = rot.asHomochoric()
Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric()) o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric()
ok = np.allclose(m,o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9
def test_Quaternion(self,default): def test_Quaternion(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asCubochoric(), m = rot.asCubochoric()
Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()) o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()
ok = np.allclose(m,o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
@pytest.mark.parametrize('conversion',[Rotation.qu2om,
Rotation.qu2eu,
Rotation.qu2ax,
Rotation.qu2ro,
Rotation.qu2ho])
def test_quaternion_vectorization(self,default,conversion):
qu = np.array([rot.asQuaternion() for rot in default])
conversion(qu.reshape(qu.shape[0]//2,-1,4))
co = conversion(qu)
for q,c in zip(qu,co):
print(q,c)
assert np.allclose(conversion(q),c)
@pytest.mark.parametrize('conversion',[Rotation.om2qu,
Rotation.om2eu,
Rotation.om2ax,
Rotation.om2ro,
Rotation.om2ho,
])
def test_matrix_vectorization(self,default,conversion):
om = np.array([rot.asMatrix() for rot in default])
conversion(om.reshape(om.shape[0]//2,-1,3,3))
co = conversion(om)
for o,c in zip(om,co):
print(o,c)
assert np.allclose(conversion(o),c)
@pytest.mark.parametrize('conversion',[Rotation.eu2qu,
Rotation.eu2om,
Rotation.eu2ax,
Rotation.eu2ro,
Rotation.eu2ho,
])
def test_Euler_vectorization(self,default,conversion):
eu = np.array([rot.asEulers() for rot in default])
conversion(eu.reshape(eu.shape[0]//2,-1,3))
co = conversion(eu)
for e,c in zip(eu,co):
print(e,c)
assert np.allclose(conversion(e),c)
@pytest.mark.parametrize('conversion',[Rotation.ax2qu,
Rotation.ax2om,
Rotation.ax2eu,
Rotation.ax2ro,
Rotation.ax2ho,
])
def test_axisAngle_vectorization(self,default,conversion):
ax = np.array([rot.asAxisAngle() for rot in default])
conversion(ax.reshape(ax.shape[0]//2,-1,4))
co = conversion(ax)
for a,c in zip(ax,co):
print(a,c)
assert np.allclose(conversion(a),c)
@pytest.mark.parametrize('conversion',[Rotation.ro2qu,
Rotation.ro2om,
Rotation.ro2eu,
Rotation.ro2ax,
Rotation.ro2ho,
])
def test_Rodrigues_vectorization(self,default,conversion):
ro = np.array([rot.asRodrigues() for rot in default])
conversion(ro.reshape(ro.shape[0]//2,-1,4))
co = conversion(ro)
for r,c in zip(ro,co):
print(r,c)
assert np.allclose(conversion(r),c)
@pytest.mark.parametrize('conversion',[Rotation.ho2qu,
Rotation.ho2om,
Rotation.ho2eu,
Rotation.ho2ax,
Rotation.ho2ro,
])
def test_homochoric_vectorization(self,default,conversion):
ho = np.array([rot.asHomochoric() for rot in default])
conversion(ho.reshape(ho.shape[0]//2,-1,3))
co = conversion(ho)
for h,c in zip(ho,co):
print(h,c)
assert np.allclose(conversion(h),c)

20
src/Lambert.f90
View File

@ -75,7 +75,7 @@ pure function Lambert_CubeToBall(cube) result(ball)
real(pReal), dimension(2) :: T real(pReal), dimension(2) :: T
real(pReal) :: c, s, q real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal real(pReal), parameter :: eps = 1.0e-8_pReal
integer, dimension(3) :: p integer, dimension(3,2) :: p
integer, dimension(2) :: order integer, dimension(2) :: order
if (maxval(abs(cube)) > AP/2.0+eps) then if (maxval(abs(cube)) > AP/2.0+eps) then
@ -89,7 +89,7 @@ pure function Lambert_CubeToBall(cube) result(ball)
else center else center
! get pyramide and scale by grid parameter ratio ! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube) p = GetPyramidOrder(cube)
XYZ = cube(p) * sc XYZ = cube(p(:,1)) * sc
! intercept all the points along the z-axis ! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then special: if (all(dEq0(XYZ(1:2)))) then
@ -112,7 +112,7 @@ pure function Lambert_CubeToBall(cube) result(ball)
endif special endif special
! reverse the coordinates back to order according to the original pyramid number ! reverse the coordinates back to order according to the original pyramid number
ball = LamXYZ(p) ball = LamXYZ(p(:,2))
endif center endif center
@ -130,7 +130,7 @@ pure function Lambert_BallToCube(xyz) result(cube)
real(pReal), dimension(3) :: cube, xyz1, xyz3 real(pReal), dimension(3) :: cube, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2 real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt real(pReal) :: rs, qxy, q2, sq2, q, tt
integer, dimension(3) :: p integer, dimension(3,2) :: p
rs = norm2(xyz) rs = norm2(xyz)
if (rs > R1) then if (rs > R1) then
@ -142,7 +142,7 @@ pure function Lambert_BallToCube(xyz) result(cube)
cube = 0.0_pReal cube = 0.0_pReal
else center else center
p = GetPyramidOrder(xyz) p = GetPyramidOrder(xyz)
xyz3 = xyz(p) xyz3 = xyz(p(:,1))
! inverse M_3 ! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
@ -166,7 +166,7 @@ pure function Lambert_BallToCube(xyz) result(cube)
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc
! reverse the coordinates back to order according to the original pyramid number ! reverse the coordinates back to order according to the original pyramid number
cube = xyz1(p) cube = xyz1(p(:,2))
endif center endif center
@ -181,17 +181,17 @@ end function Lambert_BallToCube
pure function GetPyramidOrder(xyz) pure function GetPyramidOrder(xyz)
real(pReal),intent(in),dimension(3) :: xyz real(pReal),intent(in),dimension(3) :: xyz
integer, dimension(3) :: GetPyramidOrder integer, dimension(3,2) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
GetPyramidOrder = [1,2,3] GetPyramidOrder = reshape([[1,2,3],[1,2,3]],[3,2])
else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. &
((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
GetPyramidOrder = [2,3,1] GetPyramidOrder = reshape([[2,3,1],[3,1,2]],[3,2])
else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. &
((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
GetPyramidOrder = [3,1,2] GetPyramidOrder = reshape([[3,1,2],[2,3,1]],[3,2])
else else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if end if

11
src/rotations.f90
View File

@ -432,18 +432,17 @@ pure function qu2eu(qu) result(eu)
real(pReal), intent(in), dimension(4) :: qu real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: eu real(pReal), dimension(3) :: eu
real(pReal) :: q12, q03, chi, chiInv real(pReal) :: q12, q03, chi
q03 = qu(1)**2+qu(4)**2 q03 = qu(1)**2+qu(4)**2
q12 = qu(2)**2+qu(3)**2 q12 = qu(2)**2+qu(3)**2
chi = sqrt(q03*q12) chi = sqrt(q03*q12)
degenerated: if (dEq0(chi)) then degenerated: if (dEq0(q12)) then
eu = merge([atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal], & eu = [atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal]
[atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal], & elseif (dEq0(q03)) then
dEq0(q12)) eu = [atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal]
else degenerated else degenerated
chiInv = 1.0_pReal/chi
eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), & eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
atan2( 2.0_pReal*chi, q03-q12 ), & atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]