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

transition to new orientation class 

forward-backward conversion quite stable
This commit is contained in:
Martin Diehl 2019-02-23 21:47:16 +01:00
parent affe65eb55
commit b3455c825e
5 changed files with 57 additions and 320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 35bed9722ddecc342719bafac32590e9ab94d053 Subproject commit f0090997df817f0a0b5a480a60e81929875b1010

6
processing/post/addOrientations.py
View File

@ -109,8 +109,8 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
crystalrotation = np.array(options.crystalrotation[1:4] + (options.crystalrotation[0],)) # Compatibility hack crystalrotation = np.array(options.crystalrotation[1:4] + (options.crystalrotation[0],)) # Compatibility hack
labrotation = np.array(options.labrotation[1:4] + (options.labrotation[0],)) # Compatibility hack labrotation = np.array(options.labrotation[1:4] + (options.labrotation[0],)) # Compatibility hack
r = damask.Rotation.fromAngleAxis(crystalrotation,options.degrees) # crystal frame rotation r = damask.Rotation.fromAxisAngle(crystalrotation,options.degrees) # crystal frame rotation
R = damask.Rotation.fromAngleAxis(labrotation,options.degrees) # lab frame rotation R = damask.Rotation.fromAxisAngle(labrotation,options.degrees) # lab frame rotation
# --- loop over input files ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
@ -183,7 +183,7 @@ for name in filenames:
elif output == 'rodrigues': table.data_append(o.asRodrigues()) elif output == 'rodrigues': table.data_append(o.asRodrigues())
elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees)) elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees))
elif output == 'matrix': table.data_append(o.asMatrix()) elif output == 'matrix': table.data_append(o.asMatrix())
elif output == 'angleaxis': table.data_append(o.asAngleAxis(degrees=options.degrees)) elif output == 'angleaxis': table.data_append(o.asAxisAngle(degrees=options.degrees))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

2
processing/post/rotateData.py
View File

@ -41,7 +41,7 @@ if options.data is None:
parser.error('no data column specified.') parser.error('no data column specified.')
rotation = np.array(options.rotation[1:4]+(options.rotation[0],)) # Compatibility hack rotation = np.array(options.rotation[1:4]+(options.rotation[0],)) # Compatibility hack
r = damask.Rotation.fromAngleAxis(rotation,options.degrees,normalise=True) r = damask.Rotation.fromAxisAngle(rotation,options.degrees,normalise=True)
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------

2
processing/pre/geom_addPrimitive.py
View File

@ -64,7 +64,7 @@ if options.dimension is None:
parser.error('no dimension specified.') parser.error('no dimension specified.')
if options.angleaxis is not None: if options.angleaxis is not None:
ax = np.array(options.angleaxis[1:4] + (options.angleaxis[0],)) # Compatibility hack ax = np.array(options.angleaxis[1:4] + (options.angleaxis[0],)) # Compatibility hack
rotation = damask.Rotation.fromAngleAxis(ax,options.degrees,normalise=True) rotation = damask.Rotation.fromAxisAngle(ax,options.degrees,normalise=True)
elif options.quaternion is not None: elif options.quaternion is not None:
rotation = damask.Rotation.fromQuaternion(options.quaternion) rotation = damask.Rotation.fromQuaternion(options.quaternion)
else: else:

363
python/damask/orientation.py
View File

