Coverage for src/pygnd/rotations.py: 100%

1####################################################################################################

6# Redistribution and use in source and binary forms, with or without modification, are

7# permitted provided that the following conditions are met:

9# - Redistributions of source code must retain the above copyright notice, this list

10# of conditions and the following disclaimer.

11# - Redistributions in binary form must reproduce the above copyright notice, this

12# list of conditions and the following disclaimer in the documentation and/or

13# other materials provided with the distribution.

14# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names

15# of its contributors may be used to endorse or promote products derived from

16# this software without specific prior written permission.

17#

18# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

19# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

20# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

21# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE

22# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

23# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

24# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

25# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,

26# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE

27# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

28####################################################################################################

29"""Conversions between rotation/orientation representations used for EBSD crystallography data.

31Every function follows the naming convention `x2y`, converting representation

32`x` to representation `y`. The representations, and the array shape each uses,

33are:

35- `qu` - quaternion, shape `(..., 4)`, in the form `(w, x, y, z)`.

36- `om` - rotation/orientation matrix, shape `(..., 3, 3)`.

37- `eu` - Bunge-Euler angles, shape `(..., 3)`, in the form `(phi1, Phi, phi2)`, in radians.

38- `ax` - axis-angle pair, shape `(..., 4)`: a unit rotation axis followed by the

39 rotation angle in radians, in the form `(n1, n2, n3, omega)`.

40- `ro` - Rodrigues-Frank vector, shape `(..., 4)`: a unit rotation axis followed

41 by `tan(omega / 2)` (which is `inf` for a 180 degree rotation), in the form

42 `(n1, n2, n3, tan(omega / 2))`.

43- `ho` - homochoric vector, shape `(..., 3)`.

44- `cu` - cubochoric vector, shape `(..., 3)`.

46"""

48import numpy as np

51_EPSIJK = 1 # @private - Convention for the sign of the Levi-Civita symbol used throughout these conversions

53_SC = np.pi ** (1.0 / 6.0) / 6.0 ** (1.0 / 6.0) # @private - Parameter for conversion to/from cubochoric space

54_BETA = np.pi ** (5.0 / 6.0) / 6.0 ** (1.0 / 6.0) / 2.0 # @private - Parameter for conversion to/from cubochoric space

55_R1 = (3.0 * np.pi / 4.0) ** (1.0 / 3.0) # @private - Parameter for conversion to/from cubochoric space

58def qu2om(qu):

59 """Quaternion to rotation matrix.

61 Args:

62 qu: shape (..., 4) quaternions in the form (w, x, y, z).

64 Returns:

65 np.ndarray of shape (..., 3, 3): rotation matrices.

66 """

67 qq = qu[..., 0:1] ** 2 - (qu[..., 1:2] ** 2 + qu[..., 2:3] ** 2 + qu[..., 3:4] ** 2)

68 om = np.block(

69 [

70 qq + 2.0 * qu[..., 1:2] ** 2,

71 2.0 * (qu[..., 2:3] * qu[..., 1:2] - _EPSIJK * qu[..., 0:1] * qu[..., 3:4]),

72 2.0 * (qu[..., 3:4] * qu[..., 1:2] + _EPSIJK * qu[..., 0:1] * qu[..., 2:3]),

73 2.0 * (qu[..., 1:2] * qu[..., 2:3] + _EPSIJK * qu[..., 0:1] * qu[..., 3:4]),

74 qq + 2.0 * qu[..., 2:3] ** 2,

75 2.0 * (qu[..., 3:4] * qu[..., 2:3] - _EPSIJK * qu[..., 0:1] * qu[..., 1:2]),

76 2.0 * (qu[..., 1:2] * qu[..., 3:4] - _EPSIJK * qu[..., 0:1] * qu[..., 2:3]),

77 2.0 * (qu[..., 2:3] * qu[..., 3:4] + _EPSIJK * qu[..., 0:1] * qu[..., 1:2]),

78 qq + 2.0 * qu[..., 3:4] ** 2,

79 ]

80 ).reshape(qu.shape[:-1] + (3, 3))

81 return om

84def qu2eu(qu):

85 """Quaternion to Bunge-Euler angles.

87 Args:

88 qu: shape (..., 4) quaternions in the form (w, x, y, z).

90 Returns:

91 np.ndarray of shape (..., 3): Bunge-Euler angles (phi1, Phi, phi2) in radians.

92 """

93 q02 = qu[..., 0:1] * qu[..., 2:3]

94 q13 = qu[..., 1:2] * qu[..., 3:4]

95 q01 = qu[..., 0:1] * qu[..., 1:2]

96 q23 = qu[..., 2:3] * qu[..., 3:4]

97 q03_s = qu[..., 0:1] ** 2 + qu[..., 3:4] ** 2

98 q12_s = qu[..., 1:2] ** 2 + qu[..., 2:3] ** 2

99 chi = np.sqrt(q03_s * q12_s)

100

