Skip to content

Commit 11a5225

Browse files
robertwbrobertwb
committed
Cleaner plot3d rotation around arbitrary axis
1 parent ff2f4de commit 11a5225

File tree

1 file changed

+128
-28
lines changed

1 file changed

+128
-28
lines changed

‎src/sage/plot/plot3d/transform.pyx

+128-28
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,9 @@
1-
from math import atan2, sin, cos, atan, sqrt, acos
1+
cdef extern from *:
2+
double sin(double)
3+
double cos(double)
4+
double sqrt(double)
5+
6+
#from math import atan2, sin, cos, atan, sqrt, acos
27

38
include "point_c.pxi"
49

@@ -34,23 +39,7 @@ cdef class Transformation:
3439
if rot is not None:
3540
# rotate about v by theta
3641
vx, vy, vz, theta = rot
37-
if vx == vy == 0:
38-
t = RDF(pi/2)
39-
else:
40-
t = atan2(vy,vz) + pi/2
41-
m = self.rotX(t) * m
42-
new_y = vy*cos(t) - vz*sin(t)
43-
# v = [vx, new_y, 0]
44-
if vx == new_y == 0:
45-
s = RDF(pi/2)
46-
else:
47-
s = atan2(vx,new_y) + pi/2
48-
m = self.rotZ(s) * m
49-
# v = [new_x, 0, 0]
50-
m = self.rotX(theta) * m
51-
# now put back to our former reference frame
52-
m = self.rotZ(-s) * m
53-
m = self.rotX(-t) * m
42+
m *= rotate_arbitrary((vx, vy, vz), theta)
5443

5544
self.matrix = m.augment(matrix(RDF, 3, 1, list(trans))) \
5645
.stack(matrix(RDF, 1, 4, [0,0,0,1]))
@@ -61,6 +50,8 @@ cdef class Transformation:
6150
for i from 0 <= i < 12:
6251
self._matrix_data[i] = m_data[i]
6352

53+
def get_matrix(self):
54+
return self.matrix.copy()
6455

6556
def is_skew(self, eps=1e-5):
6657
dx, dy, dz = self.matrix.submatrix(0,0,3,3).columns()
@@ -77,16 +68,6 @@ cdef class Transformation:
7768
len_a = a.dot_product(a)
7869
return max([abs(len_a - b.dot_product(b)) for b in basis]) < eps
7970