@ -261,31 +261,33 @@ class Rotation:
# convert to different orientation representations (numpy arrays) # convert to different orientation representations (numpy arrays)
def asQuaternion(self): def asQuaternion(self):
"""Unit quaternion: (q, [p_1, p_2, p_3])"""
return self.quaternion.asArray() return self.quaternion.asArray()
def asEulers(self, def asEulers(self,
degrees = False): degrees = False):
"""Bunge-Euler angles: (φ_1, ϕ, φ_2)"""
eu = qu2eu(self.quaternion.asArray()) eu = qu2eu(self.quaternion.asArray())
if degrees: eu = np.degrees(eu) if degrees: eu = np.degrees(eu)
return eu return eu
def asAngleAxis(self, def asAxisAngle(self,
degrees = False): degrees = False):
"""Axis-angle pair: ([n_1, n_2, n_3], ω)"""
ax = qu2ax(self.quaternion.asArray()) ax = qu2ax(self.quaternion.asArray())
if degrees: ax[3] = np.degrees(ax[3]) if degrees: ax[3] = np.degrees(ax[3])
return ax return ax
def asMatrix(self): def asMatrix(self):
"""Rotation matrix"""
return qu2om(self.quaternion.asArray()) return qu2om(self.quaternion.asArray())
def asRodrigues(self): def asRodrigues(self):
"""Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2))"""
return qu2ro(self.quaternion.asArray()) return qu2ro(self.quaternion.asArray())
def asHomochoric(self): def asHomochoric(self):
"""Homochoric vector: (h_1, h_2, h_3)"""
return qu2ho(self.quaternion.asArray()) return qu2ho(self.quaternion.asArray())
def asCubochoric(self): def asCubochoric(self):
@ -322,7 +324,7 @@ class Rotation:
return cls(eu2qu(eu)) return cls(eu2qu(eu))
@classmethod @classmethod
def fromAngleAxis(cls, def fromAxisAngle(cls,
angleAxis, angleAxis,
degrees = False, degrees = False,
normalise = False, normalise = False,
@ -373,6 +375,27 @@ class Rotation:
return cls(ro2qu(ro)) return cls(ro2qu(ro))
@classmethod
def fromHomochoric(cls,
homochoric,
P = -1):
ho = homochoric if isinstance(homochoric, np.ndarray) else np.array(homochoric)
if P > 0: ho *= -1 # convert from P=1 to P=-1
return cls(ho2qu(ho))
@classmethod
def fromCubochoric(cls,
cubochoric,
P = -1):
cu = cubochoric if isinstance(cubochoric, np.ndarray) else np.array(cubochoric)
ho = cu2ho(cu)
if P > 0: ho *= -1 # convert from P=1 to P=-1
return cls(ho2qu(ho))
def __mul__(self, other): def __mul__(self, other):
""" """
@ -672,28 +695,6 @@ class Quaternion:
def asRodrigues(self): def asRodrigues(self):
return np.inf*np.ones(3) if np.isclose(self.q,0.0) else self.p/self.q return np.inf*np.ones(3) if np.isclose(self.q,0.0) else self.p/self.q
def asEulers(self,
degrees = False):
"""Orientation as Bunge-Euler angles."""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
q03 = self.q**2 + self.p[2]**2
q12 = self.p[0]**2 + self.p[1]**2
chi = np.sqrt(q03*q12)
if np.isclose(chi,0.0) and np.isclose(q12,0.0):
eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0])
elif np.isclose(chi,0.0) and np.isclose(q03,0.0):
eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0])
else:
eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi),
math.atan2(2*chi,q03-q12),
math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi),
])
eulers %= 2.0*math.pi # enforce positive angles
return np.degrees(eulers) if degrees else eulers
# # Static constructors # # Static constructors
@classmethod @classmethod
@ -1506,18 +1507,8 @@ class Orientation:
'Symmetry: {}'.format(self.symmetry), 'Symmetry: {}'.format(self.symmetry),
'Quaternion: {}'.format(self.quaternion), 'Quaternion: {}'.format(self.quaternion),
'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ), 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
]) ])
def asQuaternion(self):
return self.quaternion.asList()
def asEulers(self,
degrees = False,
):
return self.quaternion.asEulers(degrees)
eulers = property(asEulers)
def asRodrigues(self): def asRodrigues(self):
return self.quaternion.asRodrigues() return self.quaternion.asRodrigues()
rodrigues = property(asRodrigues) rodrigues = property(asRodrigues)
@ -1526,7 +1517,6 @@ class Orientation:
degrees = False, degrees = False,
flat = False): flat = False):
return self.quaternion.asAngleAxis(degrees,flat) return self.quaternion.asAngleAxis(degrees,flat)
angleAxis = property(asAngleAxis)
def asMatrix(self): def asMatrix(self):
return self.quaternion.asMatrix() return self.quaternion.asMatrix()
@ -1643,251 +1633,6 @@ class Orientation:
symmetry = reference.symmetry.lattice) 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
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)
#################################################################################################### ####################################################################################################
# Code below available according to the followin conditions on https://github.com/MarDiehl/3Drotations # Code below available according to the followin conditions on https://github.com/MarDiehl/3Drotations
#################################################################################################### ####################################################################################################
@ -1920,10 +1665,10 @@ class Orientation:
#################################################################################################### ####################################################################################################
def isone(a): def isone(a):
return np.isclose(a,1.0,atol=1.0e-15,rtol=0.0) return np.isclose(a,1.0,atol=1.0e-7,rtol=0.0)
def iszero(a): def iszero(a):
return np.isclose(a,0.0,atol=1.0e-300,rtol=0.0) return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
def eu2om(eu): def eu2om(eu):
@ -2063,8 +1808,7 @@ def ho2ax(ho):
for i in range(2,16): for i in range(2,16):
hm *= hmag_squared hm *= hmag_squared
s += tfit[i] * hm s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(s)) # ToDo: Check sanity check in reference implementation ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
return ax return ax
@ -2155,32 +1899,16 @@ def qu2om(qu):
return om if P > 0.0 else om.T return om if P > 0.0 else om.T
def om2qu(om):
"""Orientation matrix to quaternion"""
s = [+om[0,0] +om[1,1] +om[2,2] +1.0,
+om[0,0] -om[1,1] -om[2,2] +1.0,
-om[0,0] +om[1,1] -om[2,2] +1.0,
-om[0,0] -om[1,1] +om[2,2] +1.0]
s = np.maximum(np.zeros(4),s)
qu = np.sqrt(s)*0.5*np.array([1.0,P,P,P])
# verify the signs (q0 always positive)
#ToDo: Here I donot understand the original shortcut from paper to implementation
qu /= np.linalg.norm(qu)
if any(isone(abs(qu))): qu[np.where(np.logical_not(isone(qu)))] = 0.0
if om[2,1] < om[1,2]: qu[1] *= -1.0
if om[0,2] < om[2,0]: qu[2] *= -1.0
if om[1,0] < om[0,1]: qu[3] *= -1.0
if any(om2ax(om)[0:3]*qu[1:4] < 0.0): print('sign problem',om2ax(om),qu) # something is wrong here
return qu
def qu2ax(qu): def qu2ax(qu):
"""Quaternion to axis angle""" """
omega = 2.0 * np.arccos(qu[0]) Quaternion to axis angle
if iszero(omega): # return axis as [001] if the angle is zero
Modified version of the original formulation, should be numerically more stable
"""
if isone(abs(qu[0])): # set axis to [001] if the angle is 0/360
ax = [ 0.0, 0.0, 1.0, 0.0 ] ax = [ 0.0, 0.0, 1.0, 0.0 ]
elif not iszero(qu[0]): elif not iszero(qu[0]):
omega = 2.0 * np.arccos(qu[0])
s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2)
ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ]
else: else:
@ -2196,14 +1924,14 @@ def qu2ro(qu):
else: else:
s = np.linalg.norm([qu[1],qu[2],qu[3]]) s = np.linalg.norm([qu[1],qu[2],qu[3]])
ro = [0.0,0.0,P,0.0] if iszero(s) else \ ro = [0.0,0.0,P,0.0] if iszero(s) else \
[ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(qu[0]))] [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] # avoid numerical difficulties
return np.array(ro) return np.array(ro)
def qu2ho(qu): def qu2ho(qu):
"""Quaternion to homochoric""" """Quaternion to homochoric"""
omega = 2.0 * np.arccos(qu[0]) omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties
if iszero(omega): if iszero(omega):
ho = np.array([ 0.0, 0.0, 0.0 ]) ho = np.array([ 0.0, 0.0, 0.0 ])
@ -2290,6 +2018,15 @@ def om2cu(om):
return ho2cu(om2ho(om)) return ho2cu(om2ho(om))
def om2qu(om):
"""
Orientation matrix to quaternion
The original formulation (direct conversion) had numerical issues
"""
return ax2qu(om2ax(om))
def ax2cu(ax): def ax2cu(ax):
"""Axis angle to cubochoric""" """Axis angle to cubochoric"""
return ho2cu(ax2ho(ax)) return ho2cu(ax2ho(ax))