101 eu = np.where(

102 np.abs(q12_s) < 1.0e-8,

103 np.block(

104 [

105 np.arctan2(

106 -_EPSIJK * 2.0 * qu[..., 0:1] * qu[..., 3:4],

107 qu[..., 0:1] ** 2 - qu[..., 3:4] ** 2,

108 ),

109 np.zeros(qu.shape[:-1] + (2,)),

110 ]

111 ),

112 np.where(

113 np.abs(q03_s) < 1.0e-8,

114 np.block(

115 [

116 np.arctan2(

117 2.0 * qu[..., 1:2] * qu[..., 2:3], qu[..., 1:2] ** 2 - qu[..., 2:3] ** 2

118 ),

119 np.broadcast_to(np.pi, qu.shape[:-1] + (1,)),

120 np.zeros(qu.shape[:-1] + (1,)),

121 ]

122 ),

123 np.block(

124 [

125 np.arctan2((-_EPSIJK * q02 + q13) * chi, (-_EPSIJK * q01 - q23) * chi),

126 np.arctan2(2.0 * chi, q03_s - q12_s),

127 np.arctan2((_EPSIJK * q02 + q13) * chi, (-_EPSIJK * q01 + q23) * chi),

128 ]

129 ),

130 ),

131 )

132 # reduce Euler angles to definition range

133 eu[np.abs(eu) < 1.0e-6] = 0.0

134 eu = np.where(

135 eu < 0, (eu + 2.0 * np.pi) % np.array([2.0 * np.pi, np.pi, 2.0 * np.pi]), eu

136 ) # needed?

137 return eu

138

139

140def qu2ax(qu):

141 """Quaternion to axis-angle pair.

142

143 Modified version of the original formulation, should be numerically more stable.

144

145 Args:

146 qu: shape (..., 4) quaternions in the form (w, x, y, z).

147

148 Returns:

149 np.ndarray of shape (..., 4): axis-angle pairs (n1, n2, n3, omega), with

150 omega in radians.

151 """

152 with np.errstate(invalid="ignore", divide="ignore"):

153 s = np.sign(qu[..., 0:1]) / np.sqrt(

154 qu[..., 1:2] ** 2 + qu[..., 2:3] ** 2 + qu[..., 3:4] ** 2

155 )

156 omega = 2.0 * np.arccos(np.clip(qu[..., 0:1], -1.0, 1.0))

157 ax = np.where(

158 np.broadcast_to(qu[..., 0:1] < 1.0e-8, qu.shape),

159 np.block([qu[..., 1:4], np.broadcast_to(np.pi, qu.shape[:-1] + (1,))]),

160 np.block([qu[..., 1:4] * s, omega]),

161 )

162 ax[np.isclose(qu[..., 0], 1.0, rtol=0.0)] = [0.0, 0.0, 1.0, 0.0]

163 return ax

164

165

166def qu2ro(qu):

167 """Quaternion to Rodrigues-Frank vector.

168

169 Args:

170 qu: shape (..., 4) quaternions in the form (w, x, y, z).

171

172 Returns:

173 np.ndarray of shape (..., 4): Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

174 """

175 with np.errstate(invalid="ignore", divide="ignore"):

176 s = np.linalg.norm(qu[..., 1:4], axis=-1, keepdims=True)

177 ro = np.where(

178 np.broadcast_to(np.abs(qu[..., 0:1]) < 1.0e-12, qu.shape),

179 np.block(

180 [

181 qu[..., 1:2],

182 qu[..., 2:3],

183 qu[..., 3:4],

184 np.broadcast_to(np.inf, qu.shape[:-1] + (1,)),

185 ]

186 ),

187 np.block(

188 [

189 qu[..., 1:2] / s,

190 qu[..., 2:3] / s,

191 qu[..., 3:4] / s,

192 np.tan(np.arccos(np.clip(qu[..., 0:1], -1.0, 1.0))),

193 ]

194 ),

195 )

196 ro[np.abs(s).squeeze(-1) < 1.0e-12] = [0.0, 0.0, _EPSIJK, 0.0]

197 return ro

198

199

200def qu2ho(qu):

201 """Quaternion to homochoric vector.

202

203 Args:

204 qu: shape (..., 4) quaternions in the form (w, x, y, z).

205

206 Returns:

207 np.ndarray of shape (..., 3): homochoric vectors.

208 """

209 with np.errstate(invalid="ignore"):

210 omega = 2.0 * np.arccos(np.clip(qu[..., 0:1], -1.0, 1.0))

211 ho = np.where(

212 np.abs(omega) < 1.0e-12,

213 np.zeros(3),

214 qu[..., 1:4]

215 / np.linalg.norm(qu[..., 1:4], axis=-1, keepdims=True)

216 * (0.75 * (omega - np.sin(omega))) ** (1.0 / 3.0),

217 )

218 return ho

219

220

221def qu2cu(qu):

222 """Quaternion to cubochoric vector.

223

224 Args:

225 qu: shape (..., 4) quaternions in the form (w, x, y, z).

226

227 Returns:

228 np.ndarray of shape (..., 3): cubochoric vectors.

229 """

230 return ho2cu(qu2ho(qu))

231

232

233# ---------- Rotation matrix ----------

234

235

236def om2qu(om):

237 """Rotation matrix to quaternion.

238

239 This formulation is from

240 www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion.

241 The original formulation had issues.

242

243 Args:

244 om: shape (..., 3, 3) rotation matrices.

245

246 Returns:

247 np.ndarray of shape (..., 4): quaternions in the form (w, x, y, z).

248 """

249 trace = om[..., 0, 0:1] + om[..., 1, 1:2] + om[..., 2, 2:3]

250

251 with np.errstate(invalid="ignore", divide="ignore"):

