From 025b2db70e8808f8fc8d98332616e6f196dd251c Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 Nov 2015 21:37:27 -0800 Subject: [PATCH 1/3] coordinate rotations --- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/coordutils.py | 80 ++++++++++++++++++++++++++++++++++ tests/utils/test_coordutils.py | 35 +++++++++++++++ 3 files changed, 116 insertions(+), 1 deletion(-) create mode 100644 SimPEG/Utils/coordutils.py create mode 100644 tests/utils/test_coordutils.py diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 6637a138..5280ae79 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,4 +7,4 @@ from ipythonutils import easyAnimate as animate from CounterUtils import * import ModelBuilder import SolverUtils - +from coordutils import * diff --git a/SimPEG/Utils/coordutils.py b/SimPEG/Utils/coordutils.py new file mode 100644 index 00000000..e26a751e --- /dev/null +++ b/SimPEG/Utils/coordutils.py @@ -0,0 +1,80 @@ +import numpy as np + +def crossProd(v0,v1): + """ + Cross product of 2 vectors + + :param numpy.array v0: vector of length 3 + :param numpy.array v1: vector of length 3 + :rtype: numpy.array + :return: cross product of v0,v1 + """ + # ensure both n0, n1 are vectors of length 1 + assert len(v0) == 3, "Length of v0 should be 3" + assert len(v1) == 3, "Length of v1 should be 3" + + v2 = np.zeros(3,dtype=float) + + v2[0] = v0[1]*v1[2] - v1[1]*v0[2] + v2[1] = v1[0]*v0[2] - v0[0]*v1[2] + v2[2] = v0[0]*v1[1] - v1[0]*v0[1] + + return v2 + +def rotationMatrixFromNormals(v0,v1,tol=1e-20): + """ + Performs the minimum number of rotations to define a rotation from the direction indicated by the vector n0 to the direction indicated by n1. + The axis of rotation is n0 x n1 + https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula + + :param numpy.array v0: vector of length 3 + :param numpy.array v1: vector of length 3 + :param tol = 1e-20: tolerance. If the norm of the cross product between the two vectors is below this, no rotation is performed + :rtype: numpy.array, 3x3 + :return: rotation matrix which rotates the frame so that n0 is aligned with n1 + + """ + + # ensure both n0, n1 are vectors of length 1 + assert len(v0) == 3, "Length of n0 should be 3" + assert len(v1) == 3, "Length of n1 should be 3" + + # ensure both are true normals + n0 = v0*1./np.linalg.norm(v0) + n1 = v1*1./np.linalg.norm(v1) + + n0dotn1 = n0.dot(n1) + + # define the rotation axis, which is the cross product of the two vectors + rotAx = crossProd(n0,n1) + + if np.linalg.norm(rotAx) < tol: + return np.eye(3,dtype=float) + + rotAx *= 1./np.linalg.norm(rotAx) + + cosT = n0dotn1/(np.linalg.norm(n0)*np.linalg.norm(n1)) + sinT = np.sqrt(1.-n0dotn1**2) + + ux = np.array([[0., -rotAx[2], rotAx[1]], [rotAx[2], 0., -rotAx[0]], [-rotAx[1], rotAx[0], 0.]],dtype=float) + + return np.eye(3,dtype=float) + sinT*ux + (1.-cosT)*(ux.dot(ux)) + + +def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]): + """ + rotates a grid so that the vector n0 is aligned with the vector n1 + + :param numpy.array n0: vector of length 3, should have norm 1 + :param numpy.array n1: vector of length 3, should have norm 1 + :param numpy.array x0: vector of length 3, point about which we perform the rotation + :rtype: numpy.array, 3x3 + :return: rotation matrix which rotates the frame so that n0 is aligned with n1 + """ + + R = rotationMatrixFromNormals(n0, n1) + + assert XYZ.shape[1] == 3, "Grid XYZ should be 3 wide" + assert len(x0) == 3, "x0 should have length 3" + + return (XYZ - x0).dot(R.T) + x0 \ No newline at end of file diff --git a/tests/utils/test_coordutils.py b/tests/utils/test_coordutils.py new file mode 100644 index 00000000..7770671d --- /dev/null +++ b/tests/utils/test_coordutils.py @@ -0,0 +1,35 @@ +import unittest, os +import numpy as np +from SimPEG import Utils + +tol = 1e-15 + +class coorUtilsTest(unittest.TestCase): + + def test_crossProd(self): + self.assertTrue(np.linalg.norm(Utils.coordutils.crossProd(np.r_[1.,0.,0.],np.r_[0.,1.,0.]) - np.r_[0.,0.,1.]) < tol) + self.assertTrue(np.linalg.norm(Utils.coordutils.crossProd(np.r_[0.,1.,0.],np.r_[0.,0.,1.]) - np.r_[1.,0.,0.]) < tol) + self.assertTrue(np.linalg.norm(Utils.coordutils.crossProd(np.r_[0.,0.,1.],np.r_[1.,0.,0.]) - np.r_[0.,1.,0.]) < tol) + + def test_rotationMatrixFromNormals(self): + v0 = np.random.rand(3) + v0 *= 1./np.linalg.norm(v0) + v1 = np.random.rand(3) + v1 *= 1./np.linalg.norm(v1) + Rf = Utils.coordutils.rotationMatrixFromNormals(v0,v1) + Ri = Utils.coordutils.rotationMatrixFromNormals(v1,v0) + + self.assertTrue(np.linalg.norm(Utils.mkvc(Rf.dot(v0) - v1)) < tol) + self.assertTrue(np.linalg.norm(Utils.mkvc(Ri.dot(v1) - v0)) < tol) + + def test_rotatePointsFromNormals(self): + v0 = np.random.rand(3) + v0*= 1./np.linalg.norm(v0) + v1 = np.random.rand(3) + v1*= 1./np.linalg.norm(v1) + + self.assertTrue(np.linalg.norm(Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1))-v1) < tol) + +if __name__ == '__main__': + unittest.main() + From 96b855d71d9d6f24a8e78dafe33f0602cec5d21c Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 5 Nov 2015 16:10:43 -0800 Subject: [PATCH 2/3] use np.cross --- SimPEG/Utils/coordutils.py | 23 +---------------------- tests/utils/test_coordutils.py | 5 ----- 2 files changed, 1 insertion(+), 27 deletions(-) diff --git a/SimPEG/Utils/coordutils.py b/SimPEG/Utils/coordutils.py index e26a751e..0329495a 100644 --- a/SimPEG/Utils/coordutils.py +++ b/SimPEG/Utils/coordutils.py @@ -1,26 +1,5 @@ import numpy as np -def crossProd(v0,v1): - """ - Cross product of 2 vectors - - :param numpy.array v0: vector of length 3 - :param numpy.array v1: vector of length 3 - :rtype: numpy.array - :return: cross product of v0,v1 - """ - # ensure both n0, n1 are vectors of length 1 - assert len(v0) == 3, "Length of v0 should be 3" - assert len(v1) == 3, "Length of v1 should be 3" - - v2 = np.zeros(3,dtype=float) - - v2[0] = v0[1]*v1[2] - v1[1]*v0[2] - v2[1] = v1[0]*v0[2] - v0[0]*v1[2] - v2[2] = v0[0]*v1[1] - v1[0]*v0[1] - - return v2 - def rotationMatrixFromNormals(v0,v1,tol=1e-20): """ Performs the minimum number of rotations to define a rotation from the direction indicated by the vector n0 to the direction indicated by n1. @@ -46,7 +25,7 @@ def rotationMatrixFromNormals(v0,v1,tol=1e-20): n0dotn1 = n0.dot(n1) # define the rotation axis, which is the cross product of the two vectors - rotAx = crossProd(n0,n1) + rotAx = np.cross(n0,n1) if np.linalg.norm(rotAx) < tol: return np.eye(3,dtype=float) diff --git a/tests/utils/test_coordutils.py b/tests/utils/test_coordutils.py index 7770671d..dc036529 100644 --- a/tests/utils/test_coordutils.py +++ b/tests/utils/test_coordutils.py @@ -6,11 +6,6 @@ tol = 1e-15 class coorUtilsTest(unittest.TestCase): - def test_crossProd(self): - self.assertTrue(np.linalg.norm(Utils.coordutils.crossProd(np.r_[1.,0.,0.],np.r_[0.,1.,0.]) - np.r_[0.,0.,1.]) < tol) - self.assertTrue(np.linalg.norm(Utils.coordutils.crossProd(np.r_[0.,1.,0.],np.r_[0.,0.,1.]) - np.r_[1.,0.,0.]) < tol) - self.assertTrue(np.linalg.norm(Utils.coordutils.crossProd(np.r_[0.,0.,1.],np.r_[1.,0.,0.]) - np.r_[0.,1.,0.]) < tol) - def test_rotationMatrixFromNormals(self): v0 = np.random.rand(3) v0 *= 1./np.linalg.norm(v0) From 7b0bea4e1d4cac9c3fd6fffa1cdf5e7f96adaa41 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 5 Nov 2015 16:22:31 -0800 Subject: [PATCH 3/3] make sure subtraction of x0 is a matrix with length of number of points, add test for giving an array XYZ --- SimPEG/Utils/coordutils.py | 5 ++++- tests/utils/test_coordutils.py | 17 ++++++++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/SimPEG/Utils/coordutils.py b/SimPEG/Utils/coordutils.py index 0329495a..260e1a3b 100644 --- a/SimPEG/Utils/coordutils.py +++ b/SimPEG/Utils/coordutils.py @@ -1,4 +1,5 @@ import numpy as np +from SimPEG.Utils import mkvc def rotationMatrixFromNormals(v0,v1,tol=1e-20): """ @@ -56,4 +57,6 @@ def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]): assert XYZ.shape[1] == 3, "Grid XYZ should be 3 wide" assert len(x0) == 3, "x0 should have length 3" - return (XYZ - x0).dot(R.T) + x0 \ No newline at end of file + X0 = np.ones([XYZ.shape[0],1])*mkvc(x0) + + return (XYZ - X0).dot(R.T) + X0 # equivalent to (R*(XYZ - X0)).T + X0 \ No newline at end of file diff --git a/tests/utils/test_coordutils.py b/tests/utils/test_coordutils.py index dc036529..b17afcf7 100644 --- a/tests/utils/test_coordutils.py +++ b/tests/utils/test_coordutils.py @@ -23,7 +23,22 @@ class coorUtilsTest(unittest.TestCase): v1 = np.random.rand(3) v1*= 1./np.linalg.norm(v1) - self.assertTrue(np.linalg.norm(Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1))-v1) < tol) + v2 = Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1)) + + self.assertTrue(np.linalg.norm(v2-v1) < tol) + + def test_rotateMatrixFromNormals(self): + n0 = np.random.rand(3) + n0*= 1./np.linalg.norm(n0) + n1 = np.random.rand(3) + n1*= 1./np.linalg.norm(n1) + + scale = np.random.rand(100,1) + XYZ0 = scale * n0 + XYZ1 = scale * n1 + + XYZ2 = Utils.coordutils.rotatePointsFromNormals(XYZ0,n0,n1) + self.assertTrue(np.linalg.norm(Utils.mkvc(XYZ1) - Utils.mkvc(XYZ2))/np.linalg.norm(Utils.mkvc(XYZ1)) < tol) if __name__ == '__main__': unittest.main()