multiple fixes and improvements.
1) asAngleAxis got “degrees” switch, fixed buggy output for negative angles, now angles are always non-negative. 2) removed rotateBys. 3) asEulers corrected for cases w==z or x==y. 4) added method symmetryQuats. 5) calculation of disorientation (finally) fixed, returns tuple (disorientation quaternion, index symA, index symB, boolean whether A—>B or B—>A).
This commit is contained in:
Philip Eisenlohr
parent
858fe7e897
commit
636eb8a087
|
|
@ -34,12 +34,19 @@ class Quaternion:
|
|||
|
||||
# w is the real part, (x, y, z) are the imaginary parts
|
||||
|
||||
def __init__(self, quatArray=[1.0,0.0,0.0,0.0]):
|
||||
# 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]):
|
||||
self.w, \
|
||||
self.x, \
|
||||
self.y, \
|
||||
self.z = quatArray
|
||||
self = self.homomorph()
|
||||
self.homomorph()
|
||||
|
||||
def __iter__(self):
|
||||
return iter([self.w,self.x,self.y,self.z])
|
||||
|
|
@ -51,7 +58,7 @@ class Quaternion:
|
|||
copy = __copy__
|
||||
|
||||
def __repr__(self):
|
||||
return 'Quaternion(real=%+.4f, imag=<%+.4f, %+.4f, %+.4f>)' % \
|
||||
return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \
|
||||
(self.w, self.x, self.y, self.z)
|
||||
|
||||
def __pow__(self, exponent):
|
||||
|
|
@ -237,18 +244,6 @@ class Quaternion:
|
|||
self.z = 0.
|
||||
return self
|
||||
|
||||
def rotateBy_angleaxis(self, angle, axis):
|
||||
self *= Quaternion.fromAngleAxis(angle, axis)
|
||||
return self
|
||||
|
||||
def rotateBy_Eulers(self, eulers):
|
||||
self *= Quaternion.fromEulers(eulers, type)
|
||||
return self
|
||||
|
||||
def rotateBy_matrix(self, m):
|
||||
self *= Quaternion.fromMatrix(m)
|
||||
return self
|
||||
|
||||
def normalize(self):
|
||||
d = self.magnitude()
|
||||
if d > 0.0:
|
||||
|
|
@ -299,7 +294,8 @@ class Quaternion:
|
|||
[ 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):
|
||||
def asAngleAxis(self,
|
||||
degrees = False):
|
||||
if self.w > 1:
|
||||
self.normalize()
|
||||
|
||||
|
|
@ -308,18 +304,22 @@ class Quaternion:
|
|||
y = 2*self.w * s
|
||||
|
||||
angle = math.atan2(y,x)
|
||||
if angle < 0.0:
|
||||
angle *= -1.
|
||||
s *= -1.
|
||||
|
||||
return angle, np.array([1.0, 0.0, 0.0] if angle < 1e-3 else [self.x / s, self.y / s, self.z / s])
|
||||
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):
|
||||
if self.w != 0.0:
|
||||
return np.array([self.x, self.y, self.z])/self.w
|
||||
else:
|
||||
return np.array([float('inf')]*3)
|
||||
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):
|
||||
def asEulers(self,
|
||||
type = 'bunge',
|
||||
degrees = False,
|
||||
standardRange = False):
|
||||
'''
|
||||
conversion taken from:
|
||||
conversion of ACTIVE rotation to Euler angles taken from:
|
||||
Melcher, A.; Unser, A.; Reichhardt, M.; Nestler, B.; Pötschke, 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
|
||||
|
|
@ -327,11 +327,11 @@ class Quaternion:
|
|||
angles = [0.0,0.0,0.0]
|
||||
|
||||
if type.lower() == 'bunge' or type.lower() == 'zxz':
|
||||
if abs(self.x - self.y) < 1e-8:
|
||||
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 - self.z) < 1e-8:
|
||||
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)
|
||||
|
|
@ -394,8 +394,8 @@ class Quaternion:
|
|||
|
||||
@classmethod
|
||||
def fromAngleAxis(cls, angle, axis):
|
||||
if not isinstance(axis, np.ndarray): axis = np.array(axis)
|
||||
axis /= np.linalg.norm(axis)
|
||||
if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d')
|
||||
axis = axis.astype(float)/np.linalg.norm(axis)
|
||||
s = math.sin(angle / 2.0)
|
||||
w = math.cos(angle / 2.0)
|
||||
x = axis[0] * s
|
||||
|
|
@ -437,8 +437,8 @@ class Quaternion:
|
|||
if m.shape != (3,3) and np.prod(m.shape) == 9:
|
||||
m = m.reshape(3,3)
|
||||
|
||||
tr=m[0,0]+m[1,1]+m[2,2]
|
||||
if tr > 0.00000001:
|
||||
tr = np.trace(m)
|
||||
if tr > 1e-8:
|
||||
s = math.sqrt(tr + 1.0)*2.0
|
||||
|
||||
return cls(
|
||||
|
|
@ -555,9 +555,9 @@ class Symmetry:
|
|||
def __cmp__(self,other):
|
||||
return cmp(Symmetry.lattices.index(self.lattice),Symmetry.lattices.index(other.lattice))
|
||||
|
||||
def equivalentQuaternions(self,quaternion):
|
||||
def symmetryQuats(self):
|
||||
'''
|
||||
List of symmetrically equivalent quaternions based on own symmetry.
|
||||
List of symmetry operations as quaternions.
|
||||
'''
|
||||
if self.lattice == 'cubic':
|
||||
symQuats = [
|
||||
|
|
@ -589,17 +589,17 @@ class Symmetry:
|
|||
elif self.lattice == 'hexagonal':
|
||||
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.5*math.sqrt(3), 0.0, 0.0, 0.5 ],
|
||||
[-0.5*math.sqrt(3), 0.0, 0.0,-0.5 ],
|
||||
[ 0.0, 0.5*math.sqrt(3), 0.5, 0.0 ],
|
||||
[ 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.5, 0.0, 0.0, 0.5*math.sqrt(3) ],
|
||||
[-0.5, 0.0, 0.0, 0.5*math.sqrt(3) ],
|
||||
[ 0.0, 0.5*math.sqrt(3), 0.5, 0.0 ],
|
||||
]
|
||||
elif self.lattice == 'tetragonal':
|
||||
symQuats = [
|
||||
|
|
@ -624,7 +624,14 @@ class Symmetry:
|
|||
[ 1.0,0.0,0.0,0.0 ],
|
||||
]
|
||||
|
||||
return [quaternion*Quaternion(q) for q in symQuats]
|
||||
return map(Quaternion,symQuats)
|
||||
|
||||
|
||||
def equivalentQuaternions(self,quaternion):
|
||||
'''
|
||||
List of symmetrically equivalent quaternions based on own symmetry.
|
||||
'''
|
||||
return [quaternion*Quaternion(q) for q in self.symmetryQuats()]
|
||||
|
||||
|
||||
def inFZ(self,R):
|
||||
|
|
@ -663,24 +670,25 @@ class Symmetry:
|
|||
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 and self.inFZ(R)
|
||||
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 and self.inFZ(R)
|
||||
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 and self.inFZ(R)
|
||||
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 and self.inFZ(R)
|
||||
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon
|
||||
|
||||
else:
|
||||
return True
|
||||
|
||||
|
||||
def inSST(self,vector,color = False):
|
||||
def inSST(self,
|
||||
vector,
|
||||
color = False):
|
||||
'''
|
||||
Check whether given vector falls into standard stereographic triangle of own symmetry.
|
||||
Return inverse pole figure color if requested.
|
||||
|
|
@ -720,7 +728,9 @@ class Symmetry:
|
|||
if np.all(basis == 0.0):
|
||||
theComponents = -np.ones(3,'d')
|
||||
else:
|
||||
theComponents = np.dot(basis,np.array([vector[0],vector[1],abs(vector[2])]))
|
||||
v = np.array(vector,dtype = float)
|
||||
v[2] = abs(v[2]) # z component projects identical for positive and negative values
|
||||
theComponents = np.dot(basis,v)
|
||||
|
||||
inSST = np.all(theComponents >= 0.0)
|
||||
|
||||
|
|
@ -799,8 +809,9 @@ class Orientation:
|
|||
return self.quaternion.asRodrigues()
|
||||
rodrigues = property(asRodrigues)
|
||||
|
||||
def asAngleAxis(self):
|
||||
return self.quaternion.asAngleAxis()
|
||||
def asAngleAxis(self,
|
||||
degrees = False):
|
||||
return self.quaternion.asAngleAxis(degrees)
|
||||
angleAxis = property(asAngleAxis)
|
||||
|
||||
def asMatrix(self):
|
||||
|
|
@ -818,7 +829,7 @@ class Orientation:
|
|||
def equivalentOrientations(self):
|
||||
return map(lambda q: Orientation(quaternion = q, symmetry = self.symmetry.lattice),
|
||||
self.equivalentQuaternions())
|
||||
equiOrientation = property(equivalentQuaternions)
|
||||
equiOrientations = property(equivalentQuaternions)
|
||||
|
||||
def reduced(self):
|
||||
'''
|
||||
|
|
@ -834,22 +845,28 @@ class Orientation:
|
|||
def disorientation(self,other):
|
||||
'''
|
||||
Disorientation between myself and given other orientation
|
||||
(either reduced according to my own symmetry or given one)
|
||||
(currently needs to be of same symmetry.
|
||||
look into A. Heinz and P. Neumann 1991 for cases with differing sym.)
|
||||
'''
|
||||
|
||||
lowerSymmetry = min(self.symmetry,other.symmetry)
|
||||
breaker = False
|
||||
if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.')
|
||||
|
||||
for me in self.symmetry.equivalentQuaternions(self.quaternion):
|
||||
me.conjugate()
|
||||
for they in other.symmetry.equivalentQuaternions(other.quaternion):
|
||||
theQ = they * me
|
||||
breaker = lowerSymmetry.inDisorientationSST(theQ.asRodrigues()) #\
|
||||
# or lowerSymmetry.inDisorientationSST(theQ.conjugated().asRodrigues())
|
||||
misQ = self.quaternion.conjugated()*other.quaternion
|
||||
|
||||
for i,sA in enumerate(self.symmetry.symmetryQuats()):
|
||||
for j,sB in enumerate(other.symmetry.symmetryQuats()):
|
||||
theQ = sA.conjugated()*misQ*sB
|
||||
for k in xrange(2):
|
||||
theQ.conjugate()
|
||||
hitSST = other.symmetry.inDisorientationSST(theQ)
|
||||
hitFZ = self.symmetry.inFZ(theQ)
|
||||
breaker = hitSST and hitFZ
|
||||
if breaker: break
|
||||
if breaker: break
|
||||
if breaker: break
|
||||
|
||||
return Orientation(quaternion=theQ,symmetry=self.symmetry.lattice) #, me.conjugated(), they
|
||||
return (Orientation(quaternion=theQ,symmetry=self.symmetry.lattice),
|
||||
i,j,k == 1) # disorientation, own sym, other sym, self-->other: True, self<--other: False
|
||||
|
||||
|
||||
def inversePole(self,axis,SST = True):
|
||||
|
|
|
|||
Loading…
Reference in New Issue