252 s = [

253 0.5 / np.sqrt(1.0 + trace),

254 2.0 * np.sqrt(1.0 + om[..., 0, 0:1] - om[..., 1, 1:2] - om[..., 2, 2:3]),

255 2.0 * np.sqrt(1.0 + om[..., 1, 1:2] - om[..., 2, 2:3] - om[..., 0, 0:1]),

256 2.0 * np.sqrt(1.0 + om[..., 2, 2:3] - om[..., 0, 0:1] - om[..., 1, 1:2]),

257 ]

258 qu = np.where(

259 trace > 0,

260 np.block(

261 [

262 0.25 / s[0],

263 (om[..., 2, 1:2] - om[..., 1, 2:3]) * s[0],

264 (om[..., 0, 2:3] - om[..., 2, 0:1]) * s[0],

265 (om[..., 1, 0:1] - om[..., 0, 1:2]) * s[0],

266 ]

267 ),

268 np.where(

269 om[..., 0, 0:1] > np.maximum(om[..., 1, 1:2], om[..., 2, 2:3]),

270 np.block(

271 [

272 (om[..., 2, 1:2] - om[..., 1, 2:3]) / s[1],

273 0.25 * s[1],

274 (om[..., 0, 1:2] + om[..., 1, 0:1]) / s[1],

275 (om[..., 0, 2:3] + om[..., 2, 0:1]) / s[1],

276 ]

277 ),

278 np.where(

279 om[..., 1, 1:2] > om[..., 2, 2:3],

280 np.block(

281 [

282 (om[..., 0, 2:3] - om[..., 2, 0:1]) / s[2],

283 (om[..., 0, 1:2] + om[..., 1, 0:1]) / s[2],

284 0.25 * s[2],

285 (om[..., 1, 2:3] + om[..., 2, 1:2]) / s[2],

286 ]

287 ),

288 np.block(

289 [

290 (om[..., 1, 0:1] - om[..., 0, 1:2]) / s[3],

291 (om[..., 0, 2:3] + om[..., 2, 0:1]) / s[3],

292 (om[..., 1, 2:3] + om[..., 2, 1:2]) / s[3],

293 0.25 * s[3],

294 ]

295 ),

296 ),

297 ),

298 ) * np.array([1, _EPSIJK, _EPSIJK, _EPSIJK])

299 qu[qu[..., 0] < 0] *= -1

300 return qu

301

302

303def om2eu(om):

304 """Rotation matrix to Bunge-Euler angles.

305

306 Args:

307 om: shape (..., 3, 3) rotation matrices.

308

309 Returns:

310 np.ndarray of shape (..., 3): Bunge-Euler angles (phi1, Phi, phi2) in radians.

311 """

312 with np.errstate(invalid="ignore", divide="ignore"):

313 zeta = 1.0 / np.sqrt(1.0 - om[..., 2, 2:3] ** 2)

314 eu = np.where(

315 np.isclose(np.abs(om[..., 2, 2:3]), 1.0, 1e-9),

316 np.block(

317 [

318 np.arctan2(om[..., 0, 1:2], om[..., 0, 0:1]),

319 np.pi * 0.5 * (1 - om[..., 2, 2:3]),

320 np.zeros(om.shape[:-2] + (1,)),

321 ]

322 ),

323 np.block(

324 [

325 np.arctan2(om[..., 2, 0:1] * zeta, -om[..., 2, 1:2] * zeta),

326 np.arccos(om[..., 2, 2:3]),

327 np.arctan2(om[..., 0, 2:3] * zeta, +om[..., 1, 2:3] * zeta),

328 ]

329 ),

330 )

331 eu[np.abs(eu) < 1.0e-8] = 0.0

332 eu = np.where(eu < 0, (eu + 2.0 * np.pi) % np.array([2.0 * np.pi, np.pi, 2.0 * np.pi]), eu)

333 return eu

334

335

336def om2ax(om):

337 """Rotation matrix to axis-angle pair.

338

339 Args:

340 om: shape (..., 3, 3) rotation matrices.

341

342 Returns:

343 np.ndarray of shape (..., 4): axis-angle pairs (n1, n2, n3, omega), with

344 omega in radians.

345 """

346 # return qu2ax(om2qu(om)) # HOTFIX

347 diag_delta = -_EPSIJK * np.block(

348 [

349 om[..., 1, 2:3] - om[..., 2, 1:2],

350 om[..., 2, 0:1] - om[..., 0, 2:3],

351 om[..., 0, 1:2] - om[..., 1, 0:1],

352 ]

353 )

354 t = 0.5 * (om.trace(axis2=-2, axis1=-1) - 1.0).reshape(om.shape[:-2] + (1,))

355 w, vr = np.linalg.eig(om)

356 # mask duplicated real eigenvalues

357 w[np.isclose(w[..., 0], 1.0 + 0.0j), 1:] = 0.0

358 w[np.isclose(w[..., 1], 1.0 + 0.0j), 2:] = 0.0

359 vr = np.swapaxes(vr, -1, -2)

360 ax = np.where(

361 np.abs(diag_delta) < 1e-12,

362 np.real(vr[np.isclose(w, 1.0 + 0.0j)]).reshape(om.shape[:-2] + (3,)),

363 np.abs(np.real(vr[np.isclose(w, 1.0 + 0.0j)]).reshape(om.shape[:-2] + (3,)))

364 * np.sign(diag_delta),

365 )

