From f0bb50b97d6cd9897cdeb261a14158164f9d0c30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 2 May 2020 15:16:26 +0200 Subject: [PATCH 01/32] vectorized rotation function --- python/damask/_rotation.py | 42 ++++++++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 716977a79..f442561cf 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -83,7 +83,7 @@ class Rotation: Todo ---- Document details active/passive) - considere rotation of (3,3,3,3)-matrix + consider rotation of (3,3,3,3)-matrix """ if self.quaternion.shape != (4,): @@ -99,9 +99,7 @@ class Rotation: elif isinstance(other, (tuple,np.ndarray)): if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:]) - B = 2.0 * ( self.quaternion[1]*other[0] - + self.quaternion[2]*other[1] - + self.quaternion[3]*other[2]) + B = 2.0 * np.dot(self.quaternion[1:],other) C = 2.0 * _P*self.quaternion[0] return np.array([ @@ -119,6 +117,42 @@ class Rotation: return NotImplemented + def __matmul__(self, other): + """ + Rotation. + + details to be discussed + """ + shape = self.quaternion.shape[:-1] + + if isinstance(other, Rotation): # rotate a rotation + q_m = self.quaternion[...,0].reshape(shape+(1,)) + p_m = self.quaternion[...,1:] + q_o = other.quaternion[...,0].reshape(shape+(1,)) + p_o = other.quaternion[...,1:] + q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(shape+(1,))) + p = q_m*p_m + q_o*p_m + _P * np.cross(p_m,p_o) + return self.__class__(np.block([q,p])).standardize() + + elif isinstance(other,np.ndarray): + if shape + (3,) == other.shape: + q_m = self.quaternion[...,0] + p_m = self.quaternion[...,1:] + A = q_m**2.0 - np.einsum('...i,...i',p_m,p_m) + B = 2.0 * np.einsum('...i,...i',p_m,p_m) + C = 2.0 * _P * q_m + return np.block([(A * other[...,i]).reshape(shape+(1,)) + + (B * p_m[...,i]).reshape(shape+(1,)) + + (C * ( p_m[...,(i+1)%3]*other[...,(i+2)%3]\ + - p_m[...,(i+2)%3]*other[...,(i+1)%3])).reshape(shape+(1,)) + for i in [0,1,2]]) + if shape + (3,3) == other.shape: + R = self.asMatrix() + return np.einsum('...im,...jn,...mn',R,R,other) + if shape + (3,3,3,3) == other.shape: + R = self.asMatrix() + return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) + def inverse(self): """In-place inverse rotation/backward rotation.""" self.quaternion[...,1:] *= -1 From ef4a4dad4af42e78f29a95cc18bcdadda6514d96 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 2 May 2020 15:50:46 +0200 Subject: [PATCH 02/32] shape property and numpy-like broadcasting this makes it easy to apply a single rotation to a field --- python/damask/_rotation.py | 39 +++++++++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index f442561cf..cb5b2fbd1 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -53,6 +53,12 @@ class Rotation: """ self.quaternion = quaternion.copy() + + @property + def shape(self): + return self.quaternion.shape[:-1] + + def __copy__(self): """Copy.""" return self.__class__(self.quaternion) @@ -123,35 +129,35 @@ class Rotation: details to be discussed """ - shape = self.quaternion.shape[:-1] - if isinstance(other, Rotation): # rotate a rotation - q_m = self.quaternion[...,0].reshape(shape+(1,)) + q_m = self.quaternion[...,0].reshape(self.shape+(1,)) p_m = self.quaternion[...,1:] - q_o = other.quaternion[...,0].reshape(shape+(1,)) + q_o = other.quaternion[...,0].reshape(self.shape+(1,)) p_o = other.quaternion[...,1:] - q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(shape+(1,))) + q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,))) p = q_m*p_m + q_o*p_m + _P * np.cross(p_m,p_o) return self.__class__(np.block([q,p])).standardize() elif isinstance(other,np.ndarray): - if shape + (3,) == other.shape: + if self.shape + (3,) == other.shape: q_m = self.quaternion[...,0] p_m = self.quaternion[...,1:] A = q_m**2.0 - np.einsum('...i,...i',p_m,p_m) B = 2.0 * np.einsum('...i,...i',p_m,p_m) C = 2.0 * _P * q_m - return np.block([(A * other[...,i]).reshape(shape+(1,)) + - (B * p_m[...,i]).reshape(shape+(1,)) + + return np.block([(A * other[...,i]).reshape(self.shape+(1,)) + + (B * p_m[...,i]).reshape(self.shape+(1,)) + (C * ( p_m[...,(i+1)%3]*other[...,(i+2)%3]\ - - p_m[...,(i+2)%3]*other[...,(i+1)%3])).reshape(shape+(1,)) + - p_m[...,(i+2)%3]*other[...,(i+1)%3])).reshape(self.shape+(1,)) for i in [0,1,2]]) - if shape + (3,3) == other.shape: + if self.shape + (3,3) == other.shape: R = self.asMatrix() return np.einsum('...im,...jn,...mn',R,R,other) - if shape + (3,3,3,3) == other.shape: + if self.shape + (3,3,3,3) == other.shape: R = self.asMatrix() return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) + else: + raise ValueError def inverse(self): """In-place inverse rotation/backward rotation.""" @@ -186,6 +192,17 @@ class Rotation: return other*self.inversed() + def broadcast_to(self,shape): + if self.shape == (): + q = np.broadcast_to(self.quaternion,shape+(4,)) + else: + q = np.block([np.broadcast_to(self.quaternion[...,0:1],shape+(1,)), + np.broadcast_to(self.quaternion[...,1:2],shape+(1,)), + np.broadcast_to(self.quaternion[...,2:3],shape+(1,)), + np.broadcast_to(self.quaternion[...,3:4],shape+(1,))]) + return self.__class__(q) + + def average(self,other): """ Calculate the average rotation. From 2dc46b783aa80812dc44957886a06f61c50660e7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 00:58:40 +0200 Subject: [PATCH 03/32] simplified and tested --- python/damask/_rotation.py | 7 ++----- python/tests/test_Rotation.py | 7 +++++++ 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 059f148a4..38294f75e 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -112,11 +112,8 @@ class Rotation: B = 2.0 * np.dot(self.quaternion[1:],other) C = 2.0 * _P*self.quaternion[0] - return np.array([ - A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]), - A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]), - A*other[2] + B*self.quaternion[3] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]), - ]) + return A*other + B*self.quaternion[1:] + C * np.cross(self.quaternion[1:],other) + elif other.shape == (3,3,): # rotate a single (3x3)-matrix return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) elif other.shape == (3,3,3,3,): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index b7442035f..ffa1548d7 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -300,3 +300,10 @@ class TestRotation: f = Rotation._get_pyramid_order(a,'forward') b = Rotation._get_pyramid_order(a,'backward') assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a) + + + @pytest.mark.parametrize('x',[np.random.rand(3), + np.random.rand(3,3)]) + def test_rotation_identity(self,x): + R = Rotation() + assert np.allclose(x,R*x) From e2ba294b751b3606821ad6da9b4e96570e11e79d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 09:40:32 +0200 Subject: [PATCH 04/32] bugfix: wrong variable --- python/damask/_rotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 38294f75e..a9f11c6c2 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -131,12 +131,12 @@ class Rotation: details to be discussed """ if isinstance(other, Rotation): # rotate a rotation - q_m = self.quaternion[...,0].reshape(self.shape+(1,)) + q_m = self.quaternion[...,0:1] p_m = self.quaternion[...,1:] - q_o = other.quaternion[...,0].reshape(self.shape+(1,)) + q_o = other.quaternion[...,0:1] p_o = other.quaternion[...,1:] q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,))) - p = q_m*p_m + q_o*p_m + _P * np.cross(p_m,p_o) + p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o) return self.__class__(np.block([q,p])).standardize() elif isinstance(other,np.ndarray): From 347e3b8c60763f0bae1d7de0cd5496b363003e66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 09:40:44 +0200 Subject: [PATCH 05/32] using new table class + vectorized rotation --- processing/post/addOrientations.py | 133 +++++++++-------------------- processing/post/rotateData.py | 13 ++- 2 files changed, 46 insertions(+), 100 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index e1752b0a4..b6292b18c 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np @@ -24,14 +25,8 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can """, version = scriptID) -representations = { - 'quaternion': ['qu',4], - 'rodrigues': ['ro',4], - 'Rodrigues': ['Ro',3], - 'eulers': ['eu',3], - 'matrix': ['om',9], - 'angleaxis': ['ax',4], - } +representations = ['quaternion', 'rodrigues', 'eulers', 'matrix', 'axisangle'] + parser.add_option('-o', '--output', @@ -93,10 +88,10 @@ parser.set_defaults(output = [], ) (options, filenames) = parser.parse_args() +if filenames == []: filenames = [None] -#options.output = list(map(lambda x: x.lower(), options.output)) if options.output == [] or (not set(options.output).issubset(set(representations))): - parser.error('output must be chosen from {}.'.format(', '.join(representations))) + parser.error('output must be chosen from {}.'.format(', '.join(representations))) input = [options.eulers is not None, options.rodrigues is not None, @@ -109,95 +104,47 @@ input = [options.eulers is not None, if np.sum(input) != 1: parser.error('needs exactly one input format.') -(label,dim,inputtype) = [(options.eulers,representations['eulers'][1],'eulers'), - (options.rodrigues,representations['rodrigues'][1],'rodrigues'), - ([options.x,options.y,options.z],[3,3,3],'frame'), - (options.matrix,representations['matrix'][1],'matrix'), - (options.quaternion,representations['quaternion'][1],'quaternion'), - ][np.where(input)[0][0]] # select input label that was requested - -r = damask.Rotation.fromAxisAngle(np.array(options.crystalrotation),options.degrees,normalise=True) -R = damask.Rotation.fromAxisAngle(np.array(options.labrotation),options.degrees,normalise=True) - - -# --- loop over input files ------------------------------------------------------------------------ - -if filenames == []: filenames = [None] +r = damask.Rotation.from_axis_angle(np.array(options.crystalrotation),options.degrees,normalise=True) +R = damask.Rotation.from_axis_angle(np.array(options.labrotation),options.degrees,normalise=True) for name in filenames: - try: - table = damask.ASCIItable(name = name) - except IOError: - continue - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) -# ------------------------------------------ read header ------------------------------------------ + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - table.head_read() + if options.eulers is not None: + label = options.eulers + print(np.max(table.get(options.eulers),axis=0)) + o = damask.Rotation.from_Eulers(table.get(options.eulers), options.degrees) + elif options.rodrigues is not None: + label = options.rodrigues + o = damask.Rotation.from_Rodrigues(table.get(options.rodrigues)) + elif options.matrix is not None: + label = options.matrix + o = damask.Rotation.from_matrix(table.get(options.matrix).reshape(-1,3,3)) + elif options.x is not None: + label = '<{},{},{}>'.format(options.x,options.y,options.z) + M = np.block([table.get(options.x),table.get(options.y),table.get(options.z)]).reshape(-1,3,3) + o = damask.Rotation.from_matrix(M/np.linalg.norm(M,axis=0)) + elif options.quaternion is not None: + label = options.quaternion + o = damask.Rotation.from_quaternion(table.get(options.quaternion)) -# ------------------------------------------ sanity checks ----------------------------------------- + o = r.broadcast_to(o.shape) @ o @ R.broadcast_to(o.shape) - errors = [] - remarks = [] - - if not np.all(table.label_dimension(label) == dim): errors.append('input {} does not have dimension {}.'.format(label,dim)) - else: column = table.label_index(label) + #if options.lattice is not None: + # o = damask.Orientation(rotation = o,lattice = options.lattice).reduced().rotation - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue -# ------------------------------------------ assemble header --------------------------------------- + if 'rodrigues' in options.output: + table.add('ro({})'.format(label),o.as_rodrigues(), scriptID+' '+' '.join(sys.argv[1:])) + if 'eulers' in options.output: + table.add('eu({})'.format(label),o.as_Eulers(options.degrees), scriptID+' '+' '.join(sys.argv[1:])) + if 'quaternion' in options.output: + table.add('qu({})'.format(label),o.as_quaternion(), scriptID+' '+' '.join(sys.argv[1:])) + if 'matrix' in options.output: + table.add('om({})'.format(label),o.as_matrix(), scriptID+' '+' '.join(sys.argv[1:])) + if 'axisangle' in options.output: + table.add('om({})'.format(label),o.as_axisangle(options.degrees), scriptID+' '+' '.join(sys.argv[1:])) - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - for output in options.output: - if output in representations: - table.labels_append(['{}_{}({})'.format(i+1,representations[output][0],label) \ - for i in range(representations[output][1])]) - table.head_write() - -# ------------------------------------------ process data ------------------------------------------ - - outputAlive = True - while outputAlive and table.data_read(): # read next data line of ASCII table - if inputtype == 'eulers': - d = representations['eulers'][1] - o = damask.Rotation.fromEulers(list(map(float,table.data[column:column+d])),options.degrees) - - elif inputtype == 'rodrigues': - d = representations['rodrigues'][1] - o = damask.Rotation.fromRodrigues(list(map(float,table.data[column:column+d]))) - - elif inputtype == 'matrix': - d = representations['matrix'][1] - o = damask.Rotation.fromMatrix(np.array(list(map(float,table.data[column:column+d]))).reshape(3,3)) - - elif inputtype == 'frame': - M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ - table.data[column[1]:column[1]+3] + \ - table.data[column[2]:column[2]+3]))).reshape(3,3).T - o = damask.Rotation.fromMatrix(M/np.linalg.norm(M,axis=0)) - - elif inputtype == 'quaternion': - d = representations['quaternion'][1] - o = damask.Rotation.fromQuaternion(list(map(float,table.data[column:column+d]))) - - o = r*o*R # apply additional lab and crystal frame rotations - - if options.lattice is not None: - o = damask.Orientation(rotation = o,lattice = options.lattice).reduced().rotation - - for output in options.output: - if output == 'quaternion': table.data_append(o.asQuaternion()) - elif output == 'rodrigues': table.data_append(o.asRodrigues()) - elif output == 'Rodrigues': table.data_append(o.asRodrigues(vector=True)) - elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees)) - elif output == 'matrix': table.data_append(o.asMatrix()) - elif output == 'angleaxis': table.data_append(o.asAxisAngle(degrees=options.degrees)) - outputAlive = table.data_write() # output processed line - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close ASCII tables + table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/rotateData.py b/processing/post/rotateData.py index d270eb79e..ee172c64d 100755 --- a/processing/post/rotateData.py +++ b/processing/post/rotateData.py @@ -37,27 +37,26 @@ parser.add_option('--degrees', parser.set_defaults(rotation = (1.,1.,1.,0), # no rotation about (1,1,1) degrees = False, ) - + (options,filenames) = parser.parse_args() if filenames == []: filenames = [None] if options.data is None: - parser.error('no data column specified.') + parser.error('no data column specified.') -r = damask.Rotation.fromAxisAngle(options.rotation,options.degrees,normalise=True) +r = damask.Rotation.from_axis_angle(options.rotation,options.degrees,normalise=True) for name in filenames: damask.util.report(scriptName,name) - + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) for data in options.data: d = table.get(data) if table.shapes[data] == (9,): d=d.reshape(-1,3,3) - for i,l in enumerate(d): - d[i] = r*l + d = r.broadcast_to(d.shape[0:1]) @ d if table.shapes[data] == (9,): d=d.reshape(-1,9) - + table.set(data,d,scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) From 235f3c0df7df26bc1d42ce50dd3131b2f03edf67 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 11:16:52 +0200 Subject: [PATCH 06/32] bugfix: copy and paste error --- python/damask/_rotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a9f11c6c2..3a046363f 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -52,7 +52,7 @@ class Rotation: Parameters ---------- quaternion : numpy.ndarray, optional - Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check. + Unit quaternion that follows the conventions. Use .from_quaternion to perform a sanity check. """ self.quaternion = quaternion.copy() @@ -144,7 +144,7 @@ class Rotation: q_m = self.quaternion[...,0] p_m = self.quaternion[...,1:] A = q_m**2.0 - np.einsum('...i,...i',p_m,p_m) - B = 2.0 * np.einsum('...i,...i',p_m,p_m) + B = 2.0 * np.einsum('...i,...i',p_m,other) C = 2.0 * _P * q_m return np.block([(A * other[...,i]).reshape(self.shape+(1,)) + (B * p_m[...,i]).reshape(self.shape+(1,)) + @@ -481,7 +481,7 @@ class Rotation: else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) - return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) + return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) @staticmethod def from_random(shape=None): From 9fa153916347d1413facbaf2d2de7f510ca0404f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 11:17:12 +0200 Subject: [PATCH 07/32] using new names --- processing/pre/geom_addPrimitive.py | 4 ++-- processing/pre/geom_fromDREAM3D.py | 4 ++-- processing/pre/geom_rotate.py | 8 ++++---- python/damask/_lattice.py | 2 +- python/damask/_orientation.py | 2 +- python/tests/test_Orientation.py | 6 +++--- 6 files changed, 13 insertions(+), 13 deletions(-) diff --git a/processing/pre/geom_addPrimitive.py b/processing/pre/geom_addPrimitive.py index f33ba27b1..c04fcf565 100755 --- a/processing/pre/geom_addPrimitive.py +++ b/processing/pre/geom_addPrimitive.py @@ -84,9 +84,9 @@ if [options.angleaxis,options.quaternion].count(None) == 0: parser.error('more than one rotation specified.') if options.angleaxis is not None: - rotation = damask.Rotation.fromAxisAngle(np.array(options.angleaxis),options.degrees,normalise=True) + rotation = damask.Rotation.from_axis_angle(np.array(options.angleaxis),options.degrees,normalise=True) elif options.quaternion is not None: - rotation = damask.Rotation.fromQuaternion(options.quaternion) + rotation = damask.Rotation.from_quaternion(options.quaternion) else: rotation = damask.Rotation() diff --git a/processing/pre/geom_fromDREAM3D.py b/processing/pre/geom_fromDREAM3D.py index 7d5b1442d..b000aa811 100755 --- a/processing/pre/geom_fromDREAM3D.py +++ b/processing/pre/geom_fromDREAM3D.py @@ -97,7 +97,7 @@ for name in filenames: dataset = os.path.join(group_pointwise,options.quaternion) try: quats = np.reshape(inFile[dataset][...],(np.product(grid),4)) - rot = [damask.Rotation.fromQuaternion(q,True,P=+1) for q in quats] + rot = [damask.Rotation.from_quaternion(q,True,P=+1) for q in quats] except KeyError: errors.append('Pointwise orientation (quaternion) data ({}) not readable'.format(dataset)) @@ -123,7 +123,7 @@ for name in filenames: dataset = os.path.join(group_average,options.quaternion) try: - rot = [damask.Rotation.fromQuaternion(q,True,P=+1) for q in inFile[dataset][...][1:]] # skip first entry (unindexed) + rot = [damask.Rotation.from_quaternion(q,True,P=+1) for q in inFile[dataset][...][1:]] # skip first entry (unindexed) except KeyError: errors.append('Average orientation data ({}) not readable'.format(dataset)) diff --git a/processing/pre/geom_rotate.py b/processing/pre/geom_rotate.py index 977e00b65..7a774200b 100755 --- a/processing/pre/geom_rotate.py +++ b/processing/pre/geom_rotate.py @@ -59,13 +59,13 @@ if [options.rotation,options.eulers,options.matrix,options.quaternion].count(Non parser.error('no rotation specified.') if options.quaternion is not None: - rot = damask.Rotation.fromQuaternion(np.array(options.quaternion)) # we might need P=+1 here, too... + rot = damask.Rotation.from_quaternion(np.array(options.quaternion)) # we might need P=+1 here, too... if options.rotation is not None: - rot = damask.Rotation.fromAxisAngle(np.array(options.rotation),degrees=options.degrees,normalise=True,P=+1) + rot = damask.Rotation.from_axis_angle(np.array(options.rotation),degrees=options.degrees,normalise=True,P=+1) if options.matrix is not None: - rot = damask.Rotation.fromMatrix(np.array(options.Matrix)) + rot = damask.Rotation.from_matrix(np.array(options.Matrix)) if options.eulers is not None: - rot = damask.Rotation.fromEulers(np.array(options.eulers),degrees=options.degrees) + rot = damask.Rotation.from_Eulers(np.array(options.eulers),degrees=options.degrees) eulers = rot.asEulers(degrees=True) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 42aa0e9bd..5c8166508 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -634,6 +634,6 @@ class Lattice: otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id]) otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]) - r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix.T,myMatrix))) + r['rotations'].append(Rotation.from_matrix(np.dot(otherMatrix.T,myMatrix))) return r diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e6c732007..c145dc42b 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -36,7 +36,7 @@ class Orientation: if isinstance(rotation, Rotation): self.rotation = rotation else: - self.rotation = Rotation.fromQuaternion(rotation) # assume quaternion + self.rotation = Rotation.from_quaternion(rotation) # assume quaternion if self.rotation.quaternion.shape != (4,): raise NotImplementedError('Support for multiple rotations missing') diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 08f060887..bb8810c49 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -14,7 +14,7 @@ n = 1000 @pytest.fixture def default(): """A set of n random rotations.""" - return [Rotation.fromRandom() for r in range(n)] + return [Rotation.from_random() for r in range(n)] @pytest.fixture def reference_dir(reference_dir_base): @@ -36,7 +36,7 @@ class TestOrientation: @pytest.mark.parametrize('lattice',Lattice.lattices) def test_IPF(self,lattice): direction = np.random.random(3)*2.0-1 - for rot in [Rotation.fromRandom() for r in range(n//100)]: + for rot in [Rotation.from_random() for r in range(n//100)]: R = damask.Orientation(rot,lattice) color = R.IPFcolor(direction) for equivalent in R.equivalentOrientations(): @@ -45,7 +45,7 @@ class TestOrientation: @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relationship_forward_backward(self,model,lattice): - ori = Orientation(Rotation.fromRandom(),lattice) + ori = Orientation(Rotation.from_random(),lattice) for i,r in enumerate(ori.relatedOrientations(model)): ori2 = r.relatedOrientations(model)[i] misorientation = ori.rotation.misorientation(ori2.rotation) From 19638168e63ce2fc23cd7f60e63e7c0f3230948b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 13:23:32 +0200 Subject: [PATCH 08/32] more sensible type checking and errors --- python/damask/_rotation.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 3a046363f..7abe06110 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -98,7 +98,7 @@ class Rotation: """ if self.quaternion.shape != (4,): raise NotImplementedError('Support for multiple rotations missing') - if isinstance(other, Rotation): # rotate a rotation + if isinstance(other, Rotation): self_q = self.quaternion[0] self_p = self.quaternion[1:] other_q = other.quaternion[0] @@ -106,22 +106,23 @@ class Rotation: R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p), self_q*other_p + other_q*self_p + _P * np.cross(self_p,other_p))) return R.standardize() - elif isinstance(other, (tuple,np.ndarray)): - if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid + elif isinstance(other, np.ndarray): + if other.shape == (3,): A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:]) B = 2.0 * np.dot(self.quaternion[1:],other) C = 2.0 * _P*self.quaternion[0] return A*other + B*self.quaternion[1:] + C * np.cross(self.quaternion[1:],other) - elif other.shape == (3,3,): # rotate a single (3x3)-matrix + elif other.shape == (3,3,): return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) elif other.shape == (3,3,3,3,): - raise NotImplementedError('Support for rotation of 4th order tensors missing') + R = self.asMatrix() + return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) else: - return NotImplemented + raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') else: - return NotImplemented + raise TypeError('Cannot rotate {}'.format(type(other))) def __matmul__(self, other): @@ -130,7 +131,7 @@ class Rotation: details to be discussed """ - if isinstance(other, Rotation): # rotate a rotation + if isinstance(other, Rotation): q_m = self.quaternion[...,0:1] p_m = self.quaternion[...,1:] q_o = other.quaternion[...,0:1] @@ -157,8 +158,10 @@ class Rotation: if self.shape + (3,3,3,3) == other.shape: R = self.asMatrix() return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) + else: + raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') else: - raise ValueError + raise TypeError('Cannot rotate {}'.format(type(other))) def inverse(self): """In-place inverse rotation/backward rotation.""" From 61ac40c259773671291f0fa7400d7eaeb73a8c45 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 13:44:07 +0200 Subject: [PATCH 09/32] more tests --- python/tests/test_Orientation.py | 13 +++------- python/tests/test_Rotation.py | 44 ++++++++++++++++++++++++++++---- 2 files changed, 43 insertions(+), 14 deletions(-) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index bb8810c49..8dc526325 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -8,13 +8,8 @@ import damask from damask import Rotation from damask import Orientation from damask import Lattice - -n = 1000 -@pytest.fixture -def default(): - """A set of n random rotations.""" - return [Rotation.from_random() for r in range(n)] +n = 1000 @pytest.fixture def reference_dir(reference_dir_base): @@ -28,10 +23,10 @@ class TestOrientation: {'label':'green','RGB':[0,1,0],'direction':[0,1,1]}, {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}]) @pytest.mark.parametrize('lattice',['fcc','bcc']) - def test_IPF_cubic(self,default,color,lattice): + def test_IPF_cubic(self,color,lattice): cube = damask.Orientation(damask.Rotation(),lattice) for direction in set(permutations(np.array(color['direction']))): - assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB'])) + assert np.allclose(cube.IPFcolor(np.array(direction)),np.array(color['RGB'])) @pytest.mark.parametrize('lattice',Lattice.lattices) def test_IPF(self,lattice): @@ -57,7 +52,7 @@ class TestOrientation: reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model)) ori = Orientation(Rotation(),lattice) eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)]) - if update: + if update: coords = np.array([(1,i+1) for i,x in enumerate(eu)]) table = damask.Table(eu,{'Eulers':(3,)}) table.add('pos',coords) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index ffa1548d7..23ab03ff6 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -177,7 +177,7 @@ class TestRotation: (Rotation.from_Rodrigues, np.array([1,0,0,-1])), (Rotation.from_Rodrigues, np.array([1,1,0,1])), (Rotation.from_homochoric, np.array([2,2,2])) ]) - def test_invalid(self,function,invalid): + def test_invalid_value(self,function,invalid): with pytest.raises(ValueError): function(invalid) @@ -302,8 +302,42 @@ class TestRotation: assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a) - @pytest.mark.parametrize('x',[np.random.rand(3), - np.random.rand(3,3)]) - def test_rotation_identity(self,x): + @pytest.mark.parametrize('data',[np.random.rand(3), + np.random.rand(3,3), + np.random.rand(3,3,3,3)]) + def test_rotate_identity(self,data): R = Rotation() - assert np.allclose(x,R*x) + assert np.allclose(data,R*data) + + @pytest.mark.parametrize('data',[np.random.rand(3), + np.random.rand(3,3), + np.random.rand(3,3,3,3)]) + def test_rotate_360deg(self,data): + phi_1 = np.random.random() * np.pi + phi_2 = 2*np.pi - phi_1 + R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.])) + R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2])) + assert np.allclose(data,R_2*(R_1*data)) + + @pytest.mark.parametrize('data',[np.random.rand(3), + np.random.rand(3,3), + np.random.rand(3,3,3,3)]) + def test_rotate_inverse(self,data): + R = Rotation.from_random() + assert np.allclose(data,R.inversed()*(R*data)) + + @pytest.mark.parametrize('data',[np.random.rand(4), + np.random.rand(3,2), + np.random.rand(3,2,3,3)]) + def test_rotate_invalid_shape(self,data): + R = Rotation.from_random() + with pytest.raises(ValueError): + R*data + + @pytest.mark.parametrize('data',['does_not_work', + (1,2), + 5]) + def test_rotate_invalid_type(self,data): + R = Rotation.from_random() + with pytest.raises(TypeError): + R*data From a90865c8777616925afedc632f0b747435415022 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 16 May 2020 21:56:30 +0200 Subject: [PATCH 10/32] non-vectorized versions not needed anymore using them only for testing purposes --- python/damask/_rotation.py | 672 +++++++++------------------- python/tests/rotation_conversion.py | 399 +++++++++++++++++ python/tests/test_Rotation.py | 127 +++--- 3 files changed, 658 insertions(+), 540 deletions(-) create mode 100644 python/tests/rotation_conversion.py diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7abe06110..5258145e0 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -549,67 +549,42 @@ class Rotation: #---------- Quaternion ---------- @staticmethod def qu2om(qu): - if len(qu.shape) == 1: - """Quaternion to rotation matrix.""" - qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) - om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) - - om[0,1] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) - om[1,0] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) - om[1,2] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) - om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) - om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) - om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) - else: - qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) - om = np.block([qq + 2.0*qu[...,1:2]**2, - 2.0*(qu[...,2:3]*qu[...,1:2]+qu[...,0:1]*qu[...,3:4]), - 2.0*(qu[...,3:4]*qu[...,1:2]-qu[...,0:1]*qu[...,2:3]), - 2.0*(qu[...,1:2]*qu[...,2:3]-qu[...,0:1]*qu[...,3:4]), - qq + 2.0*qu[...,2:3]**2, - 2.0*(qu[...,3:4]*qu[...,2:3]+qu[...,0:1]*qu[...,1:2]), - 2.0*(qu[...,1:2]*qu[...,3:4]+qu[...,0:1]*qu[...,2:3]), - 2.0*(qu[...,2:3]*qu[...,3:4]-qu[...,0:1]*qu[...,1:2]), - qq + 2.0*qu[...,3:4]**2, - ]).reshape(qu.shape[:-1]+(3,3)) + qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) + om = np.block([qq + 2.0*qu[...,1:2]**2, + 2.0*(qu[...,2:3]*qu[...,1:2]+qu[...,0:1]*qu[...,3:4]), + 2.0*(qu[...,3:4]*qu[...,1:2]-qu[...,0:1]*qu[...,2:3]), + 2.0*(qu[...,1:2]*qu[...,2:3]-qu[...,0:1]*qu[...,3:4]), + qq + 2.0*qu[...,2:3]**2, + 2.0*(qu[...,3:4]*qu[...,2:3]+qu[...,0:1]*qu[...,1:2]), + 2.0*(qu[...,1:2]*qu[...,3:4]+qu[...,0:1]*qu[...,2:3]), + 2.0*(qu[...,2:3]*qu[...,3:4]-qu[...,0:1]*qu[...,1:2]), + qq + 2.0*qu[...,3:4]**2, + ]).reshape(qu.shape[:-1]+(3,3)) return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) @staticmethod def qu2eu(qu): """Quaternion to Bunge-Euler angles.""" - if len(qu.shape) == 1: - q03 = qu[0]**2+qu[3]**2 - q12 = qu[1]**2+qu[2]**2 - chi = np.sqrt(q03*q12) - if np.abs(q12) < 1.e-8: - eu = np.array([np.arctan2(-_P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) - elif np.abs(q03) < 1.e-8: - eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) - else: - eu = np.array([np.arctan2((-_P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), - np.arctan2( 2.0*chi, q03-q12 ), - np.arctan2(( _P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) - else: - q02 = qu[...,0:1]*qu[...,2:3] - q13 = qu[...,1:2]*qu[...,3:4] - q01 = qu[...,0:1]*qu[...,1:2] - q23 = qu[...,2:3]*qu[...,3:4] - q03_s = qu[...,0:1]**2+qu[...,3:4]**2 - q12_s = qu[...,1:2]**2+qu[...,2:3]**2 - chi = np.sqrt(q03_s*q12_s) + q02 = qu[...,0:1]*qu[...,2:3] + q13 = qu[...,1:2]*qu[...,3:4] + q01 = qu[...,0:1]*qu[...,1:2] + q23 = qu[...,2:3]*qu[...,3:4] + q03_s = qu[...,0:1]**2+qu[...,3:4]**2 + q12_s = qu[...,1:2]**2+qu[...,2:3]**2 + chi = np.sqrt(q03_s*q12_s) - eu = np.where(np.abs(q12_s) < 1.0e-8, - np.block([np.arctan2(-_P*2.0*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2), - np.zeros(qu.shape[:-1]+(2,))]), - np.where(np.abs(q03_s) < 1.0e-8, - np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2), - np.broadcast_to(np.pi,qu.shape[:-1]+(1,)), - np.zeros(qu.shape[:-1]+(1,))]), - np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi), - np.arctan2( 2.0*chi, q03_s-q12_s ), - np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)]) - ) - ) + eu = np.where(np.abs(q12_s) < 1.0e-8, + np.block([np.arctan2(-_P*2.0*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2), + np.zeros(qu.shape[:-1]+(2,))]), + np.where(np.abs(q03_s) < 1.0e-8, + np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2), + np.broadcast_to(np.pi,qu.shape[:-1]+(1,)), + np.zeros(qu.shape[:-1]+(1,))]), + np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi), + np.arctan2( 2.0*chi, q03_s-q12_s ), + np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)]) + ) + ) # reduce Euler angles to definition range eu[np.abs(eu)<1.e-6] = 0.0 eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) @@ -622,65 +597,38 @@ class Rotation: Modified version of the original formulation, should be numerically more stable """ - if len(qu.shape) == 1: - if np.abs(np.sum(qu[1:4]**2)) < 1.e-6: # set axis to [001] if the angle is 0/360 - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - elif qu[0] > 1.e-6: - s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) - omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) - ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ]) - else: - ax = ax = np.array([ qu[1], qu[2], qu[3], np.pi]) - else: - with np.errstate(invalid='ignore',divide='ignore'): - s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) - omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) - ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-6,qu.shape), - np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), - np.block([qu[...,1:4]*s,omega])) - ax[np.sum(np.abs(qu[...,1:4])**2,axis=-1) < 1.0e-6,] = [0.0, 0.0, 1.0, 0.0] + with np.errstate(invalid='ignore',divide='ignore'): + s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) + omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) + ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-6,qu.shape), + np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:4]*s,omega])) + ax[np.sum(np.abs(qu[...,1:4])**2,axis=-1) < 1.0e-6,] = [0.0, 0.0, 1.0, 0.0] return ax @staticmethod def qu2ro(qu): """Quaternion to Rodrigues-Frank vector.""" - if len(qu.shape) == 1: - if iszero(qu[0]): - ro = np.array([qu[1], qu[2], qu[3], np.inf]) - else: - s = np.linalg.norm(qu[1:4]) - ro = np.array([0.0,0.0,_P,0.0] if iszero(s) else \ - [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))]) - else: - with np.errstate(invalid='ignore',divide='ignore'): - s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) - ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape), - np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]), - np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s, - np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0))) - ]) - ) - ro[np.abs(s).squeeze(-1) < 1.0e-12] = [0.0,0.0,_P,0.0] + with np.errstate(invalid='ignore',divide='ignore'): + s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) + ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape), + np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s, + np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0))) + ]) + ) + ro[np.abs(s).squeeze(-1) < 1.0e-12] = [0.0,0.0,_P,0.0] return ro @staticmethod def qu2ho(qu): """Quaternion to homochoric vector.""" - if len(qu.shape) == 1: - omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) - if np.abs(omega) < 1.0e-12: - ho = np.zeros(3) - else: - ho = np.array([qu[1], qu[2], qu[3]]) - f = 0.75 * ( omega - np.sin(omega) ) - ho = ho/np.linalg.norm(ho) * f**(1./3.) - else: - with np.errstate(invalid='ignore'): - omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) - ho = np.where(np.abs(omega) < 1.0e-12, - np.zeros(3), - qu[...,1:4]/np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) \ - * (0.75*(omega - np.sin(omega)))**(1./3.)) + with np.errstate(invalid='ignore'): + omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) + ho = np.where(np.abs(omega) < 1.0e-12, + np.zeros(3), + qu[...,1:4]/np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) \ + * (0.75*(omega - np.sin(omega)))**(1./3.)) return ho @staticmethod @@ -702,27 +650,18 @@ class Rotation: @staticmethod def om2eu(om): """Rotation matrix to Bunge-Euler angles.""" - if len(om.shape) == 2: - if not np.isclose(np.abs(om[2,2]),1.0,1.e-4): - zeta = 1.0/np.sqrt(1.0-om[2,2]**2) - eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), - np.arccos(om[2,2]), - np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) - else: - eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation - else: - with np.errstate(invalid='ignore',divide='ignore'): - zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) - eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4), - np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]), - np.pi*0.5*(1-om[...,2,2:3]), - np.zeros(om.shape[:-2]+(1,)), - ]), - np.block([np.arctan2(om[...,2,0:1]*zeta,-om[...,2,1:2]*zeta), - np.arccos(om[...,2,2:3]), - np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta) - ]) - ) + with np.errstate(invalid='ignore',divide='ignore'): + zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) + eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4), + np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]), + np.pi*0.5*(1-om[...,2,2:3]), + np.zeros(om.shape[:-2]+(1,)), + ]), + np.block([np.arctan2(om[...,2,0:1]*zeta,-om[...,2,1:2]*zeta), + np.arccos(om[...,2,2:3]), + np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta) + ]) + ) eu[np.abs(eu)<1.e-6] = 0.0 eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) return eu @@ -731,40 +670,23 @@ class Rotation: @staticmethod def om2ax(om): """Rotation matrix to axis angle pair.""" - if len(om.shape) == 2: - ax=np.empty(4) - - # first get the rotation angle - t = 0.5*(om.trace() -1.0) - ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - if np.abs(ax[3])<1.e-6: - ax = np.array([ 0.0, 0.0, 1.0, 0.0]) - else: - w,vr = np.linalg.eig(om) - # next, find the eigenvalue (1,0j) - i = np.where(np.isclose(w,1.0+0.0j))[0][0] - ax[0:3] = np.real(vr[0:3,i]) - diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - diagDelta[np.abs(diagDelta)<1.e-6] = 1.0 - ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) - else: - diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], - om[...,2,0:1]-om[...,0,2:3], - om[...,0,1:2]-om[...,1,0:1] - ]) - diag_delta[np.abs(diag_delta)<1.e-6] = 1.0 - t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,)) - w,vr = np.linalg.eig(om) - # mask duplicated real eigenvalues - w[np.isclose(w[...,0],1.0+0.0j),1:] = 0. - w[np.isclose(w[...,1],1.0+0.0j),2:] = 0. - vr = np.swapaxes(vr,-1,-2) - ax = np.where(np.abs(diag_delta)<0, - np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)), - np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \ - *np.sign(diag_delta)) - ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))]) - ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0] + diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], + om[...,2,0:1]-om[...,0,2:3], + om[...,0,1:2]-om[...,1,0:1] + ]) + diag_delta[np.abs(diag_delta)<1.e-6] = 1.0 + t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,)) + w,vr = np.linalg.eig(om) + # mask duplicated real eigenvalues + w[np.isclose(w[...,0],1.0+0.0j),1:] = 0. + w[np.isclose(w[...,1],1.0+0.0j),2:] = 0. + vr = np.swapaxes(vr,-1,-2) + ax = np.where(np.abs(diag_delta)<0, + np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)), + np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \ + *np.sign(diag_delta)) + ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))]) + ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0] return ax @@ -788,103 +710,61 @@ class Rotation: @staticmethod def eu2qu(eu): """Bunge-Euler angles to quaternion.""" - if len(eu.shape) == 1: - ee = 0.5*eu - cPhi = np.cos(ee[1]) - sPhi = np.sin(ee[1]) - qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), - -_P*sPhi*np.cos(ee[0]-ee[2]), - -_P*sPhi*np.sin(ee[0]-ee[2]), - -_P*cPhi*np.sin(ee[0]+ee[2]) ]) - if qu[0] < 0.0: qu*=-1 - else: - ee = 0.5*eu - cPhi = np.cos(ee[...,1:2]) - sPhi = np.sin(ee[...,1:2]) - qu = np.block([ cPhi*np.cos(ee[...,0:1]+ee[...,2:3]), - -_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]), - -_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]), - -_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])]) - qu[qu[...,0]<0.0]*=-1 + ee = 0.5*eu + cPhi = np.cos(ee[...,1:2]) + sPhi = np.sin(ee[...,1:2]) + qu = np.block([ cPhi*np.cos(ee[...,0:1]+ee[...,2:3]), + -_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]), + -_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]), + -_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])]) + qu[qu[...,0]<0.0]*=-1 return qu @staticmethod def eu2om(eu): """Bunge-Euler angles to rotation matrix.""" - if len(eu.shape) == 1: - c = np.cos(eu) - s = np.sin(eu) - - om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], - [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], - [+s[0]*s[1], -c[0]*s[1], +c[1] ]]) - else: - c = np.cos(eu) - s = np.sin(eu) - om = np.block([+c[...,0:1]*c[...,2:3]-s[...,0:1]*s[...,2:3]*c[...,1:2], - +s[...,0:1]*c[...,2:3]+c[...,0:1]*s[...,2:3]*c[...,1:2], - +s[...,2:3]*s[...,1:2], - -c[...,0:1]*s[...,2:3]-s[...,0:1]*c[...,2:3]*c[...,1:2], - -s[...,0:1]*s[...,2:3]+c[...,0:1]*c[...,2:3]*c[...,1:2], - +c[...,2:3]*s[...,1:2], - +s[...,0:1]*s[...,1:2], - -c[...,0:1]*s[...,1:2], - +c[...,1:2] - ]).reshape(eu.shape[:-1]+(3,3)) + c = np.cos(eu) + s = np.sin(eu) + om = np.block([+c[...,0:1]*c[...,2:3]-s[...,0:1]*s[...,2:3]*c[...,1:2], + +s[...,0:1]*c[...,2:3]+c[...,0:1]*s[...,2:3]*c[...,1:2], + +s[...,2:3]*s[...,1:2], + -c[...,0:1]*s[...,2:3]-s[...,0:1]*c[...,2:3]*c[...,1:2], + -s[...,0:1]*s[...,2:3]+c[...,0:1]*c[...,2:3]*c[...,1:2], + +c[...,2:3]*s[...,1:2], + +s[...,0:1]*s[...,1:2], + -c[...,0:1]*s[...,1:2], + +c[...,1:2] + ]).reshape(eu.shape[:-1]+(3,3)) om[np.abs(om)<1.e-12] = 0.0 return om @staticmethod def eu2ax(eu): """Bunge-Euler angles to axis angle pair.""" - if len(eu.shape) == 1: - t = np.tan(eu[1]*0.5) - sigma = 0.5*(eu[0]+eu[2]) - delta = 0.5*(eu[0]-eu[2]) - tau = np.linalg.norm([t,np.sin(sigma)]) - alpha = np.pi if iszero(np.cos(sigma)) else \ - 2.0*np.arctan(tau/np.cos(sigma)) - - if np.abs(alpha)<1.e-6: - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - else: - ax = -_P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front - ax = np.append(ax,alpha) - if alpha < 0.0: ax *= -1.0 # ensure alpha is positive - else: - t = np.tan(eu[...,1:2]*0.5) - sigma = 0.5*(eu[...,0:1]+eu[...,2:3]) - delta = 0.5*(eu[...,0:1]-eu[...,2:3]) - tau = np.linalg.norm(np.block([t,np.sin(sigma)]),axis=-1,keepdims=True) - alpha = np.where(np.abs(np.cos(sigma))<1.e-12,np.pi,2.0*np.arctan(tau/np.cos(sigma))) - with np.errstate(invalid='ignore',divide='ignore'): - ax = np.where(np.broadcast_to(np.abs(alpha)<1.0e-12,eu.shape[:-1]+(4,)), - [0.0,0.0,1.0,0.0], - np.block([-_P/tau*t*np.cos(delta), - -_P/tau*t*np.sin(delta), - -_P/tau* np.sin(sigma), - alpha - ])) - ax[(alpha<0.0).squeeze()] *=-1 + t = np.tan(eu[...,1:2]*0.5) + sigma = 0.5*(eu[...,0:1]+eu[...,2:3]) + delta = 0.5*(eu[...,0:1]-eu[...,2:3]) + tau = np.linalg.norm(np.block([t,np.sin(sigma)]),axis=-1,keepdims=True) + alpha = np.where(np.abs(np.cos(sigma))<1.e-12,np.pi,2.0*np.arctan(tau/np.cos(sigma))) + with np.errstate(invalid='ignore',divide='ignore'): + ax = np.where(np.broadcast_to(np.abs(alpha)<1.0e-12,eu.shape[:-1]+(4,)), + [0.0,0.0,1.0,0.0], + np.block([-_P/tau*t*np.cos(delta), + -_P/tau*t*np.sin(delta), + -_P/tau* np.sin(sigma), + alpha + ])) + ax[(alpha<0.0).squeeze()] *=-1 return ax @staticmethod def eu2ro(eu): """Bunge-Euler angles to Rodrigues-Frank vector.""" - if len(eu.shape) == 1: - ro = Rotation.eu2ax(eu) # convert to axis angle pair representation - if ro[3] >= np.pi: # Differs from original implementation. check convention 5 - ro[3] = np.inf - elif iszero(ro[3]): - ro = np.array([ 0.0, 0.0, _P, 0.0 ]) - else: - ro[3] = np.tan(ro[3]*0.5) - else: - ax = Rotation.eu2ax(eu) - ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)]) - ro[ax[...,3]>=np.pi,3] = np.inf - ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, _P, 0.0 ] + ax = Rotation.eu2ax(eu) + ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)]) + ro[ax[...,3]>=np.pi,3] = np.inf + ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, _P, 0.0 ] return ro @staticmethod @@ -902,45 +782,26 @@ class Rotation: @staticmethod def ax2qu(ax): """Axis angle pair to quaternion.""" - if len(ax.shape) == 1: - if np.abs(ax[3])<1.e-6: - qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) - else: - c = np.cos(ax[3]*0.5) - s = np.sin(ax[3]*0.5) - qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - else: - c = np.cos(ax[...,3:4]*.5) - s = np.sin(ax[...,3:4]*.5) - qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s])) + c = np.cos(ax[...,3:4]*.5) + s = np.sin(ax[...,3:4]*.5) + qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s])) return qu @staticmethod def ax2om(ax): """Axis angle pair to rotation matrix.""" - if len(ax.shape) == 1: - c = np.cos(ax[3]) - s = np.sin(ax[3]) - omc = 1.0-c - om=np.diag(ax[0:3]**2*omc + c) - - for idx in [[0,1,2],[1,2,0],[2,0,1]]: - q = omc*ax[idx[0]] * ax[idx[1]] - om[idx[0],idx[1]] = q + s*ax[idx[2]] - om[idx[1],idx[0]] = q - s*ax[idx[2]] - else: - c = np.cos(ax[...,3:4]) - s = np.sin(ax[...,3:4]) - omc = 1. -c - om = np.block([c+omc*ax[...,0:1]**2, - omc*ax[...,0:1]*ax[...,1:2] + s*ax[...,2:3], - omc*ax[...,0:1]*ax[...,2:3] - s*ax[...,1:2], - omc*ax[...,0:1]*ax[...,1:2] - s*ax[...,2:3], - c+omc*ax[...,1:2]**2, - omc*ax[...,1:2]*ax[...,2:3] + s*ax[...,0:1], - omc*ax[...,0:1]*ax[...,2:3] + s*ax[...,1:2], - omc*ax[...,1:2]*ax[...,2:3] - s*ax[...,0:1], - c+omc*ax[...,2:3]**2]).reshape(ax.shape[:-1]+(3,3)) + c = np.cos(ax[...,3:4]) + s = np.sin(ax[...,3:4]) + omc = 1. -c + om = np.block([c+omc*ax[...,0:1]**2, + omc*ax[...,0:1]*ax[...,1:2] + s*ax[...,2:3], + omc*ax[...,0:1]*ax[...,2:3] - s*ax[...,1:2], + omc*ax[...,0:1]*ax[...,1:2] - s*ax[...,2:3], + c+omc*ax[...,1:2]**2, + omc*ax[...,1:2]*ax[...,2:3] + s*ax[...,0:1], + omc*ax[...,0:1]*ax[...,2:3] + s*ax[...,1:2], + omc*ax[...,1:2]*ax[...,2:3] - s*ax[...,0:1], + c+omc*ax[...,2:3]**2]).reshape(ax.shape[:-1]+(3,3)) return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) @staticmethod @@ -951,33 +812,19 @@ class Rotation: @staticmethod def ax2ro(ax): """Axis angle pair to Rodrigues-Frank vector.""" - if len(ax.shape) == 1: - if np.abs(ax[3])<1.e-6: - ro = [ 0.0, 0.0, _P, 0.0 ] - else: - ro = [ax[0], ax[1], ax[2]] - # 180 degree case - ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ - [np.tan(ax[3]*0.5)] - ro = np.array(ro) - else: - ro = np.block([ax[...,:3], - np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), - np.inf, - np.tan(ax[...,3:4]*0.5)) - ]) - ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0] + ro = np.block([ax[...,:3], + np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), + np.inf, + np.tan(ax[...,3:4]*0.5)) + ]) + ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0] return ro @staticmethod def ax2ho(ax): """Axis angle pair to homochoric vector.""" - if len(ax.shape) == 1: - f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) - ho = ax[0:3] * f - else: - f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) - ho = ax[...,:3] * f + f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) + ho = ax[...,:3] * f return ho @staticmethod @@ -1005,36 +852,19 @@ class Rotation: @staticmethod def ro2ax(ro): """Rodrigues-Frank vector to axis angle pair.""" - if len(ro.shape) == 1: - if np.abs(ro[3]) < 1.e-6: - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - elif not np.isfinite(ro[3]): - ax = np.array([ ro[0], ro[1], ro[2], np.pi ]) - else: - angle = 2.0*np.arctan(ro[3]) - ta = np.linalg.norm(ro[0:3]) - ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ]) - else: - with np.errstate(invalid='ignore',divide='ignore'): - ax = np.where(np.isfinite(ro[...,3:4]), - np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]), - np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)])) - ax[np.abs(ro[...,3]) < 1.e-6] = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + with np.errstate(invalid='ignore',divide='ignore'): + ax = np.where(np.isfinite(ro[...,3:4]), + np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]), + np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)])) + ax[np.abs(ro[...,3]) < 1.e-6] = np.array([ 0.0, 0.0, 1.0, 0.0 ]) return ax @staticmethod def ro2ho(ro): """Rodrigues-Frank vector to homochoric vector.""" - if len(ro.shape) == 1: - if np.sum(ro[0:3]**2.0) < 1.e-6: - ho = np.zeros(3) - else: - f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi - ho = ro[0:3] * (0.75*f)**(1.0/3.0) - else: - f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) - ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape), - np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0)) + f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) + ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape), + np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0)) return ho @staticmethod @@ -1070,31 +900,16 @@ class Rotation: +0.0001703481934140054, -0.00012062065004116828, +0.000059719705868660826, -0.00001980756723965647, +0.000003953714684212874, -0.00000036555001439719544]) - if len(ho.shape) == 1: - # normalize h and store the magnitude - hmag_squared = np.sum(ho**2.) - if iszero(hmag_squared): - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - else: - hm = hmag_squared - - # convert the magnitude to the rotation angle - s = tfit[0] + tfit[1] * hmag_squared - for i in range(2,16): - hm *= hmag_squared - s += tfit[i] * hm - ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) - else: - hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True) - hm = hmag_squared.copy() - s = tfit[0] + tfit[1] * hmag_squared - for i in range(2,16): - hm *= hmag_squared - s += tfit[i] * hm - with np.errstate(invalid='ignore'): - ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-6,ho.shape[:-1]+(4,)), - [ 0.0, 0.0, 1.0, 0.0 ], - np.block([ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))])) + hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True) + hm = hmag_squared.copy() + s = tfit[0] + tfit[1] * hmag_squared + for i in range(2,16): + hm *= hmag_squared + s += tfit[i] * hm + with np.errstate(invalid='ignore'): + ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-6,ho.shape[:-1]+(4,)), + [ 0.0, 0.0, 1.0, 0.0 ], + np.block([ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))])) return ax @staticmethod @@ -1113,60 +928,31 @@ class Rotation: https://doi.org/10.1088/0965-0393/22/7/075013 """ - if len(ho.shape) == 1: - rs = np.linalg.norm(ho) + rs = np.linalg.norm(ho,axis=-1,keepdims=True) - if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16): - cu = np.zeros(3) - else: - xyz3 = ho[Rotation._get_pyramid_order(ho,'forward')] + xyz3 = np.take_along_axis(ho,Rotation._get_pyramid_order(ho,'forward'),-1) - # inverse M_3 - xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) + with np.errstate(invalid='ignore',divide='ignore'): + # inverse M_3 + xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) ) + qxy = np.sum(xyz2**2,axis=-1,keepdims=True) - # inverse M_2 - qxy = np.sum(xyz2**2) + q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2 + sq2 = np.sqrt(q2) + q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)) + tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\ + +np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) + T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]), + np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]), + np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q + T_inv[xyz2<0.0] *= -1.0 + T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0 + cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \ + * rs/np.sqrt(6.0/np.pi), + ])/ _sc - if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16): - Tinv = np.zeros(2) - else: - q2 = qxy + np.max(np.abs(xyz2))**2 - sq2 = np.sqrt(q2) - q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) - tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) - Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ - np.array([np.arccos(tt)/np.pi*12.0,1.0]) - Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) - - # inverse M_1 - cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc - cu = cu[Rotation._get_pyramid_order(ho,'backward')] - else: - rs = np.linalg.norm(ho,axis=-1,keepdims=True) - - xyz3 = np.take_along_axis(ho,Rotation._get_pyramid_order(ho,'forward'),-1) - - with np.errstate(invalid='ignore',divide='ignore'): - # inverse M_3 - xyz2 = xyz3[...,0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[...,2:3])) ) - qxy = np.sum(xyz2**2,axis=-1,keepdims=True) - - q2 = qxy + np.max(np.abs(xyz2),axis=-1,keepdims=True)**2 - sq2 = np.sqrt(q2) - q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)) - tt = np.clip((np.min(np.abs(xyz2),axis=-1,keepdims=True)**2\ - +np.max(np.abs(xyz2),axis=-1,keepdims=True)*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) - T_inv = np.where(np.abs(xyz2[...,1:2]) <= np.abs(xyz2[...,0:1]), - np.block([np.ones_like(tt),np.arccos(tt)/np.pi*12.0]), - np.block([np.arccos(tt)/np.pi*12.0,np.ones_like(tt)]))*q - T_inv[xyz2<0.0] *= -1.0 - T_inv[np.broadcast_to(np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-12),T_inv.shape)] = 0.0 - cu = np.block([T_inv, np.where(xyz3[...,2:3]<0.0,-np.ones_like(xyz3[...,2:3]),np.ones_like(xyz3[...,2:3])) \ - * rs/np.sqrt(6.0/np.pi), - ])/ _sc - - cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 - cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1) + cu[np.isclose(np.sum(np.abs(ho),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 + cu = np.take_along_axis(cu,Rotation._get_pyramid_order(ho,'backward'),-1) return cu @@ -1207,64 +993,34 @@ class Rotation: https://doi.org/10.1088/0965-0393/22/7/075013 """ - if len(cu.shape) == 1: - # transform to the sphere grid via the curved square, and intercept the zero point - if np.allclose(cu,0.0,rtol=0.0,atol=1.0e-16): - ho = np.zeros(3) - else: - # get pyramide and scale by grid parameter ratio - XYZ = cu[Rotation._get_pyramid_order(cu,'forward')] * _sc + with np.errstate(invalid='ignore',divide='ignore'): + # get pyramide and scale by grid parameter ratio + XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc + order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1]) + q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \ + / np.where(order,XYZ[...,0:1],XYZ[...,1:2]) + c = np.cos(q) + s = np.sin(q) + q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \ + * np.where(order,XYZ[...,0:1],XYZ[...,1:2]) - # intercept all the points along the z-axis - if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): - ho = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) - else: - order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] - q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] - c = np.cos(q) - s = np.sin(q) - q = _R1*2.0**0.25/_beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) - T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q + T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q - # transform to sphere grid (inverse Lambert) - # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero - c = np.sum(T**2) - s = c * np.pi/24.0 /XYZ[2]**2 - c = c * np.sqrt(np.pi/24.0)/XYZ[2] + # transform to sphere grid (inverse Lambert) + c = np.sum(T**2,axis=-1,keepdims=True) + s = c * np.pi/24.0 /XYZ[...,2:3]**2 + c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3] + q = np.sqrt( 1.0 - s) - q = np.sqrt( 1.0 - s ) - ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) + ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16), + np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]), + np.block([np.where(order,T[...,0:1],T[...,1:2])*q, + np.where(order,T[...,1:2],T[...,0:1])*q, + np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c]) + ) - ho = ho[Rotation._get_pyramid_order(cu,'backward')] - else: - with np.errstate(invalid='ignore',divide='ignore'): - # get pyramide and scale by grid parameter ratio - XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc - order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1]) - q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \ - / np.where(order,XYZ[...,0:1],XYZ[...,1:2]) - c = np.cos(q) - s = np.sin(q) - q = _R1*2.0**0.25/_beta/ np.sqrt(np.sqrt(2.0)-c) \ - * np.where(order,XYZ[...,0:1],XYZ[...,1:2]) - - T = np.block([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q - - # transform to sphere grid (inverse Lambert) - c = np.sum(T**2,axis=-1,keepdims=True) - s = c * np.pi/24.0 /XYZ[...,2:3]**2 - c = c * np.sqrt(np.pi/24.0)/XYZ[...,2:3] - q = np.sqrt( 1.0 - s) - - ho = np.where(np.isclose(np.sum(np.abs(XYZ[...,0:2]),axis=-1,keepdims=True),0.0,rtol=0.0,atol=1.0e-16), - np.block([np.zeros_like(XYZ[...,0:2]),np.sqrt(6.0/np.pi) *XYZ[...,2:3]]), - np.block([np.where(order,T[...,0:1],T[...,1:2])*q, - np.where(order,T[...,1:2],T[...,0:1])*q, - np.sqrt(6.0/np.pi) * XYZ[...,2:3] - c]) - ) - - ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 - ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1) + ho[np.isclose(np.sum(np.abs(cu),axis=-1),0.0,rtol=0.0,atol=1.0e-16)] = 0.0 + ho = np.take_along_axis(ho,Rotation._get_pyramid_order(cu,'backward'),-1) return ho @@ -1288,20 +1044,10 @@ class Rotation: https://doi.org/10.1088/0965-0393/22/7/075013 """ - order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]), + order = {'forward': np.array([[0,1,2],[1,2,0],[2,0,1]]), 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} - if len(xyz.shape) == 1: - if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \ - np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]: - p = 0 - elif np.maximum(abs(xyz[1]),abs(xyz[2])) <= xyz[0] or \ - np.maximum(abs(xyz[1]),abs(xyz[2])) <=-xyz[0]: - p = 1 - elif np.maximum(abs(xyz[2]),abs(xyz[0])) <= xyz[1] or \ - np.maximum(abs(xyz[2]),abs(xyz[0])) <=-xyz[1]: - p = 2 - else: - p = np.where(np.maximum(np.abs(xyz[...,0]),np.abs(xyz[...,1])) <= np.abs(xyz[...,2]),0, - np.where(np.maximum(np.abs(xyz[...,1]),np.abs(xyz[...,2])) <= np.abs(xyz[...,0]),1,2)) + + p = np.where(np.maximum(np.abs(xyz[...,0]),np.abs(xyz[...,1])) <= np.abs(xyz[...,2]),0, + np.where(np.maximum(np.abs(xyz[...,1]),np.abs(xyz[...,2])) <= np.abs(xyz[...,0]),1,2)) return order[direction][p] diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py new file mode 100644 index 000000000..5ca10f443 --- /dev/null +++ b/python/tests/rotation_conversion.py @@ -0,0 +1,399 @@ +#################################################################################################### +# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations +#################################################################################################### +# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH +# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without modification, are +# permitted provided that the following conditions are met: +# +# - Redistributions of source code must retain the above copyright notice, this list +# of conditions and the following disclaimer. +# - Redistributions in binary form must reproduce the above copyright notice, this +# list of conditions and the following disclaimer in the documentation and/or +# other materials provided with the distribution. +# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +# of its contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +#################################################################################################### + +import numpy as np + +_P = -1 + +# parameters for conversion from/to cubochoric +_sc = np.pi**(1./6.)/6.**(1./6.) +_beta = np.pi**(5./6.)/6.**(1./6.)/2. +_R1 = (3.*np.pi/4.)**(1./3.) + +def iszero(a): + return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) + +#---------- Quaternion ---------- +def qu2om(qu): + """Quaternion to rotation matrix.""" + qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) + om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) + + om[0,1] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) + om[1,0] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) + om[1,2] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) + om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) + om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) + om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) + return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + +def qu2eu(qu): + """Quaternion to Bunge-Euler angles.""" + q03 = qu[0]**2+qu[3]**2 + q12 = qu[1]**2+qu[2]**2 + chi = np.sqrt(q03*q12) + if np.abs(q12) < 1.e-8: + eu = np.array([np.arctan2(-_P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) + elif np.abs(q03) < 1.e-8: + eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) + else: + eu = np.array([np.arctan2((-_P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), + np.arctan2( 2.0*chi, q03-q12 ), + np.arctan2(( _P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) + # reduce Euler angles to definition range + eu[np.abs(eu)<1.e-6] = 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + +def qu2ax(qu): + """ + Quaternion to axis angle pair. + + Modified version of the original formulation, should be numerically more stable + """ + if np.abs(np.sum(qu[1:4]**2)) < 1.e-6: # set axis to [001] if the angle is 0/360 + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + elif qu[0] > 1.e-6: + s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) + ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ]) + else: + ax = ax = np.array([ qu[1], qu[2], qu[3], np.pi]) + return ax + +def qu2ro(qu): + """Quaternion to Rodrigues-Frank vector.""" + if iszero(qu[0]): + ro = np.array([qu[1], qu[2], qu[3], np.inf]) + else: + s = np.linalg.norm(qu[1:4]) + ro = np.array([0.0,0.0,_P,0.0] if iszero(s) else \ + [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))]) + return ro + +def qu2ho(qu): + """Quaternion to homochoric vector.""" + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) + if np.abs(omega) < 1.0e-12: + ho = np.zeros(3) + else: + ho = np.array([qu[1], qu[2], qu[3]]) + f = 0.75 * ( omega - np.sin(omega) ) + ho = ho/np.linalg.norm(ho) * f**(1./3.) + return ho + + +#---------- Rotation matrix ---------- +def om2eu(om): + """Rotation matrix to Bunge-Euler angles.""" + if not np.isclose(np.abs(om[2,2]),1.0,1.e-4): + zeta = 1.0/np.sqrt(1.0-om[2,2]**2) + eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), + np.arccos(om[2,2]), + np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) + else: + eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation + eu[np.abs(eu)<1.e-6] = 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + +def om2ax(om): + """Rotation matrix to axis angle pair.""" + ax=np.empty(4) + + # first get the rotation angle + t = 0.5*(om.trace() -1.0) + ax[3] = np.arccos(np.clip(t,-1.0,1.0)) + if np.abs(ax[3])<1.e-6: + ax = np.array([ 0.0, 0.0, 1.0, 0.0]) + else: + w,vr = np.linalg.eig(om) + # next, find the eigenvalue (1,0j) + i = np.where(np.isclose(w,1.0+0.0j))[0][0] + ax[0:3] = np.real(vr[0:3,i]) + diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) + diagDelta[np.abs(diagDelta)<1.e-6] = 1.0 + ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) + return ax + +#---------- Bunge-Euler angles ---------- +def eu2qu(eu): + """Bunge-Euler angles to quaternion.""" + ee = 0.5*eu + cPhi = np.cos(ee[1]) + sPhi = np.sin(ee[1]) + qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), + -_P*sPhi*np.cos(ee[0]-ee[2]), + -_P*sPhi*np.sin(ee[0]-ee[2]), + -_P*cPhi*np.sin(ee[0]+ee[2]) ]) + if qu[0] < 0.0: qu*=-1 + return qu + +def eu2om(eu): + """Bunge-Euler angles to rotation matrix.""" + c = np.cos(eu) + s = np.sin(eu) + + om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], + [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], + [+s[0]*s[1], -c[0]*s[1], +c[1] ]]) + om[np.abs(om)<1.e-12] = 0.0 + return om + +def eu2ax(eu): + """Bunge-Euler angles to axis angle pair.""" + t = np.tan(eu[1]*0.5) + sigma = 0.5*(eu[0]+eu[2]) + delta = 0.5*(eu[0]-eu[2]) + tau = np.linalg.norm([t,np.sin(sigma)]) + alpha = np.pi if iszero(np.cos(sigma)) else \ + 2.0*np.arctan(tau/np.cos(sigma)) + + if np.abs(alpha)<1.e-6: + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + else: + ax = -_P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front + ax = np.append(ax,alpha) + if alpha < 0.0: ax *= -1.0 # ensure alpha is positive + return ax + +def eu2ro(eu): + """Bunge-Euler angles to Rodrigues-Frank vector.""" + ro = eu2ax(eu) # convert to axis angle pair representation + if ro[3] >= np.pi: # Differs from original implementation. check convention 5 + ro[3] = np.inf + elif iszero(ro[3]): + ro = np.array([ 0.0, 0.0, _P, 0.0 ]) + else: + ro[3] = np.tan(ro[3]*0.5) + return ro + +#---------- Axis angle pair ---------- +def ax2qu(ax): + """Axis angle pair to quaternion.""" + if np.abs(ax[3])<1.e-6: + qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) + else: + c = np.cos(ax[3]*0.5) + s = np.sin(ax[3]*0.5) + qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) + return qu + +def ax2om(ax): + """Axis angle pair to rotation matrix.""" + c = np.cos(ax[3]) + s = np.sin(ax[3]) + omc = 1.0-c + om=np.diag(ax[0:3]**2*omc + c) + + for idx in [[0,1,2],[1,2,0],[2,0,1]]: + q = omc*ax[idx[0]] * ax[idx[1]] + om[idx[0],idx[1]] = q + s*ax[idx[2]] + om[idx[1],idx[0]] = q - s*ax[idx[2]] + return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + +def ax2ro(ax): + """Axis angle pair to Rodrigues-Frank vector.""" + if np.abs(ax[3])<1.e-6: + ro = [ 0.0, 0.0, _P, 0.0 ] + else: + ro = [ax[0], ax[1], ax[2]] + # 180 degree case + ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ + [np.tan(ax[3]*0.5)] + ro = np.array(ro) + return ro + +def ax2ho(ax): + """Axis angle pair to homochoric vector.""" + f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) + ho = ax[0:3] * f + return ho + + +#---------- Rodrigues-Frank vector ---------- +def ro2ax(ro): + """Rodrigues-Frank vector to axis angle pair.""" + if np.abs(ro[3]) < 1.e-6: + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + elif not np.isfinite(ro[3]): + ax = np.array([ ro[0], ro[1], ro[2], np.pi ]) + else: + angle = 2.0*np.arctan(ro[3]) + ta = np.linalg.norm(ro[0:3]) + ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ]) + return ax + +def ro2ho(ro): + """Rodrigues-Frank vector to homochoric vector.""" + if np.sum(ro[0:3]**2.0) < 1.e-6: + ho = np.zeros(3) + else: + f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi + ho = ro[0:3] * (0.75*f)**(1.0/3.0) + return ho + +#---------- Homochoric vector---------- +def ho2ax(ho): + """Homochoric vector to axis angle pair.""" + tfit = np.array([+1.0000000000018852, -0.5000000002194847, + -0.024999992127593126, -0.003928701544781374, + -0.0008152701535450438, -0.0002009500426119712, + -0.00002397986776071756, -0.00008202868926605841, + +0.00012448715042090092, -0.0001749114214822577, + +0.0001703481934140054, -0.00012062065004116828, + +0.000059719705868660826, -0.00001980756723965647, + +0.000003953714684212874, -0.00000036555001439719544]) + # normalize h and store the magnitude + hmag_squared = np.sum(ho**2.) + if iszero(hmag_squared): + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + else: + hm = hmag_squared + + # convert the magnitude to the rotation angle + s = tfit[0] + tfit[1] * hmag_squared + for i in range(2,16): + hm *= hmag_squared + s += tfit[i] * hm + ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) + return ax + +def ho2cu(ho): + """ + Homochoric vector to cubochoric vector. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + rs = np.linalg.norm(ho) + + if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16): + cu = np.zeros(3) + else: + xyz3 = ho[_get_pyramid_order(ho,'forward')] + + # inverse M_3 + xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) + + # inverse M_2 + qxy = np.sum(xyz2**2) + + if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16): + Tinv = np.zeros(2) + else: + q2 = qxy + np.max(np.abs(xyz2))**2 + sq2 = np.sqrt(q2) + q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) + tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) + Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ + np.array([np.arccos(tt)/np.pi*12.0,1.0]) + Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) + + # inverse M_1 + cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc + cu = cu[_get_pyramid_order(ho,'backward')] + return cu + +#---------- Cubochoric ---------- +def cu2ho(cu): + """ + Cubochoric vector to homochoric vector. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + # transform to the sphere grid via the curved square, and intercept the zero point + if np.allclose(cu,0.0,rtol=0.0,atol=1.0e-16): + ho = np.zeros(3) + else: + # get pyramide and scale by grid parameter ratio + XYZ = cu[_get_pyramid_order(cu,'forward')] * _sc + + # intercept all the points along the z-axis + if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): + ho = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) + else: + order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] + q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] + c = np.cos(q) + s = np.sin(q) + q = _R1*2.0**0.25/_beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) + T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q + + # transform to sphere grid (inverse Lambert) + # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero + c = np.sum(T**2) + s = c * np.pi/24.0 /XYZ[2]**2 + c = c * np.sqrt(np.pi/24.0)/XYZ[2] + + q = np.sqrt( 1.0 - s ) + ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) + + ho = ho[_get_pyramid_order(cu,'backward')] + return ho + +def _get_pyramid_order(xyz,direction=None): + """ + Get order of the coordinates. + + Depending on the pyramid in which the point is located, the order need to be adjusted. + + Parameters + ---------- + xyz : numpy.ndarray + coordinates of a point on a uniform refinable grid on a ball or + in a uniform refinable cubical grid. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]), + 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} + if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \ + np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]: + p = 0 + elif np.maximum(abs(xyz[1]),abs(xyz[2])) <= xyz[0] or \ + np.maximum(abs(xyz[1]),abs(xyz[2])) <=-xyz[0]: + p = 1 + elif np.maximum(abs(xyz[2]),abs(xyz[0])) <= xyz[1] or \ + np.maximum(abs(xyz[2]),abs(xyz[0])) <=-xyz[1]: + p = 2 + return order[direction][p] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 23ab03ff6..4a6350fda 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -4,6 +4,7 @@ import pytest import numpy as np from damask import Rotation +import rotation_conversion n = 1100 atol=1.e-4 @@ -181,111 +182,83 @@ class TestRotation: with pytest.raises(ValueError): function(invalid) - @pytest.mark.parametrize('conversion',[Rotation.qu2om, - Rotation.qu2eu, - Rotation.qu2ax, - Rotation.qu2ro, - Rotation.qu2ho, - Rotation.qu2cu - ]) - def test_quaternion_vectorization(self,default,conversion): + @pytest.mark.parametrize('vectorized, single',[(Rotation.qu2om,rotation_conversion.qu2om), + (Rotation.qu2eu,rotation_conversion.qu2eu), + (Rotation.qu2ax,rotation_conversion.qu2ax), + (Rotation.qu2ro,rotation_conversion.qu2ro), + (Rotation.qu2ho,rotation_conversion.qu2ho)]) + def test_quaternion_vectorization(self,default,vectorized,single): qu = np.array([rot.as_quaternion() for rot in default]) - conversion(qu.reshape(qu.shape[0]//2,-1,4)) - co = conversion(qu) + vectorized(qu.reshape(qu.shape[0]//2,-1,4)) + co = vectorized(qu) for q,c in zip(qu,co): print(q,c) - assert np.allclose(conversion(q),c) + assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q)) - @pytest.mark.parametrize('conversion',[Rotation.om2qu, - Rotation.om2eu, - Rotation.om2ax, - Rotation.om2ro, - Rotation.om2ho, - Rotation.om2cu - ]) - def test_matrix_vectorization(self,default,conversion): + + @pytest.mark.parametrize('vectorized, single',[(Rotation.om2eu,rotation_conversion.om2eu), + (Rotation.om2ax,rotation_conversion.om2ax)]) + def test_matrix_vectorization(self,default,vectorized,single): om = np.array([rot.as_matrix() for rot in default]) - conversion(om.reshape(om.shape[0]//2,-1,3,3)) - co = conversion(om) + vectorized(om.reshape(om.shape[0]//2,-1,3,3)) + co = vectorized(om) for o,c in zip(om,co): print(o,c) - assert np.allclose(conversion(o),c) + assert np.allclose(single(o),c) and np.allclose(single(o),vectorized(o)) - @pytest.mark.parametrize('conversion',[Rotation.eu2qu, - Rotation.eu2om, - Rotation.eu2ax, - Rotation.eu2ro, - Rotation.eu2ho, - Rotation.eu2cu - ]) - def test_Euler_vectorization(self,default,conversion): + @pytest.mark.parametrize('vectorized, single',[(Rotation.eu2qu,rotation_conversion.eu2qu), + (Rotation.eu2om,rotation_conversion.eu2om), + (Rotation.eu2ax,rotation_conversion.eu2ax), + (Rotation.eu2ro,rotation_conversion.eu2ro)]) + def test_Euler_vectorization(self,default,vectorized,single): eu = np.array([rot.as_Eulers() for rot in default]) - conversion(eu.reshape(eu.shape[0]//2,-1,3)) - co = conversion(eu) + vectorized(eu.reshape(eu.shape[0]//2,-1,3)) + co = vectorized(eu) for e,c in zip(eu,co): print(e,c) - assert np.allclose(conversion(e),c) + assert np.allclose(single(e),c) and np.allclose(single(e),vectorized(e)) - @pytest.mark.parametrize('conversion',[Rotation.ax2qu, - Rotation.ax2om, - Rotation.ax2eu, - Rotation.ax2ro, - Rotation.ax2ho, - Rotation.ax2cu - ]) - def test_axisAngle_vectorization(self,default,conversion): + @pytest.mark.parametrize('vectorized, single',[(Rotation.ax2qu,rotation_conversion.ax2qu), + (Rotation.ax2om,rotation_conversion.ax2om), + (Rotation.ax2ro,rotation_conversion.ax2ro), + (Rotation.ax2ho,rotation_conversion.ax2ho)]) + def test_axisAngle_vectorization(self,default,vectorized,single): ax = np.array([rot.as_axis_angle() for rot in default]) - conversion(ax.reshape(ax.shape[0]//2,-1,4)) - co = conversion(ax) + vectorized(ax.reshape(ax.shape[0]//2,-1,4)) + co = vectorized(ax) for a,c in zip(ax,co): print(a,c) - assert np.allclose(conversion(a),c) + assert np.allclose(single(a),c) and np.allclose(single(a),vectorized(a)) - @pytest.mark.parametrize('conversion',[Rotation.ro2qu, - Rotation.ro2om, - Rotation.ro2eu, - Rotation.ro2ax, - Rotation.ro2ho, - Rotation.ro2cu - ]) - def test_Rodrigues_vectorization(self,default,conversion): + @pytest.mark.parametrize('vectorized, single',[(Rotation.ro2ax,rotation_conversion.ro2ax), + (Rotation.ro2ho,rotation_conversion.ro2ho)]) + def test_Rodrigues_vectorization(self,default,vectorized,single): ro = np.array([rot.as_Rodrigues() for rot in default]) - conversion(ro.reshape(ro.shape[0]//2,-1,4)) - co = conversion(ro) + vectorized(ro.reshape(ro.shape[0]//2,-1,4)) + co = vectorized(ro) for r,c in zip(ro,co): print(r,c) - assert np.allclose(conversion(r),c) + assert np.allclose(single(r),c) and np.allclose(single(r),vectorized(r)) - @pytest.mark.parametrize('conversion',[Rotation.ho2qu, - Rotation.ho2om, - Rotation.ho2eu, - Rotation.ho2ax, - Rotation.ho2ro, - Rotation.ho2cu - ]) - def test_homochoric_vectorization(self,default,conversion): + @pytest.mark.parametrize('vectorized, single',[(Rotation.ho2ax,rotation_conversion.ho2ax), + (Rotation.ho2cu,rotation_conversion.ho2cu)]) + def test_homochoric_vectorization(self,default,vectorized,single): ho = np.array([rot.as_homochoric() for rot in default]) - conversion(ho.reshape(ho.shape[0]//2,-1,3)) - co = conversion(ho) + vectorized(ho.reshape(ho.shape[0]//2,-1,3)) + co = vectorized(ho) for h,c in zip(ho,co): print(h,c) - assert np.allclose(conversion(h),c) + assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) - @pytest.mark.parametrize('conversion',[Rotation.cu2qu, - Rotation.cu2om, - Rotation.cu2eu, - Rotation.cu2ax, - Rotation.cu2ro, - Rotation.cu2ho - ]) - def test_cubochoric_vectorization(self,default,conversion): + @pytest.mark.parametrize('vectorized, single',[(Rotation.cu2ho,rotation_conversion.cu2ho)]) + def test_cubochoric_vectorization(self,default,vectorized,single): cu = np.array([rot.as_cubochoric() for rot in default]) - conversion(cu.reshape(cu.shape[0]//2,-1,3)) - co = conversion(cu) + vectorized(cu.reshape(cu.shape[0]//2,-1,3)) + co = vectorized(cu) for u,c in zip(cu,co): print(u,c) - assert np.allclose(conversion(u),c) + assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u)) @pytest.mark.parametrize('direction',['forward', 'backward']) From 652ece6bb3c46084a2256125551d61f72d2655a4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 17 May 2020 07:25:17 +0200 Subject: [PATCH 11/32] fix: wrong capitalization --- processing/post/addOrientations.py | 2 +- python/damask/_rotation.py | 3 --- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index b6292b18c..706a959ad 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -137,7 +137,7 @@ for name in filenames: if 'rodrigues' in options.output: - table.add('ro({})'.format(label),o.as_rodrigues(), scriptID+' '+' '.join(sys.argv[1:])) + table.add('ro({})'.format(label),o.as_Rodrigues(), scriptID+' '+' '.join(sys.argv[1:])) if 'eulers' in options.output: table.add('eu({})'.format(label),o.as_Eulers(options.degrees), scriptID+' '+' '.join(sys.argv[1:])) if 'quaternion' in options.output: diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 5258145e0..44484ccd9 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -9,9 +9,6 @@ _sc = np.pi**(1./6.)/6.**(1./6.) _beta = np.pi**(5./6.)/6.**(1./6.)/2. _R1 = (3.*np.pi/4.)**(1./3.) -def iszero(a): - return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) - class Rotation: u""" From 743e91a78d6f9d1509238d6886758cc2d5f8cad6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 17 May 2020 08:01:34 +0200 Subject: [PATCH 12/32] cleaning --- processing/pre/geom_fromDREAM3D.py | 2 +- processing/pre/geom_fromTable.py | 2 +- processing/pre/geom_rotate.py | 2 +- python/damask/_rotation.py | 18 ++---------------- 4 files changed, 5 insertions(+), 19 deletions(-) diff --git a/processing/pre/geom_fromDREAM3D.py b/processing/pre/geom_fromDREAM3D.py index b000aa811..c48698e53 100755 --- a/processing/pre/geom_fromDREAM3D.py +++ b/processing/pre/geom_fromDREAM3D.py @@ -140,7 +140,7 @@ for name in filenames: config_header = [''] for i in range(np.nanmax(microstructure)): config_header += ['[{}{}]'.format(label,i+1), - '(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*rot[i].asEulers(degrees = True)), + '(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*rot[i].as_Eulers(degrees = True)), ] config_header += [''] for i in range(np.nanmax(microstructure)): diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index f6c9ba38e..be5b3a23c 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -89,7 +89,7 @@ for name in filenames: for i,data in enumerate(unique): ori = damask.Rotation(data[0:4]) config_header += ['[Grain{}]'.format(i+1), - '(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*ori.asEulers(degrees = True)), + '(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*ori.as_Eulers(degrees = True)), ] if options.axes is not None: config_header += ['axes\t{} {} {}'.format(*options.axes)] diff --git a/processing/pre/geom_rotate.py b/processing/pre/geom_rotate.py index 7a774200b..dc2e8e2a3 100755 --- a/processing/pre/geom_rotate.py +++ b/processing/pre/geom_rotate.py @@ -67,7 +67,7 @@ if options.matrix is not None: if options.eulers is not None: rot = damask.Rotation.from_Eulers(np.array(options.eulers),degrees=options.degrees) -eulers = rot.asEulers(degrees=True) +eulers = rot.as_Eulers(degrees=True) if filenames == []: filenames = [None] diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 44484ccd9..bf5a2a4f3 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -305,15 +305,6 @@ class Rotation: """ return np.einsum('...i,...j',self.quaternion,self.quaternion) - # for compatibility (old names do not follow convention) - asM = M - asQuaternion = as_quaternion - asEulers = as_Eulers - asAxisAngle = as_axis_angle - asMatrix = as_matrix - asRodrigues = as_Rodrigues - asHomochoric = as_homochoric - asCubochoric = as_cubochoric ################################################################################################ # Static constructors. The input data needs to follow the conventions, options allow to @@ -505,13 +496,8 @@ class Rotation: # for compatibility (old names do not follow convention) fromQuaternion = from_quaternion fromEulers = from_Eulers - fromAxisAngle = from_axis_angle - fromBasis = from_basis - fromMatrix = from_matrix - fromRodrigues = from_Rodrigues - fromHomochoric = from_homochoric - fromCubochoric = from_cubochoric - fromRandom = from_random + asAxisAngle = as_axis_angle + #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations From a25dd1c4389165a74237706b0e83b4964892fb1d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 17 May 2020 08:01:34 +0200 Subject: [PATCH 13/32] cleaning --- python/damask/_rotation.py | 20 ++++++++++---------- python/tests/test_Orientation.py | 2 +- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index bf5a2a4f3..e2d5963fe 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -35,8 +35,8 @@ class Rotation: ----- Vector "a" (defined in coordinate system "A") is passively rotated resulting in new coordinates "b" when expressed in system "B". - b = Q * a - b = np.dot(Q.asMatrix(),a) + b = Q @ a + b = np.dot(Q.as_matrix(),a) """ @@ -73,8 +73,8 @@ class Rotation: raise NotImplementedError('Support for multiple rotations missing') return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), - 'Matrix:\n{}'.format(self.asMatrix()), - 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)), + 'Matrix:\n{}'.format(self.as_matrix()), + 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)), ]) @@ -112,9 +112,9 @@ class Rotation: return A*other + B*self.quaternion[1:] + C * np.cross(self.quaternion[1:],other) elif other.shape == (3,3,): - return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) + return np.dot(self.as_matrix(),np.dot(other,self.as_matrix().T)) elif other.shape == (3,3,3,3,): - R = self.asMatrix() + R = self.as_matrix() return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) else: raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') @@ -150,10 +150,10 @@ class Rotation: - p_m[...,(i+2)%3]*other[...,(i+1)%3])).reshape(self.shape+(1,)) for i in [0,1,2]]) if self.shape + (3,3) == other.shape: - R = self.asMatrix() + R = self.as_matrix() return np.einsum('...im,...jn,...mn',R,R,other) if self.shape + (3,3,3,3) == other.shape: - R = self.asMatrix() + R = self.as_matrix() return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) else: raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') @@ -654,8 +654,8 @@ class Rotation: def om2ax(om): """Rotation matrix to axis angle pair.""" diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], - om[...,2,0:1]-om[...,0,2:3], - om[...,0,1:2]-om[...,1,0:1] + om[...,2,0:1]-om[...,0,2:3], + om[...,0,1:2]-om[...,1,0:1] ]) diag_delta[np.abs(diag_delta)<1.e-6] = 1.0 t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,)) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 8dc526325..a8b8afdac 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -51,7 +51,7 @@ class TestOrientation: def test_relationship_reference(self,update,reference_dir,model,lattice): reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model)) ori = Orientation(Rotation(),lattice) - eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)]) + eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.relatedOrientations(model)]) if update: coords = np.array([(1,i+1) for i,x in enumerate(eu)]) table = damask.Table(eu,{'Eulers':(3,)}) From 3dc90ddb943a7be8d55416031af8e0281a413d66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 May 2020 15:54:58 +0200 Subject: [PATCH 14/32] use the formula from the paper, not from the reference implementation a few multiplications should be faster than a transpose --- python/damask/_rotation.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index e2d5963fe..032cd1ac2 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -534,16 +534,16 @@ class Rotation: def qu2om(qu): qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) om = np.block([qq + 2.0*qu[...,1:2]**2, - 2.0*(qu[...,2:3]*qu[...,1:2]+qu[...,0:1]*qu[...,3:4]), - 2.0*(qu[...,3:4]*qu[...,1:2]-qu[...,0:1]*qu[...,2:3]), - 2.0*(qu[...,1:2]*qu[...,2:3]-qu[...,0:1]*qu[...,3:4]), + 2.0*(qu[...,2:3]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,3:4]), + 2.0*(qu[...,3:4]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,2:3]), + 2.0*(qu[...,1:2]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,3:4]), qq + 2.0*qu[...,2:3]**2, - 2.0*(qu[...,3:4]*qu[...,2:3]+qu[...,0:1]*qu[...,1:2]), - 2.0*(qu[...,1:2]*qu[...,3:4]+qu[...,0:1]*qu[...,2:3]), - 2.0*(qu[...,2:3]*qu[...,3:4]-qu[...,0:1]*qu[...,1:2]), + 2.0*(qu[...,3:4]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,1:2]), + 2.0*(qu[...,1:2]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,2:3]), + 2.0*(qu[...,2:3]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,1:2]), qq + 2.0*qu[...,3:4]**2, ]).reshape(qu.shape[:-1]+(3,3)) - return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + return om @staticmethod def qu2eu(qu): @@ -641,7 +641,7 @@ class Rotation: np.zeros(om.shape[:-2]+(1,)), ]), np.block([np.arctan2(om[...,2,0:1]*zeta,-om[...,2,1:2]*zeta), - np.arccos(om[...,2,2:3]), + np.arccos( om[...,2,2:3]), np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta) ]) ) From b200894a40782420b624a2058ff40cbd5097e95d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 May 2020 18:41:48 +0200 Subject: [PATCH 15/32] bugfix and further test --- python/damask/_rotation.py | 41 +++++++++++++++++++++-------- python/tests/rotation_conversion.py | 17 ++++++++++++ python/tests/test_Rotation.py | 3 ++- 3 files changed, 49 insertions(+), 12 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 032cd1ac2..5af3e1874 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -534,13 +534,13 @@ class Rotation: def qu2om(qu): qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) om = np.block([qq + 2.0*qu[...,1:2]**2, - 2.0*(qu[...,2:3]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,3:4]), - 2.0*(qu[...,3:4]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,2:3]), - 2.0*(qu[...,1:2]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,3:4]), + 2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]), + 2.0*(qu[...,3:4]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,2:3]), + 2.0*(qu[...,1:2]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,3:4]), qq + 2.0*qu[...,2:3]**2, - 2.0*(qu[...,3:4]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,1:2]), - 2.0*(qu[...,1:2]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,2:3]), - 2.0*(qu[...,2:3]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,1:2]), + 2.0*(qu[...,3:4]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,1:2]), + 2.0*(qu[...,1:2]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,2:3]), + 2.0*(qu[...,2:3]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,1:2]), qq + 2.0*qu[...,3:4]**2, ]).reshape(qu.shape[:-1]+(3,3)) return om @@ -626,9 +626,31 @@ class Rotation: """ Rotation matrix to quaternion. - The original formulation (direct conversion) had (numerical?) issues + This formulation is from www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion. + The original formulation had issues. """ - return Rotation.eu2qu(Rotation.om2eu(om)) + def x(a): + trace = a[0,0] + a[1,1] + a[2,2] + if trace > 0: + s = 0.5 / np.sqrt(trace+ 1.0) + qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s]) + else: + if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ): + s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2]) + qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s]) + elif (a[1,1] > a[2,2]): + s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2]) + qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s]) + else: + s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) + qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) + return qu + if len(om.shape) > 2: + om_ = om.reshape(-1,3,3) + qu = np.array([x(o) for o in om_]).reshape(om.shape[:-2]+(4,)) + else: + qu = x(om) + return qu*np.array([1,_P,_P,_P]) @staticmethod def om2eu(om): @@ -649,7 +671,6 @@ class Rotation: eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) return eu - @staticmethod def om2ax(om): """Rotation matrix to axis angle pair.""" @@ -672,7 +693,6 @@ class Rotation: ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0] return ax - @staticmethod def om2ro(om): """Rotation matrix to Rodrigues-Frank vector.""" @@ -703,7 +723,6 @@ class Rotation: qu[qu[...,0]<0.0]*=-1 return qu - @staticmethod def eu2om(eu): """Bunge-Euler angles to rotation matrix.""" diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index 5ca10f443..ddd13495f 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -112,6 +112,23 @@ def qu2ho(qu): #---------- Rotation matrix ---------- +def om2qu(a): + trace = a[0,0] + a[1,1] + a[2,2] + if trace > 0: + s = 0.5 / np.sqrt(trace+ 1.0) + qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s]) + else: + if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ): + s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2]) + qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s]) + elif (a[1,1] > a[2,2]): + s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2]) + qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s]) + else: + s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) + qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) + return qu + def om2eu(om): """Rotation matrix to Bunge-Euler angles.""" if not np.isclose(np.abs(om[2,2]),1.0,1.e-4): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 4a6350fda..71aff04a9 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -196,7 +196,8 @@ class TestRotation: assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.om2eu,rotation_conversion.om2eu), + @pytest.mark.parametrize('vectorized, single',[(Rotation.om2qu,rotation_conversion.om2qu), + (Rotation.om2eu,rotation_conversion.om2eu), (Rotation.om2ax,rotation_conversion.om2ax)]) def test_matrix_vectorization(self,default,vectorized,single): om = np.array([rot.as_matrix() for rot in default]) From b6eebcd7044a2867b5b96aa97114b4a2a77cfc47 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 18 May 2020 20:16:50 +0200 Subject: [PATCH 16/32] small fixes (testing P=+1) --- python/damask/_rotation.py | 10 +++++----- python/tests/rotation_conversion.py | 16 +++++++++------- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 5af3e1874..e38a2bc56 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -716,10 +716,10 @@ class Rotation: ee = 0.5*eu cPhi = np.cos(ee[...,1:2]) sPhi = np.sin(ee[...,1:2]) - qu = np.block([ cPhi*np.cos(ee[...,0:1]+ee[...,2:3]), - -_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]), - -_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]), - -_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])]) + qu = np.block([ cPhi*np.cos(ee[...,0:1]+ee[...,2:3]), + -_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]), + -_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]), + -_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])]) qu[qu[...,0]<0.0]*=-1 return qu @@ -804,7 +804,7 @@ class Rotation: omc*ax[...,0:1]*ax[...,2:3] + s*ax[...,1:2], omc*ax[...,1:2]*ax[...,2:3] - s*ax[...,0:1], c+omc*ax[...,2:3]**2]).reshape(ax.shape[:-1]+(3,3)) - return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + return om if _P < 0.0 else np.swapaxes(om,-1,-2) @staticmethod def ax2eu(ax): diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index ddd13495f..82451bd84 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -31,12 +31,14 @@ import numpy as np -_P = -1 +from damask import _rotation + +_P = _rotation._P # parameters for conversion from/to cubochoric -_sc = np.pi**(1./6.)/6.**(1./6.) -_beta = np.pi**(5./6.)/6.**(1./6.)/2. -_R1 = (3.*np.pi/4.)**(1./3.) +_sc = _rotation._sc +_beta = _rotation._beta +_R1 = _rotation._R1 def iszero(a): return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) @@ -53,7 +55,7 @@ def qu2om(qu): om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) - return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + return om if _P < 0.0 else np.swapaxes(om,-1,-2) def qu2eu(qu): """Quaternion to Bunge-Euler angles.""" @@ -127,7 +129,7 @@ def om2qu(a): else: s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) - return qu + return qu*np.array([1.,_P,_P,_P]) def om2eu(om): """Rotation matrix to Bunge-Euler angles.""" @@ -167,7 +169,7 @@ def eu2qu(eu): ee = 0.5*eu cPhi = np.cos(ee[1]) sPhi = np.sin(ee[1]) - qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), + qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), -_P*sPhi*np.cos(ee[0]-ee[2]), -_P*sPhi*np.sin(ee[0]-ee[2]), -_P*cPhi*np.sin(ee[0]+ee[2]) ]) From 1a3a4a800e6f6b4616eb8b76dfb18473f63945e0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 07:35:58 +0200 Subject: [PATCH 17/32] vectorized --- python/damask/_rotation.py | 54 ++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index e38a2bc56..09921a194 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -629,28 +629,38 @@ class Rotation: This formulation is from www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion. The original formulation had issues. """ - def x(a): - trace = a[0,0] + a[1,1] + a[2,2] - if trace > 0: - s = 0.5 / np.sqrt(trace+ 1.0) - qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s]) - else: - if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ): - s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2]) - qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s]) - elif (a[1,1] > a[2,2]): - s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2]) - qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s]) - else: - s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) - qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) - return qu - if len(om.shape) > 2: - om_ = om.reshape(-1,3,3) - qu = np.array([x(o) for o in om_]).reshape(om.shape[:-2]+(4,)) - else: - qu = x(om) - return qu*np.array([1,_P,_P,_P]) + trace = om[...,0,0:1]+om[...,1,1:2]+om[...,2,2:3] + + with np.errstate(invalid='ignore',divide='ignore'): + s = [ + 0.5 / np.sqrt( 1.0 + trace), + 2.0 * np.sqrt( 1.0 + om[...,0,0:1] - om[...,1,1:2] - om[...,2,2:3]), + 2.0 * np.sqrt( 1.0 + om[...,1,1:2] - om[...,2,2:3] - om[...,0,0:1]), + 2.0 * np.sqrt( 1.0 + om[...,2,2:3] - om[...,0,0:1] - om[...,1,1:2] ) + ] + qu= np.where(trace>0, + np.block([0.25 / s[0], + (om[...,2,1:2] - om[...,1,2:3] ) * s[0], + (om[...,0,2:3] - om[...,2,0:1] ) * s[0], + (om[...,1,0:1] - om[...,0,1:2] ) * s[0]]), + np.where(om[...,0,0:1] > np.maximum(om[...,1,1:2],om[...,2,2:3]), + np.block([(om[...,2,1:2] - om[...,1,2:3]) / s[1], + 0.25 * s[1], + (om[...,0,1:2] + om[...,1,0:1]) / s[1], + (om[...,0,2:3] + om[...,2,0:1]) / s[1]]), + np.where(om[...,1,1:2] > om[...,2,2:3], + np.block([(om[...,0,2:3] - om[...,2,0:1]) / s[2], + (om[...,0,1:2] + om[...,1,0:1]) / s[2], + 0.25 * s[2], + (om[...,1,2:3] + om[...,2,1:2]) / s[2]]), + np.block([(om[...,1,0:1] - om[...,0,1:2]) / s[3], + (om[...,0,2:3] + om[...,2,0:1]) / s[3], + (om[...,1,2:3] + om[...,2,1:2]) / s[3], + 0.25 * s[3]]), + ) + ) + )*np.array([1,_P,_P,_P]) + return qu @staticmethod def om2eu(om): From ebf05a279eb2693a37dcce13efa590a67ab5bdc9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 08:59:03 +0200 Subject: [PATCH 18/32] standard name --- python/damask/_rotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 09921a194..949a6c3ab 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -311,7 +311,7 @@ class Rotation: # relax the conventions. @staticmethod def from_quaternion(quaternion, - acceptHomomorph = False, + accept_homomorph = False, P = -1): qu = np.array(quaternion,dtype=float) @@ -319,7 +319,7 @@ class Rotation: raise ValueError('Invalid shape.') if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 - if acceptHomomorph: + if accept_homomorph: qu[qu[...,0] < 0.0] *= -1 else: if np.any(qu[...,0] < 0.0): @@ -472,7 +472,7 @@ class Rotation: else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) - return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) + return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),accept_homomorph = True) @staticmethod def from_random(shape=None): From 065c624f94acdacce2fce12be0aca72c1cbe3c30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 09:10:14 +0200 Subject: [PATCH 19/32] fix for backward compatibility --- python/damask/_rotation.py | 1 + python/tests/test_Rotation.py | 36 +++++++++++++++++++++++------------ 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 949a6c3ab..f45e04ff4 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -494,6 +494,7 @@ class Rotation: # for compatibility (old names do not follow convention) + asM = M fromQuaternion = from_quaternion fromEulers = from_Eulers asAxisAngle = as_axis_angle diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 71aff04a9..183fbd19e 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -91,20 +91,25 @@ def reference_dir(reference_dir_base): class TestRotation: - def test_Eulers(self,default): + @pytest.mark.parametrize('degrees',[True,False]) + def test_Eulers(self,default,degrees): for rot in default: m = rot.as_quaternion() - o = Rotation.from_Eulers(rot.as_Eulers()).as_quaternion() + o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion() ok = np.allclose(m,o,atol=atol) if np.isclose(rot.as_quaternion()[0],0.0,atol=atol): ok = ok or np.allclose(m*-1.,o,atol=atol) print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) - def test_AxisAngle(self,default): + @pytest.mark.parametrize('P',[1,-1]) + @pytest.mark.parametrize('normalise',[True,False]) + @pytest.mark.parametrize('degrees',[True,False]) + def test_AxisAngle(self,default,degrees,normalise,P): + c = np.array([P*-1,P*-1,P*-1,1.]) for rot in default: m = rot.as_Eulers() - o = Rotation.from_axis_angle(rot.as_axis_angle()).as_Eulers() + o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalise,P).as_Eulers() u = np.array([np.pi*2,np.pi,np.pi*2]) ok = np.allclose(m,o,atol=atol) ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol) @@ -124,36 +129,43 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 - def test_Rodrigues(self,default): + @pytest.mark.parametrize('P',[1,-1]) + @pytest.mark.parametrize('normalise',[True,False]) + def test_Rodrigues(self,default,normalise,P): + c = np.array([P*-1,P*-1,P*-1,1.]) for rot in default: m = rot.as_matrix() - o = Rotation.from_Rodrigues(rot.as_Rodrigues()).as_matrix() + o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalise,P).as_matrix() ok = np.allclose(m,o,atol=atol) print(m,o) assert ok and np.isclose(np.linalg.det(o),1.0) - def test_Homochoric(self,default): + @pytest.mark.parametrize('P',[1,-1]) + def test_Homochoric(self,default,P): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: m = rot.as_Rodrigues() - o = Rotation.from_homochoric(rot.as_homochoric()).as_Rodrigues() + o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues() ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = ok or np.isclose(m[3],0.0,atol=atol) print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) - def test_Cubochoric(self,default): + @pytest.mark.parametrize('P',[1,-1]) + def test_Cubochoric(self,default,P): for rot in default: m = rot.as_homochoric() - o = Rotation.from_cubochoric(rot.as_cubochoric()).as_homochoric() + o = Rotation.from_cubochoric(rot.as_cubochoric()*P*-1,P).as_homochoric() ok = np.allclose(m,o,atol=atol) print(m,o,rot.as_quaternion()) assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 - def test_Quaternion(self,default): + @pytest.mark.parametrize('P',[1,-1]) + def test_Quaternion(self,default,P): + c = np.array([1,P*-1,P*-1,P*-1]) for rot in default: m = rot.as_cubochoric() - o = Rotation.from_quaternion(rot.as_quaternion()).as_cubochoric() + o = Rotation.from_quaternion(rot.as_quaternion()*c,False,P).as_cubochoric() ok = np.allclose(m,o,atol=atol) print(m,o,rot.as_quaternion()) assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 From de79a5af43772629878e0bafe6f928de569e289a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 09:55:23 +0200 Subject: [PATCH 20/32] directly test the internal conversions some work to do ... --- python/damask/_rotation.py | 5 +-- python/tests/rotation_conversion.py | 1 + python/tests/test_Rotation.py | 50 +++++++++++++++++++++++++++++ 3 files changed, 54 insertions(+), 2 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index f45e04ff4..167552b68 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -468,8 +468,8 @@ class Rotation: weights = np.ones(N,dtype='i') for i,(r,n) in enumerate(zip(rotations,weights)): - M = r.asM() * n if i == 0 \ - else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa + M = r.M() * n if i == 0 \ + else M + r.M() * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),accept_homomorph = True) @@ -661,6 +661,7 @@ class Rotation: ) ) )*np.array([1,_P,_P,_P]) + qu[qu[...,0]<0] *=-1 return qu @staticmethod diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index 82451bd84..5219f3260 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -129,6 +129,7 @@ def om2qu(a): else: s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) + if qu[0]<0: qu*=-1 return qu*np.array([1.,_P,_P,_P]) def om2eu(om): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 183fbd19e..4dc98a797 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -91,6 +91,56 @@ def reference_dir(reference_dir_base): class TestRotation: + @pytest.mark.parametrize('forward,backward',[(Rotation.qu2om,Rotation.om2qu), + (Rotation.qu2eu,Rotation.eu2qu), + (Rotation.qu2ax,Rotation.ax2qu), + (Rotation.qu2ro,Rotation.ro2qu), + (Rotation.qu2ho,Rotation.ho2qu), + (Rotation.qu2cu,Rotation.cu2qu)]) + def test_quaternion_internal(self,default,forward,backward): + for rot in default: + m = rot.as_quaternion() + o = backward(forward(m)) + ok = np.allclose(m,o,atol=atol) + if np.isclose(rot.as_quaternion()[0],0.0,atol=atol): + ok = ok or np.allclose(m*-1.,o,atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and np.isclose(np.linalg.norm(o),1.0) + + @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om)]) + #(Rotation.om2eu,Rotation.eu2om), + #(Rotation.om2ax,Rotation.ax2om), + #(Rotation.om2ro,Rotation.ro2om), + #(Rotation.om2ho,Rotation.ho2om), + #(Rotation.om2cu,Rotation.cu2om)]) + def test_matrix_internal(self,default,forward,backward): + for rot in default: + m = rot.as_matrix() + o = backward(forward(m)) + ok = np.allclose(m,o,atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and np.isclose(np.linalg.det(o),1.0) + + @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu)]) + #(Rotation.eu2om,Rotation.om2eu), + #(Rotation.eu2ax,Rotation.ax2eu), + #(Rotation.eu2ro,Rotation.ro2eu), + #(Rotation.eu2ho,Rotation.ho2eu), + #(Rotation.eu2cu,Rotation.cu2eu)]) + def test_Eulers_internal(self,default,forward,backward): + for rot in default: + m = rot.as_Eulers() + o = backward(forward(m)) + u = np.array([np.pi*2,np.pi,np.pi*2]) + ok = np.allclose(m,o,atol=atol) + ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol) + if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol): + sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]]) + ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() + + @pytest.mark.parametrize('degrees',[True,False]) def test_Eulers(self,default,degrees): for rot in default: From a6b0aaffbac0d1188e85f5862b0e3c7452405dc4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 19:41:50 +0200 Subject: [PATCH 21/32] more tests Hotfix needed for axis angle to matrix (not used in DAMASK, needs further investigation) --- python/damask/_rotation.py | 11 +++++---- python/tests/rotation_conversion.py | 11 +++++---- python/tests/test_Rotation.py | 38 ++++++++++++++++++++--------- 3 files changed, 39 insertions(+), 21 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 167552b68..ee497e72b 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -669,7 +669,7 @@ class Rotation: """Rotation matrix to Bunge-Euler angles.""" with np.errstate(invalid='ignore',divide='ignore'): zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) - eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4), + eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-9), np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]), np.pi*0.5*(1-om[...,2,2:3]), np.zeros(om.shape[:-2]+(1,)), @@ -679,30 +679,30 @@ class Rotation: np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta) ]) ) - eu[np.abs(eu)<1.e-6] = 0.0 + eu[np.abs(eu)<1.e-8] = 0.0 eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) return eu @staticmethod def om2ax(om): """Rotation matrix to axis angle pair.""" + return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], om[...,2,0:1]-om[...,0,2:3], om[...,0,1:2]-om[...,1,0:1] ]) - diag_delta[np.abs(diag_delta)<1.e-6] = 1.0 t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,)) w,vr = np.linalg.eig(om) # mask duplicated real eigenvalues w[np.isclose(w[...,0],1.0+0.0j),1:] = 0. w[np.isclose(w[...,1],1.0+0.0j),2:] = 0. vr = np.swapaxes(vr,-1,-2) - ax = np.where(np.abs(diag_delta)<0, + ax = np.where(np.abs(diag_delta)<1e-12, np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)), np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \ *np.sign(diag_delta)) ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))]) - ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0] + ax[np.abs(ax[...,3])<1.e-8] = [ 0.0, 0.0, 1.0, 0.0] return ax @staticmethod @@ -804,6 +804,7 @@ class Rotation: @staticmethod def ax2om(ax): """Axis angle pair to rotation matrix.""" + return Rotation.qu2om(Rotation.ax2qu(ax)) # HOTFIX c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) omc = 1. -c diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index 5219f3260..ab503f402 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -134,25 +134,26 @@ def om2qu(a): def om2eu(om): """Rotation matrix to Bunge-Euler angles.""" - if not np.isclose(np.abs(om[2,2]),1.0,1.e-4): + if not np.isclose(np.abs(om[2,2]),1.0,1.e-9): zeta = 1.0/np.sqrt(1.0-om[2,2]**2) eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), np.arccos(om[2,2]), np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) else: eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation - eu[np.abs(eu)<1.e-6] = 0.0 + eu[np.abs(eu)<1.e-8] = 0.0 eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) return eu def om2ax(om): """Rotation matrix to axis angle pair.""" + return qu2ax(om2qu(om)) # HOTFIX ax=np.empty(4) # first get the rotation angle t = 0.5*(om.trace() -1.0) ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - if np.abs(ax[3])<1.e-6: + if np.abs(ax[3])<1.e-8: ax = np.array([ 0.0, 0.0, 1.0, 0.0]) else: w,vr = np.linalg.eig(om) @@ -160,8 +161,7 @@ def om2ax(om): i = np.where(np.isclose(w,1.0+0.0j))[0][0] ax[0:3] = np.real(vr[0:3,i]) diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - diagDelta[np.abs(diagDelta)<1.e-6] = 1.0 - ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) + ax[0:3] = np.where(np.abs(diagDelta)<1e-12, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) return ax #---------- Bunge-Euler angles ---------- @@ -229,6 +229,7 @@ def ax2qu(ax): def ax2om(ax): """Axis angle pair to rotation matrix.""" + return qu2om(ax2qu(ax)) # HOTFIX c = np.cos(ax[3]) s = np.sin(ax[3]) omc = 1.0-c diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 4dc98a797..adf507e49 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -107,10 +107,10 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om)]) - #(Rotation.om2eu,Rotation.eu2om), - #(Rotation.om2ax,Rotation.ax2om), - #(Rotation.om2ro,Rotation.ro2om), + @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om), + (Rotation.om2eu,Rotation.eu2om), + (Rotation.om2ax,Rotation.ax2om), + (Rotation.om2ro,Rotation.ro2om)]) #(Rotation.om2ho,Rotation.ho2om), #(Rotation.om2cu,Rotation.cu2om)]) def test_matrix_internal(self,default,forward,backward): @@ -121,12 +121,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.det(o),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu)]) - #(Rotation.eu2om,Rotation.om2eu), - #(Rotation.eu2ax,Rotation.ax2eu), - #(Rotation.eu2ro,Rotation.ro2eu), - #(Rotation.eu2ho,Rotation.ho2eu), - #(Rotation.eu2cu,Rotation.cu2eu)]) + @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu), + (Rotation.eu2om,Rotation.om2eu), + (Rotation.eu2ax,Rotation.ax2eu), + (Rotation.eu2ro,Rotation.ro2eu), + (Rotation.eu2ho,Rotation.ho2eu), + (Rotation.eu2cu,Rotation.cu2eu)]) def test_Eulers_internal(self,default,forward,backward): for rot in default: m = rot.as_Eulers() @@ -140,6 +140,22 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() + @pytest.mark.parametrize('forward,backward',[(Rotation.ax2qu,Rotation.qu2ax), + (Rotation.ax2om,Rotation.om2ax), + (Rotation.ax2eu,Rotation.eu2ax), + (Rotation.ax2ro,Rotation.ro2ax), + (Rotation.ax2ho,Rotation.ho2ax), + (Rotation.ax2cu,Rotation.cu2ax)]) + def test_axis_angle_internal(self,default,forward,backward): + for rot in default: + m = rot.as_axis_angle() + o = backward(forward(m)) + ok = np.allclose(m,o,atol=atol) + if np.isclose(m[3],np.pi,atol=atol): + ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7 + @pytest.mark.parametrize('degrees',[True,False]) def test_Eulers(self,default,degrees): @@ -177,7 +193,7 @@ class TestRotation: if np.isclose(m[3],np.pi,atol=atol): ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) print(m,o,rot.as_quaternion()) - assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7 @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('normalise',[True,False]) From 1c53a37de485420da8ca8fa6cbdea6612b60824d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 21:03:30 +0200 Subject: [PATCH 22/32] more tests and adjustments to tolerances --- python/damask/_rotation.py | 5 ++--- python/tests/rotation_conversion.py | 7 +++---- python/tests/test_Rotation.py | 6 +++--- 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index ee497e72b..0f0a081bc 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -584,10 +584,10 @@ class Rotation: with np.errstate(invalid='ignore',divide='ignore'): s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) - ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-6,qu.shape), + ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape), np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), np.block([qu[...,1:4]*s,omega])) - ax[np.sum(np.abs(qu[...,1:4])**2,axis=-1) < 1.0e-6,] = [0.0, 0.0, 1.0, 0.0] + ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0] return ax @staticmethod @@ -804,7 +804,6 @@ class Rotation: @staticmethod def ax2om(ax): """Axis angle pair to rotation matrix.""" - return Rotation.qu2om(Rotation.ax2qu(ax)) # HOTFIX c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) omc = 1. -c diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index ab503f402..3aca1a212 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -81,9 +81,9 @@ def qu2ax(qu): Modified version of the original formulation, should be numerically more stable """ - if np.abs(np.sum(qu[1:4]**2)) < 1.e-6: # set axis to [001] if the angle is 0/360 + if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - elif qu[0] > 1.e-6: + elif qu[0] > 1.e-8: s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ]) @@ -229,7 +229,6 @@ def ax2qu(ax): def ax2om(ax): """Axis angle pair to rotation matrix.""" - return qu2om(ax2qu(ax)) # HOTFIX c = np.cos(ax[3]) s = np.sin(ax[3]) omc = 1.0-c @@ -239,7 +238,7 @@ def ax2om(ax): q = omc*ax[idx[0]] * ax[idx[1]] om[idx[0],idx[1]] = q + s*ax[idx[2]] om[idx[1],idx[0]] = q - s*ax[idx[2]] - return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + return om if _P < 0.0 else om.T def ax2ro(ax): """Axis angle pair to Rodrigues-Frank vector.""" diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index adf507e49..f9565a5e2 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -110,9 +110,9 @@ class TestRotation: @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om), (Rotation.om2eu,Rotation.eu2om), (Rotation.om2ax,Rotation.ax2om), - (Rotation.om2ro,Rotation.ro2om)]) - #(Rotation.om2ho,Rotation.ho2om), - #(Rotation.om2cu,Rotation.cu2om)]) + (Rotation.om2ro,Rotation.ro2om), + (Rotation.om2ho,Rotation.ho2om), + (Rotation.om2cu,Rotation.cu2om)]) def test_matrix_internal(self,default,forward,backward): for rot in default: m = rot.as_matrix() From 3e002691793b36590c6bce536cffc2a6b8d01260 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 23:24:12 +0200 Subject: [PATCH 23/32] again some more tests --- python/damask/_rotation.py | 8 ++++++-- python/tests/test_Rotation.py | 33 +++++++++++++++++++++++++++++++-- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 0f0a081bc..fc124f5bb 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -112,10 +112,14 @@ class Rotation: return A*other + B*self.quaternion[1:] + C * np.cross(self.quaternion[1:],other) elif other.shape == (3,3,): - return np.dot(self.as_matrix(),np.dot(other,self.as_matrix().T)) + R = self.as_matrix() + return np.dot(R,np.dot(other,R.T)) elif other.shape == (3,3,3,3,): R = self.as_matrix() - return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) + RR = np.outer(R, R) + RRRR = np.outer(RR, RR).reshape(4 * (3,3)) + axes = ((0, 2, 4, 6), (0, 1, 2, 3)) + return np.tensordot(RRRR, other, axes) else: raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') else: diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index f9565a5e2..35df97f7b 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -154,8 +154,37 @@ class TestRotation: if np.isclose(m[3],np.pi,atol=atol): ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) print(m,o,rot.as_quaternion()) - assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7 + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-9 + @pytest.mark.parametrize('forward,backward',[(Rotation.ro2qu,Rotation.qu2ro), + #(Rotation.ro2om,Rotation.om2ro), + #(Rotation.ro2eu,Rotation.eu2ro), + (Rotation.ro2ax,Rotation.ax2ro), + (Rotation.ro2ho,Rotation.ho2ro), + (Rotation.ro2cu,Rotation.cu2ro)]) + def test_Rodrigues_internal(self,default,forward,backward): + cutoff = np.tan(np.pi*.5*(1.-1e-4)) + for rot in default: + m = rot.as_Rodrigues() + o = backward(forward(m)) + ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) + ok = ok or np.isclose(m[3],0.0,atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) + + @pytest.mark.parametrize('forward,backward',[(Rotation.ho2qu,Rotation.qu2ho), + (Rotation.ho2om,Rotation.om2ho), + #(Rotation.ho2eu,Rotation.eu2ho), + (Rotation.ho2ax,Rotation.ax2ho), + (Rotation.ho2ro,Rotation.ro2ho), + (Rotation.ho2cu,Rotation.cu2ho)]) + def test_homochoric_internal(self,default,forward,backward): + for rot in default: + m = rot.as_homochoric() + o = backward(forward(m)) + ok = np.allclose(m,o,atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 @pytest.mark.parametrize('degrees',[True,False]) def test_Eulers(self,default,degrees): @@ -193,7 +222,7 @@ class TestRotation: if np.isclose(m[3],np.pi,atol=atol): ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) print(m,o,rot.as_quaternion()) - assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7 + assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-9 @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('normalise',[True,False]) From 1fa4a07bb857fe3f6cd2da04ae88786caaa1ca81 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 19 May 2020 23:39:01 +0200 Subject: [PATCH 24/32] for backward compatibility --- python/damask/_rotation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index fc124f5bb..31732015e 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -502,6 +502,7 @@ class Rotation: fromQuaternion = from_quaternion fromEulers = from_Eulers asAxisAngle = as_axis_angle + asRodrigues = as_Rodrigues #################################################################################################### From 1d8903bb0c26f8094ac859936c9089f2985891f6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 00:17:45 +0200 Subject: [PATCH 25/32] transition code --- python/damask/_rotation.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 31732015e..d2698fd52 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -316,8 +316,11 @@ class Rotation: @staticmethod def from_quaternion(quaternion, accept_homomorph = False, - P = -1): + P = -1, + acceptHomomorph = None): + if acceptHomomorph is not None: + accept_homomorph = acceptHomomorph qu = np.array(quaternion,dtype=float) if qu.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') From 9694767997933e55f8ae41e08b5282ac838a6677 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 09:01:38 +0200 Subject: [PATCH 26/32] all seems to work now --- python/damask/_rotation.py | 8 ++++---- python/tests/rotation_conversion.py | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index d2698fd52..baa008391 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -694,7 +694,7 @@ class Rotation: @staticmethod def om2ax(om): """Rotation matrix to axis angle pair.""" - return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX + #return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], om[...,2,0:1]-om[...,0,2:3], om[...,0,1:2]-om[...,1,0:1] @@ -878,14 +878,14 @@ class Rotation: ax = np.where(np.isfinite(ro[...,3:4]), np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]), np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)])) - ax[np.abs(ro[...,3]) < 1.e-6] = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + ax[np.abs(ro[...,3]) < 1.e-8] = np.array([ 0.0, 0.0, 1.0, 0.0 ]) return ax @staticmethod def ro2ho(ro): """Rodrigues-Frank vector to homochoric vector.""" f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) - ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape), + ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape), np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0)) return ho @@ -929,7 +929,7 @@ class Rotation: hm *= hmag_squared s += tfit[i] * hm with np.errstate(invalid='ignore'): - ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-6,ho.shape[:-1]+(4,)), + ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-8,ho.shape[:-1]+(4,)), [ 0.0, 0.0, 1.0, 0.0 ], np.block([ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))])) return ax diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index 3aca1a212..58ac9329f 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -147,7 +147,7 @@ def om2eu(om): def om2ax(om): """Rotation matrix to axis angle pair.""" - return qu2ax(om2qu(om)) # HOTFIX + #return qu2ax(om2qu(om)) # HOTFIX ax=np.empty(4) # first get the rotation angle @@ -262,7 +262,7 @@ def ax2ho(ax): #---------- Rodrigues-Frank vector ---------- def ro2ax(ro): """Rodrigues-Frank vector to axis angle pair.""" - if np.abs(ro[3]) < 1.e-6: + if np.abs(ro[3]) < 1.e-8: ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) elif not np.isfinite(ro[3]): ax = np.array([ ro[0], ro[1], ro[2], np.pi ]) @@ -274,7 +274,7 @@ def ro2ax(ro): def ro2ho(ro): """Rodrigues-Frank vector to homochoric vector.""" - if np.sum(ro[0:3]**2.0) < 1.e-6: + if np.sum(ro[0:3]**2.0) < 1.e-8: ho = np.zeros(3) else: f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi From 9240dd59b22bb4bb3f7736d13ae918183b1243cd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 11:11:07 +0200 Subject: [PATCH 27/32] mark as internal all functionality (+ sanity checks) can be done with the class functionality --- python/damask/_rotation.py | 156 +++++++++++++++++----------------- python/tests/test_Rotation.py | 114 ++++++++++++------------- 2 files changed, 135 insertions(+), 135 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index baa008391..cf6eaed34 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -249,7 +249,7 @@ class Rotation: return angles in degrees. """ - eu = Rotation.qu2eu(self.quaternion) + eu = Rotation._qu2eu(self.quaternion) if degrees: eu = np.degrees(eu) return eu @@ -267,13 +267,13 @@ class Rotation: return tuple of axis and angle. """ - ax = Rotation.qu2ax(self.quaternion) + ax = Rotation._qu2ax(self.quaternion) if degrees: ax[...,3] = np.degrees(ax[...,3]) return (ax[...,:3],ax[...,3]) if pair else ax def as_matrix(self): """Rotation matrix.""" - return Rotation.qu2om(self.quaternion) + return Rotation._qu2om(self.quaternion) def as_Rodrigues(self, vector = False): @@ -286,16 +286,16 @@ class Rotation: return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). """ - ro = Rotation.qu2ro(self.quaternion) + ro = Rotation._qu2ro(self.quaternion) return ro[...,:3]*ro[...,3] if vector else ro def as_homochoric(self): """Homochoric vector: (h_1, h_2, h_3).""" - return Rotation.qu2ho(self.quaternion) + return Rotation._qu2ho(self.quaternion) def as_cubochoric(self): """Cubochoric vector: (c_1, c_2, c_3).""" - return Rotation.qu2cu(self.quaternion) + return Rotation._qu2cu(self.quaternion) def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M """ @@ -348,7 +348,7 @@ class Rotation: if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].') - return Rotation(Rotation.eu2qu(eu)) + return Rotation(Rotation._eu2qu(eu)) @staticmethod def from_axis_angle(axis_angle, @@ -368,7 +368,7 @@ class Rotation: if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)): raise ValueError('Axis angle rotation axis is not of unit length.') - return Rotation(Rotation.ax2qu(ax)) + return Rotation(Rotation._ax2qu(ax)) @staticmethod def from_basis(basis, @@ -392,7 +392,7 @@ class Rotation: or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)): raise ValueError('Orientation matrix is not orthogonal.') - return Rotation(Rotation.om2qu(om)) + return Rotation(Rotation._om2qu(om)) @staticmethod def from_matrix(om): @@ -415,7 +415,7 @@ class Rotation: if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): raise ValueError('Rodrigues vector rotation axis is not of unit length.') - return Rotation(Rotation.ro2qu(ro)) + return Rotation(Rotation._ro2qu(ro)) @staticmethod def from_homochoric(homochoric, @@ -430,7 +430,7 @@ class Rotation: if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9): raise ValueError('Homochoric coordinate outside of the sphere.') - return Rotation(Rotation.ho2qu(ho)) + return Rotation(Rotation._ho2qu(ho)) @staticmethod def from_cubochoric(cubochoric, @@ -443,10 +443,10 @@ class Rotation: if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) - ho = Rotation.cu2ho(cu) + ho = Rotation._cu2ho(cu) if P > 0: ho *= -1 # convert from P=1 to P=-1 - return Rotation(Rotation.ho2qu(ho)) + return Rotation(Rotation._ho2qu(ho)) @staticmethod @@ -540,7 +540,7 @@ class Rotation: #################################################################################################### #---------- Quaternion ---------- @staticmethod - def qu2om(qu): + def _qu2om(qu): qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2) om = np.block([qq + 2.0*qu[...,1:2]**2, 2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]), @@ -555,7 +555,7 @@ class Rotation: return om @staticmethod - def qu2eu(qu): + def _qu2eu(qu): """Quaternion to Bunge-Euler angles.""" q02 = qu[...,0:1]*qu[...,2:3] q13 = qu[...,1:2]*qu[...,3:4] @@ -583,7 +583,7 @@ class Rotation: return eu @staticmethod - def qu2ax(qu): + def _qu2ax(qu): """ Quaternion to axis angle pair. @@ -599,7 +599,7 @@ class Rotation: return ax @staticmethod - def qu2ro(qu): + def _qu2ro(qu): """Quaternion to Rodrigues-Frank vector.""" with np.errstate(invalid='ignore',divide='ignore'): s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) @@ -613,7 +613,7 @@ class Rotation: return ro @staticmethod - def qu2ho(qu): + def _qu2ho(qu): """Quaternion to homochoric vector.""" with np.errstate(invalid='ignore'): omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) @@ -624,14 +624,14 @@ class Rotation: return ho @staticmethod - def qu2cu(qu): + def _qu2cu(qu): """Quaternion to cubochoric vector.""" - return Rotation.ho2cu(Rotation.qu2ho(qu)) + return Rotation._ho2cu(Rotation._qu2ho(qu)) #---------- Rotation matrix ---------- @staticmethod - def om2qu(om): + def _om2qu(om): """ Rotation matrix to quaternion. @@ -673,7 +673,7 @@ class Rotation: return qu @staticmethod - def om2eu(om): + def _om2eu(om): """Rotation matrix to Bunge-Euler angles.""" with np.errstate(invalid='ignore',divide='ignore'): zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) @@ -692,9 +692,9 @@ class Rotation: return eu @staticmethod - def om2ax(om): + def _om2ax(om): """Rotation matrix to axis angle pair.""" - #return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX + #return Rotation._qu2ax(Rotation._om2qu(om)) # HOTFIX diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], om[...,2,0:1]-om[...,0,2:3], om[...,0,1:2]-om[...,1,0:1] @@ -714,24 +714,24 @@ class Rotation: return ax @staticmethod - def om2ro(om): + def _om2ro(om): """Rotation matrix to Rodrigues-Frank vector.""" - return Rotation.eu2ro(Rotation.om2eu(om)) + return Rotation._eu2ro(Rotation._om2eu(om)) @staticmethod - def om2ho(om): + def _om2ho(om): """Rotation matrix to homochoric vector.""" - return Rotation.ax2ho(Rotation.om2ax(om)) + return Rotation._ax2ho(Rotation._om2ax(om)) @staticmethod - def om2cu(om): + def _om2cu(om): """Rotation matrix to cubochoric vector.""" - return Rotation.ho2cu(Rotation.om2ho(om)) + return Rotation._ho2cu(Rotation._om2ho(om)) #---------- Bunge-Euler angles ---------- @staticmethod - def eu2qu(eu): + def _eu2qu(eu): """Bunge-Euler angles to quaternion.""" ee = 0.5*eu cPhi = np.cos(ee[...,1:2]) @@ -744,7 +744,7 @@ class Rotation: return qu @staticmethod - def eu2om(eu): + def _eu2om(eu): """Bunge-Euler angles to rotation matrix.""" c = np.cos(eu) s = np.sin(eu) @@ -762,7 +762,7 @@ class Rotation: return om @staticmethod - def eu2ax(eu): + def _eu2ax(eu): """Bunge-Euler angles to axis angle pair.""" t = np.tan(eu[...,1:2]*0.5) sigma = 0.5*(eu[...,0:1]+eu[...,2:3]) @@ -781,28 +781,28 @@ class Rotation: return ax @staticmethod - def eu2ro(eu): + def _eu2ro(eu): """Bunge-Euler angles to Rodrigues-Frank vector.""" - ax = Rotation.eu2ax(eu) + ax = Rotation._eu2ax(eu) ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)]) ro[ax[...,3]>=np.pi,3] = np.inf ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, _P, 0.0 ] return ro @staticmethod - def eu2ho(eu): + def _eu2ho(eu): """Bunge-Euler angles to homochoric vector.""" - return Rotation.ax2ho(Rotation.eu2ax(eu)) + return Rotation._ax2ho(Rotation._eu2ax(eu)) @staticmethod - def eu2cu(eu): + def _eu2cu(eu): """Bunge-Euler angles to cubochoric vector.""" - return Rotation.ho2cu(Rotation.eu2ho(eu)) + return Rotation._ho2cu(Rotation._eu2ho(eu)) #---------- Axis angle pair ---------- @staticmethod - def ax2qu(ax): + def _ax2qu(ax): """Axis angle pair to quaternion.""" c = np.cos(ax[...,3:4]*.5) s = np.sin(ax[...,3:4]*.5) @@ -810,7 +810,7 @@ class Rotation: return qu @staticmethod - def ax2om(ax): + def _ax2om(ax): """Axis angle pair to rotation matrix.""" c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) @@ -827,12 +827,12 @@ class Rotation: return om if _P < 0.0 else np.swapaxes(om,-1,-2) @staticmethod - def ax2eu(ax): + def _ax2eu(ax): """Rotation matrix to Bunge Euler angles.""" - return Rotation.om2eu(Rotation.ax2om(ax)) + return Rotation._om2eu(Rotation._ax2om(ax)) @staticmethod - def ax2ro(ax): + def _ax2ro(ax): """Axis angle pair to Rodrigues-Frank vector.""" ro = np.block([ax[...,:3], np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), @@ -843,36 +843,36 @@ class Rotation: return ro @staticmethod - def ax2ho(ax): + def _ax2ho(ax): """Axis angle pair to homochoric vector.""" f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) ho = ax[...,:3] * f return ho @staticmethod - def ax2cu(ax): + def _ax2cu(ax): """Axis angle pair to cubochoric vector.""" - return Rotation.ho2cu(Rotation.ax2ho(ax)) + return Rotation._ho2cu(Rotation._ax2ho(ax)) #---------- Rodrigues-Frank vector ---------- @staticmethod - def ro2qu(ro): + def _ro2qu(ro): """Rodrigues-Frank vector to quaternion.""" - return Rotation.ax2qu(Rotation.ro2ax(ro)) + return Rotation._ax2qu(Rotation._ro2ax(ro)) @staticmethod - def ro2om(ro): + def _ro2om(ro): """Rodgrigues-Frank vector to rotation matrix.""" - return Rotation.ax2om(Rotation.ro2ax(ro)) + return Rotation._ax2om(Rotation._ro2ax(ro)) @staticmethod - def ro2eu(ro): + def _ro2eu(ro): """Rodrigues-Frank vector to Bunge-Euler angles.""" - return Rotation.om2eu(Rotation.ro2om(ro)) + return Rotation._om2eu(Rotation._ro2om(ro)) @staticmethod - def ro2ax(ro): + def _ro2ax(ro): """Rodrigues-Frank vector to axis angle pair.""" with np.errstate(invalid='ignore',divide='ignore'): ax = np.where(np.isfinite(ro[...,3:4]), @@ -882,7 +882,7 @@ class Rotation: return ax @staticmethod - def ro2ho(ro): + def _ro2ho(ro): """Rodrigues-Frank vector to homochoric vector.""" f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-8,ro[...,0:3].shape), @@ -890,29 +890,29 @@ class Rotation: return ho @staticmethod - def ro2cu(ro): + def _ro2cu(ro): """Rodrigues-Frank vector to cubochoric vector.""" - return Rotation.ho2cu(Rotation.ro2ho(ro)) + return Rotation._ho2cu(Rotation._ro2ho(ro)) #---------- Homochoric vector---------- @staticmethod - def ho2qu(ho): + def _ho2qu(ho): """Homochoric vector to quaternion.""" - return Rotation.ax2qu(Rotation.ho2ax(ho)) + return Rotation._ax2qu(Rotation._ho2ax(ho)) @staticmethod - def ho2om(ho): + def _ho2om(ho): """Homochoric vector to rotation matrix.""" - return Rotation.ax2om(Rotation.ho2ax(ho)) + return Rotation._ax2om(Rotation._ho2ax(ho)) @staticmethod - def ho2eu(ho): + def _ho2eu(ho): """Homochoric vector to Bunge-Euler angles.""" - return Rotation.ax2eu(Rotation.ho2ax(ho)) + return Rotation._ax2eu(Rotation._ho2ax(ho)) @staticmethod - def ho2ax(ho): + def _ho2ax(ho): """Homochoric vector to axis angle pair.""" tfit = np.array([+1.0000000000018852, -0.5000000002194847, -0.024999992127593126, -0.003928701544781374, @@ -935,12 +935,12 @@ class Rotation: return ax @staticmethod - def ho2ro(ho): + def _ho2ro(ho): """Axis angle pair to Rodrigues-Frank vector.""" - return Rotation.ax2ro(Rotation.ho2ax(ho)) + return Rotation._ax2ro(Rotation._ho2ax(ho)) @staticmethod - def ho2cu(ho): + def _ho2cu(ho): """ Homochoric vector to cubochoric vector. @@ -980,32 +980,32 @@ class Rotation: #---------- Cubochoric ---------- @staticmethod - def cu2qu(cu): + def _cu2qu(cu): """Cubochoric vector to quaternion.""" - return Rotation.ho2qu(Rotation.cu2ho(cu)) + return Rotation._ho2qu(Rotation._cu2ho(cu)) @staticmethod - def cu2om(cu): + def _cu2om(cu): """Cubochoric vector to rotation matrix.""" - return Rotation.ho2om(Rotation.cu2ho(cu)) + return Rotation._ho2om(Rotation._cu2ho(cu)) @staticmethod - def cu2eu(cu): + def _cu2eu(cu): """Cubochoric vector to Bunge-Euler angles.""" - return Rotation.ho2eu(Rotation.cu2ho(cu)) + return Rotation._ho2eu(Rotation._cu2ho(cu)) @staticmethod - def cu2ax(cu): + def _cu2ax(cu): """Cubochoric vector to axis angle pair.""" - return Rotation.ho2ax(Rotation.cu2ho(cu)) + return Rotation._ho2ax(Rotation._cu2ho(cu)) @staticmethod - def cu2ro(cu): + def _cu2ro(cu): """Cubochoric vector to Rodrigues-Frank vector.""" - return Rotation.ho2ro(Rotation.cu2ho(cu)) + return Rotation._ho2ro(Rotation._cu2ho(cu)) @staticmethod - def cu2ho(cu): + def _cu2ho(cu): """ Cubochoric vector to homochoric vector. diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 35df97f7b..4e48fbb7c 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -91,12 +91,12 @@ def reference_dir(reference_dir_base): class TestRotation: - @pytest.mark.parametrize('forward,backward',[(Rotation.qu2om,Rotation.om2qu), - (Rotation.qu2eu,Rotation.eu2qu), - (Rotation.qu2ax,Rotation.ax2qu), - (Rotation.qu2ro,Rotation.ro2qu), - (Rotation.qu2ho,Rotation.ho2qu), - (Rotation.qu2cu,Rotation.cu2qu)]) + @pytest.mark.parametrize('forward,backward',[(Rotation._qu2om,Rotation._om2qu), + (Rotation._qu2eu,Rotation._eu2qu), + (Rotation._qu2ax,Rotation._ax2qu), + (Rotation._qu2ro,Rotation._ro2qu), + (Rotation._qu2ho,Rotation._ho2qu), + (Rotation._qu2cu,Rotation._cu2qu)]) def test_quaternion_internal(self,default,forward,backward): for rot in default: m = rot.as_quaternion() @@ -107,12 +107,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om), - (Rotation.om2eu,Rotation.eu2om), - (Rotation.om2ax,Rotation.ax2om), - (Rotation.om2ro,Rotation.ro2om), - (Rotation.om2ho,Rotation.ho2om), - (Rotation.om2cu,Rotation.cu2om)]) + @pytest.mark.parametrize('forward,backward',[(Rotation._om2qu,Rotation._qu2om), + (Rotation._om2eu,Rotation._eu2om), + (Rotation._om2ax,Rotation._ax2om), + (Rotation._om2ro,Rotation._ro2om), + (Rotation._om2ho,Rotation._ho2om), + (Rotation._om2cu,Rotation._cu2om)]) def test_matrix_internal(self,default,forward,backward): for rot in default: m = rot.as_matrix() @@ -121,12 +121,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.det(o),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu), - (Rotation.eu2om,Rotation.om2eu), - (Rotation.eu2ax,Rotation.ax2eu), - (Rotation.eu2ro,Rotation.ro2eu), - (Rotation.eu2ho,Rotation.ho2eu), - (Rotation.eu2cu,Rotation.cu2eu)]) + @pytest.mark.parametrize('forward,backward',[(Rotation._eu2qu,Rotation._qu2eu), + (Rotation._eu2om,Rotation._om2eu), + (Rotation._eu2ax,Rotation._ax2eu), + (Rotation._eu2ro,Rotation._ro2eu), + (Rotation._eu2ho,Rotation._ho2eu), + (Rotation._eu2cu,Rotation._cu2eu)]) def test_Eulers_internal(self,default,forward,backward): for rot in default: m = rot.as_Eulers() @@ -140,12 +140,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() - @pytest.mark.parametrize('forward,backward',[(Rotation.ax2qu,Rotation.qu2ax), - (Rotation.ax2om,Rotation.om2ax), - (Rotation.ax2eu,Rotation.eu2ax), - (Rotation.ax2ro,Rotation.ro2ax), - (Rotation.ax2ho,Rotation.ho2ax), - (Rotation.ax2cu,Rotation.cu2ax)]) + @pytest.mark.parametrize('forward,backward',[(Rotation._ax2qu,Rotation._qu2ax), + (Rotation._ax2om,Rotation._om2ax), + (Rotation._ax2eu,Rotation._eu2ax), + (Rotation._ax2ro,Rotation._ro2ax), + (Rotation._ax2ho,Rotation._ho2ax), + (Rotation._ax2cu,Rotation._cu2ax)]) def test_axis_angle_internal(self,default,forward,backward): for rot in default: m = rot.as_axis_angle() @@ -156,12 +156,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-9 - @pytest.mark.parametrize('forward,backward',[(Rotation.ro2qu,Rotation.qu2ro), - #(Rotation.ro2om,Rotation.om2ro), - #(Rotation.ro2eu,Rotation.eu2ro), - (Rotation.ro2ax,Rotation.ax2ro), - (Rotation.ro2ho,Rotation.ho2ro), - (Rotation.ro2cu,Rotation.cu2ro)]) + @pytest.mark.parametrize('forward,backward',[(Rotation._ro2qu,Rotation._qu2ro), + #(Rotation._ro2om,Rotation._om2ro), + #(Rotation._ro2eu,Rotation._eu2ro), + (Rotation._ro2ax,Rotation._ax2ro), + (Rotation._ro2ho,Rotation._ho2ro), + (Rotation._ro2cu,Rotation._cu2ro)]) def test_Rodrigues_internal(self,default,forward,backward): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: @@ -172,12 +172,12 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) - @pytest.mark.parametrize('forward,backward',[(Rotation.ho2qu,Rotation.qu2ho), - (Rotation.ho2om,Rotation.om2ho), - #(Rotation.ho2eu,Rotation.eu2ho), - (Rotation.ho2ax,Rotation.ax2ho), - (Rotation.ho2ro,Rotation.ro2ho), - (Rotation.ho2cu,Rotation.cu2ho)]) + @pytest.mark.parametrize('forward,backward',[(Rotation._ho2qu,Rotation._qu2ho), + (Rotation._ho2om,Rotation._om2ho), + #(Rotation._ho2eu,Rotation._eu2ho), + (Rotation._ho2ax,Rotation._ax2ho), + (Rotation._ho2ro,Rotation._ro2ho), + (Rotation._ho2cu,Rotation._cu2ho)]) def test_homochoric_internal(self,default,forward,backward): for rot in default: m = rot.as_homochoric() @@ -289,11 +289,11 @@ class TestRotation: with pytest.raises(ValueError): function(invalid) - @pytest.mark.parametrize('vectorized, single',[(Rotation.qu2om,rotation_conversion.qu2om), - (Rotation.qu2eu,rotation_conversion.qu2eu), - (Rotation.qu2ax,rotation_conversion.qu2ax), - (Rotation.qu2ro,rotation_conversion.qu2ro), - (Rotation.qu2ho,rotation_conversion.qu2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,rotation_conversion.qu2om), + (Rotation._qu2eu,rotation_conversion.qu2eu), + (Rotation._qu2ax,rotation_conversion.qu2ax), + (Rotation._qu2ro,rotation_conversion.qu2ro), + (Rotation._qu2ho,rotation_conversion.qu2ho)]) def test_quaternion_vectorization(self,default,vectorized,single): qu = np.array([rot.as_quaternion() for rot in default]) vectorized(qu.reshape(qu.shape[0]//2,-1,4)) @@ -303,9 +303,9 @@ class TestRotation: assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.om2qu,rotation_conversion.om2qu), - (Rotation.om2eu,rotation_conversion.om2eu), - (Rotation.om2ax,rotation_conversion.om2ax)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,rotation_conversion.om2qu), + (Rotation._om2eu,rotation_conversion.om2eu), + (Rotation._om2ax,rotation_conversion.om2ax)]) def test_matrix_vectorization(self,default,vectorized,single): om = np.array([rot.as_matrix() for rot in default]) vectorized(om.reshape(om.shape[0]//2,-1,3,3)) @@ -314,10 +314,10 @@ class TestRotation: print(o,c) assert np.allclose(single(o),c) and np.allclose(single(o),vectorized(o)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.eu2qu,rotation_conversion.eu2qu), - (Rotation.eu2om,rotation_conversion.eu2om), - (Rotation.eu2ax,rotation_conversion.eu2ax), - (Rotation.eu2ro,rotation_conversion.eu2ro)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._eu2qu,rotation_conversion.eu2qu), + (Rotation._eu2om,rotation_conversion.eu2om), + (Rotation._eu2ax,rotation_conversion.eu2ax), + (Rotation._eu2ro,rotation_conversion.eu2ro)]) def test_Euler_vectorization(self,default,vectorized,single): eu = np.array([rot.as_Eulers() for rot in default]) vectorized(eu.reshape(eu.shape[0]//2,-1,3)) @@ -326,10 +326,10 @@ class TestRotation: print(e,c) assert np.allclose(single(e),c) and np.allclose(single(e),vectorized(e)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.ax2qu,rotation_conversion.ax2qu), - (Rotation.ax2om,rotation_conversion.ax2om), - (Rotation.ax2ro,rotation_conversion.ax2ro), - (Rotation.ax2ho,rotation_conversion.ax2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._ax2qu,rotation_conversion.ax2qu), + (Rotation._ax2om,rotation_conversion.ax2om), + (Rotation._ax2ro,rotation_conversion.ax2ro), + (Rotation._ax2ho,rotation_conversion.ax2ho)]) def test_axisAngle_vectorization(self,default,vectorized,single): ax = np.array([rot.as_axis_angle() for rot in default]) vectorized(ax.reshape(ax.shape[0]//2,-1,4)) @@ -339,8 +339,8 @@ class TestRotation: assert np.allclose(single(a),c) and np.allclose(single(a),vectorized(a)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.ro2ax,rotation_conversion.ro2ax), - (Rotation.ro2ho,rotation_conversion.ro2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,rotation_conversion.ro2ax), + (Rotation._ro2ho,rotation_conversion.ro2ho)]) def test_Rodrigues_vectorization(self,default,vectorized,single): ro = np.array([rot.as_Rodrigues() for rot in default]) vectorized(ro.reshape(ro.shape[0]//2,-1,4)) @@ -349,8 +349,8 @@ class TestRotation: print(r,c) assert np.allclose(single(r),c) and np.allclose(single(r),vectorized(r)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.ho2ax,rotation_conversion.ho2ax), - (Rotation.ho2cu,rotation_conversion.ho2cu)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,rotation_conversion.ho2ax), + (Rotation._ho2cu,rotation_conversion.ho2cu)]) def test_homochoric_vectorization(self,default,vectorized,single): ho = np.array([rot.as_homochoric() for rot in default]) vectorized(ho.reshape(ho.shape[0]//2,-1,3)) @@ -359,7 +359,7 @@ class TestRotation: print(h,c) assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) - @pytest.mark.parametrize('vectorized, single',[(Rotation.cu2ho,rotation_conversion.cu2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,rotation_conversion.cu2ho)]) def test_cubochoric_vectorization(self,default,vectorized,single): cu = np.array([rot.as_cubochoric() for rot in default]) vectorized(cu.reshape(cu.shape[0]//2,-1,3)) From b59d773689170a14571321322ff364fbfd51baec Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 15:26:49 +0200 Subject: [PATCH 28/32] store all test data together --- python/tests/rotation_conversion.py | 419 -------------------- python/tests/test_Rotation.py | 590 ++++++++++++++++++++++++---- 2 files changed, 506 insertions(+), 503 deletions(-) delete mode 100644 python/tests/rotation_conversion.py diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py deleted file mode 100644 index 58ac9329f..000000000 --- a/python/tests/rotation_conversion.py +++ /dev/null @@ -1,419 +0,0 @@ -#################################################################################################### -# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations -#################################################################################################### -# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH -# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University -# All rights reserved. -# -# Redistribution and use in source and binary forms, with or without modification, are -# permitted provided that the following conditions are met: -# -# - Redistributions of source code must retain the above copyright notice, this list -# of conditions and the following disclaimer. -# - Redistributions in binary form must reproduce the above copyright notice, this -# list of conditions and the following disclaimer in the documentation and/or -# other materials provided with the distribution. -# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -# of its contributors may be used to endorse or promote products derived from -# this software without specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE -# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -#################################################################################################### - -import numpy as np - -from damask import _rotation - -_P = _rotation._P - -# parameters for conversion from/to cubochoric -_sc = _rotation._sc -_beta = _rotation._beta -_R1 = _rotation._R1 - -def iszero(a): - return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) - -#---------- Quaternion ---------- -def qu2om(qu): - """Quaternion to rotation matrix.""" - qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) - om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) - - om[0,1] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) - om[1,0] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) - om[1,2] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) - om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) - om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) - om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) - return om if _P < 0.0 else np.swapaxes(om,-1,-2) - -def qu2eu(qu): - """Quaternion to Bunge-Euler angles.""" - q03 = qu[0]**2+qu[3]**2 - q12 = qu[1]**2+qu[2]**2 - chi = np.sqrt(q03*q12) - if np.abs(q12) < 1.e-8: - eu = np.array([np.arctan2(-_P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) - elif np.abs(q03) < 1.e-8: - eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) - else: - eu = np.array([np.arctan2((-_P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), - np.arctan2( 2.0*chi, q03-q12 ), - np.arctan2(( _P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) - # reduce Euler angles to definition range - eu[np.abs(eu)<1.e-6] = 0.0 - eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) - return eu - -def qu2ax(qu): - """ - Quaternion to axis angle pair. - - Modified version of the original formulation, should be numerically more stable - """ - if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - elif qu[0] > 1.e-8: - s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) - omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) - ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ]) - else: - ax = ax = np.array([ qu[1], qu[2], qu[3], np.pi]) - return ax - -def qu2ro(qu): - """Quaternion to Rodrigues-Frank vector.""" - if iszero(qu[0]): - ro = np.array([qu[1], qu[2], qu[3], np.inf]) - else: - s = np.linalg.norm(qu[1:4]) - ro = np.array([0.0,0.0,_P,0.0] if iszero(s) else \ - [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))]) - return ro - -def qu2ho(qu): - """Quaternion to homochoric vector.""" - omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) - if np.abs(omega) < 1.0e-12: - ho = np.zeros(3) - else: - ho = np.array([qu[1], qu[2], qu[3]]) - f = 0.75 * ( omega - np.sin(omega) ) - ho = ho/np.linalg.norm(ho) * f**(1./3.) - return ho - - -#---------- Rotation matrix ---------- -def om2qu(a): - trace = a[0,0] + a[1,1] + a[2,2] - if trace > 0: - s = 0.5 / np.sqrt(trace+ 1.0) - qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s]) - else: - if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ): - s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2]) - qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s]) - elif (a[1,1] > a[2,2]): - s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2]) - qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s]) - else: - s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) - qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) - if qu[0]<0: qu*=-1 - return qu*np.array([1.,_P,_P,_P]) - -def om2eu(om): - """Rotation matrix to Bunge-Euler angles.""" - if not np.isclose(np.abs(om[2,2]),1.0,1.e-9): - zeta = 1.0/np.sqrt(1.0-om[2,2]**2) - eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), - np.arccos(om[2,2]), - np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) - else: - eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation - eu[np.abs(eu)<1.e-8] = 0.0 - eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) - return eu - -def om2ax(om): - """Rotation matrix to axis angle pair.""" - #return qu2ax(om2qu(om)) # HOTFIX - ax=np.empty(4) - - # first get the rotation angle - t = 0.5*(om.trace() -1.0) - ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - if np.abs(ax[3])<1.e-8: - ax = np.array([ 0.0, 0.0, 1.0, 0.0]) - else: - w,vr = np.linalg.eig(om) - # next, find the eigenvalue (1,0j) - i = np.where(np.isclose(w,1.0+0.0j))[0][0] - ax[0:3] = np.real(vr[0:3,i]) - diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - ax[0:3] = np.where(np.abs(diagDelta)<1e-12, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) - return ax - -#---------- Bunge-Euler angles ---------- -def eu2qu(eu): - """Bunge-Euler angles to quaternion.""" - ee = 0.5*eu - cPhi = np.cos(ee[1]) - sPhi = np.sin(ee[1]) - qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), - -_P*sPhi*np.cos(ee[0]-ee[2]), - -_P*sPhi*np.sin(ee[0]-ee[2]), - -_P*cPhi*np.sin(ee[0]+ee[2]) ]) - if qu[0] < 0.0: qu*=-1 - return qu - -def eu2om(eu): - """Bunge-Euler angles to rotation matrix.""" - c = np.cos(eu) - s = np.sin(eu) - - om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], - [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], - [+s[0]*s[1], -c[0]*s[1], +c[1] ]]) - om[np.abs(om)<1.e-12] = 0.0 - return om - -def eu2ax(eu): - """Bunge-Euler angles to axis angle pair.""" - t = np.tan(eu[1]*0.5) - sigma = 0.5*(eu[0]+eu[2]) - delta = 0.5*(eu[0]-eu[2]) - tau = np.linalg.norm([t,np.sin(sigma)]) - alpha = np.pi if iszero(np.cos(sigma)) else \ - 2.0*np.arctan(tau/np.cos(sigma)) - - if np.abs(alpha)<1.e-6: - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - else: - ax = -_P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front - ax = np.append(ax,alpha) - if alpha < 0.0: ax *= -1.0 # ensure alpha is positive - return ax - -def eu2ro(eu): - """Bunge-Euler angles to Rodrigues-Frank vector.""" - ro = eu2ax(eu) # convert to axis angle pair representation - if ro[3] >= np.pi: # Differs from original implementation. check convention 5 - ro[3] = np.inf - elif iszero(ro[3]): - ro = np.array([ 0.0, 0.0, _P, 0.0 ]) - else: - ro[3] = np.tan(ro[3]*0.5) - return ro - -#---------- Axis angle pair ---------- -def ax2qu(ax): - """Axis angle pair to quaternion.""" - if np.abs(ax[3])<1.e-6: - qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) - else: - c = np.cos(ax[3]*0.5) - s = np.sin(ax[3]*0.5) - qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - return qu - -def ax2om(ax): - """Axis angle pair to rotation matrix.""" - c = np.cos(ax[3]) - s = np.sin(ax[3]) - omc = 1.0-c - om=np.diag(ax[0:3]**2*omc + c) - - for idx in [[0,1,2],[1,2,0],[2,0,1]]: - q = omc*ax[idx[0]] * ax[idx[1]] - om[idx[0],idx[1]] = q + s*ax[idx[2]] - om[idx[1],idx[0]] = q - s*ax[idx[2]] - return om if _P < 0.0 else om.T - -def ax2ro(ax): - """Axis angle pair to Rodrigues-Frank vector.""" - if np.abs(ax[3])<1.e-6: - ro = [ 0.0, 0.0, _P, 0.0 ] - else: - ro = [ax[0], ax[1], ax[2]] - # 180 degree case - ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ - [np.tan(ax[3]*0.5)] - ro = np.array(ro) - return ro - -def ax2ho(ax): - """Axis angle pair to homochoric vector.""" - f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) - ho = ax[0:3] * f - return ho - - -#---------- Rodrigues-Frank vector ---------- -def ro2ax(ro): - """Rodrigues-Frank vector to axis angle pair.""" - if np.abs(ro[3]) < 1.e-8: - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - elif not np.isfinite(ro[3]): - ax = np.array([ ro[0], ro[1], ro[2], np.pi ]) - else: - angle = 2.0*np.arctan(ro[3]) - ta = np.linalg.norm(ro[0:3]) - ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ]) - return ax - -def ro2ho(ro): - """Rodrigues-Frank vector to homochoric vector.""" - if np.sum(ro[0:3]**2.0) < 1.e-8: - ho = np.zeros(3) - else: - f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi - ho = ro[0:3] * (0.75*f)**(1.0/3.0) - return ho - -#---------- Homochoric vector---------- -def ho2ax(ho): - """Homochoric vector to axis angle pair.""" - tfit = np.array([+1.0000000000018852, -0.5000000002194847, - -0.024999992127593126, -0.003928701544781374, - -0.0008152701535450438, -0.0002009500426119712, - -0.00002397986776071756, -0.00008202868926605841, - +0.00012448715042090092, -0.0001749114214822577, - +0.0001703481934140054, -0.00012062065004116828, - +0.000059719705868660826, -0.00001980756723965647, - +0.000003953714684212874, -0.00000036555001439719544]) - # normalize h and store the magnitude - hmag_squared = np.sum(ho**2.) - if iszero(hmag_squared): - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - else: - hm = hmag_squared - - # convert the magnitude to the rotation angle - s = tfit[0] + tfit[1] * hmag_squared - for i in range(2,16): - hm *= hmag_squared - s += tfit[i] * hm - ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) - return ax - -def ho2cu(ho): - """ - Homochoric vector to cubochoric vector. - - References - ---------- - D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 - https://doi.org/10.1088/0965-0393/22/7/075013 - - """ - rs = np.linalg.norm(ho) - - if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16): - cu = np.zeros(3) - else: - xyz3 = ho[_get_pyramid_order(ho,'forward')] - - # inverse M_3 - xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) - - # inverse M_2 - qxy = np.sum(xyz2**2) - - if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16): - Tinv = np.zeros(2) - else: - q2 = qxy + np.max(np.abs(xyz2))**2 - sq2 = np.sqrt(q2) - q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) - tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) - Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ - np.array([np.arccos(tt)/np.pi*12.0,1.0]) - Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) - - # inverse M_1 - cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc - cu = cu[_get_pyramid_order(ho,'backward')] - return cu - -#---------- Cubochoric ---------- -def cu2ho(cu): - """ - Cubochoric vector to homochoric vector. - - References - ---------- - D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 - https://doi.org/10.1088/0965-0393/22/7/075013 - - """ - # transform to the sphere grid via the curved square, and intercept the zero point - if np.allclose(cu,0.0,rtol=0.0,atol=1.0e-16): - ho = np.zeros(3) - else: - # get pyramide and scale by grid parameter ratio - XYZ = cu[_get_pyramid_order(cu,'forward')] * _sc - - # intercept all the points along the z-axis - if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): - ho = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) - else: - order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] - q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] - c = np.cos(q) - s = np.sin(q) - q = _R1*2.0**0.25/_beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) - T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q - - # transform to sphere grid (inverse Lambert) - # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero - c = np.sum(T**2) - s = c * np.pi/24.0 /XYZ[2]**2 - c = c * np.sqrt(np.pi/24.0)/XYZ[2] - - q = np.sqrt( 1.0 - s ) - ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) - - ho = ho[_get_pyramid_order(cu,'backward')] - return ho - -def _get_pyramid_order(xyz,direction=None): - """ - Get order of the coordinates. - - Depending on the pyramid in which the point is located, the order need to be adjusted. - - Parameters - ---------- - xyz : numpy.ndarray - coordinates of a point on a uniform refinable grid on a ball or - in a uniform refinable cubical grid. - - References - ---------- - D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 - https://doi.org/10.1088/0965-0393/22/7/075013 - - """ - order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]), - 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} - if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \ - np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]: - p = 0 - elif np.maximum(abs(xyz[1]),abs(xyz[2])) <= xyz[0] or \ - np.maximum(abs(xyz[1]),abs(xyz[2])) <=-xyz[0]: - p = 1 - elif np.maximum(abs(xyz[2]),abs(xyz[0])) <= xyz[1] or \ - np.maximum(abs(xyz[2]),abs(xyz[0])) <=-xyz[1]: - p = 2 - return order[direction][p] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 4e48fbb7c..92f39b4fb 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -4,7 +4,9 @@ import pytest import numpy as np from damask import Rotation -import rotation_conversion +from damask import _rotation + + n = 1100 atol=1.e-4 @@ -13,68 +15,69 @@ scatter=1.e-2 @pytest.fixture def default(): """A set of n random rotations.""" - specials = np.array( - [np.array([ 1.0, 0.0, 0.0, 0.0]), - #----------------------------------------------- - np.array([0.0, 1.0, 0.0, 0.0]), - np.array([0.0, 0.0, 1.0, 0.0]), - np.array([0.0, 0.0, 0.0, 1.0]), - np.array([0.0,-1.0, 0.0, 0.0]), - np.array([0.0, 0.0,-1.0, 0.0]), - np.array([0.0, 0.0, 0.0,-1.0]), - #----------------------------------------------- - np.array([1.0, 1.0, 0.0, 0.0])/np.sqrt(2.), - np.array([1.0, 0.0, 1.0, 0.0])/np.sqrt(2.), - np.array([1.0, 0.0, 0.0, 1.0])/np.sqrt(2.), - np.array([0.0, 1.0, 1.0, 0.0])/np.sqrt(2.), - np.array([0.0, 1.0, 0.0, 1.0])/np.sqrt(2.), - np.array([0.0, 0.0, 1.0, 1.0])/np.sqrt(2.), - #----------------------------------------------- - np.array([1.0,-1.0, 0.0, 0.0])/np.sqrt(2.), - np.array([1.0, 0.0,-1.0, 0.0])/np.sqrt(2.), - np.array([1.0, 0.0, 0.0,-1.0])/np.sqrt(2.), - np.array([0.0, 1.0,-1.0, 0.0])/np.sqrt(2.), - np.array([0.0, 1.0, 0.0,-1.0])/np.sqrt(2.), - np.array([0.0, 0.0, 1.0,-1.0])/np.sqrt(2.), - #----------------------------------------------- - np.array([0.0, 1.0,-1.0, 0.0])/np.sqrt(2.), - np.array([0.0, 1.0, 0.0,-1.0])/np.sqrt(2.), - np.array([0.0, 0.0, 1.0,-1.0])/np.sqrt(2.), - #----------------------------------------------- - np.array([0.0,-1.0,-1.0, 0.0])/np.sqrt(2.), - np.array([0.0,-1.0, 0.0,-1.0])/np.sqrt(2.), - np.array([0.0, 0.0,-1.0,-1.0])/np.sqrt(2.), - #----------------------------------------------- - np.array([1.0, 1.0, 1.0, 0.0])/np.sqrt(3.), - np.array([1.0, 1.0, 0.0, 1.0])/np.sqrt(3.), - np.array([1.0, 0.0, 1.0, 1.0])/np.sqrt(3.), - np.array([1.0,-1.0, 1.0, 0.0])/np.sqrt(3.), - np.array([1.0,-1.0, 0.0, 1.0])/np.sqrt(3.), - np.array([1.0, 0.0,-1.0, 1.0])/np.sqrt(3.), - np.array([1.0, 1.0,-1.0, 0.0])/np.sqrt(3.), - np.array([1.0, 1.0, 0.0,-1.0])/np.sqrt(3.), - np.array([1.0, 0.0, 1.0,-1.0])/np.sqrt(3.), - np.array([1.0,-1.0,-1.0, 0.0])/np.sqrt(3.), - np.array([1.0,-1.0, 0.0,-1.0])/np.sqrt(3.), - np.array([1.0, 0.0,-1.0,-1.0])/np.sqrt(3.), - #----------------------------------------------- - np.array([0.0, 1.0, 1.0, 1.0])/np.sqrt(3.), - np.array([0.0, 1.0,-1.0, 1.0])/np.sqrt(3.), - np.array([0.0, 1.0, 1.0,-1.0])/np.sqrt(3.), - np.array([0.0,-1.0, 1.0, 1.0])/np.sqrt(3.), - np.array([0.0,-1.0,-1.0, 1.0])/np.sqrt(3.), - np.array([0.0,-1.0, 1.0,-1.0])/np.sqrt(3.), - np.array([0.0,-1.0,-1.0,-1.0])/np.sqrt(3.), - #----------------------------------------------- - np.array([1.0, 1.0, 1.0, 1.0])/2., - np.array([1.0,-1.0, 1.0, 1.0])/2., - np.array([1.0, 1.0,-1.0, 1.0])/2., - np.array([1.0, 1.0, 1.0,-1.0])/2., - np.array([1.0,-1.0,-1.0, 1.0])/2., - np.array([1.0,-1.0, 1.0,-1.0])/2., - np.array([1.0, 1.0,-1.0,-1.0])/2., - np.array([1.0,-1.0,-1.0,-1.0])/2., - ]) + specials = np.array([ + [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,-1.0, 0.0, 0.0], + [0.0, 0.0,-1.0, 0.0], + [0.0, 0.0, 0.0,-1.0], + #---------------------- + [1.0, 1.0, 0.0, 0.0], + [1.0, 0.0, 1.0, 0.0], + [1.0, 0.0, 0.0, 1.0], + [0.0, 1.0, 1.0, 0.0], + [0.0, 1.0, 0.0, 1.0], + [0.0, 0.0, 1.0, 1.0], + #---------------------- + [1.0,-1.0, 0.0, 0.0], + [1.0, 0.0,-1.0, 0.0], + [1.0, 0.0, 0.0,-1.0], + [0.0, 1.0,-1.0, 0.0], + [0.0, 1.0, 0.0,-1.0], + [0.0, 0.0, 1.0,-1.0], + #---------------------- + [0.0, 1.0,-1.0, 0.0], + [0.0, 1.0, 0.0,-1.0], + [0.0, 0.0, 1.0,-1.0], + #---------------------- + [0.0,-1.0,-1.0, 0.0], + [0.0,-1.0, 0.0,-1.0], + [0.0, 0.0,-1.0,-1.0], + #---------------------- + [1.0, 1.0, 1.0, 0.0], + [1.0, 1.0, 0.0, 1.0], + [1.0, 0.0, 1.0, 1.0], + [1.0,-1.0, 1.0, 0.0], + [1.0,-1.0, 0.0, 1.0], + [1.0, 0.0,-1.0, 1.0], + [1.0, 1.0,-1.0, 0.0], + [1.0, 1.0, 0.0,-1.0], + [1.0, 0.0, 1.0,-1.0], + [1.0,-1.0,-1.0, 0.0], + [1.0,-1.0, 0.0,-1.0], + [1.0, 0.0,-1.0,-1.0], + #---------------------- + [0.0, 1.0, 1.0, 1.0], + [0.0, 1.0,-1.0, 1.0], + [0.0, 1.0, 1.0,-1.0], + [0.0,-1.0, 1.0, 1.0], + [0.0,-1.0,-1.0, 1.0], + [0.0,-1.0, 1.0,-1.0], + [0.0,-1.0,-1.0,-1.0], + #---------------------- + [1.0, 1.0, 1.0, 1.0], + [1.0,-1.0, 1.0, 1.0], + [1.0, 1.0,-1.0, 1.0], + [1.0, 1.0, 1.0,-1.0], + [1.0,-1.0,-1.0, 1.0], + [1.0,-1.0, 1.0,-1.0], + [1.0, 1.0,-1.0,-1.0], + [1.0,-1.0,-1.0,-1.0], + ]) + specials /= np.linalg.norm(specials,axis=1).reshape(-1,1) specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape) specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1) specials_scatter[specials_scatter[:,0]<0]*=-1 @@ -89,6 +92,425 @@ def reference_dir(reference_dir_base): return os.path.join(reference_dir_base,'Rotation') +#################################################################################################### +# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations +#################################################################################################### +# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH +# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without modification, are +# permitted provided that the following conditions are met: +# +# - Redistributions of source code must retain the above copyright notice, this list +# of conditions and the following disclaimer. +# - Redistributions in binary form must reproduce the above copyright notice, this +# list of conditions and the following disclaimer in the documentation and/or +# other materials provided with the distribution. +# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +# of its contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +#################################################################################################### +_P = _rotation._P + +# parameters for conversion from/to cubochoric +_sc = _rotation._sc +_beta = _rotation._beta +_R1 = _rotation._R1 + +def iszero(a): + return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) + +#---------- Quaternion ---------- +def qu2om(qu): + """Quaternion to rotation matrix.""" + qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) + om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) + + om[0,1] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) + om[1,0] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) + om[1,2] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) + om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) + om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) + om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) + return om if _P < 0.0 else np.swapaxes(om,-1,-2) + +def qu2eu(qu): + """Quaternion to Bunge-Euler angles.""" + q03 = qu[0]**2+qu[3]**2 + q12 = qu[1]**2+qu[2]**2 + chi = np.sqrt(q03*q12) + if np.abs(q12) < 1.e-8: + eu = np.array([np.arctan2(-_P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) + elif np.abs(q03) < 1.e-8: + eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) + else: + eu = np.array([np.arctan2((-_P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), + np.arctan2( 2.0*chi, q03-q12 ), + np.arctan2(( _P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) + # reduce Euler angles to definition range + eu[np.abs(eu)<1.e-6] = 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + +def qu2ax(qu): + """ + Quaternion to axis angle pair. + + Modified version of the original formulation, should be numerically more stable + """ + if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + elif qu[0] > 1.e-8: + s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) + ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ]) + else: + ax = ax = np.array([ qu[1], qu[2], qu[3], np.pi]) + return ax + +def qu2ro(qu): + """Quaternion to Rodrigues-Frank vector.""" + if iszero(qu[0]): + ro = np.array([qu[1], qu[2], qu[3], np.inf]) + else: + s = np.linalg.norm(qu[1:4]) + ro = np.array([0.0,0.0,_P,0.0] if iszero(s) else \ + [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))]) + return ro + +def qu2ho(qu): + """Quaternion to homochoric vector.""" + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) + if np.abs(omega) < 1.0e-12: + ho = np.zeros(3) + else: + ho = np.array([qu[1], qu[2], qu[3]]) + f = 0.75 * ( omega - np.sin(omega) ) + ho = ho/np.linalg.norm(ho) * f**(1./3.) + return ho + + +#---------- Rotation matrix ---------- +def om2qu(a): + trace = a[0,0] + a[1,1] + a[2,2] + if trace > 0: + s = 0.5 / np.sqrt(trace+ 1.0) + qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s]) + else: + if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ): + s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2]) + qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s]) + elif (a[1,1] > a[2,2]): + s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2]) + qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s]) + else: + s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] ) + qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s]) + if qu[0]<0: qu*=-1 + return qu*np.array([1.,_P,_P,_P]) + +def om2eu(om): + """Rotation matrix to Bunge-Euler angles.""" + if not np.isclose(np.abs(om[2,2]),1.0,1.e-9): + zeta = 1.0/np.sqrt(1.0-om[2,2]**2) + eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), + np.arccos(om[2,2]), + np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) + else: + eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation + eu[np.abs(eu)<1.e-8] = 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + +def om2ax(om): + """Rotation matrix to axis angle pair.""" + #return qu2ax(om2qu(om)) # HOTFIX + ax=np.empty(4) + + # first get the rotation angle + t = 0.5*(om.trace() -1.0) + ax[3] = np.arccos(np.clip(t,-1.0,1.0)) + if np.abs(ax[3])<1.e-8: + ax = np.array([ 0.0, 0.0, 1.0, 0.0]) + else: + w,vr = np.linalg.eig(om) + # next, find the eigenvalue (1,0j) + i = np.where(np.isclose(w,1.0+0.0j))[0][0] + ax[0:3] = np.real(vr[0:3,i]) + diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) + ax[0:3] = np.where(np.abs(diagDelta)<1e-12, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) + return ax + +#---------- Bunge-Euler angles ---------- +def eu2qu(eu): + """Bunge-Euler angles to quaternion.""" + ee = 0.5*eu + cPhi = np.cos(ee[1]) + sPhi = np.sin(ee[1]) + qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), + -_P*sPhi*np.cos(ee[0]-ee[2]), + -_P*sPhi*np.sin(ee[0]-ee[2]), + -_P*cPhi*np.sin(ee[0]+ee[2]) ]) + if qu[0] < 0.0: qu*=-1 + return qu + +def eu2om(eu): + """Bunge-Euler angles to rotation matrix.""" + c = np.cos(eu) + s = np.sin(eu) + + om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], + [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], + [+s[0]*s[1], -c[0]*s[1], +c[1] ]]) + om[np.abs(om)<1.e-12] = 0.0 + return om + +def eu2ax(eu): + """Bunge-Euler angles to axis angle pair.""" + t = np.tan(eu[1]*0.5) + sigma = 0.5*(eu[0]+eu[2]) + delta = 0.5*(eu[0]-eu[2]) + tau = np.linalg.norm([t,np.sin(sigma)]) + alpha = np.pi if iszero(np.cos(sigma)) else \ + 2.0*np.arctan(tau/np.cos(sigma)) + + if np.abs(alpha)<1.e-6: + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + else: + ax = -_P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front + ax = np.append(ax,alpha) + if alpha < 0.0: ax *= -1.0 # ensure alpha is positive + return ax + +def eu2ro(eu): + """Bunge-Euler angles to Rodrigues-Frank vector.""" + ro = eu2ax(eu) # convert to axis angle pair representation + if ro[3] >= np.pi: # Differs from original implementation. check convention 5 + ro[3] = np.inf + elif iszero(ro[3]): + ro = np.array([ 0.0, 0.0, _P, 0.0 ]) + else: + ro[3] = np.tan(ro[3]*0.5) + return ro + +#---------- Axis angle pair ---------- +def ax2qu(ax): + """Axis angle pair to quaternion.""" + if np.abs(ax[3])<1.e-6: + qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) + else: + c = np.cos(ax[3]*0.5) + s = np.sin(ax[3]*0.5) + qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) + return qu + +def ax2om(ax): + """Axis angle pair to rotation matrix.""" + c = np.cos(ax[3]) + s = np.sin(ax[3]) + omc = 1.0-c + om=np.diag(ax[0:3]**2*omc + c) + + for idx in [[0,1,2],[1,2,0],[2,0,1]]: + q = omc*ax[idx[0]] * ax[idx[1]] + om[idx[0],idx[1]] = q + s*ax[idx[2]] + om[idx[1],idx[0]] = q - s*ax[idx[2]] + return om if _P < 0.0 else om.T + +def ax2ro(ax): + """Axis angle pair to Rodrigues-Frank vector.""" + if np.abs(ax[3])<1.e-6: + ro = [ 0.0, 0.0, _P, 0.0 ] + else: + ro = [ax[0], ax[1], ax[2]] + # 180 degree case + ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ + [np.tan(ax[3]*0.5)] + ro = np.array(ro) + return ro + +def ax2ho(ax): + """Axis angle pair to homochoric vector.""" + f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) + ho = ax[0:3] * f + return ho + + +#---------- Rodrigues-Frank vector ---------- +def ro2ax(ro): + """Rodrigues-Frank vector to axis angle pair.""" + if np.abs(ro[3]) < 1.e-8: + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + elif not np.isfinite(ro[3]): + ax = np.array([ ro[0], ro[1], ro[2], np.pi ]) + else: + angle = 2.0*np.arctan(ro[3]) + ta = np.linalg.norm(ro[0:3]) + ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ]) + return ax + +def ro2ho(ro): + """Rodrigues-Frank vector to homochoric vector.""" + if np.sum(ro[0:3]**2.0) < 1.e-8: + ho = np.zeros(3) + else: + f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi + ho = ro[0:3] * (0.75*f)**(1.0/3.0) + return ho + +#---------- Homochoric vector---------- +def ho2ax(ho): + """Homochoric vector to axis angle pair.""" + tfit = np.array([+1.0000000000018852, -0.5000000002194847, + -0.024999992127593126, -0.003928701544781374, + -0.0008152701535450438, -0.0002009500426119712, + -0.00002397986776071756, -0.00008202868926605841, + +0.00012448715042090092, -0.0001749114214822577, + +0.0001703481934140054, -0.00012062065004116828, + +0.000059719705868660826, -0.00001980756723965647, + +0.000003953714684212874, -0.00000036555001439719544]) + # normalize h and store the magnitude + hmag_squared = np.sum(ho**2.) + if iszero(hmag_squared): + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + else: + hm = hmag_squared + + # convert the magnitude to the rotation angle + s = tfit[0] + tfit[1] * hmag_squared + for i in range(2,16): + hm *= hmag_squared + s += tfit[i] * hm + ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) + return ax + +def ho2cu(ho): + """ + Homochoric vector to cubochoric vector. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + rs = np.linalg.norm(ho) + + if np.allclose(ho,0.0,rtol=0.0,atol=1.0e-16): + cu = np.zeros(3) + else: + xyz3 = ho[_get_pyramid_order(ho,'forward')] + + # inverse M_3 + xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) + + # inverse M_2 + qxy = np.sum(xyz2**2) + + if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16): + Tinv = np.zeros(2) + else: + q2 = qxy + np.max(np.abs(xyz2))**2 + sq2 = np.sqrt(q2) + q = (_beta/np.sqrt(2.0)/_R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) + tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) + Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ + np.array([np.arccos(tt)/np.pi*12.0,1.0]) + Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) + + # inverse M_1 + cu = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /_sc + cu = cu[_get_pyramid_order(ho,'backward')] + return cu + +#---------- Cubochoric ---------- +def cu2ho(cu): + """ + Cubochoric vector to homochoric vector. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + # transform to the sphere grid via the curved square, and intercept the zero point + if np.allclose(cu,0.0,rtol=0.0,atol=1.0e-16): + ho = np.zeros(3) + else: + # get pyramide and scale by grid parameter ratio + XYZ = cu[_get_pyramid_order(cu,'forward')] * _sc + + # intercept all the points along the z-axis + if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16): + ho = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) + else: + order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] + q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] + c = np.cos(q) + s = np.sin(q) + q = _R1*2.0**0.25/_beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) + T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q + + # transform to sphere grid (inverse Lambert) + # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero + c = np.sum(T**2) + s = c * np.pi/24.0 /XYZ[2]**2 + c = c * np.sqrt(np.pi/24.0)/XYZ[2] + + q = np.sqrt( 1.0 - s ) + ho = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) + + ho = ho[_get_pyramid_order(cu,'backward')] + return ho + +def _get_pyramid_order(xyz,direction=None): + """ + Get order of the coordinates. + + Depending on the pyramid in which the point is located, the order need to be adjusted. + + Parameters + ---------- + xyz : numpy.ndarray + coordinates of a point on a uniform refinable grid on a ball or + in a uniform refinable cubical grid. + + References + ---------- + D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 + https://doi.org/10.1088/0965-0393/22/7/075013 + + """ + order = {'forward':np.array([[0,1,2],[1,2,0],[2,0,1]]), + 'backward':np.array([[0,1,2],[2,0,1],[1,2,0]])} + if np.maximum(abs(xyz[0]),abs(xyz[1])) <= xyz[2] or \ + np.maximum(abs(xyz[0]),abs(xyz[1])) <=-xyz[2]: + p = 0 + elif np.maximum(abs(xyz[1]),abs(xyz[2])) <= xyz[0] or \ + np.maximum(abs(xyz[1]),abs(xyz[2])) <=-xyz[0]: + p = 1 + elif np.maximum(abs(xyz[2]),abs(xyz[0])) <= xyz[1] or \ + np.maximum(abs(xyz[2]),abs(xyz[0])) <=-xyz[1]: + p = 2 + return order[direction][p] + +#################################################################################################### +#################################################################################################### + + class TestRotation: @pytest.mark.parametrize('forward,backward',[(Rotation._qu2om,Rotation._om2qu), @@ -289,11 +711,11 @@ class TestRotation: with pytest.raises(ValueError): function(invalid) - @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,rotation_conversion.qu2om), - (Rotation._qu2eu,rotation_conversion.qu2eu), - (Rotation._qu2ax,rotation_conversion.qu2ax), - (Rotation._qu2ro,rotation_conversion.qu2ro), - (Rotation._qu2ho,rotation_conversion.qu2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,qu2om), + (Rotation._qu2eu,qu2eu), + (Rotation._qu2ax,qu2ax), + (Rotation._qu2ro,qu2ro), + (Rotation._qu2ho,qu2ho)]) def test_quaternion_vectorization(self,default,vectorized,single): qu = np.array([rot.as_quaternion() for rot in default]) vectorized(qu.reshape(qu.shape[0]//2,-1,4)) @@ -303,9 +725,9 @@ class TestRotation: assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q)) - @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,rotation_conversion.om2qu), - (Rotation._om2eu,rotation_conversion.om2eu), - (Rotation._om2ax,rotation_conversion.om2ax)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu), + (Rotation._om2eu,om2eu), + (Rotation._om2ax,om2ax)]) def test_matrix_vectorization(self,default,vectorized,single): om = np.array([rot.as_matrix() for rot in default]) vectorized(om.reshape(om.shape[0]//2,-1,3,3)) @@ -314,10 +736,10 @@ class TestRotation: print(o,c) assert np.allclose(single(o),c) and np.allclose(single(o),vectorized(o)) - @pytest.mark.parametrize('vectorized, single',[(Rotation._eu2qu,rotation_conversion.eu2qu), - (Rotation._eu2om,rotation_conversion.eu2om), - (Rotation._eu2ax,rotation_conversion.eu2ax), - (Rotation._eu2ro,rotation_conversion.eu2ro)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._eu2qu,eu2qu), + (Rotation._eu2om,eu2om), + (Rotation._eu2ax,eu2ax), + (Rotation._eu2ro,eu2ro)]) def test_Euler_vectorization(self,default,vectorized,single): eu = np.array([rot.as_Eulers() for rot in default]) vectorized(eu.reshape(eu.shape[0]//2,-1,3)) @@ -326,10 +748,10 @@ class TestRotation: print(e,c) assert np.allclose(single(e),c) and np.allclose(single(e),vectorized(e)) - @pytest.mark.parametrize('vectorized, single',[(Rotation._ax2qu,rotation_conversion.ax2qu), - (Rotation._ax2om,rotation_conversion.ax2om), - (Rotation._ax2ro,rotation_conversion.ax2ro), - (Rotation._ax2ho,rotation_conversion.ax2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._ax2qu,ax2qu), + (Rotation._ax2om,ax2om), + (Rotation._ax2ro,ax2ro), + (Rotation._ax2ho,ax2ho)]) def test_axisAngle_vectorization(self,default,vectorized,single): ax = np.array([rot.as_axis_angle() for rot in default]) vectorized(ax.reshape(ax.shape[0]//2,-1,4)) @@ -339,8 +761,8 @@ class TestRotation: assert np.allclose(single(a),c) and np.allclose(single(a),vectorized(a)) - @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,rotation_conversion.ro2ax), - (Rotation._ro2ho,rotation_conversion.ro2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax), + (Rotation._ro2ho,ro2ho)]) def test_Rodrigues_vectorization(self,default,vectorized,single): ro = np.array([rot.as_Rodrigues() for rot in default]) vectorized(ro.reshape(ro.shape[0]//2,-1,4)) @@ -349,8 +771,8 @@ class TestRotation: print(r,c) assert np.allclose(single(r),c) and np.allclose(single(r),vectorized(r)) - @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,rotation_conversion.ho2ax), - (Rotation._ho2cu,rotation_conversion.ho2cu)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax), + (Rotation._ho2cu,ho2cu)]) def test_homochoric_vectorization(self,default,vectorized,single): ho = np.array([rot.as_homochoric() for rot in default]) vectorized(ho.reshape(ho.shape[0]//2,-1,3)) @@ -359,7 +781,7 @@ class TestRotation: print(h,c) assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) - @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,rotation_conversion.cu2ho)]) + @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)]) def test_cubochoric_vectorization(self,default,vectorized,single): cu = np.array([rot.as_cubochoric() for rot in default]) vectorized(cu.reshape(cu.shape[0]//2,-1,3)) From bb419d49dfc81bb70a0a5c42a95ef587a9378461 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 16:19:12 +0200 Subject: [PATCH 29/32] polishing --- python/damask/_rotation.py | 3 +- python/tests/test_Rotation.py | 198 +++++++++++++++++++--------------- 2 files changed, 114 insertions(+), 87 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index cf6eaed34..4081f2844 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -9,7 +9,6 @@ _sc = np.pi**(1./6.)/6.**(1./6.) _beta = np.pi**(5./6.)/6.**(1./6.)/2. _R1 = (3.*np.pi/4.)**(1./3.) - class Rotation: u""" Orientation stored with functionality for conversion to different representations. @@ -427,7 +426,7 @@ class Rotation: if P > 0: ho *= -1 # convert from P=1 to P=-1 - if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9): + if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9): raise ValueError('Homochoric coordinate outside of the sphere.') return Rotation(Rotation._ho2qu(ho)) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 92f39b4fb..f424fd9ab 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -520,6 +520,7 @@ class TestRotation: (Rotation._qu2ho,Rotation._ho2qu), (Rotation._qu2cu,Rotation._cu2qu)]) def test_quaternion_internal(self,default,forward,backward): + """Ensure invariance of conversion from quaternion and back.""" for rot in default: m = rot.as_quaternion() o = backward(forward(m)) @@ -536,6 +537,7 @@ class TestRotation: (Rotation._om2ho,Rotation._ho2om), (Rotation._om2cu,Rotation._cu2om)]) def test_matrix_internal(self,default,forward,backward): + """Ensure invariance of conversion from rotation matrix and back.""" for rot in default: m = rot.as_matrix() o = backward(forward(m)) @@ -550,6 +552,7 @@ class TestRotation: (Rotation._eu2ho,Rotation._ho2eu), (Rotation._eu2cu,Rotation._cu2eu)]) def test_Eulers_internal(self,default,forward,backward): + """Ensure invariance of conversion from Euler angles and back.""" for rot in default: m = rot.as_Eulers() o = backward(forward(m)) @@ -569,6 +572,7 @@ class TestRotation: (Rotation._ax2ho,Rotation._ho2ax), (Rotation._ax2cu,Rotation._cu2ax)]) def test_axis_angle_internal(self,default,forward,backward): + """Ensure invariance of conversion from axis angle angles pair and back.""" for rot in default: m = rot.as_axis_angle() o = backward(forward(m)) @@ -585,6 +589,7 @@ class TestRotation: (Rotation._ro2ho,Rotation._ho2ro), (Rotation._ro2cu,Rotation._cu2ro)]) def test_Rodrigues_internal(self,default,forward,backward): + """Ensure invariance of conversion from Rodrigues-Frank vector and back.""" cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: m = rot.as_Rodrigues() @@ -601,12 +606,114 @@ class TestRotation: (Rotation._ho2ro,Rotation._ro2ho), (Rotation._ho2cu,Rotation._cu2ho)]) def test_homochoric_internal(self,default,forward,backward): + """Ensure invariance of conversion from homochoric vector and back.""" for rot in default: m = rot.as_homochoric() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) print(m,o,rot.as_quaternion()) - assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 + assert ok and np.linalg.norm(o) < _R1 + 1.e-9 + + @pytest.mark.parametrize('forward,backward',[(Rotation._cu2qu,Rotation._qu2cu), + (Rotation._cu2om,Rotation._om2cu), + #(Rotation._cu2eu,Rotation._eu2cu), + (Rotation._cu2ax,Rotation._ax2cu), + (Rotation._cu2ro,Rotation._ro2cu), + (Rotation._cu2ho,Rotation._ho2cu)]) + def test_cubochoric_internal(self,default,forward,backward): + """Ensure invariance of conversion from cubochoric vector and back.""" + for rot in default: + m = rot.as_cubochoric() + o = backward(forward(m)) + ok = np.allclose(m,o,atol=atol) + print(m,o,rot.as_quaternion()) + assert ok and np.max(np.abs(o)) < np.pi**(2./3.) * 0.5 + 1.e-9 + + @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,qu2om), + (Rotation._qu2eu,qu2eu), + (Rotation._qu2ax,qu2ax), + (Rotation._qu2ro,qu2ro), + (Rotation._qu2ho,qu2ho)]) + def test_quaternion_vectorization(self,default,vectorized,single): + """Check vectorized implementation for quaternion against single point calculation.""" + qu = np.array([rot.as_quaternion() for rot in default]) + vectorized(qu.reshape(qu.shape[0]//2,-1,4)) + co = vectorized(qu) + for q,c in zip(qu,co): + print(q,c) + assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q)) + + + @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu), + (Rotation._om2eu,om2eu), + (Rotation._om2ax,om2ax)]) + def test_matrix_vectorization(self,default,vectorized,single): + """Check vectorized implementation for rotation matrix against single point calculation.""" + om = np.array([rot.as_matrix() for rot in default]) + vectorized(om.reshape(om.shape[0]//2,-1,3,3)) + co = vectorized(om) + for o,c in zip(om,co): + print(o,c) + assert np.allclose(single(o),c) and np.allclose(single(o),vectorized(o)) + + @pytest.mark.parametrize('vectorized, single',[(Rotation._eu2qu,eu2qu), + (Rotation._eu2om,eu2om), + (Rotation._eu2ax,eu2ax), + (Rotation._eu2ro,eu2ro)]) + def test_Eulers_vectorization(self,default,vectorized,single): + """Check vectorized implementation for Euler angles against single point calculation.""" + eu = np.array([rot.as_Eulers() for rot in default]) + vectorized(eu.reshape(eu.shape[0]//2,-1,3)) + co = vectorized(eu) + for e,c in zip(eu,co): + print(e,c) + assert np.allclose(single(e),c) and np.allclose(single(e),vectorized(e)) + + @pytest.mark.parametrize('vectorized, single',[(Rotation._ax2qu,ax2qu), + (Rotation._ax2om,ax2om), + (Rotation._ax2ro,ax2ro), + (Rotation._ax2ho,ax2ho)]) + def test_axis_angle_vectorization(self,default,vectorized,single): + """Check vectorized implementation for axis angle pair against single point calculation.""" + ax = np.array([rot.as_axis_angle() for rot in default]) + vectorized(ax.reshape(ax.shape[0]//2,-1,4)) + co = vectorized(ax) + for a,c in zip(ax,co): + print(a,c) + assert np.allclose(single(a),c) and np.allclose(single(a),vectorized(a)) + + + @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax), + (Rotation._ro2ho,ro2ho)]) + def test_Rodrigues_vectorization(self,default,vectorized,single): + """Check vectorized implementation for Rodrigues-Frank vector against single point calculation.""" + ro = np.array([rot.as_Rodrigues() for rot in default]) + vectorized(ro.reshape(ro.shape[0]//2,-1,4)) + co = vectorized(ro) + for r,c in zip(ro,co): + print(r,c) + assert np.allclose(single(r),c) and np.allclose(single(r),vectorized(r)) + + @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax), + (Rotation._ho2cu,ho2cu)]) + def test_homochoric_vectorization(self,default,vectorized,single): + """Check vectorized implementation for homochoric vector against single point calculation.""" + ho = np.array([rot.as_homochoric() for rot in default]) + vectorized(ho.reshape(ho.shape[0]//2,-1,3)) + co = vectorized(ho) + for h,c in zip(ho,co): + print(h,c) + assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) + + @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)]) + def test_cubochoric_vectorization(self,default,vectorized,single): + """Check vectorized implementation for cubochoric vector against single point calculation.""" + cu = np.array([rot.as_cubochoric() for rot in default]) + vectorized(cu.reshape(cu.shape[0]//2,-1,3)) + co = vectorized(cu) + for u,c in zip(cu,co): + print(u,c) + assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u)) @pytest.mark.parametrize('degrees',[True,False]) def test_Eulers(self,default,degrees): @@ -622,7 +729,7 @@ class TestRotation: @pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('normalise',[True,False]) @pytest.mark.parametrize('degrees',[True,False]) - def test_AxisAngle(self,default,degrees,normalise,P): + def test_axis_angle(self,default,degrees,normalise,P): c = np.array([P*-1,P*-1,P*-1,1.]) for rot in default: m = rot.as_Eulers() @@ -636,7 +743,7 @@ class TestRotation: print(m,o,rot.as_quaternion()) assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() - def test_Matrix(self,default): + def test_matrix(self,default): for rot in default: m = rot.as_axis_angle() o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle() @@ -658,7 +765,7 @@ class TestRotation: assert ok and np.isclose(np.linalg.det(o),1.0) @pytest.mark.parametrize('P',[1,-1]) - def test_Homochoric(self,default,P): + def test_homochoric(self,default,P): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in default: m = rot.as_Rodrigues() @@ -669,7 +776,7 @@ class TestRotation: assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) @pytest.mark.parametrize('P',[1,-1]) - def test_Cubochoric(self,default,P): + def test_cubochoric(self,default,P): for rot in default: m = rot.as_homochoric() o = Rotation.from_cubochoric(rot.as_cubochoric()*P*-1,P).as_homochoric() @@ -678,7 +785,7 @@ class TestRotation: assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 @pytest.mark.parametrize('P',[1,-1]) - def test_Quaternion(self,default,P): + def test_quaternion(self,default,P): c = np.array([1,P*-1,P*-1,P*-1]) for rot in default: m = rot.as_cubochoric() @@ -711,85 +818,6 @@ class TestRotation: with pytest.raises(ValueError): function(invalid) - @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,qu2om), - (Rotation._qu2eu,qu2eu), - (Rotation._qu2ax,qu2ax), - (Rotation._qu2ro,qu2ro), - (Rotation._qu2ho,qu2ho)]) - def test_quaternion_vectorization(self,default,vectorized,single): - qu = np.array([rot.as_quaternion() for rot in default]) - vectorized(qu.reshape(qu.shape[0]//2,-1,4)) - co = vectorized(qu) - for q,c in zip(qu,co): - print(q,c) - assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q)) - - - @pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu), - (Rotation._om2eu,om2eu), - (Rotation._om2ax,om2ax)]) - def test_matrix_vectorization(self,default,vectorized,single): - om = np.array([rot.as_matrix() for rot in default]) - vectorized(om.reshape(om.shape[0]//2,-1,3,3)) - co = vectorized(om) - for o,c in zip(om,co): - print(o,c) - assert np.allclose(single(o),c) and np.allclose(single(o),vectorized(o)) - - @pytest.mark.parametrize('vectorized, single',[(Rotation._eu2qu,eu2qu), - (Rotation._eu2om,eu2om), - (Rotation._eu2ax,eu2ax), - (Rotation._eu2ro,eu2ro)]) - def test_Euler_vectorization(self,default,vectorized,single): - eu = np.array([rot.as_Eulers() for rot in default]) - vectorized(eu.reshape(eu.shape[0]//2,-1,3)) - co = vectorized(eu) - for e,c in zip(eu,co): - print(e,c) - assert np.allclose(single(e),c) and np.allclose(single(e),vectorized(e)) - - @pytest.mark.parametrize('vectorized, single',[(Rotation._ax2qu,ax2qu), - (Rotation._ax2om,ax2om), - (Rotation._ax2ro,ax2ro), - (Rotation._ax2ho,ax2ho)]) - def test_axisAngle_vectorization(self,default,vectorized,single): - ax = np.array([rot.as_axis_angle() for rot in default]) - vectorized(ax.reshape(ax.shape[0]//2,-1,4)) - co = vectorized(ax) - for a,c in zip(ax,co): - print(a,c) - assert np.allclose(single(a),c) and np.allclose(single(a),vectorized(a)) - - - @pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax), - (Rotation._ro2ho,ro2ho)]) - def test_Rodrigues_vectorization(self,default,vectorized,single): - ro = np.array([rot.as_Rodrigues() for rot in default]) - vectorized(ro.reshape(ro.shape[0]//2,-1,4)) - co = vectorized(ro) - for r,c in zip(ro,co): - print(r,c) - assert np.allclose(single(r),c) and np.allclose(single(r),vectorized(r)) - - @pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax), - (Rotation._ho2cu,ho2cu)]) - def test_homochoric_vectorization(self,default,vectorized,single): - ho = np.array([rot.as_homochoric() for rot in default]) - vectorized(ho.reshape(ho.shape[0]//2,-1,3)) - co = vectorized(ho) - for h,c in zip(ho,co): - print(h,c) - assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h)) - - @pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)]) - def test_cubochoric_vectorization(self,default,vectorized,single): - cu = np.array([rot.as_cubochoric() for rot in default]) - vectorized(cu.reshape(cu.shape[0]//2,-1,3)) - co = vectorized(cu) - for u,c in zip(cu,co): - print(u,c) - assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u)) - @pytest.mark.parametrize('direction',['forward', 'backward']) def test_pyramid_vectorization(self,direction): From 128a96e7f6f44c5ff6dafa26ba7536b055dad3f9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 18:13:51 +0200 Subject: [PATCH 30/32] vectorized formula is enough --- python/damask/_rotation.py | 47 +++----------------------------------- 1 file changed, 3 insertions(+), 44 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 4081f2844..3f0cc7107 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -77,9 +77,9 @@ class Rotation: ]) - def __mul__(self, other): + def __matmul__(self, other): """ - Multiplication. + Rotation of vector, second or fourth order tensor, or rotation object. Parameters ---------- @@ -89,47 +89,6 @@ class Rotation: Todo ---- Document details active/passive) - consider rotation of (3,3,3,3)-matrix - - """ - if self.quaternion.shape != (4,): - raise NotImplementedError('Support for multiple rotations missing') - if isinstance(other, Rotation): - self_q = self.quaternion[0] - self_p = self.quaternion[1:] - other_q = other.quaternion[0] - other_p = other.quaternion[1:] - R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p), - self_q*other_p + other_q*self_p + _P * np.cross(self_p,other_p))) - return R.standardize() - elif isinstance(other, np.ndarray): - if other.shape == (3,): - A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:]) - B = 2.0 * np.dot(self.quaternion[1:],other) - C = 2.0 * _P*self.quaternion[0] - - return A*other + B*self.quaternion[1:] + C * np.cross(self.quaternion[1:],other) - - elif other.shape == (3,3,): - R = self.as_matrix() - return np.dot(R,np.dot(other,R.T)) - elif other.shape == (3,3,3,3,): - R = self.as_matrix() - RR = np.outer(R, R) - RRRR = np.outer(RR, RR).reshape(4 * (3,3)) - axes = ((0, 2, 4, 6), (0, 1, 2, 3)) - return np.tensordot(RRRR, other, axes) - else: - raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') - else: - raise TypeError('Cannot rotate {}'.format(type(other))) - - - def __matmul__(self, other): - """ - Rotation. - - details to be discussed """ if isinstance(other, Rotation): q_m = self.quaternion[...,0:1] @@ -505,7 +464,7 @@ class Rotation: fromEulers = from_Eulers asAxisAngle = as_axis_angle asRodrigues = as_Rodrigues - + __mul__ = __matmul__ #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations From 353fd3ceb602f52c027b4b9359dac24dcca761f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 19:10:16 +0200 Subject: [PATCH 31/32] more tests now 95% test coverage of Rotation class --- python/damask/_rotation.py | 3 +- python/tests/test_Rotation.py | 85 ++++++++++++++++++++++++++++++++--- 2 files changed, 81 insertions(+), 7 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 3f0cc7107..cac0cec3e 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -88,7 +88,8 @@ class Rotation: Todo ---- - Document details active/passive) + Check rotation of 4th order tensor + """ if isinstance(other, Rotation): q_m = self.quaternion[...,0:1] diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index f424fd9ab..a9768afb8 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -510,6 +510,54 @@ def _get_pyramid_order(xyz,direction=None): #################################################################################################### #################################################################################################### +def mul(me, other): + """ + Multiplication. + + Parameters + ---------- + other : numpy.ndarray or Rotation + Vector, second or fourth order tensor, or rotation object that is rotated. + + Todo + ---- + Document details active/passive) + consider rotation of (3,3,3,3)-matrix + + """ + if me.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') + if isinstance(other, Rotation): + me_q = me.quaternion[0] + me_p = me.quaternion[1:] + other_q = other.quaternion[0] + other_p = other.quaternion[1:] + R = me.__class__(np.append(me_q*other_q - np.dot(me_p,other_p), + me_q*other_p + other_q*me_p + _P * np.cross(me_p,other_p))) + return R.standardize() + elif isinstance(other, np.ndarray): + if other.shape == (3,): + A = me.quaternion[0]**2.0 - np.dot(me.quaternion[1:],me.quaternion[1:]) + B = 2.0 * np.dot(me.quaternion[1:],other) + C = 2.0 * _P*me.quaternion[0] + + return A*other + B*me.quaternion[1:] + C * np.cross(me.quaternion[1:],other) + + elif other.shape == (3,3,): + R = me.as_matrix() + return np.dot(R,np.dot(other,R.T)) + elif other.shape == (3,3,3,3,): + R = me.as_matrix() + return np.einsum('ia,jb,kc,ld,abcd->ijkl',R,R,R,R,other) + RR = np.outer(R, R) + RRRR = np.outer(RR, RR).reshape(4 * (3,3)) + axes = ((0, 2, 4, 6), (0, 1, 2, 3)) + return np.tensordot(RRRR, other, axes) + else: + raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') + else: + raise TypeError('Cannot rotate {}'.format(type(other))) + class TestRotation: @@ -616,7 +664,7 @@ class TestRotation: @pytest.mark.parametrize('forward,backward',[(Rotation._cu2qu,Rotation._qu2cu), (Rotation._cu2om,Rotation._om2cu), - #(Rotation._cu2eu,Rotation._eu2cu), + (Rotation._cu2eu,Rotation._eu2cu), (Rotation._cu2ax,Rotation._ax2cu), (Rotation._cu2ro,Rotation._ro2cu), (Rotation._cu2ho,Rotation._ho2cu)]) @@ -626,6 +674,8 @@ class TestRotation: m = rot.as_cubochoric() o = backward(forward(m)) ok = np.allclose(m,o,atol=atol) + if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)): + ok = ok or np.allclose(m*-1.,o,atol=atol) print(m,o,rot.as_quaternion()) assert ok and np.max(np.abs(o)) < np.pi**(2./3.) * 0.5 + 1.e-9 @@ -785,21 +835,32 @@ class TestRotation: assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9 @pytest.mark.parametrize('P',[1,-1]) - def test_quaternion(self,default,P): - c = np.array([1,P*-1,P*-1,P*-1]) + @pytest.mark.parametrize('accept_homomorph',[True,False]) + def test_quaternion(self,default,P,accept_homomorph): + c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1) for rot in default: m = rot.as_cubochoric() - o = Rotation.from_quaternion(rot.as_quaternion()*c,False,P).as_cubochoric() + o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric() ok = np.allclose(m,o,atol=atol) + if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)): + ok = ok or np.allclose(m*-1.,o,atol=atol) print(m,o,rot.as_quaternion()) assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 + @pytest.mark.parametrize('reciprocal',[True,False]) + def test_basis(self,default,reciprocal): + for rot in default: + om = rot.as_matrix() + 0.1*np.eye(3) + rot = Rotation.from_basis(om,False,reciprocal=reciprocal) + assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) + @pytest.mark.parametrize('function',[Rotation.from_quaternion, Rotation.from_Eulers, Rotation.from_axis_angle, Rotation.from_matrix, Rotation.from_Rodrigues, - Rotation.from_homochoric]) + Rotation.from_homochoric, + Rotation.from_cubochoric]) def test_invalid_shape(self,function): invalid_shape = np.random.random(np.random.randint(8,32,(3))) with pytest.raises(ValueError): @@ -813,7 +874,8 @@ class TestRotation: (Rotation.from_matrix, np.random.rand(3,3)), (Rotation.from_Rodrigues, np.array([1,0,0,-1])), (Rotation.from_Rodrigues, np.array([1,1,0,1])), - (Rotation.from_homochoric, np.array([2,2,2])) ]) + (Rotation.from_homochoric, np.array([2,2,2])), + (Rotation.from_cubochoric, np.array([1.1,0,0])) ]) def test_invalid_value(self,function,invalid): with pytest.raises(ValueError): function(invalid) @@ -833,6 +895,17 @@ class TestRotation: assert np.all(np.take_along_axis(np.take_along_axis(a,f,-1),b,-1) == a) + @pytest.mark.parametrize('data',[np.random.rand(5,3), + np.random.rand(5,3,3), + np.random.rand(5,3,3,3,3)]) + def test_rotate_vectorization(self,default,data): + for rot in default: + v = rot.broadcast_to((5,)) @ data + for i in range(data.shape[0]): + print(i-data[i]) + assert np.allclose(mul(rot,data[i]),v[i]) + + @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) From 16f30a10d0f85bdd277b372f2cfbee47e69d3e2a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 May 2020 23:50:08 +0200 Subject: [PATCH 32/32] polishing --- python/damask/_orientation.py | 8 ++++---- python/damask/_rotation.py | 3 +-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index c145dc42b..7c3c614fb 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -65,8 +65,8 @@ class Orientation: r = b*aInv for k in range(2): r.inverse() - breaker = self.lattice.symmetry.inFZ(r.asRodrigues(vector=True)) \ - and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues(vector=True))) + breaker = self.lattice.symmetry.inFZ(r.as_Rodrigues(vector=True)) \ + and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True))) if breaker: break if breaker: break if breaker: break @@ -75,7 +75,7 @@ class Orientation: # ... own sym, other sym, # self-->other: True, self<--other: False def inFZ(self): - return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True)) + return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True)) def equivalentOrientations(self,members=[]): @@ -97,7 +97,7 @@ class Orientation: def reduced(self): """Transform orientation to fall into fundamental zone according to symmetry.""" for me in self.equivalentOrientations(): - if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break + if self.lattice.symmetry.inFZ(me.rotation.as_Rodrigues(vector=True)): break return self.__class__(me.rotation,self.lattice) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index cac0cec3e..387ee84cb 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -119,7 +119,7 @@ class Rotation: R = self.as_matrix() return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other) else: - raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') + raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors') else: raise TypeError('Cannot rotate {}'.format(type(other))) @@ -464,7 +464,6 @@ class Rotation: fromQuaternion = from_quaternion fromEulers = from_Eulers asAxisAngle = as_axis_angle - asRodrigues = as_Rodrigues __mul__ = __matmul__ ####################################################################################################