From e931b716fd6347f64161d50704c0c4a8184b8720 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 11 Feb 2019 23:11:11 +0100 Subject: [PATCH] conversion routines from Marc de Graefs 3D rotation repository Python version available on https://github.com/MarDiehl/3Drotations --- python/damask/Lambert.py | 122 ++++++++++ python/damask/orientation.py | 443 +++++++++++++++++++++++++++++++++++ 2 files changed, 565 insertions(+) create mode 100644 python/damask/Lambert.py diff --git a/python/damask/Lambert.py b/python/damask/Lambert.py new file mode 100644 index 000000000..9972b7965 --- /dev/null +++ b/python/damask/Lambert.py @@ -0,0 +1,122 @@ +#################################################################################################### +# Code below available according to below 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 + +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 CubeToBall(cube): + + if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5: + raise ValueError + + # transform to the sphere grid via the curved square, and intercept the zero point + if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300): + ball = np.zeros(3) + else: + # get pyramide and scale by grid parameter ratio + p = GetPyramidOrder(cube) + XYZ = cube[p] * sc + + # intercept all the points along the z-axis + if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300): + ball = 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 ) + ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) + + # reverse the coordinates back to the regular order according to the original pyramid number + ball = ball[p] + + return ball + + +def BallToCube(ball): + + rs = np.linalg.norm(ball) + if rs > R1: + raise ValueError + + if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300): + cube = np.zeros(3) + else: + p = GetPyramidOrder(ball) + xyz3 = ball[p] + + # 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-300): + 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 + cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc + + # reverst the coordinates back to the regular order according to the original pyramid number + cube = cube[p] + + return cube + +def GetPyramidOrder(xyz): + + if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \ + (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]): + return [0,1,2] + elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \ + (abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]): + return [1,2,0] + elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \ + (abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]): + return [2,0,1] diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 1bc850734..73b7620ca 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -6,6 +6,9 @@ import math,os import numpy as np +from . import Lambert + +P = -1 # ****************************************************************************************** class Quaternion: @@ -1093,3 +1096,443 @@ class Orientation: rot=np.dot(otherMatrix,myMatrix.T) return Orientation(matrix=np.dot(rot,self.asMatrix()),symmetry=targetSymmetry) + +#################################################################################################### +# Code below available according to below 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. +#################################################################################################### + +def isone(a): + return np.isclose(a,1.0,atol=1.0e-15,rtol=0.0) + +def iszero(a): + return np.isclose(a,0.0,atol=1.0e-300,rtol=0.0) + + +def eu2om(eu): + """Euler angles to orientation 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.where(iszero(om))] = 0.0 + return om + + +def eu2ax(eu): + """Euler angles to axis angle""" + 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 iszero(alpha): + 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): + """Euler angles to Rodrigues vector""" + ro = eu2ax(eu) # convert to axis angle 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 + + +def eu2qu(eu): + """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.homomorph() !ToDo: Check with original + return qu + + +def om2eu(om): + """Euler angles to orientation matrix""" + if isone(om[2,2]**2): + 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: + 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)]) + + # reduce Euler angles to definition range, i.e a lower limit of 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 ax2om(ax): + """Axis angle to orientation 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 qu2eu(qu): + """Quaternion to Euler angles""" + q03 = qu[0]**2+qu[3]**2 + q12 = qu[1]**2+qu[2]**2 + chi = np.sqrt(q03*q12) + + if iszero(chi): + eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \ + np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) + else: + #chiInv = 1.0/chi ToDo: needed for what? + 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, i.e a lower limit of 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 ax2ho(ax): + """Axis angle to homochoric""" + f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) + ho = ax[0:3] * f + return ho + + +def ho2ax(ho): + """Homochoric to axis angle""" + 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(s)) # ToDo: Check sanity check in reference implementation + + return ax + + +def om2ax(om): + """Orientation matrix to axis angle""" + 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 iszero(ax[3]): + ax = [ 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 = 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(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) + + return np.array(ax) + + +def ro2ax(ro): + """Rodrigues vector to axis angle""" + ta = ro[3] + + if iszero(ta): + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elif not np.isfinite(ta): + ax = [ ro[0], ro[1], ro[2], np.pi ] + else: + angle = 2.0*np.arctan(ta) + ta = 1.0/np.linalg.norm(ro[0:3]) + ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] + + return np.array(ax) + + +def ax2ro(ax): + """Axis angle to Rodrigues vector""" + if iszero(ax[3]): + 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)] + + return np.array(ro) + + +def ax2qu(ax): + """Axis angle to quaternion""" + if iszero(ax[3]): + 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 ro2ho(ro): + """Rodrigues vector to homochoric""" + if iszero(np.sum(ro[0:3]**2.0)): + ho = [ 0.0, 0.0, 0.0 ] + 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 np.array(ho) + + +def qu2om(qu): + """Quaternion to orientation 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[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) + om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) + om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) + om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) + om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) + om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) + return om if P > 0.0 else om.T + + +def om2qu(om): + """Orientation matrix to quaternion""" + s = [+om[0,0] +om[1,1] +om[2,2] +1.0, + +om[0,0] -om[1,1] -om[2,2] +1.0, + -om[0,0] +om[1,1] -om[2,2] +1.0, + -om[0,0] -om[1,1] +om[2,2] +1.0] + s = np.maximum(np.zeros(4),s) + qu = np.sqrt(s)*0.5*np.array([1.0,P,P,P]) + # verify the signs (q0 always positive) + #ToDo: Here I donot understand the original shortcut from paper to implementation + + qu /= np.linalg.norm(qu) + if any(isone(abs(qu))): qu[np.where(np.logical_not(isone(qu)))] = 0.0 + if om[2,1] < om[1,2]: qu[1] *= -1.0 + if om[0,2] < om[2,0]: qu[2] *= -1.0 + if om[1,0] < om[0,1]: qu[3] *= -1.0 + if any(om2ax(om)[0:3]*qu[1:4] < 0.0): print(om2ax(om),qu) # something is wrong here + return qu + +def qu2ax(qu): + """Quaternion to axis angle""" + omega = 2.0 * np.arccos(qu[0]) + if iszero(omega): # return axis as [001] if the angle is zero + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elif not iszero(qu[0]): + s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) + ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] + else: + ax = [ qu[1], qu[2], qu[3], np.pi] + + return np.array(ax) + + +def qu2ro(qu): + """Quaternion to Rodrigues vector""" + if iszero(qu[0]): + ro = [qu[1], qu[2], qu[3], np.inf] + else: + s = np.linalg.norm([qu[1],qu[2],qu[3]]) + ro = [0.0,0.0,P,0.0] if iszero(s) else \ + [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(qu[0]))] + + return np.array(ro) + + +def qu2ho(qu): + """Quaternion to homochoric""" + omega = 2.0 * np.arccos(qu[0]) + + if iszero(omega): + ho = np.array([ 0.0, 0.0, 0.0 ]) + 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 + + +def ho2cu(ho): + """Homochoric to cubochoric""" + return Lambert.BallToCube(ho) + + +def cu2ho(cu): + """Cubochoric to homochoric""" + return Lambert.CubeToBall(cu) + + +def ro2eu(ro): + """Rodrigues vector to orientation matrix""" + return om2eu(ro2om(ro)) + + +def eu2ho(eu): + """Euler angles to homochoric""" + return ax2ho(eu2ax(eu)) + + +def om2ro(om): + """Orientation matrix to Rodriques vector""" + return eu2ro(om2eu(om)) + + +def om2ho(om): + """Orientation matrix to homochoric""" + return ax2ho(om2ax(om)) + + +def ax2eu(ax): + """Orientation matrix to Euler angles""" + return om2eu(ax2om(ax)) + + +def ro2om(ro): + """Rodgrigues vector to orientation matrix""" + return ax2om(ro2ax(ro)) + + +def ro2qu(ro): + """Rodrigues vector to quaternion""" + return ax2qu(ro2ax(ro)) + + +def ho2eu(ho): + """Homochoric to Euler angles""" + return ax2eu(ho2ax(ho)) + + +def ho2om(ho): + """Homochoric to orientation matrix""" + return ax2om(ho2ax(ho)) + + +def ho2ro(ho): + """Axis angle to Rodriques vector""" + return ax2ro(ho2ax(ho)) + + +def ho2qu(ho): + """Homochoric to quaternion""" + return ax2qu(ho2ax(ho)) + + +def eu2cu(eu): + """Euler angles to cubochoric""" + return ho2cu(eu2ho(eu)) + + +def om2cu(om): + """Orientation matrix to cubochoric""" + return ho2cu(om2ho(om)) + + +def ax2cu(ax): + """Axis angle to cubochoric""" + return ho2cu(ax2ho(ax)) + + +def ro2cu(ro): + """Rodrigues vector to cubochoric""" + return ho2cu(ro2ho(ro)) + + +def qu2cu(qu): + """Quaternion to cubochoric""" + return ho2cu(qu2ho(qu)) + + +def cu2eu(cu): + """Cubochoric to Euler angles""" + return ho2eu(cu2ho(cu)) + + +def cu2om(cu): + """Cubochoric to orientation matrix""" + return ho2om(cu2ho(cu)) + + +def cu2ax(cu): + """Cubochoric to axis angle""" + return ho2ax(cu2ho(cu)) + + +def cu2ro(cu): + """Cubochoric to Rodrigues vector""" + return ho2ro(cu2ho(cu)) + + +def cu2qu(cu): + """Cubochoric to quaternion""" + return ho2qu(cu2ho(cu))