366 ax = np.block([ax, np.arccos(np.clip(t, -1.0, 1.0))])

367 ax[np.abs(ax[..., 3]) < 1.0e-8] = [0.0, 0.0, 1.0, 0.0]

368 return ax

369

370

371def om2ro(om):

372 """Rotation matrix to Rodrigues-Frank vector.

373

374 Args:

375 om: shape (..., 3, 3) rotation matrices.

376

377 Returns:

378 np.ndarray of shape (..., 4): Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

379 """

380 return eu2ro(om2eu(om))

381

382

383def om2ho(om):

384 """Rotation matrix to homochoric vector.

385

386 Args:

387 om: shape (..., 3, 3) rotation matrices.

388

389 Returns:

390 np.ndarray of shape (..., 3): homochoric vectors.

391 """

392 return ax2ho(om2ax(om))

393

394

395def om2cu(om):

396 """Rotation matrix to cubochoric vector.

397

398 Args:

399 om: shape (..., 3, 3) rotation matrices.

400

401 Returns:

402 np.ndarray of shape (..., 3): cubochoric vectors.

403 """

404 return ho2cu(om2ho(om))

405

406

407# ---------- Bunge-Euler angles ----------

408

409

410def eu2qu(eu):

411 """Bunge-Euler angles to quaternion.

412

413 Args:

414 eu: shape (..., 3) Bunge-Euler angles (phi1, Phi, phi2) in radians.

415

416 Returns:

417 np.ndarray of shape (..., 4): quaternions in the form (w, x, y, z).

418 """

419 ee = 0.5 * eu

420 cPhi = np.cos(ee[..., 1:2])

421 sPhi = np.sin(ee[..., 1:2])

422 qu = np.block(

423 [

424 cPhi * np.cos(ee[..., 0:1] + ee[..., 2:3]),

425 -_EPSIJK * sPhi * np.cos(ee[..., 0:1] - ee[..., 2:3]),

426 -_EPSIJK * sPhi * np.sin(ee[..., 0:1] - ee[..., 2:3]),

427 -_EPSIJK * cPhi * np.sin(ee[..., 0:1] + ee[..., 2:3]),

428 ]

429 )

430 qu[qu[..., 0] < 0.0] *= -1

431 return qu

432

433

434def eu2om(eu):

435 """Bunge-Euler angles to rotation matrix.

436

437 Args:

438 eu: shape (..., 3) Bunge-Euler angles (phi1, Phi, phi2) in radians.

439

440 Returns:

441 np.ndarray of shape (..., 3, 3): rotation matrices.

442 """

443 c = np.cos(eu)

444 s = np.sin(eu)

445 om = np.block(

446 [

447 +c[..., 0:1] * c[..., 2:3] - s[..., 0:1] * s[..., 2:3] * c[..., 1:2],

448 +s[..., 0:1] * c[..., 2:3] + c[..., 0:1] * s[..., 2:3] * c[..., 1:2],

449 +s[..., 2:3] * s[..., 1:2],

450 -c[..., 0:1] * s[..., 2:3] - s[..., 0:1] * c[..., 2:3] * c[..., 1:2],

451 -s[..., 0:1] * s[..., 2:3] + c[..., 0:1] * c[..., 2:3] * c[..., 1:2],

452 +c[..., 2:3] * s[..., 1:2],

453 +s[..., 0:1] * s[..., 1:2],

454 -c[..., 0:1] * s[..., 1:2],

455 +c[..., 1:2],

456 ]

457 ).reshape(eu.shape[:-1] + (3, 3))

458 om[np.abs(om) < 1.0e-12] = 0.0

459 return om

460

461

462def eu2ax(eu):

463 """Bunge-Euler angles to axis-angle pair.

464

465 Args:

466 eu: shape (..., 3) Bunge-Euler angles (phi1, Phi, phi2) in radians.

467

468 Returns:

469 np.ndarray of shape (..., 4): axis-angle pairs (n1, n2, n3, omega), with

470 omega in radians.

471 """

472 t = np.tan(eu[..., 1:2] * 0.5)

473 sigma = 0.5 * (eu[..., 0:1] + eu[..., 2:3])

474 delta = 0.5 * (eu[..., 0:1] - eu[..., 2:3])

475 tau = np.linalg.norm(np.block([t, np.sin(sigma)]), axis=-1, keepdims=True)

476 alpha = np.where(np.abs(np.cos(sigma)) < 1.0e-12, np.pi, 2.0 * np.arctan(tau / np.cos(sigma)))

477 with np.errstate(invalid="ignore", divide="ignore"):

478 ax = np.where(

479 np.broadcast_to(np.abs(alpha) < 1.0e-12, eu.shape[:-1] + (4,)),

480 [0.0, 0.0, 1.0, 0.0],

481 np.block(

482 [

483 -_EPSIJK / tau * t * np.cos(delta),

484 -_EPSIJK / tau * t * np.sin(delta),

485 -_EPSIJK / tau * np.sin(sigma),

486 alpha,

487 ]

488 ),

489 )

490 ax[(alpha < 0.0).squeeze()] *= -1

491 return ax

492

493

494def eu2ro(eu):

495 """Bunge-Euler angles to Rodrigues-Frank vector.

496

497 Args:

498 eu: shape (..., 3) Bunge-Euler angles (phi1, Phi, phi2) in radians.

499

500 Returns:

501 np.ndarray of shape (..., 4): Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

502 """

