Merge branch 'no-python-quaternion' into 'development'
No python quaternion See merge request damask/DAMASK!122
This commit is contained in:
commit
87767f0cc2
|
@ -3,9 +3,8 @@ import math
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
|
||||||
from . import Lambert
|
from . import Lambert
|
||||||
from .quaternion import Quaternion
|
|
||||||
from .quaternion import P
|
|
||||||
|
|
||||||
|
P = -1
|
||||||
|
|
||||||
####################################################################################################
|
####################################################################################################
|
||||||
class Rotation:
|
class Rotation:
|
||||||
|
@ -50,10 +49,7 @@ class Rotation:
|
||||||
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
|
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
if isinstance(quaternion,Quaternion):
|
self.quaternion = quaternion.copy()
|
||||||
self.quaternion = quaternion.copy()
|
|
||||||
else:
|
|
||||||
self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4])
|
|
||||||
|
|
||||||
def __copy__(self):
|
def __copy__(self):
|
||||||
"""Copy."""
|
"""Copy."""
|
||||||
|
@ -65,9 +61,9 @@ class Rotation:
|
||||||
def __repr__(self):
|
def __repr__(self):
|
||||||
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
|
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
|
||||||
return '\n'.join([
|
return '\n'.join([
|
||||||
'{}'.format(self.quaternion),
|
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
|
||||||
'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
|
'Matrix:\n{}'.format(self.asMatrix()),
|
||||||
'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
|
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)),
|
||||||
])
|
])
|
||||||
|
|
||||||
def __mul__(self, other):
|
def __mul__(self, other):
|
||||||
|
@ -86,19 +82,25 @@ class Rotation:
|
||||||
|
|
||||||
"""
|
"""
|
||||||
if isinstance(other, Rotation): # rotate a rotation
|
if isinstance(other, Rotation): # rotate a rotation
|
||||||
return self.__class__(self.quaternion * other.quaternion).standardize()
|
self_q = self.quaternion[0]
|
||||||
elif isinstance(other, np.ndarray):
|
self_p = self.quaternion[1:]
|
||||||
if other.shape == (3,): # rotate a single (3)-vector
|
other_q = other.quaternion[0]
|
||||||
( x, y, z) = self.quaternion.p
|
other_p = other.quaternion[1:]
|
||||||
(Vx,Vy,Vz) = other[0:3]
|
R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
|
||||||
A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p)
|
self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p)))
|
||||||
B = 2.0 * (x*Vx + y*Vy + z*Vz)
|
return R.standardize()
|
||||||
C = 2.0 * P*self.quaternion.q
|
elif isinstance(other, (tuple,np.ndarray)):
|
||||||
|
if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid
|
||||||
|
A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
|
||||||
|
B = 2.0 * ( self.quaternion[1]*other[0]
|
||||||
|
+ self.quaternion[2]*other[1]
|
||||||
|
+ self.quaternion[3]*other[2])
|
||||||
|
C = 2.0 * P*self.quaternion[0]
|
||||||
|
|
||||||
return np.array([
|
return np.array([
|
||||||
A*Vx + B*x + C*(y*Vz - z*Vy),
|
A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
|
||||||
A*Vy + B*y + C*(z*Vx - x*Vz),
|
A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
|
||||||
A*Vz + B*z + C*(x*Vy - y*Vx),
|
A*other[2] + B*self.quaternion[3] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
|
||||||
])
|
])
|
||||||
elif other.shape == (3,3,): # rotate a single (3x3)-matrix
|
elif other.shape == (3,3,): # rotate a single (3x3)-matrix
|
||||||
return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
|
return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
|
||||||
|
@ -106,25 +108,13 @@ class Rotation:
|
||||||
raise NotImplementedError
|
raise NotImplementedError
|
||||||
else:
|
else:
|
||||||
return NotImplemented
|
return NotImplemented
|
||||||
elif isinstance(other, tuple): # used to rotate a meshgrid-tuple
|
|
||||||
( x, y, z) = self.quaternion.p
|
|
||||||
(Vx,Vy,Vz) = other[0:3]
|
|
||||||
A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p)
|
|
||||||
B = 2.0 * (x*Vx + y*Vy + z*Vz)
|
|
||||||
C = 2.0 * P*self.quaternion.q
|
|
||||||
|
|
||||||
return np.array([
|
|
||||||
A*Vx + B*x + C*(y*Vz - z*Vy),
|
|
||||||
A*Vy + B*y + C*(z*Vx - x*Vz),
|
|
||||||
A*Vz + B*z + C*(x*Vy - y*Vx),
|
|
||||||
])
|
|
||||||
else:
|
else:
|
||||||
return NotImplemented
|
return NotImplemented
|
||||||
|
|
||||||
|
|
||||||
def inverse(self):
|
def inverse(self):
|
||||||
"""In-place inverse rotation/backward rotation."""
|
"""In-place inverse rotation/backward rotation."""
|
||||||
self.quaternion.conjugate()
|
self.quaternion[1:] *= -1
|
||||||
return self
|
return self
|
||||||
|
|
||||||
def inversed(self):
|
def inversed(self):
|
||||||
|
@ -134,7 +124,7 @@ class Rotation:
|
||||||
|
|
||||||
def standardize(self):
|
def standardize(self):
|
||||||
"""In-place quaternion representation with positive q."""
|
"""In-place quaternion representation with positive q."""
|
||||||
if self.quaternion.q < 0.0: self.quaternion.homomorph()
|
if self.quaternion[0] < 0.0: self.quaternion*=-1
|
||||||
return self
|
return self
|
||||||
|
|
||||||
def standardized(self):
|
def standardized(self):
|
||||||
|
@ -171,8 +161,7 @@ class Rotation:
|
||||||
################################################################################################
|
################################################################################################
|
||||||
# convert to different orientation representations (numpy arrays)
|
# convert to different orientation representations (numpy arrays)
|
||||||
|
|
||||||
def asQuaternion(self,
|
def asQuaternion(self):
|
||||||
quaternion = False):
|
|
||||||
"""
|
"""
|
||||||
Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
|
Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
|
||||||
|
|
||||||
|
@ -182,7 +171,7 @@ class Rotation:
|
||||||
return quaternion as DAMASK object.
|
return quaternion as DAMASK object.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return self.quaternion if quaternion else self.quaternion.asArray()
|
return self.quaternion
|
||||||
|
|
||||||
def asEulers(self,
|
def asEulers(self,
|
||||||
degrees = False):
|
degrees = False):
|
||||||
|
@ -195,7 +184,7 @@ class Rotation:
|
||||||
return angles in degrees.
|
return angles in degrees.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
eu = qu2eu(self.quaternion.asArray())
|
eu = qu2eu(self.quaternion)
|
||||||
if degrees: eu = np.degrees(eu)
|
if degrees: eu = np.degrees(eu)
|
||||||
return eu
|
return eu
|
||||||
|
|
||||||
|
@ -213,13 +202,13 @@ class Rotation:
|
||||||
return tuple of axis and angle.
|
return tuple of axis and angle.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
ax = qu2ax(self.quaternion.asArray())
|
ax = qu2ax(self.quaternion)
|
||||||
if degrees: ax[3] = np.degrees(ax[3])
|
if degrees: ax[3] = np.degrees(ax[3])
|
||||||
return (ax[:3],np.degrees(ax[3])) if pair else ax
|
return (ax[:3],np.degrees(ax[3])) if pair else ax
|
||||||
|
|
||||||
def asMatrix(self):
|
def asMatrix(self):
|
||||||
"""Rotation matrix."""
|
"""Rotation matrix."""
|
||||||
return qu2om(self.quaternion.asArray())
|
return qu2om(self.quaternion)
|
||||||
|
|
||||||
def asRodrigues(self,
|
def asRodrigues(self,
|
||||||
vector = False):
|
vector = False):
|
||||||
|
@ -232,16 +221,16 @@ class Rotation:
|
||||||
return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
|
return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
|
||||||
|
|
||||||
"""
|
"""
|
||||||
ro = qu2ro(self.quaternion.asArray())
|
ro = qu2ro(self.quaternion)
|
||||||
return ro[:3]*ro[3] if vector else ro
|
return ro[:3]*ro[3] if vector else ro
|
||||||
|
|
||||||
def asHomochoric(self):
|
def asHomochoric(self):
|
||||||
"""Homochoric vector: (h_1, h_2, h_3)."""
|
"""Homochoric vector: (h_1, h_2, h_3)."""
|
||||||
return qu2ho(self.quaternion.asArray())
|
return qu2ho(self.quaternion)
|
||||||
|
|
||||||
def asCubochoric(self):
|
def asCubochoric(self):
|
||||||
"""Cubochoric vector: (c_1, c_2, c_3)."""
|
"""Cubochoric vector: (c_1, c_2, c_3)."""
|
||||||
return qu2cu(self.quaternion.asArray())
|
return qu2cu(self.quaternion)
|
||||||
|
|
||||||
def asM(self):
|
def asM(self):
|
||||||
"""
|
"""
|
||||||
|
@ -253,19 +242,18 @@ class Rotation:
|
||||||
https://doi.org/10.2514/1.28949
|
https://doi.org/10.2514/1.28949
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return self.quaternion.asM()
|
return np.outer(self.quaternion,self.quaternion)
|
||||||
|
|
||||||
|
|
||||||
################################################################################################
|
################################################################################################
|
||||||
# static constructors. The input data needs to follow the convention, options allow to
|
# static constructors. The input data needs to follow the convention, options allow to
|
||||||
# relax these convections
|
# relax these convections
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromQuaternion(cls,
|
def fromQuaternion(quaternion,
|
||||||
quaternion,
|
|
||||||
acceptHomomorph = False,
|
acceptHomomorph = False,
|
||||||
P = -1):
|
P = -1):
|
||||||
|
|
||||||
qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \
|
qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \
|
||||||
else np.array(quaternion,dtype=float)
|
else np.array(quaternion,dtype=float)
|
||||||
if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1
|
if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1
|
||||||
if qu[0] < 0.0:
|
if qu[0] < 0.0:
|
||||||
|
@ -276,11 +264,10 @@ class Rotation:
|
||||||
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.\n{} {} {} {}'.format(*qu))
|
||||||
|
|
||||||
return cls(qu)
|
return Rotation(qu)
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromEulers(cls,
|
def fromEulers(eulers,
|
||||||
eulers,
|
|
||||||
degrees = False):
|
degrees = False):
|
||||||
|
|
||||||
eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
|
eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
|
||||||
|
@ -289,11 +276,10 @@ class Rotation:
|
||||||
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π].\n{} {} {}.'.format(*eu))
|
||||||
|
|
||||||
return cls(eu2qu(eu))
|
return Rotation(eu2qu(eu))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromAxisAngle(cls,
|
def fromAxisAngle(angleAxis,
|
||||||
angleAxis,
|
|
||||||
degrees = False,
|
degrees = False,
|
||||||
normalise = False,
|
normalise = False,
|
||||||
P = -1):
|
P = -1):
|
||||||
|
@ -308,11 +294,10 @@ class Rotation:
|
||||||
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.\n{} {} {}'.format(*ax[0:3]))
|
||||||
|
|
||||||
return cls(ax2qu(ax))
|
return Rotation(ax2qu(ax))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromBasis(cls,
|
def fromBasis(basis,
|
||||||
basis,
|
|
||||||
orthonormal = True,
|
orthonormal = True,
|
||||||
reciprocal = False,
|
reciprocal = False,
|
||||||
):
|
):
|
||||||
|
@ -331,18 +316,16 @@ class Rotation:
|
||||||
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.\n{}'.format(om))
|
||||||
|
|
||||||
return cls(om2qu(om))
|
return Rotation(om2qu(om))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromMatrix(cls,
|
def fromMatrix(om,
|
||||||
om,
|
|
||||||
):
|
):
|
||||||
|
|
||||||
return cls.fromBasis(om)
|
return Rotation.fromBasis(om)
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromRodrigues(cls,
|
def fromRodrigues(rodrigues,
|
||||||
rodrigues,
|
|
||||||
normalise = False,
|
normalise = False,
|
||||||
P = -1):
|
P = -1):
|
||||||
|
|
||||||
|
@ -355,35 +338,32 @@ class Rotation:
|
||||||
if ro[3] < 0.0:
|
if ro[3] < 0.0:
|
||||||
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
|
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
|
||||||
|
|
||||||
return cls(ro2qu(ro))
|
return Rotation(ro2qu(ro))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromHomochoric(cls,
|
def fromHomochoric(homochoric,
|
||||||
homochoric,
|
P = -1):
|
||||||
P = -1):
|
|
||||||
|
|
||||||
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
|
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
|
||||||
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
|
||||||
|
|
||||||
return cls(ho2qu(ho))
|
return Rotation(ho2qu(ho))
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromCubochoric(cls,
|
def fromCubochoric(cubochoric,
|
||||||
cubochoric,
|
P = -1):
|
||||||
P = -1):
|
|
||||||
|
|
||||||
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)
|
||||||
ho = cu2ho(cu)
|
ho = 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
|
||||||
|
|
||||||
return cls(ho2qu(ho))
|
return Rotation(ho2qu(ho))
|
||||||
|
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromAverage(cls,
|
def fromAverage(rotations,
|
||||||
rotations,
|
|
||||||
weights = []):
|
weights = []):
|
||||||
"""
|
"""
|
||||||
Average rotation.
|
Average rotation.
|
||||||
|
@ -413,19 +393,18 @@ class Rotation:
|
||||||
else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa
|
else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa
|
||||||
eig, vec = np.linalg.eig(M/N)
|
eig, vec = np.linalg.eig(M/N)
|
||||||
|
|
||||||
return cls.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
|
return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
|
||||||
|
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromRandom(cls):
|
def fromRandom():
|
||||||
r = np.random.random(3)
|
r = np.random.random(3)
|
||||||
A = np.sqrt(r[2])
|
A = np.sqrt(r[2])
|
||||||
B = np.sqrt(1.0-r[2])
|
B = np.sqrt(1.0-r[2])
|
||||||
w = np.cos(2.0*np.pi*r[0])*A
|
return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
|
||||||
x = np.sin(2.0*np.pi*r[1])*B
|
np.sin(2.0*np.pi*r[1])*B,
|
||||||
y = np.cos(2.0*np.pi*r[1])*B
|
np.cos(2.0*np.pi*r[1])*B,
|
||||||
z = np.sin(2.0*np.pi*r[0])*A
|
np.sin(2.0*np.pi*r[0])*A])).standardize()
|
||||||
return cls.fromQuaternion([w,x,y,z],acceptHomomorph=True)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
@ -1194,9 +1173,8 @@ class Orientation:
|
||||||
return color
|
return color
|
||||||
|
|
||||||
|
|
||||||
@classmethod
|
@staticmethod
|
||||||
def fromAverage(cls,
|
def fromAverage(orientations,
|
||||||
orientations,
|
|
||||||
weights = []):
|
weights = []):
|
||||||
"""Create orientation from average of list of orientations."""
|
"""Create orientation from average of list of orientations."""
|
||||||
if not all(isinstance(item, Orientation) for item in orientations):
|
if not all(isinstance(item, Orientation) for item in orientations):
|
||||||
|
|
|
@ -1,212 +0,0 @@
|
||||||
import numpy as np
|
|
||||||
|
|
||||||
P = -1 # convention (see DOI:10.1088/0965-0393/23/8/083501)
|
|
||||||
|
|
||||||
####################################################################################################
|
|
||||||
class Quaternion:
|
|
||||||
u"""
|
|
||||||
Quaternion with basic operations
|
|
||||||
|
|
||||||
q is the real part, p = (x, y, z) are the imaginary parts.
|
|
||||||
Definition of multiplication depends on variable P, P ∈ {-1,1}.
|
|
||||||
"""
|
|
||||||
|
|
||||||
def __init__(self,
|
|
||||||
q = 0.0,
|
|
||||||
p = np.zeros(3,dtype=float)):
|
|
||||||
"""Initializes to identity unless specified"""
|
|
||||||
self.q = float(q)
|
|
||||||
self.p = np.array(p,dtype=float)
|
|
||||||
|
|
||||||
|
|
||||||
def __copy__(self):
|
|
||||||
"""Copy"""
|
|
||||||
return self.__class__(q=self.q,
|
|
||||||
p=self.p.copy())
|
|
||||||
|
|
||||||
copy = __copy__
|
|
||||||
|
|
||||||
|
|
||||||
def __iter__(self):
|
|
||||||
"""Components"""
|
|
||||||
return iter(self.asList())
|
|
||||||
|
|
||||||
def __repr__(self):
|
|
||||||
"""Readable string"""
|
|
||||||
return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
|
|
||||||
|
|
||||||
|
|
||||||
def __add__(self, other):
|
|
||||||
"""Addition"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
return self.__class__(q=self.q + other.q,
|
|
||||||
p=self.p + other.p)
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __iadd__(self, other):
|
|
||||||
"""In-place addition"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
self.q += other.q
|
|
||||||
self.p += other.p
|
|
||||||
return self
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __pos__(self):
|
|
||||||
"""Unary positive operator"""
|
|
||||||
return self
|
|
||||||
|
|
||||||
|
|
||||||
def __sub__(self, other):
|
|
||||||
"""Subtraction"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
return self.__class__(q=self.q - other.q,
|
|
||||||
p=self.p - other.p)
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __isub__(self, other):
|
|
||||||
"""In-place subtraction"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
self.q -= other.q
|
|
||||||
self.p -= other.p
|
|
||||||
return self
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __neg__(self):
|
|
||||||
"""Unary negative operator"""
|
|
||||||
return self * -1.0
|
|
||||||
|
|
||||||
|
|
||||||
def __mul__(self, other):
|
|
||||||
"""Multiplication with quaternion or scalar"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
|
|
||||||
p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
|
|
||||||
elif isinstance(other, (int, float)):
|
|
||||||
return self.__class__(q=self.q*other,
|
|
||||||
p=self.p*other)
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __imul__(self, other):
|
|
||||||
"""In-place multiplication with quaternion or scalar"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
self.q = self.q*other.q - np.dot(self.p,other.p)
|
|
||||||
self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
|
|
||||||
return self
|
|
||||||
elif isinstance(other, (int, float)):
|
|
||||||
self.q *= other
|
|
||||||
self.p *= other
|
|
||||||
return self
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
|
|
||||||
def __truediv__(self, other):
|
|
||||||
"""Divsion with quaternion or scalar"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
s = other.conjugate()/abs(other)**2.
|
|
||||||
return self.__class__(q=self.q * s,
|
|
||||||
p=self.p * s)
|
|
||||||
elif isinstance(other, (int, float)):
|
|
||||||
return self.__class__(q=self.q / other,
|
|
||||||
p=self.p / other)
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __itruediv__(self, other):
|
|
||||||
"""In-place divsion with quaternion or scalar"""
|
|
||||||
if isinstance(other, Quaternion):
|
|
||||||
s = other.conjugate()/abs(other)**2.
|
|
||||||
self *= s
|
|
||||||
return self
|
|
||||||
elif isinstance(other, (int, float)):
|
|
||||||
self.q /= other
|
|
||||||
self.p /= other
|
|
||||||
return self
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
|
|
||||||
def __pow__(self, exponent):
|
|
||||||
"""Power"""
|
|
||||||
if isinstance(exponent, (int, float)):
|
|
||||||
omega = np.acos(self.q)
|
|
||||||
return self.__class__(q= np.cos(exponent*omega),
|
|
||||||
p=self.p * np.sin(exponent*omega)/np.sin(omega))
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
def __ipow__(self, exponent):
|
|
||||||
"""In-place power"""
|
|
||||||
if isinstance(exponent, (int, float)):
|
|
||||||
omega = np.acos(self.q)
|
|
||||||
self.q = np.cos(exponent*omega)
|
|
||||||
self.p *= np.sin(exponent*omega)/np.sin(omega)
|
|
||||||
return self
|
|
||||||
else:
|
|
||||||
return NotImplemented
|
|
||||||
|
|
||||||
|
|
||||||
def __abs__(self):
|
|
||||||
"""Norm"""
|
|
||||||
return np.sqrt(self.q ** 2 + np.dot(self.p,self.p))
|
|
||||||
|
|
||||||
magnitude = __abs__
|
|
||||||
|
|
||||||
|
|
||||||
def __eq__(self,other):
|
|
||||||
"""Equal (sufficiently close) to each other"""
|
|
||||||
return np.isclose(( self-other).magnitude(),0.0) \
|
|
||||||
or np.isclose((-self-other).magnitude(),0.0)
|
|
||||||
|
|
||||||
def __ne__(self,other):
|
|
||||||
"""Not equal (sufficiently close) to each other"""
|
|
||||||
return not self.__eq__(other)
|
|
||||||
|
|
||||||
|
|
||||||
def asM(self):
|
|
||||||
"""Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)"""
|
|
||||||
return np.outer(self.asArray(),self.asArray())
|
|
||||||
|
|
||||||
def asArray(self):
|
|
||||||
"""As numpy array"""
|
|
||||||
return np.array((self.q,self.p[0],self.p[1],self.p[2]))
|
|
||||||
|
|
||||||
def asList(self):
|
|
||||||
"""As list"""
|
|
||||||
return [self.q]+list(self.p)
|
|
||||||
|
|
||||||
|
|
||||||
def normalize(self):
|
|
||||||
"""Normalizes in-place"""
|
|
||||||
d = self.magnitude()
|
|
||||||
if d > 0.0: self /= d
|
|
||||||
return self
|
|
||||||
|
|
||||||
def normalized(self):
|
|
||||||
"""Returns normalized copy"""
|
|
||||||
return self.copy().normalize()
|
|
||||||
|
|
||||||
|
|
||||||
def conjugate(self):
|
|
||||||
"""Conjugates in-place"""
|
|
||||||
self.p = -self.p
|
|
||||||
return self
|
|
||||||
|
|
||||||
def conjugated(self):
|
|
||||||
"""Returns conjugated copy"""
|
|
||||||
return self.copy().conjugate()
|
|
||||||
|
|
||||||
|
|
||||||
def homomorph(self):
|
|
||||||
"""Homomorphs in-place"""
|
|
||||||
self *= -1.0
|
|
||||||
return self
|
|
||||||
|
|
||||||
def homomorphed(self):
|
|
||||||
"""Returns homomorphed copy"""
|
|
||||||
return self.copy().homomorph()
|
|
|
@ -1,4 +1,5 @@
|
||||||
import os
|
import os
|
||||||
|
from itertools import permutations
|
||||||
|
|
||||||
import pytest
|
import pytest
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
@ -6,6 +7,7 @@ import numpy as np
|
||||||
import damask
|
import damask
|
||||||
from damask import Rotation
|
from damask import Rotation
|
||||||
from damask import Orientation
|
from damask import Orientation
|
||||||
|
from damask import Lattice
|
||||||
|
|
||||||
n = 1000
|
n = 1000
|
||||||
|
|
||||||
|
@ -45,19 +47,37 @@ class TestRotation:
|
||||||
def test_Homochoric(self,default):
|
def test_Homochoric(self,default):
|
||||||
for rot in default:
|
for rot in default:
|
||||||
assert np.allclose(rot.asRodrigues(),
|
assert np.allclose(rot.asRodrigues(),
|
||||||
Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=5.e-5)
|
Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=1.e-4)
|
||||||
|
|
||||||
def test_Cubochoric(self,default):
|
def test_Cubochoric(self,default):
|
||||||
for rot in default:
|
for rot in default:
|
||||||
assert np.allclose(rot.asHomochoric(),
|
assert np.allclose(rot.asHomochoric(),
|
||||||
Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric(),rtol=5.e-5)
|
Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric())
|
||||||
|
|
||||||
def test_Quaternion(self,default):
|
def test_Quaternion(self,default):
|
||||||
for rot in default:
|
for rot in default:
|
||||||
assert np.allclose(rot.asCubochoric(),
|
assert np.allclose(rot.asCubochoric(),
|
||||||
Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric(),rtol=5.e-5)
|
Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric())
|
||||||
|
|
||||||
|
|
||||||
|
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]},
|
||||||
|
{'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
|
||||||
|
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
|
||||||
|
@pytest.mark.parametrize('lattice',['fcc','bcc'])
|
||||||
|
def test_IPF_cubic(self,default,color,lattice):
|
||||||
|
cube = damask.Orientation(damask.Rotation(),lattice)
|
||||||
|
for direction in set(permutations(np.array(color['direction']))):
|
||||||
|
assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
|
||||||
|
|
||||||
|
@pytest.mark.parametrize('lattice',Lattice.lattices)
|
||||||
|
def test_IPF(self,lattice):
|
||||||
|
direction = np.random.random(3)*2.0-1
|
||||||
|
for rot in [Rotation.fromRandom() for r in range(n//100)]:
|
||||||
|
R = damask.Orientation(rot,lattice)
|
||||||
|
color = R.IPFcolor(direction)
|
||||||
|
for equivalent in R.equivalentOrientations():
|
||||||
|
assert np.allclose(color,R.IPFcolor(direction))
|
||||||
|
|
||||||
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
|
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
|
||||||
@pytest.mark.parametrize('lattice',['fcc','bcc'])
|
@pytest.mark.parametrize('lattice',['fcc','bcc'])
|
||||||
def test_relationship_forward_backward(self,model,lattice):
|
def test_relationship_forward_backward(self,model,lattice):
|
||||||
|
|
Loading…
Reference in New Issue