undoes commit 6c3b5f17eb

tests fail, needs closer look. Changes incorporated into
10-consistent-orientation-conversions-in-the-damask-core-and-the-python-module

lib/damask/orientation.py

@@ -27,22 +27,15 @@ class Rodrigues:
 
 # ******************************************************************************************
 class Quaternion:
-    u"""
+    """
     Orientation represented as unit quaternion.
     
-    All methods and naming conventions based on Rowenhorst_etal2015
+    All methods and naming conventions based on http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions.
-    Convention 1: coordinate frames are right-handed
-    Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
-                  when viewing from the end point of the rotation axis unit vector towards the origin
-    Convention 3: rotations will be interpreted in the passive sense
-    Convention 4: Euler angle triplets are implemented using the Bunge convention,
-                  with the angular ranges as [0, 2π],[0, π],[0, 2π]
-    Convention 5: the rotation angle ω is limited to the interval [0, π]
 
     w is the real part, (x, y, z) are the imaginary parts.
-
+    Representation of rotation is in ACTIVE form!
-    Vector "a" (defined in coordinate system "A") is passively rotated
+    (Derived directly or through angleAxis, Euler angles, or active matrix)
-               resulting in new coordinates "b" when expressed in system "B".
+    Vector "a" (defined in coordinate system "A") is actively rotated to new coordinates "b".
     b = Q * a
     b = np.dot(Q.asMatrix(),a)
     """
@@ -316,12 +309,10 @@ class Quaternion:
       return np.outer([i for i in self],[i for i in self])
 
     def asMatrix(self):
-      qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2)
+      return np.array(
-      return 2.0*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)],
-        [[      qbarhalf + self.x**2    , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y],
+         [    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)],
-         [ self.x*self.y + self.w*self.z,      qbarhalf + self.y**2    , self.y*self.z - self.w*self.x],
+         [    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)]])
-         [ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x,      qbarhalf + self.z**2    ],
-        ])
 
     def asAngleAxis(self,
                     degrees = False):
@@ -344,23 +335,52 @@ class Quaternion:
       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,
-                 degrees = False):
+                 type = "bunge",
-      """Orientation as Bunge-Euler angles."""
+                 degrees = False,
-      q03 = self.w**2+self.z**2
+                 standardRange = False):
-      q12 = self.x**2+self.y**2
+      """
-      chi = np.sqrt(q03*q12)
+      Orientation as Bunge-Euler angles.
-      
-      if abs(chi) < 1e-10 and abs(q12) < 1e-10:
-        eulers = np.array([math.atan2(-2*self.w*self.z,self.w**2-self.z**2),0,0])
-      elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
-        eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0])
-      else:
-        eulers = np.array([math.atan2((self.x*self.z-self.w*self.y)/chi,(-self.w*self.x-self.y*self.z)/chi),
-                           math.atan2(2*chi,q03-q12),
-                           math.atan2((self.w*self.y+self.x*self.z)/chi,( self.y*self.z-self.w*self.x)/chi),
-                          ])
 
-      return np.degrees(eulers) if degrees else eulers
+      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
@@ -388,7 +408,7 @@ class Quaternion:
       halfangle = math.atan(np.linalg.norm(rodrigues))
       c = math.cos(halfangle)
       w = c
-      x,y,z = rodrigues/c
+      x,y,z = c*rodrigues
       return cls([w,x,y,z])
 
 
@@ -411,19 +431,24 @@ class Quaternion:
     @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
 
-      sigma = 0.5*(eulers[0]+eulers[2])
+      c = np.cos(0.5 * eulers)
-      delta = 0.5*(eulers[0]-eulers[2])
+      s = np.sin(0.5 * eulers)
-      c = np.cos(0.5*eulers[1])
-      s = np.sin(0.5*eulers[1])
 
-      w =  c * np.cos(sigma) 
+      if type.lower() == 'bunge' or type.lower() == 'zxz':
-      x = -s * np.cos(delta) 
+        w =   c[0] * c[1] * c[2] - s[0] * c[1] * s[2]
-      y = -s * np.sin(delta)
+        x =   c[0] * s[1] * c[2] + s[0] * s[1] * s[2]
-      z = -c * np.sin(sigma) 
+        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])
 
 
@@ -435,16 +460,49 @@ class Quaternion:
       if m.shape != (3,3) and np.prod(m.shape) == 9:
         m = m.reshape(3,3)
 
-      w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
+      tr = np.trace(m)
-      x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
+      if tr > 1e-8:
-      y = 0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2])
+        s = math.sqrt(tr + 1.0)*2.0
-      z = 0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2])
 
-      x *= -1 if m[2,1] < m[1,2] else 1
+        return cls(
-      y *= -1 if m[0,2] < m[2,0] else 1
+          [ s*0.25,
-      z *= -1 if m[1,0] < m[0,1] else 1
+            (m[2,1] - m[1,2])/s,
+            (m[0,2] - m[2,0])/s,
+            (m[1,0] - m[0,1])/s,
+          ])
 
-      return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**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
+        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
@@ -771,7 +829,7 @@ class Orientation:
       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,degrees=degrees)
+      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
@@ -797,15 +855,16 @@ class Orientation:
     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(degrees=True))) )
+           '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(degrees)
+    return self.quaternion.asEulers(type, degrees, standardRange)
   eulers = property(asEulers)
 
   def asRodrigues(self):