503 ax = eu2ax(eu)

504 ro = np.block([ax[..., :3], np.tan(ax[..., 3:4] * 0.5)])

505 ro[ax[..., 3] >= np.pi, 3] = np.inf

506 ro[np.abs(ax[..., 3]) < 1.0e-16] = [0.0, 0.0, _EPSIJK, 0.0]

507 return ro

508

509

510def eu2ho(eu):

511 """Bunge-Euler angles to homochoric vector.

512

513 Args:

514 eu: shape (..., 3) Bunge-Euler angles (phi1, Phi, phi2) in radians.

515

516 Returns:

517 np.ndarray of shape (..., 3): homochoric vectors.

518 """

519 return ax2ho(eu2ax(eu))

520

521

522def eu2cu(eu):

523 """Bunge-Euler angles to cubochoric vector.

524

525 Args:

526 eu: shape (..., 3) Bunge-Euler angles (phi1, Phi, phi2) in radians.

527

528 Returns:

529 np.ndarray of shape (..., 3): cubochoric vectors.

530 """

531 return ho2cu(eu2ho(eu))

532

533

534# ---------- Axis angle pair ----------

535

536

537def ax2qu(ax):

538 """Axis-angle pair to quaternion.

539

540 Args:

541 ax: shape (..., 4) axis-angle pairs (n1, n2, n3, omega), with omega in radians.

542

543 Returns:

544 np.ndarray of shape (..., 4): quaternions in the form (w, x, y, z).

545 """

546 c = np.cos(ax[..., 3:4] * 0.5)

547 s = np.sin(ax[..., 3:4] * 0.5)

548 qu = np.where(

549 np.abs(ax[..., 3:4]) < 1.0e-6, [1.0, 0.0, 0.0, 0.0], np.block([c, ax[..., :3] * s])

550 )

551 return qu

552

553

554def ax2om(ax):

555 """Axis-angle pair to rotation matrix.

556

557 Args:

558 ax: shape (..., 4) axis-angle pairs (n1, n2, n3, omega), with omega in radians.

559

560 Returns:

561 np.ndarray of shape (..., 3, 3): rotation matrices.

562 """

563 c = np.cos(ax[..., 3:4])

564 s = np.sin(ax[..., 3:4])

565 omc = 1.0 - c

566 om = np.block(

567 [

568 c + omc * ax[..., 0:1] ** 2,

569 omc * ax[..., 0:1] * ax[..., 1:2] + s * ax[..., 2:3],

570 omc * ax[..., 0:1] * ax[..., 2:3] - s * ax[..., 1:2],

571 omc * ax[..., 0:1] * ax[..., 1:2] - s * ax[..., 2:3],

572 c + omc * ax[..., 1:2] ** 2,

573 omc * ax[..., 1:2] * ax[..., 2:3] + s * ax[..., 0:1],

574 omc * ax[..., 0:1] * ax[..., 2:3] + s * ax[..., 1:2],

575 omc * ax[..., 1:2] * ax[..., 2:3] - s * ax[..., 0:1],

576 c + omc * ax[..., 2:3] ** 2,

577 ]

578 ).reshape(ax.shape[:-1] + (3, 3))

579 return om if _EPSIJK < 0.0 else np.swapaxes(om, -1, -2)

580

581

582def ax2eu(ax):

583 """Axis-angle pair to Bunge-Euler angles.

584

585 Args:

586 ax: shape (..., 4) axis-angle pairs (n1, n2, n3, omega), with omega in radians.

587

588 Returns:

589 np.ndarray of shape (..., 3): Bunge-Euler angles (phi1, Phi, phi2) in radians.

590 """

591 return om2eu(ax2om(ax))

592

593

594def ax2ro(ax):

595 """Axis-angle pair to Rodrigues-Frank vector.

596

597 Args:

598 ax: shape (..., 4) axis-angle pairs (n1, n2, n3, omega), with omega in radians.

599

600 Returns:

601 np.ndarray of shape (..., 4): Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

602 """

603 ro = np.block(

604 [

605 ax[..., :3],

606 np.where(

607 np.isclose(ax[..., 3:4], np.pi, atol=1.0e-15, rtol=0.0),

608 np.inf,

609 np.tan(ax[..., 3:4] * 0.5),

610 ),

611 ]

612 )

613 ro[np.abs(ax[..., 3]) < 1.0e-6] = [0.0, 0.0, _EPSIJK, 0.0]

614 return ro

615

616

617def ax2ho(ax):

618 """Axis-angle pair to homochoric vector.

619

620 Args:

621 ax: shape (..., 4) axis-angle pairs (n1, n2, n3, omega), with omega in radians.

622

623 Returns:

624 np.ndarray of shape (..., 3): homochoric vectors.

625 """

626 f = (0.75 * (ax[..., 3:4] - np.sin(ax[..., 3:4]))) ** (1.0 / 3.0)

627 ho = ax[..., :3] * f

628 return ho

629

630

631def ax2cu(ax):

632 """Axis-angle pair to cubochoric vector.

633

634 Args:

635 ax: shape (..., 4) axis-angle pairs (n1, n2, n3, omega), with omega in radians.

636

637 Returns:

638 np.ndarray of shape (..., 3): cubochoric vectors.

639 """

640 return ho2cu(ax2ho(ax))

641

642

643# ---------- Rodrigues-Frank vector ----------

