1 // clang-format off
2 /*Copyright (c) 2016 PM Larsen
3
4 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
5
6 The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
7
8 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
9 */
10
11 #include "ptm_quat.h"
12 #include <cstring>
13 #include <cmath>
14 #include <algorithm>
15
16
17 namespace ptm {
18
19 #define SIGN(x) (x >= 0 ? 1 : -1)
20
21
22 const double generator_cubic[24][4] = {
23 { 1, 0, 0, 0 },
24 { sqrt(2)/2, sqrt(2)/2, 0, 0 },
25 { sqrt(2)/2, 0, sqrt(2)/2, 0 },
26 { sqrt(2)/2, 0, 0, sqrt(2)/2 },
27 { sqrt(2)/2, 0, 0, -sqrt(2)/2 },
28 { sqrt(2)/2, 0, -sqrt(2)/2, 0 },
29 { sqrt(2)/2, -sqrt(2)/2, -0, -0 },
30 { 0.5, 0.5, 0.5, 0.5 },
31 { 0.5, 0.5, 0.5, -0.5 },
32 { 0.5, 0.5, -0.5, 0.5 },
33 { 0.5, 0.5, -0.5, -0.5 },
34 { 0.5, -0.5, 0.5, 0.5 },
35 { 0.5, -0.5, 0.5, -0.5 },
36 { 0.5, -0.5, -0.5, 0.5 },
37 { 0.5, -0.5, -0.5, -0.5 },
38 { 0, 1, 0, 0 },
39 { 0, sqrt(2)/2, sqrt(2)/2, 0 },
40 { 0, sqrt(2)/2, 0, sqrt(2)/2 },
41 { 0, sqrt(2)/2, 0, -sqrt(2)/2 },
42 { 0, sqrt(2)/2, -sqrt(2)/2, 0 },
43 { 0, 0, 1, 0 },
44 { 0, 0, sqrt(2)/2, sqrt(2)/2 },
45 { 0, 0, sqrt(2)/2, -sqrt(2)/2 },
46 { 0, 0, 0, 1 },
47 };
48
49 const double generator_diamond_cubic[12][4] = {
50 { 1, 0, 0, 0 },
51 { 0.5, 0.5, 0.5, 0.5 },
52 { 0.5, 0.5, 0.5, -0.5 },
53 { 0.5, 0.5, -0.5, 0.5 },
54 { 0.5, 0.5, -0.5, -0.5 },
55 { 0.5, -0.5, 0.5, 0.5 },
56 { 0.5, -0.5, 0.5, -0.5 },
57 { 0.5, -0.5, -0.5, 0.5 },
58 { 0.5, -0.5, -0.5, -0.5 },
59 { 0, 1, 0, 0 },
60 { 0, 0, 1, 0 },
61 { 0, 0, 0, 1 },
62 };
63
64 const double generator_hcp[6][4] = {
65 { 1, 0, 0, 0 },
66 { 0.5, 0, 0, sqrt(3)/2 },
67 { 0.5, 0, 0, -sqrt(3)/2 },
68 { 0, sqrt(3)/2, 0.5, 0 },
69 { 0, sqrt(3)/2, -0.5, 0 },
70 { 0, 0, 1, 0 },
71 };
72
73
74 const double generator_hcp_conventional[12][4] = {
75 { 1, 0, 0, 0 },
76 { sqrt(3)/2, 0, 0, 0.5 },
77 { sqrt(3)/2, 0, 0, -0.5 },
78 { 0.5, 0, 0, sqrt(3)/2 },
79 { 0.5, 0, 0, -sqrt(3)/2 },
80 { 0, 1, 0, 0 },
81 { 0, sqrt(3)/2, 0.5, 0 },
82 { 0, sqrt(3)/2, -0.5, 0 },
83 { 0, 0.5, sqrt(3)/2, 0 },
84 { 0, 0.5, -sqrt(3)/2, 0 },
85 { 0, 0, 1, 0 },
86 { 0, 0, 0, 1 },
87 };
88
89 const double generator_diamond_hexagonal[3][4] = {
90 { 1, 0, 0, 0 },
91 { 0.5, 0, 0, sqrt(3)/2 },
92 { 0.5, 0, 0, -sqrt(3)/2 },
93 };
94
95 const double generator_icosahedral[60][4] = {
96 { 1, 0, 0, 0 },
97 { (1+sqrt(5))/4, 0.5, sqrt(25-10*sqrt(5))/10, sqrt(50-10*sqrt(5))/20 },
98 { (1+sqrt(5))/4, 0.5, -sqrt(25-10*sqrt(5))/10, -sqrt(50-10*sqrt(5))/20 },
99 { (1+sqrt(5))/4, 1/(1+sqrt(5)), sqrt(10*sqrt(5)+50)/20, -sqrt(50-10*sqrt(5))/20 },
100 { (1+sqrt(5))/4, 1/(1+sqrt(5)), -sqrt(10*sqrt(5)+50)/20, sqrt(50-10*sqrt(5))/20 },
101 { (1+sqrt(5))/4, 0, sqrt(50-10*sqrt(5))/10, sqrt(50-10*sqrt(5))/20 },
102 { (1+sqrt(5))/4, 0, 0, sqrt(5./8-sqrt(5)/8) },
103 { (1+sqrt(5))/4, 0, 0, -sqrt(5./8-sqrt(5)/8) },
104 { (1+sqrt(5))/4, 0, -sqrt(50-10*sqrt(5))/10, -sqrt(50-10*sqrt(5))/20 },
105 { (1+sqrt(5))/4, -1/(1+sqrt(5)), sqrt(10*sqrt(5)+50)/20, -sqrt(50-10*sqrt(5))/20 },
106 { (1+sqrt(5))/4, -1/(1+sqrt(5)), -sqrt(10*sqrt(5)+50)/20, sqrt(50-10*sqrt(5))/20 },
107 { (1+sqrt(5))/4, -0.5, sqrt(25-10*sqrt(5))/10, sqrt(50-10*sqrt(5))/20 },
108 { (1+sqrt(5))/4, -0.5, -sqrt(25-10*sqrt(5))/10, -sqrt(50-10*sqrt(5))/20 },
109 { 0.5, (1+sqrt(5))/4, sqrt(50-10*sqrt(5))/20, -sqrt(25-10*sqrt(5))/10 },
110 { 0.5, (1+sqrt(5))/4, -sqrt(50-10*sqrt(5))/20, sqrt(25-10*sqrt(5))/10 },
111 { 0.5, 0.5, sqrt((5+2*sqrt(5))/20), sqrt(25-10*sqrt(5))/10 },
112 { 0.5, 0.5, sqrt(25-10*sqrt(5))/10, -sqrt((5+2*sqrt(5))/20) },
113 { 0.5, 0.5, -sqrt(25-10*sqrt(5))/10, sqrt((5+2*sqrt(5))/20) },
114 { 0.5, 0.5, -sqrt((5+2*sqrt(5))/20), -sqrt(25-10*sqrt(5))/10 },
115 { 0.5, 1/(1+sqrt(5)), sqrt(10*sqrt(5)+50)/20, sqrt((5+2*sqrt(5))/20) },
116 { 0.5, 1/(1+sqrt(5)), -sqrt(10*sqrt(5)+50)/20, -sqrt((5+2*sqrt(5))/20) },
117 { 0.5, 0, sqrt((5+sqrt(5))/10), -sqrt(25-10*sqrt(5))/10 },
118 { 0.5, 0, sqrt(50-10*sqrt(5))/10, -sqrt((5+2*sqrt(5))/20) },
119 { 0.5, 0, -sqrt(50-10*sqrt(5))/10, sqrt((5+2*sqrt(5))/20) },
120 { 0.5, 0, -sqrt((5+sqrt(5))/10), sqrt(25-10*sqrt(5))/10 },
121 { 0.5, -1/(1+sqrt(5)), sqrt(10*sqrt(5)+50)/20, sqrt((5+2*sqrt(5))/20) },
122 { 0.5, -1/(1+sqrt(5)), -sqrt(10*sqrt(5)+50)/20, -sqrt((5+2*sqrt(5))/20) },
123 { 0.5, -0.5, sqrt((5+2*sqrt(5))/20), sqrt(25-10*sqrt(5))/10 },
124 { 0.5, -0.5, sqrt(25-10*sqrt(5))/10, -sqrt((5+2*sqrt(5))/20) },
125 { 0.5, -0.5, -sqrt(25-10*sqrt(5))/10, sqrt((5+2*sqrt(5))/20) },
126 { 0.5, -0.5, -sqrt((5+2*sqrt(5))/20), -sqrt(25-10*sqrt(5))/10 },
127 { 0.5, -(1+sqrt(5))/4, sqrt(50-10*sqrt(5))/20, -sqrt(25-10*sqrt(5))/10 },
128 { 0.5, -(1+sqrt(5))/4, -sqrt(50-10*sqrt(5))/20, sqrt(25-10*sqrt(5))/10 },
129 { 1/(1+sqrt(5)), (1+sqrt(5))/4, sqrt(50-10*sqrt(5))/20, sqrt(10*sqrt(5)+50)/20 },
130 { 1/(1+sqrt(5)), (1+sqrt(5))/4, -sqrt(50-10*sqrt(5))/20, -sqrt(10*sqrt(5)+50)/20 },
131 { 1/(1+sqrt(5)), 0.5, sqrt((5+2*sqrt(5))/20), -sqrt(10*sqrt(5)+50)/20 },
132 { 1/(1+sqrt(5)), 0.5, -sqrt((5+2*sqrt(5))/20), sqrt(10*sqrt(5)+50)/20 },
133 { 1/(1+sqrt(5)), 0, sqrt((5+sqrt(5))/10), sqrt(10*sqrt(5)+50)/20 },
134 { 1/(1+sqrt(5)), 0, 0, sqrt(1-1/(2*sqrt(5)+6)) },
135 { 1/(1+sqrt(5)), 0, 0, -sqrt(1-1/(2*sqrt(5)+6)) },
136 { 1/(1+sqrt(5)), 0, -sqrt((5+sqrt(5))/10), -sqrt(10*sqrt(5)+50)/20 },
137 { 1/(1+sqrt(5)), -0.5, sqrt((5+2*sqrt(5))/20), -sqrt(10*sqrt(5)+50)/20 },
138 { 1/(1+sqrt(5)), -0.5, -sqrt((5+2*sqrt(5))/20), sqrt(10*sqrt(5)+50)/20 },
139 { 1/(1+sqrt(5)), -(1+sqrt(5))/4, sqrt(50-10*sqrt(5))/20, sqrt(10*sqrt(5)+50)/20 },
140 { 1/(1+sqrt(5)), -(1+sqrt(5))/4, -sqrt(50-10*sqrt(5))/20, -sqrt(10*sqrt(5)+50)/20 },
141 { 0, 1, 0, 0 },
142 { 0, (1+sqrt(5))/4, sqrt(5./8-sqrt(5)/8), 0 },
143 { 0, (1+sqrt(5))/4, sqrt(50-10*sqrt(5))/20, -sqrt(50-10*sqrt(5))/10 },
144 { 0, (1+sqrt(5))/4, -sqrt(50-10*sqrt(5))/20, sqrt(50-10*sqrt(5))/10 },
145 { 0, (1+sqrt(5))/4, -sqrt(5./8-sqrt(5)/8), 0 },
146 { 0, 0.5, sqrt((5+2*sqrt(5))/20), sqrt(50-10*sqrt(5))/10 },
147 { 0, 0.5, sqrt(25-10*sqrt(5))/10, sqrt((5+sqrt(5))/10) },
148 { 0, 0.5, -sqrt(25-10*sqrt(5))/10, -sqrt((5+sqrt(5))/10) },
149 { 0, 0.5, -sqrt((5+2*sqrt(5))/20), -sqrt(50-10*sqrt(5))/10 },
150 { 0, 1/(1+sqrt(5)), sqrt(1-1/(2*sqrt(5)+6)), 0 },
151 { 0, 1/(1+sqrt(5)), sqrt(10*sqrt(5)+50)/20, -sqrt((5+sqrt(5))/10) },
152 { 0, 1/(1+sqrt(5)), -sqrt(10*sqrt(5)+50)/20, sqrt((5+sqrt(5))/10) },
153 { 0, 1/(1+sqrt(5)), -sqrt(1-1/(2*sqrt(5)+6)), 0 },
154 { 0, 0, sqrt((5+sqrt(5))/10), -sqrt(50-10*sqrt(5))/10 },
155 { 0, 0, sqrt(50-10*sqrt(5))/10, sqrt((5+sqrt(5))/10) },
156 };
157
158
quat_rot(double * r,double * a,double * b)159 static void quat_rot(double* r, double* a, double* b)
160 {
161 b[0] = (r[0] * a[0] - r[1] * a[1] - r[2] * a[2] - r[3] * a[3]);
162 b[1] = (r[0] * a[1] + r[1] * a[0] + r[2] * a[3] - r[3] * a[2]);
163 b[2] = (r[0] * a[2] - r[1] * a[3] + r[2] * a[0] + r[3] * a[1]);
164 b[3] = (r[0] * a[3] + r[1] * a[2] - r[2] * a[1] + r[3] * a[0]);
165 }
166
rotate_quaternion_into_fundamental_zone(int num_generators,const double (* generator)[4],double * q)167 static int rotate_quaternion_into_fundamental_zone(int num_generators, const double (*generator)[4], double* q)
168 {
169 double max = 0.0;
170 int i = 0, bi = -1;
171 for (i=0;i<num_generators;i++)
172 {
173 const double* g = generator[i];
174 double t = fabs(q[0] * g[0] - q[1] * g[1] - q[2] * g[2] - q[3] * g[3]);
175 if (t > max)
176 {
177 max = t;
178 bi = i;
179 }
180 }
181
182 double f[4];
183 quat_rot(q, (double*)generator[bi], f);
184 memcpy(q, &f, 4 * sizeof(double));
185 if (q[0] < 0)
186 {
187 q[0] = -q[0];
188 q[1] = -q[1];
189 q[2] = -q[2];
190 q[3] = -q[3];
191 }
192
193 return bi;
194 }
195
rotate_quaternion_into_cubic_fundamental_zone(double * q)196 int rotate_quaternion_into_cubic_fundamental_zone(double* q)
197 {
198 return rotate_quaternion_into_fundamental_zone(24, generator_cubic, q);
199 }
200
rotate_quaternion_into_diamond_cubic_fundamental_zone(double * q)201 int rotate_quaternion_into_diamond_cubic_fundamental_zone(double* q)
202 {
203 return rotate_quaternion_into_fundamental_zone(12, generator_diamond_cubic, q);
204 }
205
rotate_quaternion_into_icosahedral_fundamental_zone(double * q)206 int rotate_quaternion_into_icosahedral_fundamental_zone(double* q)
207 {
208 return rotate_quaternion_into_fundamental_zone(60, generator_icosahedral, q);
209 }
210
rotate_quaternion_into_hcp_fundamental_zone(double * q)211 int rotate_quaternion_into_hcp_fundamental_zone(double* q)
212 {
213 return rotate_quaternion_into_fundamental_zone(6, generator_hcp, q);
214 }
215
rotate_quaternion_into_hcp_conventional_fundamental_zone(double * q)216 int rotate_quaternion_into_hcp_conventional_fundamental_zone(double* q)
217 {
218 return rotate_quaternion_into_fundamental_zone(12, generator_hcp_conventional, q);
219 }
220
rotate_quaternion_into_diamond_hexagonal_fundamental_zone(double * q)221 int rotate_quaternion_into_diamond_hexagonal_fundamental_zone(double* q)
222 {
223 return rotate_quaternion_into_fundamental_zone(3, generator_diamond_hexagonal, q);
224 }
225
quat_dot(double * a,double * b)226 double quat_dot(double* a, double* b)
227 {
228 return a[0] * b[0]
229 + a[1] * b[1]
230 + a[2] * b[2]
231 + a[3] * b[3];
232 }
233
quat_size(double * q)234 double quat_size(double* q)
235 {
236 return sqrt(quat_dot(q, q));
237 }
238
normalize_quaternion(double * q)239 void normalize_quaternion(double* q)
240 {
241 double size = quat_size(q);
242
243 q[0] /= size;
244 q[1] /= size;
245 q[2] /= size;
246 q[3] /= size;
247 }
248
quaternion_to_rotation_matrix(double * q,double * u)249 void quaternion_to_rotation_matrix(double* q, double* u)
250 {
251 double a = q[0];
252 double b = q[1];
253 double c = q[2];
254 double d = q[3];
255
256 u[0] = a*a + b*b - c*c - d*d;
257 u[1] = 2*b*c - 2*a*d;
258 u[2] = 2*b*d + 2*a*c;
259
260 u[3] = 2*b*c + 2*a*d;
261 u[4] = a*a - b*b + c*c - d*d;
262 u[5] = 2*c*d - 2*a*b;
263
264 u[6] = 2*b*d - 2*a*c;
265 u[7] = 2*c*d + 2*a*b;
266 u[8] = a*a - b*b - c*c + d*d;
267 }
268
rotation_matrix_to_quaternion(double * u,double * q)269 void rotation_matrix_to_quaternion(double* u, double* q)
270 {
271 double r11 = u[0];
272 double r12 = u[1];
273 double r13 = u[2];
274 double r21 = u[3];
275 double r22 = u[4];
276 double r23 = u[5];
277 double r31 = u[6];
278 double r32 = u[7];
279 double r33 = u[8];
280
281 q[0] = (1.0 + r11 + r22 + r33) / 4.0;
282 q[1] = (1.0 + r11 - r22 - r33) / 4.0;
283 q[2] = (1.0 - r11 + r22 - r33) / 4.0;
284 q[3] = (1.0 - r11 - r22 + r33) / 4.0;
285
286 q[0] = sqrt(std::max(0., q[0]));
287 q[1] = sqrt(std::max(0., q[1]));
288 q[2] = sqrt(std::max(0., q[2]));
289 q[3] = sqrt(std::max(0., q[3]));
290
291 double m0 = std::max(q[0], q[1]);
292 double m1 = std::max(q[2], q[3]);
293 double max = std::max(m0, m1);
294
295 int i = 0;
296 for (i=0;i<4;i++)
297 if (q[i] == max)
298 break;
299
300 if (i == 0)
301 {
302 q[1] *= SIGN(r32 - r23);
303 q[2] *= SIGN(r13 - r31);
304 q[3] *= SIGN(r21 - r12);
305 }
306 else if (i == 1)
307 {
308 q[0] *= SIGN(r32 - r23);
309 q[2] *= SIGN(r21 + r12);
310 q[3] *= SIGN(r13 + r31);
311 }
312 else if (i == 2)
313 {
314 q[0] *= SIGN(r13 - r31);
315 q[1] *= SIGN(r21 + r12);
316 q[3] *= SIGN(r32 + r23);
317 }
318 else if (i == 3)
319 {
320 q[0] *= SIGN(r21 - r12);
321 q[1] *= SIGN(r31 + r13);
322 q[2] *= SIGN(r32 + r23);
323 }
324
325 normalize_quaternion(q);
326 }
327
quat_quick_misorientation(double * q1,double * q2)328 double quat_quick_misorientation(double* q1, double* q2)
329 {
330 double t = quat_dot(q1, q2);
331 t = std::min(1., std::max(-1., t));
332 return 2 * t * t - 1;
333 }
334
quat_misorientation(double * q1,double * q2)335 double quat_misorientation(double* q1, double* q2)
336 {
337 return acos(quat_quick_misorientation(q1, q2));
338 }
339
340 }
341
342