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

add getAverageOrientation() to Orientation class as class method.

This commit is contained in:
Chen Zhang 2015-04-02 19:15:09 +00:00
parent 86f39de462
commit 43b665ba48
1 changed files with 50 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
lib/damask/orientation.py
View File

@ -1,5 +1,9 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
###################################################
# NOTE: everything here needs to be a numpy array #
###################################################
import numpy,math,random import numpy,math,random
# ****************************************************************************************** # ******************************************************************************************
@ -8,7 +12,7 @@ class Rodrigues:
def __init__(self, vector = numpy.zeros(3)): def __init__(self, vector = numpy.zeros(3)):
self.vector = vector self.vector = vector
def asQuaternion(self): def asQuaternion(self):
norm = numpy.linalg.norm(self.vector) norm = numpy.linalg.norm(self.vector)
halfAngle = numpy.arctan(norm) halfAngle = numpy.arctan(norm)
@ -24,7 +28,7 @@ class Rodrigues:
# ****************************************************************************************** # ******************************************************************************************
class Quaternion: class Quaternion:
# ****************************************************************************************** # ******************************************************************************************
# All methods and naming conventions based off # All methods and naming conventions based off
# http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions # http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions
# w is the real part, (x, y, z) are the imaginary parts # w is the real part, (x, y, z) are the imaginary parts
@ -58,7 +62,7 @@ class Quaternion:
Q.z = self.z * vRescale Q.z = self.z * vRescale
Q.w = math.cos(exponent*omega) Q.w = math.cos(exponent*omega)
return Q return Q
def __ipow__(self, exponent): def __ipow__(self, exponent):
omega = math.acos(self.w) omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega) vRescale = math.sin(exponent*omega)/math.sin(omega)
@ -67,7 +71,7 @@ class Quaternion:
self.z *= vRescale self.z *= vRescale
self.w = numpy.cos(exponent*omega) self.w = numpy.cos(exponent*omega)
return self return self
def __mul__(self, other): def __mul__(self, other):
try: # quaternion try: # quaternion
Ax = self.x Ax = self.x
@ -125,13 +129,13 @@ class Quaternion:
By = other.y By = other.y
Bz = other.z Bz = other.z
Bw = other.w Bw = other.w
self.x = Ax * Bw + Ay * Bz - Az * By + Aw * Bx self.x = Ax * Bw + Ay * Bz - Az * By + Aw * Bx
self.y = -Ax * Bz + Ay * Bw + Az * Bx + Aw * By self.y = -Ax * Bz + Ay * Bw + Az * Bx + Aw * By
self.z = Ax * By - Ay * Bx + Az * Bw + Aw * Bz self.z = Ax * By - Ay * Bx + Az * Bw + Aw * Bz
self.w = -Ax * Bx - Ay * By - Az * Bz + Aw * Bw self.w = -Ax * Bx - Ay * By - Az * Bz + Aw * Bw
except: pass except: pass
return self return self
def __div__(self, other): def __div__(self, other):
if isinstance(other, (int,float,long)): if isinstance(other, (int,float,long)):
w = self.w / other w = self.w / other
@ -211,7 +215,7 @@ class Quaternion:
(abs(-self.w-other.w) < 1e-8 and \ (abs(-self.w-other.w) < 1e-8 and \
abs(-self.x-other.x) < 1e-8 and \ abs(-self.x-other.x) < 1e-8 and \
abs(-self.y-other.y) < 1e-8 and \ abs(-self.y-other.y) < 1e-8 and \
abs(-self.z-other.z) < 1e-8) abs(-self.z-other.z) < 1e-8)
def __ne__(self,other): def __ne__(self,other):
return not __eq__(self,other) return not __eq__(self,other)
@ -223,7 +227,7 @@ class Quaternion:
return self.w ** 2 + \ return self.w ** 2 + \
self.x ** 2 + \ self.x ** 2 + \
self.y ** 2 + \ self.y ** 2 + \
self.z ** 2 self.z ** 2
def identity(self): def identity(self):
self.w = 1. self.w = 1.
@ -293,7 +297,7 @@ class Quaternion:
return numpy.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)], return numpy.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.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)]]) [ 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): def asAngleAxis(self):
if self.w > 1: if self.w > 1:
self.normalize() self.normalize()
@ -301,9 +305,9 @@ class Quaternion:
s = math.sqrt(1. - self.w**2) s = math.sqrt(1. - self.w**2)
x = 2*self.w**2 - 1. x = 2*self.w**2 - 1.
y = 2*self.w * s y = 2*self.w * s
angle = math.atan2(y,x) angle = math.atan2(y,x)
if angle < 1e-3: if angle < 1e-3:
return angle, numpy.array([1.0, 0.0, 0.0]) return angle, numpy.array([1.0, 0.0, 0.0])
else: else:
@ -356,7 +360,7 @@ class Quaternion:
# angles[2] *= -1 # angles[2] *= -1
# if angles[2] < 0.0: # if angles[2] < 0.0:
# angles[2] += 2*math.pi # angles[2] += 2*math.pi
return angles return angles
@ -435,40 +439,40 @@ class Quaternion:
(m[2,0] - m[0,2])*s, (m[2,0] - m[0,2])*s,
(m[0,1] - m[1,0])*s, (m[0,1] - m[1,0])*s,
]) ])
elif m[0,0] > m[1,1] and m[0,0] > m[2,2]: 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 t = m[0,0] - m[1,1] - m[2,2] + 1.0
s = 0.5/math.sqrt(t) s = 0.5/math.sqrt(t)
return cls( return cls(
[ (m[1,2] - m[2,1])*s, [ (m[1,2] - m[2,1])*s,
s*t, s*t,
(m[0,1] + m[1,0])*s, (m[0,1] + m[1,0])*s,
(m[2,0] + m[0,2])*s, (m[2,0] + m[0,2])*s,
]) ])
elif m[1,1] > m[2,2]: elif m[1,1] > m[2,2]:
t = -m[0,0] + m[1,1] - m[2,2] + 1.0 t = -m[0,0] + m[1,1] - m[2,2] + 1.0
s = 0.5/math.sqrt(t) s = 0.5/math.sqrt(t)
return cls( return cls(
[ (m[2,0] - m[0,2])*s, [ (m[2,0] - m[0,2])*s,
(m[0,1] + m[1,0])*s, (m[0,1] + m[1,0])*s,
s*t, s*t,
(m[1,2] + m[2,1])*s, (m[1,2] + m[2,1])*s,
]) ])
else: else:
t = -m[0,0] - m[1,1] + m[2,2] + 1.0 t = -m[0,0] - m[1,1] + m[2,2] + 1.0
s = 0.5/math.sqrt(t) s = 0.5/math.sqrt(t)
return cls( return cls(
[ (m[0,1] - m[1,0])*s, [ (m[0,1] - m[1,0])*s,
(m[2,0] + m[0,2])*s, (m[2,0] + m[0,2])*s,
(m[1,2] + m[2,1])*s, (m[1,2] + m[2,1])*s,
s*t, s*t,
]) ])
@classmethod @classmethod
def new_interpolate(cls, q1, q2, t): def new_interpolate(cls, q1, q2, t):
@ -515,7 +519,7 @@ class Symmetry:
# ****************************************************************************************** # ******************************************************************************************
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',] lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
def __init__(self, symmetry = None): def __init__(self, symmetry = None):
if isinstance(symmetry, basestring) and symmetry.lower() in Symmetry.lattices: if isinstance(symmetry, basestring) and symmetry.lower() in Symmetry.lattices:
self.lattice = symmetry.lower() self.lattice = symmetry.lower()
@ -635,7 +639,7 @@ class Symmetry:
return 1.0 >= R[0] and 1.0 >= R[1] \ 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[0] + R[1] \
and math.sqrt(2.0) >= R[2] + 1.0 and math.sqrt(2.0) >= R[2] + 1.0
elif self.lattice == 'orthorhombic': elif self.lattice == 'orthorhombic':
return 1.0 >= R[0] and 1.0 >= R[1] and 1.0 >= R[2] return 1.0 >= R[0] and 1.0 >= R[1] and 1.0 >= R[2]
else: else:
return True return True
@ -651,7 +655,7 @@ class Symmetry:
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
epsilon = 0.0 epsilon = 0.0
if self.lattice == 'cubic': if self.lattice == 'cubic':
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon and self.inFZ(R) return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon and self.inFZ(R)
@ -709,7 +713,7 @@ class Symmetry:
theComponents = -numpy.ones(3,'d') theComponents = -numpy.ones(3,'d')
else: else:
theComponents = numpy.dot(basis,numpy.array([vector[0],vector[1],abs(vector[2])])) theComponents = numpy.dot(basis,numpy.array([vector[0],vector[1],abs(vector[2])]))
inSST = numpy.all(theComponents >= 0.0) inSST = numpy.all(theComponents >= 0.0)
if color: # have to return color array if color: # have to return color array
@ -733,8 +737,8 @@ class Orientation:
# ****************************************************************************************** # ******************************************************************************************
__slots__ = ['quaternion','symmetry'] __slots__ = ['quaternion','symmetry']
def __init__(self, def __init__(self,
quaternion = Quaternion.fromIdentity(), quaternion = Quaternion.fromIdentity(),
Rodrigues = None, Rodrigues = None,
angleAxis = None, angleAxis = None,
@ -826,7 +830,7 @@ class Orientation:
''' '''
axis rotated according to orientation (using crystal symmetry to ensure location falls into SST) axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)
''' '''
for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent orientations for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent orientations
if SST: # pole requested to be within SST if SST: # pole requested to be within SST
pole = q.conjugated()*axis # align crystal direction to axis pole = q.conjugated()*axis # align crystal direction to axis
@ -849,3 +853,23 @@ class Orientation:
if inSST: break if inSST: break
return color return color
@classmethod
def getAverageOrientation(cls, orientationList):
"""RETURN THE 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.getAverageOrientation([a,b])"""
if not all(isinstance(item, Orientation) for item in orientationList):
raise TypeError("Only instances of Orientation can be averaged.")
n = len(orientationList)
tmp_m = orientationList.pop(0).quaternion.asM()
for tmp_o in orientationList:
tmp_m += tmp_o.quaternion.asM()
eig, vec = numpy.linalg.eig(tmp_m/n)
return Orientation( quaternion=Quaternion(quatArray=vec.T[eig.argmax()]) )