644

645

646def ro2qu(ro):

647 """Rodrigues-Frank vector to quaternion.

648

649 Args:

650 ro: shape (..., 4) Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

651

652 Returns:

653 np.ndarray of shape (..., 4): quaternions in the form (w, x, y, z).

654 """

655 return ax2qu(ro2ax(ro))

656

657

658def ro2om(ro):

659 """Rodrigues-Frank vector to rotation matrix.

660

661 Args:

662 ro: shape (..., 4) Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

663

664 Returns:

665 np.ndarray of shape (..., 3, 3): rotation matrices.

666 """

667 return ax2om(ro2ax(ro))

668

669

670def ro2eu(ro):

671 """Rodrigues-Frank vector to Bunge-Euler angles.

672

673 Args:

674 ro: shape (..., 4) Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

675

676 Returns:

677 np.ndarray of shape (..., 3): Bunge-Euler angles (phi1, Phi, phi2) in radians.

678 """

679 return om2eu(ro2om(ro))

680

681

682def ro2ax(ro):

683 """Rodrigues-Frank vector to axis-angle pair.

684

685 Args:

686 ro: shape (..., 4) Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

687

688 Returns:

689 np.ndarray of shape (..., 4): axis-angle pairs (n1, n2, n3, omega), with

690 omega in radians.

691 """

692 with np.errstate(invalid="ignore", divide="ignore"):

693 ax = np.where(

694 np.isfinite(ro[..., 3:4]),

695 np.block(

696 [

697 ro[..., 0:3] / np.linalg.norm(ro[..., 0:3], axis=-1, keepdims=True),

698 2.0 * np.arctan(ro[..., 3:4]),

699 ]

700 ),

701 np.block([ro[..., 0:3], np.broadcast_to(np.pi, ro[..., 3:4].shape)]),

702 )

703 ax[np.abs(ro[..., 3]) < 1.0e-8] = np.array([0.0, 0.0, 1.0, 0.0])

704 return ax

705

706

707def ro2ho(ro):

708 """Rodrigues-Frank vector to homochoric vector.

709

710 Args:

711 ro: shape (..., 4) Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

712

713 Returns:

714 np.ndarray of shape (..., 3): homochoric vectors.

715 """

716 f = np.where(

717 np.isfinite(ro[..., 3:4]),

718 2.0 * np.arctan(ro[..., 3:4]) - np.sin(2.0 * np.arctan(ro[..., 3:4])),

719 np.pi,

720 )

721 ho = np.where(

722 np.broadcast_to(

723 np.sum(ro[..., 0:3] ** 2.0, axis=-1, keepdims=True) < 1.0e-8, ro[..., 0:3].shape

724 ),

725 np.zeros(3),

726 ro[..., 0:3] * (0.75 * f) ** (1.0 / 3.0),

727 )

728 return ho

729

730

731def ro2cu(ro):

732 """Rodrigues-Frank vector to cubochoric vector.

733

734 Args:

735 ro: shape (..., 4) Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

736

737 Returns:

738 np.ndarray of shape (..., 3): cubochoric vectors.

739 """

740 return ho2cu(ro2ho(ro))

741

742

743# ---------- Homochoric vector----------

744

745

746def ho2qu(ho):

747 """Homochoric vector to quaternion.

748

749 Args:

750 ho: shape (..., 3) homochoric vectors.

751

752 Returns:

753 np.ndarray of shape (..., 4): quaternions in the form (w, x, y, z).

754 """

755 return ax2qu(ho2ax(ho))

756

757

758def ho2om(ho):

759 """Homochoric vector to rotation matrix.

760

761 Args:

762 ho: shape (..., 3) homochoric vectors.

763

764 Returns:

765 np.ndarray of shape (..., 3, 3): rotation matrices.

766 """

767 return ax2om(ho2ax(ho))

768

769

770def ho2eu(ho):

771 """Homochoric vector to Bunge-Euler angles.

772

773 Args:

774 ho: shape (..., 3) homochoric vectors.

775

776 Returns:

777 np.ndarray of shape (..., 3): Bunge-Euler angles (phi1, Phi, phi2) in radians.

778 """

779 return ax2eu(ho2ax(ho))

780

781

782def ho2ax(ho):

783 """Homochoric vector to axis-angle pair.

784

785 Args:

786 ho: shape (..., 3) homochoric vectors.

787

788 Returns:

789 np.ndarray of shape (..., 4): axis-angle pairs (n1, n2, n3, omega), with

790 omega in radians.

791 """

792 tfit = np.array(

793 [

794 +1.0000000000018852,

795 -0.5000000002194847,

796 -0.024999992127593126,

797 -0.003928701544781374,

798 -0.0008152701535450438,

799 -0.0002009500426119712,

800 -0.00002397986776071756,

801 -0.00008202868926605841,

802 +0.00012448715042090092,

803 -0.0001749114214822577,

804 +0.0001703481934140054,

805 -0.00012062065004116828,

806 +0.000059719705868660826,

807 -0.00001980756723965647,

808 +0.000003953714684212874,

809 -0.00000036555001439719544,

810 ]

811 )

812 hmag_squared = np.sum(ho**2.0, axis=-1, keepdims=True)

813 hm = hmag_squared.copy()

814 s = tfit[0] + tfit[1] * hmag_squared

815 for i in range(2, 16):

816 hm *= hmag_squared

