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'])