80-
def rotX(self, theta):
81-
return matrix(RDF, 3, 3, [1, 0, 0,
82-
0, cos(theta), -sin(theta),
83-
0, sin(theta), cos(theta)])
84-
85-
def rotZ(self, theta):
86-
return matrix(RDF, 3, 3, [cos(theta), -sin(theta), 0,
87-
sin(theta), cos(theta), 0,
88-
0, 0, 1])
89-
9071
def transform_point(self, x):
9172
Tx = self.matrix * vector(RDF, [x[0], x[1], x[2], 1])
9273
return (Tx[0], Tx[1], Tx[2])
@@ -112,3 +93,122 @@ cdef class Transformation:
11293
if self._svd is None:
11394
self._svd = self.matrix.submatrix(0,0,3,3).SVD()
11495
return self._svd[1][0,0]
96+
97+
98+
def rotate_arbitrary(v, double theta):
99+
"""
100+
Return a matrix that rotates the coordinate space about
101+
the axis v by the angle theta.
102+
103+
EXAMPLES:
104+
sage: from sage.plot.plot3d.transform import rotate_arbitrary
105+
106+
Try rotating about the axes:
107+
sage: rotate_arbitrary((1,0,0), 1)
108+
[ 1.0 0.0 0.0]
109+
[ 0.0 0.540302305868 0.841470984808]
110+
[ 0.0 -0.841470984808 0.540302305868]
111+
sage: rotate_arbitrary((0,1,0), 1)
112+
[ 0.540302305868 0.0 -0.841470984808]
113+
[ 0.0 1.0 0.0]
114+
[ 0.841470984808 0.0 0.540302305868]
115+
sage: rotate_arbitrary((0,0,1), 1)
116+
[ 0.540302305868 0.841470984808 0.0]
117+
[-0.841470984808 0.540302305868 0.0]
118+
[ 0.0 0.0 1.0]
119+
120+
These next two should be the same (up to machine epsilon)
121+
sage: rotate_arbitrary((1,1,1), 1)
122+
[ 0.693534870579 0.639056064305 -0.332590934883]
123+
[-0.332590934883 0.693534870579 0.639056064305]
124+
[ 0.639056064305 -0.332590934883 0.693534870579]
125+
sage: rotate_arbitrary((1,1,1), -1)^(-1)
126+
[ 0.693534870579 0.639056064305 -0.332590934883]
127+
[-0.332590934883 0.693534870579 0.639056064305]
128+
[ 0.639056064305 -0.332590934883 0.693534870579]
129+
130+
Make sure it does the right thing...
131+
sage: rotate_arbitrary((1,2,3), -1).det()
132+
1.0
133+
sage: rotate_arbitrary((1,1,1), 2*pi/3) * vector(RDF, (1,2,3))
134+
(2.0, 3.0, 1.0)
135+
sage: rotate_arbitrary((1,2,3), 5) * vector(RDF, (1,2,3))
136+
(1.0, 2.0, 3.0)
137+
sage: rotate_arbitrary((1,1,1), pi/7)^7
138+
[-0.333333333333 0.666666666667 0.666666666667]
139+
[ 0.666666666667 -0.333333333333 0.666666666667]
140+
[ 0.666666666667 0.666666666667 -0.333333333333]
141+
142+
143+
AUTHORS:
144+
-- Robert Bradshaw
145+
146+
147+
ALGORITHM:
148+
There is a formula. Where did it come from? Lets take
149+
a quick jaunt into SAGE's calculus package...
150+
151+
Setup some variables
152+
sage: vx,vy,vz,theta = var('x y z theta')
153+
154+
Symbolic rotation matrices about X and Y axis:
155+
sage: def rotX(theta): return matrix(SR, 3, 3, [1, 0, 0, 0, cos(theta), -sin(theta), 0, sin(theta), cos(theta)])
156+
sage: def rotZ(theta): return matrix(SR, 3, 3, [cos(theta), -sin(theta), 0, sin(theta), cos(theta), 0, 0, 0, 1])
157+
158+
Normalizing $y$ so that $|v|=1$. Perhaps there is a better
159+
way to tell maxima that $x^2+y^2+z^2=1$ which would make for
160+
a much cleaner calculation.
161+
sage: vy = sqrt(1-vx^2-vz^2)
162+
163+
Now we rotate about the $x$-axis so $v$ is in the $xy$-plane.
164+
sage: t = atan(vy/vz)+pi/2
165+
sage: m = rotX(t)
166+
sage: new_y = vy*cos(t) - vz*sin(t)
167+
168+
And rotate about the $z$ axis so $v$ lies on the $x$ axis.
169+
sage: s = atan(vx/new_y) + pi/2
170+
sage: m = rotZ(s) * m
171+
172+
Rotating about $v$ in our old system is the same as rotating
173+
about the $x$-axis in the new.
174+
sage: m = rotX(theta) * m
175+
176+
Do some simplfying here to avoid blow-up.
177+
sage: ix = [(i,j) for i in range(3) for j in range(3)]
178+
sage: for ij in ix: m[ij] = m[ij].simplify_rational()
179+
180+
Now go back to the original coordinate system.
181+
sage: m = rotZ(-s) * m
182+
sage: m = rotX(-t) * m
183+
sage: for ij in ix: m[ij] = m[ij].simplify_rational()
184+
sage: m # or show(m)
185+
[ (1 - cos(theta))*x^2 + cos(theta) -(sin(theta)*abs(z)^3 + (cos(theta) - 1)*x*z^2*sqrt(-z^2 - x^2 + 1))/z^2 (sin(theta)*sqrt(-z^2 - x^2 + 1)*abs(z)^3 + (1 - cos(theta))*x*z^4)/z^3]
186+
[ sin(theta)*abs(z) + (1 - cos(theta))*x*sqrt(-z^2 - x^2 + 1) (cos(theta) - 1)*z^2 + (cos(theta) - 1)*x^2 + 1 -(sin(theta)*x*abs(z) + (cos(theta) - 1)*z^2*sqrt(-z^2 - x^2 + 1))/z]
187+
[ -(sin(theta)*sqrt(-z^2 - x^2 + 1)*abs(z) + (cos(theta) - 1)*x*z^2)/z -((cos(theta) - 1)*z^2*sqrt(-z^2 - x^2 + 1) - sin(theta)*x*abs(z))/z (1 - cos(theta))*z^2 + cos(theta)]
188+
189+
Re-expressing some entries in terms of y and resolving the absolute
190+
values introduced by eliminating y, we get the desired result.
191+
"""
192+
cdef double x,y,z, len_v
193+
x,y,z = v
194+
len_v = sqrt(x*x+y*y+z*z)
195+
# normalize for an easier formula
196+
x /= len_v
197+
y /= len_v
198+
z /= len_v
199+
cdef double cos_t = cos(theta), sin_t = sin(theta)
200+
201+
entries = [
202+
(1 - cos_t)*x*x + cos_t,
203+
sin_t*z - (cos_t - 1)*x*y,
204+
-sin_t*y + (1 - cos_t)*x*z,
205+
206+
-sin_t*z + (1 - cos_t)*x*y,
207+
(1 - cos_t)*y*y + cos_t,
208+
sin_t*x - (cos_t - 1)*z*y,
209+
210+
sin_t*y - (cos_t - 1)*x*z,
211+
-(cos_t - 1)*z*y - sin_t*x,
212+
(1 - cos_t)*z*z + cos_t ]
213+
214+
return matrix(RDF, 3, 3, entries)

0 commit comments

Comments
 (0)