817 s += tfit[i] * hm

818 with np.errstate(invalid="ignore"):

819 ax = np.where(

820 np.broadcast_to(np.abs(hmag_squared) < 1.0e-8, ho.shape[:-1] + (4,)),

821 [0.0, 0.0, 1.0, 0.0],

822 np.block([ho / np.sqrt(hmag_squared), 2.0 * np.arccos(np.clip(s, -1.0, 1.0))]),

823 )

824 return ax

825

826

827def ho2ro(ho):

828 """Homochoric vector to Rodrigues-Frank vector.

829

830 Args:

831 ho: shape (..., 3) homochoric vectors.

832

833 Returns:

834 np.ndarray of shape (..., 4): Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

835 """

836 return ax2ro(ho2ax(ho))

837

838

839def ho2cu(ho):

840 """Homochoric vector to cubochoric vector.

841

842 Args:

843 ho: shape (..., 3) homochoric vectors.

844

845 Returns:

846 np.ndarray of shape (..., 3): cubochoric vectors.

847

848 References:

849 D. Roşca et al., Modelling and Simulation in Materials Science and

850 Engineering 22:075013, 2014. https://doi.org/10.1088/0965-0393/22/7/075013

851 """

852 rs = np.linalg.norm(ho, axis=-1, keepdims=True)

853

854 xyz3 = np.take_along_axis(ho, _get_pyramid_order(ho, "forward"), -1)

855

856 with np.errstate(invalid="ignore", divide="ignore"):

857 # inverse M_3

858 xyz2 = xyz3[..., 0:2] * np.sqrt(2.0 * rs / (rs + np.abs(xyz3[..., 2:3])))

859 qxy = np.sum(xyz2**2, axis=-1, keepdims=True)

860

861 q2 = qxy + np.max(np.abs(xyz2), axis=-1, keepdims=True) ** 2

862 sq2 = np.sqrt(q2)

863 q = (_BETA / np.sqrt(2.0) / _R1) * np.sqrt(

864 q2 * qxy / (q2 - np.max(np.abs(xyz2), axis=-1, keepdims=True) * sq2)

865 )

866 tt = np.clip(

867 (

868 np.min(np.abs(xyz2), axis=-1, keepdims=True) ** 2

869 + np.max(np.abs(xyz2), axis=-1, keepdims=True) * sq2

870 )

871 / np.sqrt(2.0)

872 / qxy,

873 -1.0,

874 1.0,

875 )

876 T_inv = (

877 np.where(

878 np.abs(xyz2[..., 1:2]) <= np.abs(xyz2[..., 0:1]),

879 np.block([np.ones_like(tt), np.arccos(tt) / np.pi * 12.0]),

880 np.block([np.arccos(tt) / np.pi * 12.0, np.ones_like(tt)]),

881 )

882 * q

883 )

884 T_inv[xyz2 < 0.0] *= -1.0

885 T_inv[np.broadcast_to(np.isclose(qxy, 0.0, rtol=0.0, atol=1.0e-12), T_inv.shape)] = 0.0

886 cu = (

887 np.block(

888 [

889 T_inv,

890 np.where(

891 xyz3[..., 2:3] < 0.0,

892 -np.ones_like(xyz3[..., 2:3]),

893 np.ones_like(xyz3[..., 2:3]),

894 )

895 * rs

896 / np.sqrt(6.0 / np.pi),

897 ]

898 )

899 / _SC

900 )

901

902 cu[np.isclose(np.sum(np.abs(ho), axis=-1), 0.0, rtol=0.0, atol=1.0e-16)] = 0.0

903 cu = np.take_along_axis(cu, _get_pyramid_order(ho, "backward"), -1)

904

905 return cu

906

907

908# ---------- Cubochoric ----------

909

910

911def cu2qu(cu):

912 """Cubochoric vector to quaternion.

913

914 Args:

915 cu: shape (..., 3) cubochoric vectors.

916

917 Returns:

918 np.ndarray of shape (..., 4): quaternions in the form (w, x, y, z).

919 """

920 return ho2qu(cu2ho(cu))

921

922

923def cu2om(cu):

924 """Cubochoric vector to rotation matrix.

925

926 Args:

927 cu: shape (..., 3) cubochoric vectors.

928

929 Returns:

930 np.ndarray of shape (..., 3, 3): rotation matrices.

931 """

932 return ho2om(cu2ho(cu))

933

934

935def cu2eu(cu):

936 """Cubochoric vector to Bunge-Euler angles.

937

938 Args:

939 cu: shape (..., 3) cubochoric vectors.

940

941 Returns:

942 np.ndarray of shape (..., 3): Bunge-Euler angles (phi1, Phi, phi2) in radians.

943 """

944 return ho2eu(cu2ho(cu))

945

946

947def cu2ax(cu):

948 """Cubochoric vector to axis-angle pair.

949

950 Args:

951 cu: shape (..., 3) cubochoric vectors.

952

953 Returns:

954 np.ndarray of shape (..., 4): axis-angle pairs (n1, n2, n3, omega), with

955 omega in radians.

956 """

957 return ho2ax(cu2ho(cu))

958

959

960def cu2ro(cu):

961 """Cubochoric vector to Rodrigues-Frank vector.

962

963 Args:

964 cu: shape (..., 3) cubochoric vectors.

965

966 Returns:

967 np.ndarray of shape (..., 4): Rodrigues-Frank vectors (n1, n2, n3, tan(omega / 2)).

968 """

