# -*- coding: UTF-8 no BOM -*-
###################################################
# NOTE: everything here needs to be a np array #
###################################################
import math,os
import numpy as np
# ******************************************************************************************
class Rodrigues:
def __init__(self, vector = np.zeros(3)):
self.vector = vector
def asQuaternion(self):
norm = np.linalg.norm(self.vector)
halfAngle = np.arctan(norm)
return Quaternion(np.cos(halfAngle),np.sin(halfAngle)*self.vector/norm)
def asAngleAxis(self):
norm = np.linalg.norm(self.vector)
halfAngle = np.arctan(norm)
return (2.0*halfAngle,self.vector/norm)
# ******************************************************************************************
class Quaternion:
"""
Orientation represented as unit quaternion.
All methods and naming conventions based on http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions.
w is the real part, (x, y, z) are the imaginary parts.
Representation of rotation is in ACTIVE form!
(Derived directly or through angleAxis, Euler angles, or active matrix)
Vector "a" (defined in coordinate system "A") is actively rotated to new coordinates "b".
b = Q * a
b = np.dot(Q.asMatrix(),a)
"""
def __init__(self,
quatArray = [1.0,0.0,0.0,0.0]):
"""Initializes to identity if not given"""
self.w, \
self.x, \
self.y, \
self.z = quatArray
self.homomorph()
def __iter__(self):
"""Components"""
return iter([self.w,self.x,self.y,self.z])
def __copy__(self):
"""Create copy"""
Q = Quaternion([self.w,self.x,self.y,self.z])
return Q
copy = __copy__
def __repr__(self):
"""Readbable string"""
return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \
(self.w, self.x, self.y, self.z)
def __pow__(self, exponent):
"""Power"""
omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega)
Q = Quaternion()
Q.w = math.cos(exponent*omega)
Q.x = self.x * vRescale
Q.y = self.y * vRescale
Q.z = self.z * vRescale
return Q
def __ipow__(self, exponent):
"""In-place power"""
omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega)
self.w = np.cos(exponent*omega)
self.x *= vRescale
self.y *= vRescale
self.z *= vRescale
return self
def __mul__(self, other):
"""Multiplication"""
try: # quaternion
Aw = self.w
Ax = self.x
Ay = self.y
Az = self.z
Bw = other.w
Bx = other.x
By = other.y
Bz = other.z
Q = Quaternion()
Q.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw
Q.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx
Q.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By
Q.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz
return Q
except: pass
try: # vector (perform active rotation, i.e. q*v*q.conjugated)
w = self.w
x = self.x
y = self.y
z = self.z
Vx = other[0]
Vy = other[1]
Vz = other[2]
return np.array([\
w * w * Vx + 2 * y * w * Vz - 2 * z * w * Vy + \
x * x * Vx + 2 * y * x * Vy + 2 * z * x * Vz - \
z * z * Vx - y * y * Vx,
2 * x * y * Vx + y * y * Vy + 2 * z * y * Vz + \
2 * w * z * Vx - z * z * Vy + w * w * Vy - \
2 * x * w * Vz - x * x * Vy,
2 * x * z * Vx + 2 * y * z * Vy + \
z * z * Vz - 2 * w * y * Vx - y * y * Vz + \
2 * w * x * Vy - x * x * Vz + w * w * Vz ])
except: pass
try: # scalar
Q = self.copy()
Q.w *= other
Q.x *= other
Q.y *= other
Q.z *= other
return Q
except:
return self.copy()
def __imul__(self, other):
"""In-place multiplication"""
try: # Quaternion
Aw = self.w
Ax = self.x
Ay = self.y
Az = self.z
Bw = other.w
Bx = other.x
By = other.y
Bz = other.z
self.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw
self.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx
self.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By
self.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz
except: pass
return self
def __div__(self, other):
"""Division"""
if isinstance(other, (int,float)):
w = self.w / other
x = self.x / other
y = self.y / other
z = self.z / other
return self.__class__([w,x,y,z])
else:
return NotImplemented
def __idiv__(self, other):
"""In-place division"""
if isinstance(other, (int,float)):
self.w /= other
self.x /= other
self.y /= other
self.z /= other
return self
def __add__(self, other):
"""Addition"""
if isinstance(other, Quaternion):
w = self.w + other.w
x = self.x + other.x
y = self.y + other.y
z = self.z + other.z
return self.__class__([w,x,y,z])
else:
return NotImplemented
def __iadd__(self, other):
"""In-place addition"""
if isinstance(other, Quaternion):
self.w += other.w
self.x += other.x
self.y += other.y
self.z += other.z
return self
def __sub__(self, other):
"""Subtraction"""
if isinstance(other, Quaternion):
Q = self.copy()
Q.w -= other.w
Q.x -= other.x
Q.y -= other.y
Q.z -= other.z
return Q
else:
return self.copy()
def __isub__(self, other):
"""In-place subtraction"""
if isinstance(other, Quaternion):
self.w -= other.w
self.x -= other.x
self.y -= other.y
self.z -= other.z
return self
def __neg__(self):
"""Additive inverse"""
self.w = -self.w
self.x = -self.x
self.y = -self.y
self.z = -self.z
return self
def __abs__(self):
"""Norm"""
return math.sqrt(self.w ** 2 + \
self.x ** 2 + \
self.y ** 2 + \
self.z ** 2)
magnitude = __abs__
def __eq__(self,other):
"""Equal at e-8 precision"""
return (abs(self.w-other.w) < 1e-8 and \
abs(self.x-other.x) < 1e-8 and \
abs(self.y-other.y) < 1e-8 and \
abs(self.z-other.z) < 1e-8) \
or \
(abs(-self.w-other.w) < 1e-8 and \
abs(-self.x-other.x) < 1e-8 and \
abs(-self.y-other.y) < 1e-8 and \
abs(-self.z-other.z) < 1e-8)
def __ne__(self,other):
"""Not equal at e-8 precision"""
return not self.__eq__(self,other)
def __cmp__(self,other):
"""Linear ordering"""
return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues()<other.Rodrigues())
def magnitude_squared(self):
return self.w ** 2 + \
self.x ** 2 + \
self.y ** 2 + \
self.z ** 2
def identity(self):
self.w = 1.
self.x = 0.
self.y = 0.
self.z = 0.
return self
def normalize(self):
d = self.magnitude()
if d > 0.0:
self /= d
return self
def conjugate(self):
self.x = -self.x
self.y = -self.y
self.z = -self.z
return self
def inverse(self):
d = self.magnitude()
if d > 0.0:
self.conjugate()
self /= d
return self
def homomorph(self):
if self.w < 0.0:
self.w = -self.w
self.x = -self.x
self.y = -self.y
self.z = -self.z
return self
def normalized(self):
return self.copy().normalize()
def conjugated(self):
return self.copy().conjugate()
def inversed(self):
return self.copy().inverse()
def homomorphed(self):
return self.copy().homomorph()
def asList(self):
return [i for i in self]
def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.)
return np.outer([i for i in self],[i for i in self])
def asMatrix(self):
return np.array(
[[1.0-2.0*(self.y*self.y+self.z*self.z), 2.0*(self.x*self.y-self.z*self.w), 2.0*(self.x*self.z+self.y*self.w)],
[ 2.0*(self.x*self.y+self.z*self.w), 1.0-2.0*(self.x*self.x+self.z*self.z), 2.0*(self.y*self.z-self.x*self.w)],
[ 2.0*(self.x*self.z-self.y*self.w), 2.0*(self.x*self.w+self.y*self.z), 1.0-2.0*(self.x*self.x+self.y*self.y)]])
def asAngleAxis(self,
degrees = False):
if self.w > 1:
self.normalize()
s = math.sqrt(1. - self.w**2)
x = 2*self.w**2 - 1.
y = 2*self.w * s
angle = math.atan2(y,x)
if angle < 0.0:
angle *= -1.
s *= -1.
return (np.degrees(angle) if degrees else angle,
np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else [self.x / s, self.y / s, self.z / s]))
def asRodrigues(self):
return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w
def asEulers(self,
type = "bunge",
degrees = False,
standardRange = False):
"""
Orientation as Bunge-Euler angles.
Conversion of ACTIVE rotation to Euler angles taken from:
Melcher, A.; Unser, A.; Reichhardt, M.; Nestler, B.; Poetschke, M.; Selzer, M.
Conversion of EBSD data by a quaternion based algorithm to be used for grain structure simulations
Technische Mechanik 30 (2010) pp 401--413.
"""
angles = [0.0,0.0,0.0]
if type.lower() == 'bunge' or type.lower() == 'zxz':
if abs(self.x) < 1e-4 and abs(self.y) < 1e-4:
x = self.w**2 - self.z**2
y = 2.*self.w*self.z
angles[0] = math.atan2(y,x)
elif abs(self.w) < 1e-4 and abs(self.z) < 1e-4:
x = self.x**2 - self.y**2
y = 2.*self.x*self.y
angles[0] = math.atan2(y,x)
angles[1] = math.pi
else:
chi = math.sqrt((self.w**2 + self.z**2)*(self.x**2 + self.y**2))
x = (self.w * self.x - self.y * self.z)/2./chi
y = (self.w * self.y + self.x * self.z)/2./chi
angles[0] = math.atan2(y,x)
x = self.w**2 + self.z**2 - (self.x**2 + self.y**2)
y = 2.*chi
angles[1] = math.atan2(y,x)
x = (self.w * self.x + self.y * self.z)/2./chi
y = (self.z * self.x - self.y * self.w)/2./chi
angles[2] = math.atan2(y,x)
if standardRange:
angles[0] %= 2*math.pi
if angles[1] < 0.0:
angles[1] += math.pi
angles[2] *= -1.0
angles[2] %= 2*math.pi
return np.degrees(angles) if degrees else angles
# # Static constructors
@classmethod
def fromIdentity(cls):
return cls()
@classmethod
def fromRandom(cls,randomSeed = None):
if randomSeed is None:
randomSeed = int(os.urandom(4).encode('hex'), 16)
np.random.seed(randomSeed)
r = np.random.random(3)
w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2])
x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2])
return cls([w,x,y,z])
@classmethod
def fromRodrigues(cls, rodrigues):
if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues)
halfangle = math.atan(np.linalg.norm(rodrigues))
c = math.cos(halfangle)
w = c
x,y,z = c*rodrigues
return cls([w,x,y,z])
@classmethod
def fromAngleAxis(cls,
angle,
axis,
degrees = False):
if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d')
axis = axis.astype(float)/np.linalg.norm(axis)
angle = np.radians(angle) if degrees else angle
s = math.sin(0.5 * angle)
w = math.cos(0.5 * angle)
x = axis[0] * s
y = axis[1] * s
z = axis[2] * s
return cls([w,x,y,z])
@classmethod
def fromEulers(cls,
eulers,
type = 'Bunge',
degrees = False):
if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d')
eulers = np.radians(eulers) if degrees else eulers
c = np.cos(0.5 * eulers)
s = np.sin(0.5 * eulers)
if type.lower() == 'bunge' or type.lower() == 'zxz':
w = c[0] * c[1] * c[2] - s[0] * c[1] * s[2]
x = c[0] * s[1] * c[2] + s[0] * s[1] * s[2]
y = - c[0] * s[1] * s[2] + s[0] * s[1] * c[2]
z = c[0] * c[1] * s[2] + s[0] * c[1] * c[2]
else:
w = c[0] * c[1] * c[2] - s[0] * s[1] * s[2]
x = s[0] * s[1] * c[2] + c[0] * c[1] * s[2]
y = s[0] * c[1] * c[2] + c[0] * s[1] * s[2]
z = c[0] * s[1] * c[2] - s[0] * c[1] * s[2]
return cls([w,x,y,z])
# Modified Method to calculate Quaternion from Orientation Matrix,
# Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
@classmethod
def fromMatrix(cls, m):
if m.shape != (3,3) and np.prod(m.shape) == 9:
m = m.reshape(3,3)
tr = np.trace(m)
if tr > 1e-8:
s = math.sqrt(tr + 1.0)*2.0
return cls(
[ s*0.25,
(m[2,1] - m[1,2])/s,
(m[0,2] - m[2,0])/s,
(m[1,0] - m[0,1])/s,
])
elif m[0,0] > m[1,1] and m[0,0] > m[2,2]:
t = m[0,0] - m[1,1] - m[2,2] + 1.0
s = 2.0*math.sqrt(t)
return cls(
[ (m[2,1] - m[1,2])/s,
s*0.25,
(m[0,1] + m[1,0])/s,
(m[2,0] + m[0,2])/s,
])
elif m[1,1] > m[2,2]:
t = -m[0,0] + m[1,1] - m[2,2] + 1.0
s = 2.0*math.sqrt(t)
return cls(
[ (m[0,2] - m[2,0])/s,
(m[0,1] + m[1,0])/s,
s*0.25,
(m[1,2] + m[2,1])/s,
])
else:
t = -m[0,0] - m[1,1] + m[2,2] + 1.0
s = 2.0*math.sqrt(t)
return cls(
[ (m[1,0] - m[0,1])/s,
(m[2,0] + m[0,2])/s,
(m[1,2] + m[2,1])/s,
s*0.25,
])
@classmethod
def new_interpolate(cls, q1, q2, t):
"""
Interpolation
See http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf
for (another?) way to interpolate quaternions.
"""
assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion)
Q = cls()
costheta = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z
if costheta < 0.:
costheta = -costheta
q1 = q1.conjugated()
elif costheta > 1:
costheta = 1
theta = math.acos(costheta)
if abs(theta) < 0.01:
Q.w = q2.w
Q.x = q2.x
Q.y = q2.y
Q.z = q2.z
return Q
sintheta = math.sqrt(1.0 - costheta * costheta)
if abs(sintheta) < 0.01:
Q.w = (q1.w + q2.w) * 0.5
Q.x = (q1.x + q2.x) * 0.5
Q.y = (q1.y + q2.y) * 0.5
Q.z = (q1.z + q2.z) * 0.5
return Q
ratio1 = math.sin((1 - t) * theta) / sintheta
ratio2 = math.sin(t * theta) / sintheta
Q.w = q1.w * ratio1 + q2.w * ratio2
Q.x = q1.x * ratio1 + q2.x * ratio2
Q.y = q1.y * ratio1 + q2.y * ratio2
Q.z = q1.z * ratio1 + q2.z * ratio2
return Q
# ******************************************************************************************
class Symmetry:
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
def __init__(self, symmetry = None):
"""Lattice with given symmetry, defaults to None"""
if isinstance(symmetry, str) and symmetry.lower() in Symmetry.lattices:
self.lattice = symmetry.lower()
else:
self.lattice = None
def __copy__(self):
"""Copy"""
return self.__class__(self.lattice)
copy = __copy__
def __repr__(self):
"""Readbable string"""
return '%s' % (self.lattice)
def __eq__(self, other):
"""Equal"""
return self.lattice == other.lattice
def __neq__(self, other):
"""Not equal"""
return not self.__eq__(other)
def __cmp__(self,other):
"""Linear ordering"""
myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice)
return (myOrder > otherOrder) - (myOrder < otherOrder)
def symmetryQuats(self,who = []):
"""List of symmetry operations as quaternions."""
if self.lattice == 'cubic':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2) ],
[ 0.0, 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2) ],
[ 0.0, 0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2) ],
[ 0.0, 0.5*math.sqrt(2), 0.0, -0.5*math.sqrt(2) ],
[ 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ],
[ 0.0, -0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ],
[ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ],
[-0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2), 0.0 ],
[-0.5*math.sqrt(2), 0.0, -0.5*math.sqrt(2), 0.0 ],
[-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0, 0.0 ],
[-0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0, 0.0 ],
]
elif self.lattice == 'hexagonal':
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
[-0.5*math.sqrt(3), 0.0, 0.0,-0.5 ],
[ 0.5, 0.0, 0.0, 0.5*math.sqrt(3) ],
[ 0.0,0.0,0.0,1.0 ],
[-0.5, 0.0, 0.0, 0.5*math.sqrt(3) ],
[-0.5*math.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,-0.5*math.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5,-0.5*math.sqrt(3), 0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,-0.5,-0.5*math.sqrt(3), 0.0 ],
[ 0.0, 0.5*math.sqrt(3), 0.5, 0.0 ],
]
elif self.lattice == 'tetragonal':
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
[ 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ],
[ 0.0,-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ],
[ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ],
[-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ],
]
elif self.lattice == 'orthorhombic':
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
]
else:
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
]
return list(map(Quaternion,
np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))]))
def equivalentQuaternions(self,
quaternion,
who = []):
"""List of symmetrically equivalent quaternions based on own symmetry."""
return [quaternion*q for q in self.symmetryQuats(who)]
def inFZ(self,R):
"""Check whether given Rodrigues vector falls into fundamental zone of own symmetry."""
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
# fundamental zone in Rodrigues space is point symmetric around origin
R = abs(R)
if self.lattice == 'cubic':
return math.sqrt(2.0)-1.0 >= R[0] \
and math.sqrt(2.0)-1.0 >= R[1] \
and math.sqrt(2.0)-1.0 >= R[2] \
and 1.0 >= R[0] + R[1] + R[2]
elif self.lattice == 'hexagonal':
return 1.0 >= R[0] and 1.0 >= R[1] and 1.0 >= R[2] \
and 2.0 >= math.sqrt(3)*R[0] + R[1] \
and 2.0 >= math.sqrt(3)*R[1] + R[0] \
and 2.0 >= math.sqrt(3) + R[2]
elif self.lattice == 'tetragonal':
return 1.0 >= R[0] and 1.0 >= R[1] \
and math.sqrt(2.0) >= R[0] + R[1] \
and math.sqrt(2.0) >= R[2] + 1.0
elif self.lattice == 'orthorhombic':
return 1.0 >= R[0] and 1.0 >= R[1] and 1.0 >= R[2]
else:
return True
def inDisorientationSST(self,R):
"""
Check whether given Rodrigues vector (of misorientation) falls into standard stereographic triangle of own symmetry.
Determination of disorientations follow the work of A. Heinz and P. Neumann:
Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals
Acta Cryst. (1991). A47, 780-789
"""
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
epsilon = 0.0
if self.lattice == 'cubic':
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon
elif self.lattice == 'hexagonal':
return R[0] >= math.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'tetragonal':
return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'orthorhombic':
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon
else:
return True
def inSST(self,
vector,
proper = False,
color = False):
"""
Check whether given vector falls into standard stereographic triangle of own symmetry.
proper considers only vectors with z >= 0, hence uses two neighboring SSTs.
Return inverse pole figure color if requested.
"""
# basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
# [1.,0.,1.]/np.sqrt(2.), # direction of green
# [1.,1.,1.]/np.sqrt(3.)]).transpose()), # direction of blue
# 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
# [1.,0.,0.], # direction of green
# [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).transpose()), # direction of blue
# 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
# [1.,0.,0.], # direction of green
# [1.,1.,0.]/np.sqrt(2.)]).transpose()), # direction of blue
# 'orthorhombic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
# [1.,0.,0.], # direction of green
# [0.,1.,0.]]).transpose()), # direction of blue
# }
if self.lattice == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.lattice == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.lattice == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.lattice == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else: # direct exit for unspecified symmetry
if color:
return (True,np.zeros(3,'d'))
else:
return True
v = np.array(vector,dtype = float)
if proper: # check both improper ...
theComponents = np.dot(basis['improper'],v)
inSST = np.all(theComponents >= 0.0)
if not inSST: # ... and proper SST
theComponents = np.dot(basis['proper'],v)
inSST = np.all(theComponents >= 0.0)
else:
v[2] = abs(v[2]) # z component projects identical
theComponents = np.dot(basis['improper'],v) # for positive and negative values
inSST = np.all(theComponents >= 0.0)
if color: # have to return color array
if inSST:
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps
rgb = np.minimum(np.ones(3,'d'),rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1
else:
rgb = np.zeros(3,'d')
return (inSST,rgb)
else:
return inSST
# code derived from https://github.com/ezag/pyeuclid
# suggested reading: http://web.mit.edu/2.998/www/QuaternionReport1.pdf
# ******************************************************************************************
class Orientation:
__slots__ = ['quaternion','symmetry']
def __init__(self,
quaternion = Quaternion.fromIdentity(),
Rodrigues = None,
angleAxis = None,
matrix = None,
Eulers = None,
random = False, # integer to have a fixed seed or True for real random
symmetry = None,
degrees = False,
):
if random: # produce random orientation
if isinstance(random, bool ):
self.quaternion = Quaternion.fromRandom()
else:
self.quaternion = Quaternion.fromRandom(randomSeed=random)
elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,): # based on given Euler angles
self.quaternion = Quaternion.fromEulers(Eulers,type='bunge',degrees=degrees)
elif isinstance(matrix, np.ndarray) : # based on given rotation matrix
self.quaternion = Quaternion.fromMatrix(matrix)
elif isinstance(angleAxis, np.ndarray) and angleAxis.shape == (4,): # based on given angle and rotation axis
self.quaternion = Quaternion.fromAngleAxis(angleAxis[0],angleAxis[1:4],degrees=degrees)
elif isinstance(Rodrigues, np.ndarray) and Rodrigues.shape == (3,): # based on given Rodrigues vector
self.quaternion = Quaternion.fromRodrigues(Rodrigues)
elif isinstance(quaternion, Quaternion): # based on given quaternion
self.quaternion = quaternion.homomorphed()
elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion-like array
self.quaternion = Quaternion(quaternion).homomorphed()
self.symmetry = Symmetry(symmetry)
def __copy__(self):
"""Copy"""
return self.__class__(quaternion=self.quaternion,symmetry=self.symmetry.lattice)
copy = __copy__
def __repr__(self):
"""Value as all implemented representations"""
return 'Symmetry: %s\n' % (self.symmetry) + \
'Quaternion: %s\n' % (self.quaternion) + \
'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \
'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers('bunge',degrees=True))) )
def asQuaternion(self):
return self.quaternion.asList()
def asEulers(self,
type = 'bunge',
degrees = False,
standardRange = False):
return self.quaternion.asEulers(type, degrees, standardRange)
eulers = property(asEulers)
def asRodrigues(self):
return self.quaternion.asRodrigues()
rodrigues = property(asRodrigues)
def asAngleAxis(self,
degrees = False):
return self.quaternion.asAngleAxis(degrees)
angleAxis = property(asAngleAxis)
def asMatrix(self):
return self.quaternion.asMatrix()
matrix = property(asMatrix)
def inFZ(self):
return self.symmetry.inFZ(self.quaternion.asRodrigues())
infz = property(inFZ)
def equivalentQuaternions(self,
who = []):
return self.symmetry.equivalentQuaternions(self.quaternion,who)
def equivalentOrientations(self,
who = []):
return [Orientation(quaternion = q, symmetry = self.symmetry.lattice) for q in self.equivalentQuaternions(who)]
def reduced(self):
"""Transform orientation to fall into fundamental zone according to symmetry"""
for me in self.symmetry.equivalentQuaternions(self.quaternion):
if self.symmetry.inFZ(me.asRodrigues()): break
return Orientation(quaternion=me,symmetry=self.symmetry.lattice)
def disorientation(self,
other,
SST = True):
"""
Disorientation between myself and given other orientation.
Rotation axis falls into SST if SST == True.
(Currently requires same symmetry for both orientations.
Look into A. Heinz and P. Neumann 1991 for cases with differing sym.)
"""
if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.')
misQ = self.quaternion.conjugated()*other.quaternion
mySymQs = self.symmetry.symmetryQuats() if SST else self.symmetry.symmetryQuats()[:1] # take all or only first sym operation
otherSymQs = other.symmetry.symmetryQuats()
for i,sA in enumerate(mySymQs):
for j,sB in enumerate(otherSymQs):
theQ = sA.conjugated()*misQ*sB
for k in range(2):
theQ.conjugate()
breaker = self.symmetry.inFZ(theQ) \
and (not SST or other.symmetry.inDisorientationSST(theQ))
if breaker: break
if breaker: break
if breaker: break
# disorientation, own sym, other sym, self-->other: True, self<--other: False
return (Orientation(quaternion = theQ,symmetry = self.symmetry.lattice),
i,j,k == 1)
def inversePole(self,
axis,
proper = False,
SST = True):
"""Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)"""
if SST: # pole requested to be within SST
for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent quaternions
pole = q.conjugated()*axis # align crystal direction to axis
if self.symmetry.inSST(pole,proper): break # found SST version
else:
pole = self.quaternion.conjugated()*axis # align crystal direction to axis
return (pole,i if SST else 0)
def IPFcolor(self,axis):
"""TSL color of inverse pole figure for given axis"""
color = np.zeros(3,'d')
for q in self.symmetry.equivalentQuaternions(self.quaternion):
pole = q.conjugated()*axis # align crystal direction to axis
inSST,color = self.symmetry.inSST(pole,color=True)
if inSST: break
return color
@classmethod
def average(cls,
orientations,
multiplicity = []):
"""
Average orientation
ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman.
Averaging Quaternions,
Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197.
doi: 10.2514/1.28949
usage:
a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal')
b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal')
avg = Orientation.average([a,b])
"""
if not all(isinstance(item, Orientation) for item in orientations):
raise TypeError("Only instances of Orientation can be averaged.")
N = len(orientations)
if multiplicity == [] or not multiplicity:
multiplicity = np.ones(N,dtype='i')
reference = orientations[0] # take first as reference
for i,(o,n) in enumerate(zip(orientations,multiplicity)):
closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0] # select sym orientation with lowest misorientation
M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa
eig, vec = np.linalg.eig(M/N)
return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])),
symmetry = reference.symmetry.lattice)
def related(self,
relationModel,
direction,
targetSymmetry = 'cubic'):
"""
Orientation relationship
positive number: fcc --> bcc
negative number: bcc --> fcc
"""
if relationModel not in ['KS','GT','GTdash','NW','Pitsch','Bain']: return None
if int(direction) == 0: return None
# KS from S. Morito et al./Journal of Alloys and Compounds 5775 (2013) S587-S592
# for KS rotation matrices also check K. Kitahara et al./Acta Materialia 54 (2006) 1279-1288
# GT from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81
# GT' from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81
# NW from H. Kitahara et al./Materials Characterization 54 (2005) 378-386
# Pitsch from Y. He et al./Acta Materialia 53 (2005) 1179-1190
# Bain from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81
variant = int(abs(direction))-1
(me,other) = (0,1) if direction > 0 else (1,0)
planes = {'KS': \
np.array([[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]]]),
'GT': \
np.array([[[ 1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ -1, 0, 1]],
[[ -1, -1, 1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, 1]],
[[ -1, 1, 1],[ -1, 0, 1]],
[[ -1, 1, 1],[ -1, 1, 0]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 1, 0, 1]],
[[ 1, -1, 1],[ 1, -1, 0]],
[[ 1, -1, 1],[ 0, -1, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 1, 0, 1]],
[[ -1, -1, 1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, 1]],
[[ -1, -1, 1],[ -1, 0, 1]],
[[ -1, 1, 1],[ -1, 1, 0]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ -1, 0, 1]],
[[ 1, -1, 1],[ 1, -1, 0]],
[[ 1, -1, 1],[ 0, -1, 1]],
[[ 1, -1, 1],[ 1, 0, 1]]]),
'GTdash': \
np.array([[[ 7, 17, 17],[ 12, 5, 17]],
[[ 17, 7, 17],[ 17, 12, 5]],
[[ 17, 17, 7],[ 5, 17, 12]],
[[ -7,-17, 17],[-12, -5, 17]],
[[-17, -7, 17],[-17,-12, 5]],
[[-17,-17, 7],[ -5,-17, 12]],
[[ 7,-17,-17],[ 12, -5,-17]],
[[ 17, -7,-17],[ 17,-12, -5]],
[[ 17,-17, -7],[ 5,-17,-12]],
[[ -7, 17,-17],[-12, 5,-17]],
[[-17, 7,-17],[-17, 12, -5]],
[[-17, 17, -7],[ -5, 17,-12]],
[[ 7, 17, 17],[ 12, 17, 5]],
[[ 17, 7, 17],[ 5, 12, 17]],
[[ 17, 17, 7],[ 17, 5, 12]],
[[ -7,-17, 17],[-12,-17, 5]],
[[-17, -7, 17],[ -5,-12, 17]],
[[-17,-17, 7],[-17, -5, 12]],
[[ 7,-17,-17],[ 12,-17, -5]],
[[ 17, -7,-17],[ 5, -12,-17]],
[[ 17,-17, 7],[ 17, -5,-12]],
[[ -7, 17,-17],[-12, 17, -5]],
[[-17, 7,-17],[ -5, 12,-17]],
[[-17, 17, -7],[-17, 5,-12]]]),
'NW': \
np.array([[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]]]),
'Pitsch': \
np.array([[[ 0, 1, 0],[ -1, 0, 1]],
[[ 0, 0, 1],[ 1, -1, 0]],
[[ 1, 0, 0],[ 0, 1, -1]],
[[ 1, 0, 0],[ 0, -1, -1]],
[[ 0, 1, 0],[ -1, 0, -1]],
[[ 0, 0, 1],[ -1, -1, 0]],
[[ 0, 1, 0],[ -1, 0, -1]],
[[ 0, 0, 1],[ -1, -1, 0]],
[[ 1, 0, 0],[ 0, -1, -1]],
[[ 1, 0, 0],[ 0, -1, 1]],
[[ 0, 1, 0],[ 1, 0, -1]],
[[ 0, 0, 1],[ -1, 1, 0]]]),
'Bain': \
np.array([[[ 1, 0, 0],[ 1, 0, 0]],
[[ 0, 1, 0],[ 0, 1, 0]],
[[ 0, 0, 1],[ 0, 0, 1]]]),
}
normals = {'KS': \
np.array([[[ -1, 0, 1],[ -1, -1, 1]],
[[ -1, 0, 1],[ -1, 1, -1]],
[[ 0, 1, -1],[ -1, -1, 1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ 1, -1, 0],[ -1, -1, 1]],
[[ 1, -1, 0],[ -1, 1, -1]],
[[ 1, 0, -1],[ -1, -1, 1]],
[[ 1, 0, -1],[ -1, 1, -1]],
[[ -1, -1, 0],[ -1, -1, 1]],
[[ -1, -1, 0],[ -1, 1, -1]],
[[ 0, 1, 1],[ -1, -1, 1]],
[[ 0, 1, 1],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ 0, -1, 1],[ -1, 1, -1]],
[[ -1, 0, -1],[ -1, -1, 1]],
[[ -1, 0, -1],[ -1, 1, -1]],
[[ 1, 1, 0],[ -1, -1, 1]],
[[ 1, 1, 0],[ -1, 1, -1]],
[[ -1, 1, 0],[ -1, -1, 1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, -1],[ -1, -1, 1]],
[[ 0, -1, -1],[ -1, 1, -1]],
[[ 1, 0, 1],[ -1, -1, 1]],
[[ 1, 0, 1],[ -1, 1, -1]]]),
'GT': \
np.array([[[ -5,-12, 17],[-17, -7, 17]],
[[ 17, -5,-12],[ 17,-17, -7]],
[[-12, 17, -5],[ -7, 17,-17]],
[[ 5, 12, 17],[ 17, 7, 17]],
[[-17, 5,-12],[-17, 17, -7]],
[[ 12,-17, -5],[ 7,-17,-17]],
[[ -5, 12,-17],[-17, 7,-17]],
[[ 17, 5, 12],[ 17, 17, 7]],
[[-12,-17, 5],[ -7,-17, 17]],
[[ 5,-12,-17],[ 17, -7,-17]],
[[-17, -5, 12],[-17,-17, 7]],
[[ 12, 17, 5],[ 7, 17, 17]],
[[ -5, 17,-12],[-17, 17, -7]],
[[-12, -5, 17],[ -7,-17, 17]],
[[ 17,-12, -5],[ 17, -7,-17]],
[[ 5,-17,-12],[ 17,-17, -7]],
[[ 12, 5, 17],[ 7, 17, 17]],
[[-17, 12, -5],[-17, 7,-17]],
[[ -5,-17, 12],[-17,-17, 7]],
[[-12, 5,-17],[ -7, 17,-17]],
[[ 17, 12, 5],[ 17, 7, 17]],
[[ 5, 17, 12],[ 17, 17, 7]],
[[ 12, -5,-17],[ 7,-17,-17]],
[[-17,-12, 5],[-17, 7, 17]]]),
'GTdash': \
np.array([[[ 0, 1, -1],[ 1, 1, -1]],
[[ -1, 0, 1],[ -1, 1, 1]],
[[ 1, -1, 0],[ 1, -1, 1]],
[[ 0, -1, -1],[ -1, -1, -1]],
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, -1, 0],[ 1, -1, -1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ 1, 0, 1],[ 1, 1, 1]],
[[ -1, -1, 0],[ -1, -1, 1]],
[[ 0, -1, -1],[ 1, -1, -1]],
[[ -1, 0, 1],[ -1, -1, 1]],
[[ -1, -1, 0],[ -1, -1, -1]],
[[ 0, -1, 1],[ 1, -1, 1]],
[[ 1, 0, -1],[ 1, 1, -1]],
[[ -1, 1, 0],[ -1, 1, 1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ -1, 0, -1],[ -1, -1, -1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ -1, 0, -1],[ -1, 1, -1]],
[[ 1, 1, 0],[ 1, 1, 1]],
[[ 0, 1, 1],[ 1, 1, 1]],
[[ 1, 0, -1],[ 1, -1, -1]],
[[ 1, 1, 0],[ 1, 1, -1]]]),
'NW': \
np.array([[[ 2, -1, -1],[ 0, -1, 1]],
[[ -1, 2, -1],[ 0, -1, 1]],
[[ -1, -1, 2],[ 0, -1, 1]],
[[ -2, -1, -1],[ 0, -1, 1]],
[[ 1, 2, -1],[ 0, -1, 1]],
[[ 1, -1, 2],[ 0, -1, 1]],
[[ 2, 1, -1],[ 0, -1, 1]],
[[ -1, -2, -1],[ 0, -1, 1]],
[[ -1, 1, 2],[ 0, -1, 1]],
[[ -1, 2, 1],[ 0, -1, 1]],
[[ -1, 2, 1],[ 0, -1, 1]],
[[ -1, -1, -2],[ 0, -1, 1]]]),
'Pitsch': \
np.array([[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, 1, 0],[ 1, 1, -1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ -1, 0, 1],[ -1, -1, 1]],
[[ 1, -1, 0],[ 1, -1, -1]],
[[ 1, 0, -1],[ 1, -1, -1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, 1, 0],[ 1, 1, -1]]]),
'Bain': \
np.array([[[ 0, 1, 0],[ 0, 1, 1]],
[[ 0, 0, 1],[ 1, 0, 1]],
[[ 1, 0, 0],[ 1, 1, 0]]]),
}
myPlane = [float(i) for i in planes[relationModel][variant,me]] # map(float, planes[...]) does not work in python 3
myPlane /= np.linalg.norm(myPlane)
myNormal = [float(i) for i in normals[relationModel][variant,me]] # map(float, planes[...]) does not work in python 3
myNormal /= np.linalg.norm(myNormal)
myMatrix = np.array([myNormal,np.cross(myPlane,myNormal),myPlane]).T
otherPlane = [float(i) for i in planes[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3
otherPlane /= np.linalg.norm(otherPlane)
otherNormal = [float(i) for i in normals[relationModel][variant,other]] # map(float, planes[...]) does not work in python 3
otherNormal /= np.linalg.norm(otherNormal)
otherMatrix = np.array([otherNormal,np.cross(otherPlane,otherNormal),otherPlane]).T
rot=np.dot(otherMatrix,myMatrix.T)
return Orientation(matrix=np.dot(rot,self.asMatrix()),symmetry=targetSymmetry)