vectorizing cubochoric conversions
This commit is contained in:
Martin Diehl
parent
7d1e0850ab
commit
5f3f87cd68
|
|
@ -241,7 +241,7 @@ class Rotation:
|
|||
"""Homochoric vector: (h_1, h_2, h_3)."""
|
||||
return Rotation.qu2ho(self.quaternion)
|
||||
|
||||
def asCubochoric(self):
|
||||
def as_cubochoric(self):
|
||||
"""Cubochoric vector: (c_1, c_2, c_3)."""
|
||||
return Rotation.qu2cu(self.quaternion)
|
||||
|
||||
|
|
@ -265,6 +265,7 @@ class Rotation:
|
|||
asMatrix = as_matrix
|
||||
asRodrigues = as_Rodrigues
|
||||
asHomochoric = as_homochoric
|
||||
asCubochoric = as_cubochoric
|
||||
|
||||
################################################################################################
|
||||
# Static constructors. The input data needs to follow the conventions, options allow to
|
||||
|
|
@ -386,7 +387,7 @@ class Rotation:
|
|||
return Rotation(Rotation.ho2qu(ho))
|
||||
|
||||
@staticmethod
|
||||
def fromCubochoric(cubochoric,
|
||||
def from_cubochoric(cubochoric,
|
||||
P = -1):
|
||||
|
||||
cu = np.array(cubochoric,dtype=float)
|
||||
|
|
@ -461,6 +462,7 @@ class Rotation:
|
|||
fromMatrix = from_matrix
|
||||
fromRodrigues = from_Rodrigues
|
||||
fromHomochoric = from_homochoric
|
||||
fromCubochoric = from_cubochoric
|
||||
fromRandom = from_random
|
||||
|
||||
####################################################################################################
|
||||
|
|
@ -1066,7 +1068,7 @@ class Rotation:
|
|||
if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16):
|
||||
cu = np.zeros(3)
|
||||
else:
|
||||
xyz3 = ho[Rotation._get_order(ho,'forward')]
|
||||
xyz3 = ho[Rotation._get_pyramid_order(ho,'forward')]
|
||||
|
||||
# inverse M_3
|
||||
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
|
||||
|
|
@ -1087,20 +1089,35 @@ class Rotation:
|
|||
|
||||
# inverse M_1
|
||||
cu = 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
|
||||
cu = cu[Rotation._get_order(ho,'backward')]
|
||||
cu = cu[Rotation._get_pyramid_order(ho,'backward')]
|
||||
else:
|
||||
rs = np.linalg.norm(ho,axis=-1,keepdims=True)
|
||||
|
||||
xyz3 = np.take_along_axis(ho,Rotation._get_pyramid_order(ho,'forward'),-1)
|
||||
|
||||
with np.errstate(invalid='ignore',divide='ignore'):
|
||||
# inverse M_3
|
||||
xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) )
|
||||
qxy = np.sum(xyz2**2,axis=-1,keepdims=True)
|
||||
|
||||
q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2
|
||||
sq2 = np.sqrt(q2)
|
||||
q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2))
|
||||
tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\
|
||||
+np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
|
||||
T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]),
|
||||
np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]),
|
||||
np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q
|
||||
T_inv[xyz2<0.0] *= -1.0
|
||||
T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0
|
||||
cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \
|
||||
* rs/np.sqrt(6.0/np.pi),
|
||||
])/ _sc
|
||||
|
||||
cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0
|
||||
cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1)
|
||||
|
||||
return cu
|
||||
else:
|
||||
q2 = qxy + np.max(np.abs(xyz2))**2
|
||||
sq2 = np.sqrt(q2)
|
||||
q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2))
|
||||
tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
|
||||
Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \
|
||||
np.array([np.arccos(tt)/np.pi*12.0,1.0])
|
||||
Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv)
|
||||
raise NotImplementedError('Support for multiple rotations missing')
|
||||
|
||||
|
||||
#---------- Cubochoric ----------
|
||||
@staticmethod
|
||||
|
|
@ -1145,7 +1162,7 @@ class Rotation:
|
|||
ho = np.zeros(3)
|
||||
else:
|
||||
# get pyramide and scale by grid parameter ratio
|
||||
XYZ = cu[Rotation._get_order(cu,'forward')] * _sc
|
||||
XYZ = cu[Rotation._get_pyramid_order(cu,'forward')] * _sc
|
||||
|
||||
# intercept all the points along the z-axis
|
||||
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16):
|
||||
|
|
@ -1163,31 +1180,46 @@ class Rotation:
|
|||
c = np.sum(T**2)
|
||||
s = c * np.pi/24.0 /XYZ[2]**2
|
||||
c = c * np.sqrt(np.pi/24.0)/XYZ[2]
|
||||
|
||||
q = np.sqrt( 1.0 - s )
|
||||
ho = 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
|
||||
ho = ho[Rotation._get_order(cu,'backward')]
|
||||
|
||||
return ho
|
||||
ho = ho[Rotation._get_pyramid_order(cu,'backward')]
|
||||
else:
|
||||
with np.errstate(invalid='ignore',divide='ignore'):
|
||||
# get pyramide and scale by grid parameter ratio
|
||||
XYZ = cu[Rotation._get_order(cu,'forward')] * _sc
|
||||
order = np.where(np.abs(XYZ[...,1]) <= np.abs(XYZ[...,0]),
|
||||
np.broadcast_to([1,0],XYZ.shape[:-1]+(2,)),
|
||||
np.broadcast_to([0,1],XYZ.shape[:-1]+(2,)))
|
||||
q = np.pi/12.0 * XYZ[order[...,0]]/XYZ[order[...,1]]
|
||||
XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc
|
||||
order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1])
|
||||
q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \
|
||||
/ np.where(order,XYZ[...,0:1],XYZ[...,1:2])
|
||||
c = np.cos(q)
|
||||
s = np.sin(q)
|
||||
q = _R1*2.0**0.25/_beta * XYZ[order[...,1]] / np.sqrt(np.sqrt(2.0)-c)
|
||||
q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \
|
||||
* np.where(order,XYZ[...,0:1],XYZ[...,1:2])
|
||||
|
||||
T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
|
||||
s = c * np.pi/24.0 /XYZ[...,2]**2
|
||||
c = c * np.sqrt(np.pi/24.0)/XYZ[...,2]
|
||||
|
||||
# transform to sphere grid (inverse Lambert)
|
||||
c = np.sum(T**2,axis=-1,keepdims=True)
|
||||
s = c * np.pi/24.0 /XYZ[...,2:3]**2
|
||||
c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3]
|
||||
q = np.sqrt( 1.0 - s)
|
||||
raise NotImplementedError('Support for multiple rotations missing')
|
||||
|
||||
ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16),
|
||||
np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]),
|
||||
np.block([np.where(order,T[...,0:1],T[...,1:2])*q,
|
||||
np.where(order,T[...,1:2],T[...,0:1])*q,
|
||||
np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c])
|
||||
)
|
||||
|
||||
ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0
|
||||
ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1)
|
||||
|
||||
return ho
|
||||
|
||||
|
||||
@staticmethod
|
||||
def _get_order(xyz,direction=None):
|
||||
def _get_pyramid_order(xyz,direction=None):
|
||||
"""
|
||||
Get order of the coordinates.
|
||||
|
||||
|
|
|
|||
|
|
@ -78,9 +78,9 @@ def default():
|
|||
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))]
|
||||
return [Rotation.from_quaternion(s) for s in specials] + \
|
||||
[Rotation.from_quaternion(s) for s in specials_scatter] + \
|
||||
[Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))]
|
||||
|
||||
@pytest.fixture
|
||||
def reference_dir(reference_dir_base):
|
||||
|
|
@ -92,41 +92,41 @@ class TestRotation:
|
|||
|
||||
def test_Eulers(self,default):
|
||||
for rot in default:
|
||||
m = rot.asQuaternion()
|
||||
o = Rotation.fromEulers(rot.asEulers()).asQuaternion()
|
||||
m = rot.as_quaternion()
|
||||
o = Rotation.from_Eulers(rot.as_Eulers()).as_quaternion()
|
||||
ok = np.allclose(m,o,atol=atol)
|
||||
if np.isclose(rot.asQuaternion()[0],0.0,atol=atol):
|
||||
if np.isclose(rot.as_quaternion()[0],0.0,atol=atol):
|
||||
ok = ok or np.allclose(m*-1.,o,atol=atol)
|
||||
print(m,o,rot.asQuaternion())
|
||||
print(m,o,rot.as_quaternion())
|
||||
assert ok and np.isclose(np.linalg.norm(o),1.0)
|
||||
|
||||
def test_AxisAngle(self,default):
|
||||
for rot in default:
|
||||
m = rot.asEulers()
|
||||
o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers()
|
||||
m = rot.as_Eulers()
|
||||
o = Rotation.from_axis_angle(rot.as_axis_angle()).as_Eulers()
|
||||
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())
|
||||
print(m,o,rot.as_quaternion())
|
||||
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):
|
||||
for rot in default:
|
||||
m = rot.asAxisAngle()
|
||||
o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle()
|
||||
m = rot.as_axis_angle()
|
||||
o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle()
|
||||
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())
|
||||
print(m,o,rot.as_quaternion())
|
||||
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9
|
||||
|
||||
def test_Rodrigues(self,default):
|
||||
for rot in default:
|
||||
m = rot.asMatrix()
|
||||
o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix()
|
||||
m = rot.as_matrix()
|
||||
o = Rotation.from_Rodrigues(rot.as_Rodrigues()).as_matrix()
|
||||
ok = np.allclose(m,o,atol=atol)
|
||||
print(m,o)
|
||||
assert ok and np.isclose(np.linalg.det(o),1.0)
|
||||
|
|
@ -134,27 +134,27 @@ class TestRotation:
|
|||
def test_Homochoric(self,default):
|
||||
cutoff = np.tan(np.pi*.5*(1.-1e-4))
|
||||
for rot in default:
|
||||
m = rot.asRodrigues()
|
||||
o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues()
|
||||
m = rot.as_Rodrigues()
|
||||
o = Rotation.from_homochoric(rot.as_homochoric()).as_Rodrigues()
|
||||
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())
|
||||
print(m,o,rot.as_quaternion())
|
||||
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0)
|
||||
|
||||
def test_Cubochoric(self,default):
|
||||
for rot in default:
|
||||
m = rot.asHomochoric()
|
||||
o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric()
|
||||
m = rot.as_homochoric()
|
||||
o = Rotation.from_cubochoric(rot.as_cubochoric()).as_homochoric()
|
||||
ok = np.allclose(m,o,atol=atol)
|
||||
print(m,o,rot.asQuaternion())
|
||||
print(m,o,rot.as_quaternion())
|
||||
assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9
|
||||
|
||||
def test_Quaternion(self,default):
|
||||
for rot in default:
|
||||
m = rot.asCubochoric()
|
||||
o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()
|
||||
m = rot.as_cubochoric()
|
||||
o = Rotation.from_quaternion(rot.as_quaternion()).as_cubochoric()
|
||||
ok = np.allclose(m,o,atol=atol)
|
||||
print(m,o,rot.asQuaternion())
|
||||
print(m,o,rot.as_quaternion())
|
||||
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
|
||||
|
||||
@pytest.mark.parametrize('function',[Rotation.from_quaternion,
|
||||
|
|
@ -185,9 +185,11 @@ class TestRotation:
|
|||
Rotation.qu2eu,
|
||||
Rotation.qu2ax,
|
||||
Rotation.qu2ro,
|
||||
Rotation.qu2ho])
|
||||
Rotation.qu2ho,
|
||||
Rotation.qu2cu
|
||||
])
|
||||
def test_quaternion_vectorization(self,default,conversion):
|
||||
qu = np.array([rot.asQuaternion() for rot in default])
|
||||
qu = np.array([rot.as_quaternion() for rot in default])
|
||||
conversion(qu.reshape(qu.shape[0]//2,-1,4))
|
||||
co = conversion(qu)
|
||||
for q,c in zip(qu,co):
|
||||
|
|
@ -199,9 +201,10 @@ class TestRotation:
|
|||
Rotation.om2ax,
|
||||
Rotation.om2ro,
|
||||
Rotation.om2ho,
|
||||
Rotation.om2cu
|
||||
])
|
||||
def test_matrix_vectorization(self,default,conversion):
|
||||
om = np.array([rot.asMatrix() for rot in default])
|
||||
om = np.array([rot.as_matrix() for rot in default])
|
||||
conversion(om.reshape(om.shape[0]//2,-1,3,3))
|
||||
co = conversion(om)
|
||||
for o,c in zip(om,co):
|
||||
|
|
@ -213,9 +216,10 @@ class TestRotation:
|
|||
Rotation.eu2ax,
|
||||
Rotation.eu2ro,
|
||||
Rotation.eu2ho,
|
||||
Rotation.eu2cu
|
||||
])
|
||||
def test_Euler_vectorization(self,default,conversion):
|
||||
eu = np.array([rot.asEulers() for rot in default])
|
||||
eu = np.array([rot.as_Eulers() for rot in default])
|
||||
conversion(eu.reshape(eu.shape[0]//2,-1,3))
|
||||
co = conversion(eu)
|
||||
for e,c in zip(eu,co):
|
||||
|
|
@ -227,9 +231,10 @@ class TestRotation:
|
|||
Rotation.ax2eu,
|
||||
Rotation.ax2ro,
|
||||
Rotation.ax2ho,
|
||||
Rotation.ax2cu
|
||||
])
|
||||
def test_axisAngle_vectorization(self,default,conversion):
|
||||
ax = np.array([rot.asAxisAngle() for rot in default])
|
||||
ax = np.array([rot.as_axis_angle() for rot in default])
|
||||
conversion(ax.reshape(ax.shape[0]//2,-1,4))
|
||||
co = conversion(ax)
|
||||
for a,c in zip(ax,co):
|
||||
|
|
@ -242,9 +247,10 @@ class TestRotation:
|
|||
Rotation.ro2eu,
|
||||
Rotation.ro2ax,
|
||||
Rotation.ro2ho,
|
||||
Rotation.ro2cu
|
||||
])
|
||||
def test_Rodrigues_vectorization(self,default,conversion):
|
||||
ro = np.array([rot.asRodrigues() for rot in default])
|
||||
ro = np.array([rot.as_Rodrigues() for rot in default])
|
||||
conversion(ro.reshape(ro.shape[0]//2,-1,4))
|
||||
co = conversion(ro)
|
||||
for r,c in zip(ro,co):
|
||||
|
|
@ -256,26 +262,41 @@ class TestRotation:
|
|||
Rotation.ho2eu,
|
||||
Rotation.ho2ax,
|
||||
Rotation.ho2ro,
|
||||
Rotation.ho2cu
|
||||
])
|
||||
def test_homochoric_vectorization(self,default,conversion):
|
||||
ho = np.array([rot.asHomochoric() for rot in default])
|
||||
ho = np.array([rot.as_homochoric() 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)
|
||||
|
||||
@pytest.mark.parametrize('conversion',[Rotation.cu2qu,
|
||||
Rotation.cu2om,
|
||||
Rotation.cu2eu,
|
||||
Rotation.cu2ax,
|
||||
Rotation.cu2ro,
|
||||
Rotation.cu2ho
|
||||
])
|
||||
def test_cubochoric_vectorization(self,default,conversion):
|
||||
cu = np.array([rot.as_cubochoric() for rot in default])
|
||||
conversion(cu.reshape(cu.shape[0]//2,-1,3))
|
||||
co = conversion(cu)
|
||||
for u,c in zip(cu,co):
|
||||
print(u,c)
|
||||
assert np.allclose(conversion(u),c)
|
||||
|
||||
@pytest.mark.parametrize('direction',['forward',
|
||||
'backward'])
|
||||
def test_pyramid_vectorization(self,direction):
|
||||
p = np.random.rand(n,3)
|
||||
o = Rotation._get_order(p,direction)
|
||||
o = Rotation._get_pyramid_order(p,direction)
|
||||
for i,o_i in enumerate(o):
|
||||
assert np.allclose(o_i,Rotation._get_order(p[i],direction))
|
||||
assert np.all(o_i==Rotation._get_pyramid_order(p[i],direction))
|
||||
|
||||
def test_pyramid_invariant(self):
|
||||
a = np.random.rand(n,3)
|
||||
f = damask.Rotation._get_order(a,'forward')
|
||||
b = damask.Rotation._get_order(a,'backward')
|
||||
f = Rotation._get_pyramid_order(a,'forward')
|
||||
b = Rotation._get_pyramid_order(a,'backward')
|
||||
assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a)
|
||||
|
|
|
|||
Loading…
Reference in New Issue