969 return ho2ro(cu2ho(cu))

970

971

972def cu2ho(cu):

973 """Cubochoric vector to homochoric vector.

974

975 Args:

976 cu: shape (..., 3) cubochoric vectors.

977

978 Returns:

979 np.ndarray of shape (..., 3): homochoric vectors.

980

981 References:

982 D. Roşca et al., Modelling and Simulation in Materials Science and

983 Engineering 22:075013, 2014. https://doi.org/10.1088/0965-0393/22/7/075013

984 """

985 with np.errstate(invalid="ignore", divide="ignore"):

986 # get pyramide and scale by grid parameter ratio

987 XYZ = np.take_along_axis(cu, _get_pyramid_order(cu, "forward"), -1) * _SC

988 order = np.abs(XYZ[..., 1:2]) <= np.abs(XYZ[..., 0:1])

989 q = (

990 np.pi

991 / 12.0

992 * np.where(order, XYZ[..., 1:2], XYZ[..., 0:1])

993 / np.where(order, XYZ[..., 0:1], XYZ[..., 1:2])

994 )

995 c = np.cos(q)

996 s = np.sin(q)

997 q = (

998 _R1

999 * 2.0**0.25

1000 / _BETA

1001 / np.sqrt(np.sqrt(2.0) - c)

1002 * np.where(order, XYZ[..., 0:1], XYZ[..., 1:2])

1003 )

1004

1005 T = np.block([(np.sqrt(2.0) * c - 1.0), np.sqrt(2.0) * s]) * q

1006

1007 # transform to sphere grid (inverse Lambert)

1008 c = np.sum(T**2, axis=-1, keepdims=True)

1009 s = c * np.pi / 24.0 / XYZ[..., 2:3] ** 2

1010 c = c * np.sqrt(np.pi / 24.0) / XYZ[..., 2:3]

1011 q = np.sqrt(1.0 - s)

1012

1013 ho = np.where(

1014 np.isclose(

1015 np.sum(np.abs(XYZ[..., 0:2]), axis=-1, keepdims=True), 0.0, rtol=0.0, atol=1.0e-16

1016 ),

1017 np.block([np.zeros_like(XYZ[..., 0:2]), np.sqrt(6.0 / np.pi) * XYZ[..., 2:3]]),

1018 np.block(

1019 [

1020 np.where(order, T[..., 0:1], T[..., 1:2]) * q,

1021 np.where(order, T[..., 1:2], T[..., 0:1]) * q,

1022 np.sqrt(6.0 / np.pi) * XYZ[..., 2:3] - c,

1023 ]

1024 ),

1025 )

1026

1027 ho[np.isclose(np.sum(np.abs(cu), axis=-1), 0.0, rtol=0.0, atol=1.0e-16)] = 0.0

1028 ho = np.take_along_axis(ho, _get_pyramid_order(cu, "backward"), -1)

1029

1030 return ho

1031

1032

1033def _get_pyramid_order(xyz, direction=None):

1034 """Get the coordinate order for the cubochoric/homochoric pyramid transform.

1035

1036 Depending on which of the six pyramids a point falls into, the coordinate

1037 order needs to be cyclically permuted before (`"forward"`) and after

1038 (`"backward"`) the transform.

1039

1040 Args:

1041 xyz: shape (..., 3) coordinates of a point on a uniform refinable grid

1042 on a ball, or in a uniform refinable cubical grid.

1043 direction: either `"forward"` or `"backward"`.

1044

1045 Returns:

1046 np.ndarray of shape (..., 3): the index order to apply via `np.take_along_axis`.

1047

1048 References:

1049 D. Roşca et al., Modelling and Simulation in Materials Science and

1050 Engineering 22:075013, 2014. https://doi.org/10.1088/0965-0393/22/7/075013

1051 """

1052 order = {

1053 "forward": np.array([[0, 1, 2], [1, 2, 0], [2, 0, 1]]),

1054 "backward": np.array([[0, 1, 2], [2, 0, 1], [1, 2, 0]]),

1055 }

1056

1057 p = np.where(

1058 np.maximum(np.abs(xyz[..., 0]), np.abs(xyz[..., 1])) <= np.abs(xyz[..., 2]),

1059 0,

1060 np.where(np.maximum(np.abs(xyz[..., 1]), np.abs(xyz[..., 2])) <= np.abs(xyz[..., 0]), 1, 2),

1061 )

1062

1063 return order[direction][p]

1064

1065

1066if __name__ == "__main__":

1067

1068 def qu_norm(qu: np.ndarray) -> np.ndarray:

1069 """

1070 Normalize quaternions to unit norm.

1071

1072 Args:

1073 qu: shape (..., 4) quaternions in form (w, x, y, z)

1074

1075 Returns:

1076 np.ndarray of normalized quaternions.

1077 """

1078 norms = np.linalg.norm(qu, axis=-1, keepdims=True)

1079 with np.errstate(divide="ignore", invalid="ignore"):

1080 return np.where(norms > 0, qu / norms, 0)

1081

1082 qu = np.random.rand(100, 4)

1083 qu = qu_norm(qu)

1084 ax = qu2ax(qu)

1085 qub = ax2qu(ax)

1086 print("Max difference quaternion-axisquaternion:", np.max(np.abs(qu - qub)))