1 /**
2  * Represents a 3D rotation of a given radians about an axis represented as an
3  * unit 3D vector.
4  * <p>
5  * This class uses double-precision components.
6  * 
7  * @author Kai Burjack
8  */
9 module doml.axis_angle_4d;
10 
11 import Math = doml.math;
12 
13 import doml.matrix_3d;
14 import doml.matrix_4d;
15 
16 
17 import doml.vector_3d;
18 import doml.vector_4d;
19 
20 import doml.quaternion_d;
21 
22 /*
23  * The MIT License
24  *
25  * Copyright (c) 2015-2021 Kai Burjack
26  %$%@^ Translated by jordan4ibanez
27  *
28  * Permission is hereby granted, free of charge, to any person obtaining a copy
29  * of this software and associated documentation files (the "Software"), to deal
30  * in the Software without restriction, including without limitation the rights
31  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
32  * copies of the Software, and to permit persons to whom the Software is
33  * furnished to do so, subject to the following conditions:
34  *
35  * The above copyright notice and this permission notice shall be included in
36  * all copies or substantial portions of the Software.
37  *
38  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
39  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
40  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
41  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
42  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
43  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
44  * THE SOFTWARE.
45  */
46 
47 /**
48  * Represents a 3D rotation of a given radians about an axis represented as an
49  * unit 3D vector.
50  * <p>
51  * This class uses double-precision components.
52  * 
53  * @author Kai Burjack
54  */
55 struct AxisAngle4d {
56     /**
57      * The angle in radians.
58      */
59     double angle = 0.0;
60     /**
61      * The x-component of the rotation axis.
62      */
63     double x = 0.0;
64     /**
65      * The y-component of the rotation axis.
66      */
67     double y = 0.0;
68     /**
69      * The z-component of the rotation axis.
70      */
71     double z = 1.0;
72     /**
73      * Create a new {@link AxisAngle4d} with the same values of <code>a</code>.
74      * 
75      * @param a
76      *            the AngleAxis4d to copy the values from
77      */
78     this(AxisAngle4d a) {
79         x = a.x;
80         y = a.y;
81         z = a.z;
82         angle = (a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI);
83     }
84 
85     /**
86      * Create a new {@link AxisAngle4d} from the given {@link Quaterniond}.
87      * <p>
88      * Reference: <a href=
89      * "http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToAngle/"
90      * >http://www.euclideanspace.com</a>
91      * 
92      * @param q
93      *            the quaternion from which to create the new AngleAxis4d
94      */
95     this(Quaterniond q) {
96         double acos = Math.safeAcos(q.w);
97         double invSqrt = Math.invsqrt(1.0 - q.w * q.w);
98         if (Math.isInfinite(invSqrt)) {
99             this.x = 0.0;
100             this.y = 0.0;
101             this.z = 1.0;
102         } else {
103             this.x = q.x * invSqrt;
104             this.y = q.y * invSqrt;
105             this.z = q.z * invSqrt;
106         }
107         this.angle = acos + acos;
108     }
109 
110     /**
111      * Create a new {@link AxisAngle4d} with the given values.
112      *
113      * @param angle
114      *            the angle in radians
115      * @param x
116      *            the x-coordinate of the rotation axis
117      * @param y
118      *            the y-coordinate of the rotation axis
119      * @param z
120      *            the z-coordinate of the rotation axis
121      */
122     this(double angle, double x, double y, double z) {
123         this.x = x;
124         this.y = y;
125         this.z = z;
126         this.angle = (angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI);
127     }
128 
129     /**
130      * Create a new {@link AxisAngle4d} with the given values.
131      *
132      * @param angle the angle in radians
133      * @param v     the rotation axis as a {@link Vector3d}
134      */
135     this(double angle, Vector3d v) {
136         this(angle, v.x, v.y, v.z);
137     }
138 
139 
140     /**
141      * Set this {@link AxisAngle4d} to the given values.
142      * 
143      * @param angle
144      *            the angle in radians
145      * @param x
146      *            the x-coordinate of the rotation axis
147      * @param y
148      *            the y-coordinate of the rotation axis
149      * @param z
150      *            the z-coordinate of the rotation axis
151      * @return this
152      */
153     ref public AxisAngle4d set(double angle, double x, double y, double z) return {
154         this.x = x;
155         this.y = y;
156         this.z = z;
157         this.angle = (angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI);
158         return this;
159     }
160 
161     /**
162      * Set this {@link AxisAngle4d} to the given values.
163      *
164      * @param angle
165      *            the angle in radians
166      * @param v    
167      *            the rotation axis as a {@link Vector3d}
168      * @return this
169      */
170     ref public AxisAngle4d set(double angle, Vector3d v) return {
171         return set(angle, v.x, v.y, v.z);
172     }
173 
174     /**
175      * Set this {@link AxisAngle4d} to be equivalent to the given
176      * {@link Quaterniond}.
177      * 
178      * @param q
179      *            the quaternion to set this AngleAxis4d from
180      * @return this
181      */
182     ref public AxisAngle4d set(Quaterniond q) return {
183         double acos = Math.safeAcos(q.w);
184         double invSqrt = Math.invsqrt(1.0f - q.w * q.w);
185         if (Math.isInfinite(invSqrt)) {
186             this.x = 0.0;
187             this.y = 0.0;
188             this.z = 1.0;
189         } else {
190             this.x = q.x * invSqrt;
191             this.y = q.y * invSqrt;
192             this.z = q.z * invSqrt;
193         }
194         this.angle = acos + acos;
195         return this;
196     }
197 
198     /**
199      * Set this {@link AxisAngle4d} to be equivalent to the rotation 
200      * of the given {@link Matrix3d}.
201      * <p>
202      * Reference: <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToAngle/">http://www.euclideanspace.com</a>
203      * 
204      * @param m
205      *            the Matrix3d to set this AngleAxis4d from
206      * @return this
207      */
208     ref public AxisAngle4d set(Matrix3d m) return {
209         double nm00 = m.m00, nm01 = m.m01, nm02 = m.m02;
210         double nm10 = m.m10, nm11 = m.m11, nm12 = m.m12;
211         double nm20 = m.m20, nm21 = m.m21, nm22 = m.m22;
212         double lenX = Math.invsqrt(m.m00 * m.m00 + m.m01 * m.m01 + m.m02 * m.m02);
213         double lenY = Math.invsqrt(m.m10 * m.m10 + m.m11 * m.m11 + m.m12 * m.m12);
214         double lenZ = Math.invsqrt(m.m20 * m.m20 + m.m21 * m.m21 + m.m22 * m.m22);
215         nm00 *= lenX; nm01 *= lenX; nm02 *= lenX;
216         nm10 *= lenY; nm11 *= lenY; nm12 *= lenY;
217         nm20 *= lenZ; nm21 *= lenZ; nm22 *= lenZ;
218         double epsilon = 1E-4, epsilon2 = 1E-3;
219         if (Math.abs(nm10 - nm01) < epsilon && Math.abs(nm20 - nm02) < epsilon && Math.abs(nm21 - nm12) < epsilon) {
220             if (Math.abs(nm10 + nm01) < epsilon2 && Math.abs(nm20 + nm02) < epsilon2 && Math.abs(nm21 + nm12) < epsilon2
221                     && Math.abs(nm00 + nm11 + nm22 - 3) < epsilon2) {
222                 x = 0;
223                 y = 0;
224                 z = 1;
225                 angle = 0;
226                 return this;
227             }
228             angle = Math.PI;
229             double xx = (nm00 + 1) / 2;
230             double yy = (nm11 + 1) / 2;
231             double zz = (nm22 + 1) / 2;
232             double xy = (nm10 + nm01) / 4;
233             double xz = (nm20 + nm02) / 4;
234             double yz = (nm21 + nm12) / 4;
235             if ((xx > yy) && (xx > zz)) {
236                 x = Math.sqrt(xx);
237                 y = xy / x;
238                 z = xz / x;
239             } else if (yy > zz) {
240                 y = Math.sqrt(yy);
241                 x = xy / y;
242                 z = yz / y;
243             } else {
244                 z = Math.sqrt(zz);
245                 x = xz / z;
246                 y = yz / z;
247             }
248             return this;
249         }
250         double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
251         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
252         x = (nm12 - nm21) / s;
253         y = (nm20 - nm02) / s;
254         z = (nm01 - nm10) / s;
255         return this;
256     }
257 
258     /**
259      * Set this {@link AxisAngle4d} to be equivalent to the rotational component 
260      * of the given {@link Matrix4d}.
261      * <p>
262      * Reference: <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToAngle/">http://www.euclideanspace.com</a>
263      * 
264      * @param m
265      *            the Matrix4d to set this AngleAxis4d from
266      * @return this
267      */
268     ref public AxisAngle4d set(Matrix4d m) return {
269         double nm00 = m.m00, nm01 = m.m01, nm02 = m.m02;
270         double nm10 = m.m10, nm11 = m.m11, nm12 = m.m12;
271         double nm20 = m.m20, nm21 = m.m21, nm22 = m.m22;
272         double lenX = Math.invsqrt(m.m00 * m.m00 + m.m01 * m.m01 + m.m02 * m.m02);
273         double lenY = Math.invsqrt(m.m10 * m.m10 + m.m11 * m.m11 + m.m12 * m.m12);
274         double lenZ = Math.invsqrt(m.m20 * m.m20 + m.m21 * m.m21 + m.m22 * m.m22);
275         nm00 *= lenX; nm01 *= lenX; nm02 *= lenX;
276         nm10 *= lenY; nm11 *= lenY; nm12 *= lenY;
277         nm20 *= lenZ; nm21 *= lenZ; nm22 *= lenZ;
278         double epsilon = 1E-4, epsilon2 = 1E-3;
279         if (Math.abs(nm10 - nm01) < epsilon && Math.abs(nm20 - nm02) < epsilon && Math.abs(nm21 - nm12) < epsilon) {
280             if (Math.abs(nm10 + nm01) < epsilon2 && Math.abs(nm20 + nm02) < epsilon2 && Math.abs(nm21 + nm12) < epsilon2
281                     && Math.abs(nm00 + nm11 + nm22 - 3) < epsilon2) {
282                 x = 0;
283                 y = 0;
284                 z = 1;
285                 angle = 0;
286                 return this;
287             }
288             angle = Math.PI;
289             double xx = (nm00 + 1) / 2;
290             double yy = (nm11 + 1) / 2;
291             double zz = (nm22 + 1) / 2;
292             double xy = (nm10 + nm01) / 4;
293             double xz = (nm20 + nm02) / 4;
294             double yz = (nm21 + nm12) / 4;
295             if ((xx > yy) && (xx > zz)) {
296                 x = Math.sqrt(xx);
297                 y = xy / x;
298                 z = xz / x;
299             } else if (yy > zz) {
300                 y = Math.sqrt(yy);
301                 x = xy / y;
302                 z = yz / y;
303             } else {
304                 z = Math.sqrt(zz);
305                 x = xz / z;
306                 y = yz / z;
307             }
308             return this;
309         }
310         double s = Math.sqrt((nm12 - nm21) * (nm12 - nm21) + (nm20 - nm02) * (nm20 - nm02) + (nm01 - nm10) * (nm01 - nm10));
311         angle = Math.safeAcos((nm00 + nm11 + nm22 - 1) / 2);
312         x = (nm12 - nm21) / s;
313         y = (nm20 - nm02) / s;
314         z = (nm01 - nm10) / s;
315         return this;
316     }
317 
318     /**
319      * Set this {@link AxisAngle4d} to the values of <code>a</code>.
320      * 
321      * @param a
322      *            the AngleAxis4f to copy the values from
323      * @return this
324      */
325     ref public AxisAngle4d set(AxisAngle4d a) return {
326         x = a.x;
327         y = a.y;
328         z = a.z;
329         angle = (a.angle < 0.0 ? Math.PI + Math.PI + a.angle % (Math.PI + Math.PI) : a.angle) % (Math.PI + Math.PI);
330         return this;
331     }
332 
333     /**
334      * Set the given {@link Quaterniond} to be equivalent to this {@link AxisAngle4d} rotation.
335      * 
336      * @see Quaterniond#set(AxisAngle4d)
337      * 
338      * @param q
339      *          the quaternion to set
340      * @return q
341      */
342     public Quaterniond get(Quaterniond q) {
343         return q.set(this);
344     }
345 
346     /**
347      * Set the given {@link Matrix4d} to a rotation transformation equivalent to this {@link AxisAngle4d}.
348      * 
349      * @see Matrix4f#set(AxisAngle4d)
350      * 
351      * @param m
352      *          the matrix to set
353      * @return m
354      */
355     public Matrix4d get(Matrix4d m) {
356         return m.set(this);
357     }
358 
359     /**
360      * Set the given {@link Matrix3d} to a rotation transformation equivalent to this {@link AxisAngle4d}.
361      * 
362      * @see Matrix3f#set(AxisAngle4d)
363      * 
364      * @param m
365      *          the matrix to set
366      * @return m
367      */
368     public Matrix3d get(Matrix3d m) {
369         return m.set(this);
370     }
371 
372     /**
373      * Set the given {@link AxisAngle4d} to this {@link AxisAngle4d}.
374      * 
375      * @param dest
376      *          will hold the result
377      * @return dest
378      */
379     public AxisAngle4d get(AxisAngle4d dest) {
380         return dest.set(this);
381     }
382 
383     /**
384      * Normalize the axis vector.
385      * 
386      * @return this
387      */
388     ref public AxisAngle4d normalize() return {
389         double invLength = Math.invsqrt(x * x + y * y + z * z);
390         x *= invLength;
391         y *= invLength;
392         z *= invLength;
393         return this;
394     }
395 
396     /**
397      * Increase the rotation angle by the given amount.
398      * <p>
399      * This method also takes care of wrapping around.
400      * 
401      * @param ang
402      *          the angle increase
403      * @return this
404      */
405     ref public AxisAngle4d rotate(double ang) return {
406         angle += ang;
407         angle = (angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI);
408         return this;
409     }
410 
411     /**
412      * Transform the given vector by the rotation transformation described by this {@link AxisAngle4d}.
413      * 
414      * @param v
415      *          the vector to transform
416      * @return v
417      */
418     public Vector3d transform(Vector3d v) {
419         return transform(v, v);
420     }
421 
422     /**
423      * Transform the given vector by the rotation transformation described by this {@link AxisAngle4d}
424      * and store the result in <code>dest</code>.
425      * 
426      * @param v
427      *          the vector to transform
428      * @param dest
429      *          will hold the result
430      * @return dest
431      */
432     public Vector3d transform(Vector3d v, Vector3d dest) {
433         double sin = Math.sin(angle);
434         double cos = Math.cosFromSin(sin, angle);
435         double dot = x * v.x + y * v.y + z * v.z;
436         dest.set(v.x * cos + sin * (y * v.z - z * v.y) + (1.0 - cos) * dot * x,
437                  v.y * cos + sin * (z * v.x - x * v.z) + (1.0 - cos) * dot * y,
438                  v.z * cos + sin * (x * v.y - y * v.x) + (1.0 - cos) * dot * z);
439         return dest;
440     }
441 
442     /**
443      * Transform the given vector by the rotation transformation described by this {@link AxisAngle4d}.
444      * 
445      * @param v
446      *          the vector to transform
447      * @return v
448      */
449     public Vector4d transform(Vector4d v) {
450         return transform(v, v);
451     }
452 
453     /**
454      * Transform the given vector by the rotation transformation described by this {@link AxisAngle4d}
455      * and store the result in <code>dest</code>.
456      * 
457      * @param v
458      *          the vector to transform
459      * @param dest
460      *          will hold the result
461      * @return dest
462      */
463     public Vector4d transform(Vector4d v, Vector4d dest) {
464         double sin = Math.sin(angle);
465         double cos = Math.cosFromSin(sin, angle);
466         double dot = x * v.x + y * v.y + z * v.z;
467         dest.set(v.x * cos + sin * (y * v.z - z * v.y) + (1.0 - cos) * dot * x,
468                  v.y * cos + sin * (z * v.x - x * v.z) + (1.0 - cos) * dot * y,
469                  v.z * cos + sin * (x * v.y - y * v.x) + (1.0 - cos) * dot * z,
470                  dest.w);
471         return dest;
472     }
473 
474     public int hashCode() {
475         immutable int prime = 31;
476         int result = 1;
477         long temp;
478         temp = Math.doubleToLongBits((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI));
479         result = prime * result + cast(int) (temp ^ (temp >>> 32));
480         temp = Math.doubleToLongBits(x);
481         result = prime * result + cast(int) (temp ^ (temp >>> 32));
482         temp = Math.doubleToLongBits(y);
483         result = prime * result + cast(int) (temp ^ (temp >>> 32));
484         temp = Math.doubleToLongBits(z);
485         result = prime * result + cast(int) (temp ^ (temp >>> 32));
486         return result;
487     }
488 
489     public bool equals(AxisAngle4d other) {
490         if (this == other)
491             return true;
492         if (Math.doubleToLongBits((angle < 0.0 ? Math.PI + Math.PI + angle % (Math.PI + Math.PI) : angle) % (Math.PI + Math.PI)) != 
493                 Math.doubleToLongBits((other.angle < 0.0 ? Math.PI + Math.PI + other.angle % (Math.PI + Math.PI) : other.angle) % (Math.PI + Math.PI)))
494             return false;
495         if (Math.doubleToLongBits(x) != Math.doubleToLongBits(other.x))
496             return false;
497         if (Math.doubleToLongBits(y) != Math.doubleToLongBits(other.y))
498             return false;
499         if (Math.doubleToLongBits(z) != Math.doubleToLongBits(other.z))
500             return false;
501         return true;
502     }
503 }