|
| 1 | +//! Example of rotation with Euler angles. |
| 2 | +//! |
| 3 | +//! This is an example of some coordinate transformations (using Euler angles) |
| 4 | +//! for illustrative purposes. Note that other crates such as |
| 5 | +//! [`cgmath`](https://crates.io/crates/cgmath) or |
| 6 | +//! [`nalgebra`](https://crates.io/crates/nalgebra) may be better-suited if |
| 7 | +//! most of your work is coordinate transformations, since they have built-in |
| 8 | +//! geometry primitives. |
| 9 | +//! |
| 10 | +//! This is the original Python program: |
| 11 | +//! |
| 12 | +//! ```python |
| 13 | +//! # Euler angles (rows) for four coordinate systems (columns). |
| 14 | +//! nelems = 4 |
| 15 | +//! bunge = np.ones((3, nelems)) |
| 16 | +//! |
| 17 | +//! # Precompute sines and cosines |
| 18 | +//! s1 = np.sin(bunge[0, :]) |
| 19 | +//! c1 = np.cos(bunge[0, :]) |
| 20 | +//! s2 = np.sin(bunge[1, :]) |
| 21 | +//! c2 = np.cos(bunge[1, :]) |
| 22 | +//! s3 = np.sin(bunge[2, :]) |
| 23 | +//! c3 = np.cos(bunge[2, :]) |
| 24 | +//! |
| 25 | +//! # Rotation matrices. |
| 26 | +//! rmat = np.zeros((3, 3, nelems), order='F') |
| 27 | +//! for i in range(nelems): |
| 28 | +//! rmat[0, 0, i] = c1[i] * c3[i] - s1[i] * s3[i] * c2[i] |
| 29 | +//! rmat[0, 1, i] = -c1[i] * s3[i] - s1[i] * c2[i] * c3[i] |
| 30 | +//! rmat[0, 2, i] = s1[i] * s2[i] |
| 31 | +//! |
| 32 | +//! rmat[1, 0, i] = s1[i] * c3[i] + c1[i] * c2[i] * s3[i] |
| 33 | +//! rmat[1, 1, i] = -s1[i] * s3[i] + c1[i] * c2[i] * c3[i] |
| 34 | +//! rmat[1, 2, i] = -c1[i] * s2[i] |
| 35 | +//! |
| 36 | +//! rmat[2, 0, i] = s2[i] * s3[i] |
| 37 | +//! rmat[2, 1, i] = s2[i] * c3[i] |
| 38 | +//! rmat[2, 2, i] = c2[i] |
| 39 | +//! |
| 40 | +//! # Unit vectors of coordinate systems to rotate. |
| 41 | +//! eye2d = np.eye(3) |
| 42 | +//! |
| 43 | +//! # Unit vectors after rotation. |
| 44 | +//! rotated = np.zeros((3, 3, nelems), order='F') |
| 45 | +//! for i in range(nelems): |
| 46 | +//! rotated[:,:,i] = rmat[:,:,i].dot(eye2d) |
| 47 | +//! ``` |
| 48 | +//! |
| 49 | +//! This is a direct translation to `ndarray`: |
| 50 | +//! |
| 51 | +//! ``` |
| 52 | +//! #[macro_use] |
| 53 | +//! extern crate ndarray; |
| 54 | +//! |
| 55 | +//! use ndarray::prelude::*; |
| 56 | +//! |
| 57 | +//! fn main() { |
| 58 | +//! let nelems = 4; |
| 59 | +//! let bunge = Array::ones((3, nelems)); |
| 60 | +//! |
| 61 | +//! let s1 = bunge.slice(s![0, ..]).mapv(f64::sin); |
| 62 | +//! let c1 = bunge.slice(s![0, ..]).mapv(f64::cos); |
| 63 | +//! let s2 = bunge.slice(s![1, ..]).mapv(f64::sin); |
| 64 | +//! let c2 = bunge.slice(s![1, ..]).mapv(f64::cos); |
| 65 | +//! let s3 = bunge.slice(s![2, ..]).mapv(f64::sin); |
| 66 | +//! let c3 = bunge.slice(s![2, ..]).mapv(f64::cos); |
| 67 | +//! |
| 68 | +//! let mut rmat = Array::zeros((3, 3, nelems).f()); |
| 69 | +//! for i in 0..nelems { |
| 70 | +//! rmat[[0, 0, i]] = c1[i] * c3[i] - s1[i] * s3[i] * c2[i]; |
| 71 | +//! rmat[[0, 1, i]] = -c1[i] * s3[i] - s1[i] * c2[i] * c3[i]; |
| 72 | +//! rmat[[0, 2, i]] = s1[i] * s2[i]; |
| 73 | +//! |
| 74 | +//! rmat[[1, 0, i]] = s1[i] * c3[i] + c1[i] * c2[i] * s3[i]; |
| 75 | +//! rmat[[1, 1, i]] = -s1[i] * s3[i] + c1[i] * c2[i] * c3[i]; |
| 76 | +//! rmat[[1, 2, i]] = -c1[i] * s2[i]; |
| 77 | +//! |
| 78 | +//! rmat[[2, 0, i]] = s2[i] * s3[i]; |
| 79 | +//! rmat[[2, 1, i]] = s2[i] * c3[i]; |
| 80 | +//! rmat[[2, 2, i]] = c2[i]; |
| 81 | +//! } |
| 82 | +//! |
| 83 | +//! let eye2d = Array::eye(3); |
| 84 | +//! |
| 85 | +//! let mut rotated = Array::zeros((3, 3, nelems).f()); |
| 86 | +//! for i in 0..nelems { |
| 87 | +//! rotated |
| 88 | +//! .slice_mut(s![.., .., i]) |
| 89 | +//! .assign({ &rmat.slice(s![.., .., i]).dot(&eye2d) }); |
| 90 | +//! } |
| 91 | +//! } |
| 92 | +//! ``` |
| 93 | +//! |
| 94 | +//! Instead of looping over indices, a cleaner (and usually faster) option is |
| 95 | +//! to zip arrays together. It's also possible to avoid some of the temporary |
| 96 | +//! memory allocations in the original program. The improved version looks like |
| 97 | +//! this: |
| 98 | +//! |
| 99 | +//! ``` |
| 100 | +//! #[macro_use] |
| 101 | +//! extern crate ndarray; |
| 102 | +//! |
| 103 | +//! use ndarray::prelude::*; |
| 104 | +//! |
| 105 | +//! fn main() { |
| 106 | +//! let nelems = 4; |
| 107 | +//! let bunge = Array2::<f64>::ones((3, nelems)); |
| 108 | +//! |
| 109 | +//! let mut rmat = Array::zeros((3, 3, nelems).f()); |
| 110 | +//! azip!(mut rmat (rmat.axis_iter_mut(Axis(2))), ref bunge (bunge.axis_iter(Axis(1))) in { |
| 111 | +//! let s1 = bunge[0].sin(); |
| 112 | +//! let c1 = bunge[0].cos(); |
| 113 | +//! let s2 = bunge[1].sin(); |
| 114 | +//! let c2 = bunge[1].cos(); |
| 115 | +//! let s3 = bunge[2].sin(); |
| 116 | +//! let c3 = bunge[2].cos(); |
| 117 | +//! |
| 118 | +//! rmat[[0, 0]] = c1 * c3 - s1 * s3 * c2; |
| 119 | +//! rmat[[0, 1]] = -c1 * s3 - s1 * c2 * c3; |
| 120 | +//! rmat[[0, 2]] = s1 * s2; |
| 121 | +//! |
| 122 | +//! rmat[[1, 0]] = s1 * c3 + c1 * c2 * s3; |
| 123 | +//! rmat[[1, 1]] = -s1 * s3 + c1 * c2 * c3; |
| 124 | +//! rmat[[1, 2]] = -c1 * s2; |
| 125 | +//! |
| 126 | +//! rmat[[2, 0]] = s2 * s3; |
| 127 | +//! rmat[[2, 1]] = s2 * c3; |
| 128 | +//! rmat[[2, 2]] = c2; |
| 129 | +//! }); |
| 130 | +//! |
| 131 | +//! let eye2d = Array2::<f64>::eye(3); |
| 132 | +//! |
| 133 | +//! let mut rotated = Array3::<f64>::zeros((3, 3, nelems).f()); |
| 134 | +//! azip!(mut rotated (rotated.axis_iter_mut(Axis(2)), rmat (rmat.axis_iter(Axis(2)))) in { |
| 135 | +//! rotated.assign({ &rmat.dot(&eye2d) }); |
| 136 | +//! }); |
| 137 | +//! } |
| 138 | +//! ``` |
0 commit comments