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