0
# Numpy Quaternion
1

2
A Python library that extends NumPy with a quaternion data type, providing comprehensive quaternion arithmetic and mathematical operations optimized for scientific computing. It implements fast ufuncs for quaternion operations, supports seamless conversion between quaternions and various rotation representations, and includes advanced functionality for interpolation, integration, and statistical operations on quaternion arrays.
3

4
## Package Information
5

6
- **Package Name**: numpy-quaternion
7
- **Language**: Python  
8
- **Installation**: `pip install numpy-quaternion`
9
- **Requirements**: Python >=3.10, numpy >=1.25, scipy >=1.5 (optional)
10

11
## Core Imports
12

13
```python
14
import quaternion
15
import numpy as np
16
```
17

18
For direct access to specific functions:
19

20
```python
21
from quaternion import (
22
    as_rotation_matrix, from_rotation_matrix,
23
    as_euler_angles, from_euler_angles, 
24
    slerp, squad, integrate_angular_velocity
25
)
26
```
27

28
## Basic Usage
29

30
```python
31
import quaternion
32
import numpy as np
33

34
# Create quaternions
35
q1 = quaternion.quaternion(1, 0, 0, 0)  # Identity quaternion
36
q2 = quaternion.quaternion(0, 1, 0, 0)  # Pure quaternion (i)
37

38
# Use predefined constants
39
identity = quaternion.one
40
x_rotation = quaternion.x
41
y_rotation = quaternion.y
42
z_rotation = quaternion.z
43

44
# Create quaternion arrays
45
q_array = np.array([q1, q2, quaternion.y, quaternion.z])
46

47
# Basic arithmetic
48
result = q1 * q2  # Quaternion multiplication
49
conjugate = q1.conjugate()  # Quaternion conjugate
50
magnitude = abs(q1)  # Quaternion norm
51

52
# Convert to rotation matrix
53
rotation_matrix = quaternion.as_rotation_matrix(q1)
54

55
# Create from Euler angles (radians)
56
alpha, beta, gamma = 0.1, 0.2, 0.3
57
q_from_euler = quaternion.from_euler_angles(alpha, beta, gamma)
58

59
# Interpolate between quaternions
60
t_in = np.linspace(0, 1, 5)
61
t_out = np.linspace(0, 1, 10)
62
R_in = np.array([quaternion.one, quaternion.x, quaternion.y])
63
R_out = quaternion.slerp(R_in, t_in, t_out)
64
```
65

66
## Architecture
67

68
The package is built around three main components:
69

70
- **C Extension Core**: Implements the fundamental quaternion dtype with optimized ufuncs for arithmetic, mathematical operations, and comparisons
71
- **Conversion Layer**: Provides seamless conversion between quaternions and rotation matrices, Euler angles, axis-angle representations, and spherical coordinates  
72
- **Advanced Operations**: Time series interpolation, calculus operations, and statistical functions for working with quaternion arrays
73

74
The quaternion dtype integrates directly with NumPy's broadcasting and array operations, enabling efficient vectorized computation on quaternion arrays.
75

76
## Capabilities
77

78
### Core Quaternion Operations
79

80
Fundamental quaternion type with arithmetic operations, mathematical functions, and array conversions. Includes predefined quaternion constants and efficient array manipulation functions.
81

82
```python { .api }
83
# Core quaternion type and constants
84
quaternion(w, x, y, z)  # Constructor
85
quaternion.zero         # quaternion(0, 0, 0, 0)
86
quaternion.one          # quaternion(1, 0, 0, 0)  
87
quaternion.x            # quaternion(0, 1, 0, 0)
88
quaternion.y            # quaternion(0, 0, 1, 0)
89
quaternion.z            # quaternion(0, 0, 0, 1)
90

91
# Array conversion functions
92
def as_float_array(a): ...          # -> ndarray with shape (..., 4)
93
def as_quat_array(a): ...           # -> quaternion array
94
def as_vector_part(q): ...          # -> ndarray with shape (..., 3)
95
def from_vector_part(v, vector_axis=-1): ...  # -> quaternion array
96
def as_spinor_array(a): ...         # -> complex array (..., 2)
97
```
98

99
[Core Operations](./core-operations.md)
100

101
### Rotation Conversions
102

103
Convert between quaternions and various rotation representations including rotation matrices, axis-angle vectors, Euler angles, and spherical coordinates.
104

105
```python { .api }
106
def as_rotation_matrix(q): ...      # -> ndarray (..., 3, 3)
107
def from_rotation_matrix(rot, nonorthogonal=True): ...  # -> quaternion array
108
def as_rotation_vector(q): ...      # -> ndarray (..., 3)
109
def from_rotation_vector(rot): ...  # -> quaternion array
110
def as_euler_angles(q): ...         # -> ndarray (..., 3)
111
def from_euler_angles(alpha_beta_gamma, beta=None, gamma=None): ...  # -> quaternion array
112
def as_spherical_coords(q): ...     # -> ndarray (..., 2)  
113
def from_spherical_coords(theta_phi, phi=None): ...  # -> quaternion array
114
def rotate_vectors(R, v, axis=-1): ...  # -> ndarray
115
```
116

117
[Rotation Conversions](./rotation-conversions.md)
118

119
### Time Series and Interpolation
120

121
Spherical interpolation methods for smooth quaternion time series, angular velocity integration, and time series analysis tools for rotational data.
122

123
```python { .api }
124
def slerp(R1, R2, t1, t2, t_out): ...  # -> quaternion array
125
def squad(R_in, t_in, t_out, unflip_input_rotors=False): ...  # -> quaternion array
126
def integrate_angular_velocity(Omega, t0, t1, R0=None, tolerance=1e-12): ...  # -> quaternion array
127
def unflip_rotors(q, axis=-1, inplace=False): ...  # -> quaternion array
128
def angular_velocity(R, t): ...     # -> ndarray (..., 3)
129
def minimal_rotation(R, t, iterations=2): ...  # -> quaternion array
130

131
# Low-level evaluation functions
132
def slerp_evaluate(R1, R2, t):
133
    """
134
    Low-level spherical linear interpolation evaluation.
135
    
136
    Args:
137
        R1 (quaternion): Starting quaternion (must be normalized)
138
        R2 (quaternion): Ending quaternion (must be normalized)  
139
        t (float or array): Interpolation parameter(s) in [0, 1]
140
        
141
    Returns:
142
        quaternion or quaternion array: Interpolated quaternion(s)
143
        
144
    Notes:
145
        Direct evaluation function used internally by slerp().
146
        Assumes inputs are properly normalized unit quaternions.
147
    """
148

149
def squad_evaluate(R0, R1, R2, R3, t):
150
    """
151
    Low-level spherical quadrangle interpolation evaluation.
152
    
153
    Args:
154
        R0 (quaternion): Control point before start (normalized)
155
        R1 (quaternion): Start quaternion (normalized)
156
        R2 (quaternion): End quaternion (normalized)
157
        R3 (quaternion): Control point after end (normalized)
158
        t (float or array): Interpolation parameter(s) in [0, 1]
159
        
160
    Returns:
161
        quaternion or quaternion array: Interpolated quaternion(s)
162
        
163
    Notes:
164
        Direct evaluation function used internally by squad().
165
        Provides C1 continuous interpolation using spherical Bézier curves.
166
    """
167
```
168

169
[Time Series and Interpolation](./time-series.md)
170

171
### Calculus Operations
172

173
Numerical differentiation and integration methods for quaternion arrays with multiple implementation options.
174

175
```python { .api }
176
def derivative(f, t):
177
    """
178
    Primary derivative function using best available method.
179
    
180
    Args:
181
        f (array_like): Function values to differentiate
182
        t (array_like): Time values
183
        
184
    Returns:
185
        ndarray: Derivative values with same shape as f
186
        
187
    Notes:
188
        Uses spline_derivative if SciPy available, otherwise finite difference.
189
        Automatically selects most accurate method for the given data.
190
    """
191

192
def definite_integral(f, t):
193
    """
194
    Primary definite integral function using best available method.
195
    
196
    Args:
197
        f (array_like): Function values to integrate
198
        t (array_like): Time values
199
        
200
    Returns:
201
        ndarray: Definite integral values
202
        
203
    Notes:
204
        Uses spline integration if SciPy available, otherwise trapezoidal rule.
205
        Automatically selects most accurate method.
206
    """
207

208
def indefinite_integral(f, t):
209
    """
210
    Primary indefinite integral function using best available method.
211
    
212
    Args:
213
        f (array_like): Function values to integrate
214
        t (array_like): Time values
215
        
216
    Returns:
217
        ndarray: Indefinite integral values with same shape as f
218
        
219
    Notes:
220
        Uses spline integration if SciPy available, otherwise trapezoidal rule.
221
        Automatically selects most accurate method.
222
    """
223

224
# Finite difference implementations
225
def fd_derivative(f, t):
226
    """
227
    Fourth-order finite difference derivative for non-uniform time steps.
228
    
229
    Args:
230
        f (array_like): Function values (1D, 2D, or 3D arrays)
231
        t (array_like): Time values (must be 1D)
232
        
233
    Returns:
234
        ndarray: Derivative values with same shape as f
235
        
236
    Notes:
237
        Uses fourth-order finite difference formula from Bowen & Smith.
238
        Requires at least 5 points for full accuracy.
239
    """
240

241
def fd_definite_integral(f, t):
242
    """
243
    Finite difference definite integral using trapezoidal rule.
244
    
245
    Args:
246
        f (array_like): Function values to integrate
247
        t (array_like): Time values
248
        
249
    Returns:
250
        ndarray: Definite integral values
251
    """
252

253
def fd_indefinite_integral(f, t):
254
    """
255
    Finite difference indefinite integral (cumulative integration).
256
    
257
    Args:
258
        f (array_like): Function values to integrate
259
        t (array_like): Time values
260
        
261
    Returns:
262
        ndarray: Indefinite integral values with same shape as f
263
    """
264

265
# Spline-based implementations (require SciPy)
266
def spline_derivative(f, t):
267
    """
268
    Spline-based derivative using cubic spline interpolation.
269
    
270
    Args:
271
        f (array_like): Function values
272
        t (array_like): Time values
273
        
274
    Returns:
275
        ndarray: Derivative values computed from spline
276
        
277
    Notes:
278
        More accurate than finite differences for smooth data.
279
        Requires SciPy. Falls back to finite differences if unavailable.
280
    """
281

282
def spline_definite_integral(f, t):
283
    """
284
    Spline-based definite integral using cubic spline interpolation.
285
    
286
    Args:
287
        f (array_like): Function values to integrate
288
        t (array_like): Time values
289
        
290
    Returns:
291
        ndarray: Definite integral computed from spline
292
        
293
    Notes:
294
        More accurate than trapezoidal rule for smooth data.
295
    """
296

297
def spline_indefinite_integral(f, t):
298
    """
299
    Spline-based indefinite integral using cubic spline interpolation.
300
    
301
    Args:
302
        f (array_like): Function values to integrate
303
        t (array_like): Time values
304
        
305
    Returns:
306
        ndarray: Indefinite integral values computed from spline
307
    """
308

309
def spline(f, t, t_out=None):
310
    """
311
    Create cubic spline interpolation of quaternion time series.
312
    
313
    Args:
314
        f (array_like): Function values to interpolate
315
        t (array_like): Time values corresponding to f
316
        t_out (array_like, optional): Output time values. If None, returns spline object.
317
        
318
    Returns:
319
        ndarray or spline object: Interpolated values or spline interpolator
320
        
321
    Notes:
322
        Available only when SciPy is installed.
323
        Returns InterpolatedUnivariateSpline object if t_out is None.
324
    """
325
```
326

327
### Distance and Comparison Functions
328

329
Comparison functions with configurable tolerances and distance metrics for quaternions and rotations.
330

331
```python { .api }
332
def allclose(a, b, rtol=4*np.finfo(float).eps, atol=0.0, equal_nan=False, verbose=False):
333
    """
334
    Test if all quaternions in arrays are equal within tolerance.
335
    
336
    Args:
337
        a (array_like): First quaternion array
338
        b (array_like): Second quaternion array
339
        rtol (float): Relative tolerance parameter (default: 4*machine epsilon)  
340
        atol (float): Absolute tolerance parameter (default: 0.0)
341
        equal_nan (bool): Whether to consider NaN values as equal (default: False)
342
        verbose (bool): Print non-close values if False is returned (default: False)
343
        
344
    Returns:
345
        bool: True if all elements are close within tolerance, False otherwise
346
    """
347

348
def rotor_intrinsic_distance(a, b):
349
    """
350
    Intrinsic (geodesic) distance between unit quaternions.
351
    
352
    Args:
353
        a (quaternion array): First quaternion array (must be normalized)
354
        b (quaternion array): Second quaternion array (must be normalized)
355
        
356
    Returns:
357
        ndarray: Intrinsic distances with same broadcasting shape as inputs
358
        
359
    Notes:
360
        Computes geodesic distance on quaternion manifold (4D unit sphere).
361
        Distance = 2 * arccos(|⟨a,b⟩|) where ⟨⟩ is quaternion inner product.
362
    """
363

364
def rotor_chordal_distance(a, b):
365
    """
366
    Chordal (Euclidean) distance between unit quaternions.
367
    
368
    Args:
369
        a (quaternion array): First quaternion array (must be normalized)
370
        b (quaternion array): Second quaternion array (must be normalized)
371
        
372
    Returns:
373
        ndarray: Chordal distances with same broadcasting shape as inputs
374
        
375
    Notes:
376
        Standard Euclidean distance in 4D space: |a - b|.
377
        Simpler to compute than intrinsic distance.
378
    """
379

380
def rotation_intrinsic_distance(a, b):
381
    """
382
    Intrinsic distance between rotations represented by quaternions.
383
    
384
    Args:
385
        a (quaternion array): First quaternion array (must be normalized)
386
        b (quaternion array): Second quaternion array (must be normalized)
387
        
388
    Returns:
389
        ndarray: Rotation distances with same broadcasting shape as inputs
390
        
391
    Notes:
392
        Accounts for quaternion double-cover: ±q represent same rotation.
393
        Distance = min(rotor_intrinsic_distance(a,b), rotor_intrinsic_distance(a,-b)).
394
    """
395

396
def rotation_chordal_distance(a, b):
397
    """
398
    Chordal distance between rotations represented by quaternions.
399
    
400
    Args:
401
        a (quaternion array): First quaternion array (must be normalized)
402
        b (quaternion array): Second quaternion array (must be normalized)
403
        
404
    Returns:
405
        ndarray: Rotation chordal distances with same broadcasting shape as inputs
406
        
407
    Notes:
408
        Accounts for quaternion double-cover: ±q represent same rotation.
409
        Distance = min(|a - b|, |a + b|).
410
    """
411
```
412

413
### Statistical Analysis Functions
414

415
Advanced statistical operations for quaternion arrays including mean computation and optimal alignment.
416

417
```python { .api }
418
def mean_rotor_in_chordal_metric(R, t=None):
419
    """
420
    Compute mean quaternion using chordal metric (least-squares sense).
421
    
422
    Args:
423
        R (quaternion array): Input quaternions (should be normalized)
424
        t (array_like, optional): Time weights for integral. If None, uses simple sum.
425
        
426
    Returns:
427
        quaternion: Normalized mean quaternion
428
        
429
    Raises:
430
        ValueError: If t is provided but len(t) < 4 or len(R) < 4
431
        
432
    Notes:
433
        Uses spline integration when t is provided (requires ≥4 points).
434
        Without t, computes simple arithmetic mean and normalizes.
435
    """
436

437
def optimal_alignment_in_chordal_metric(Ra, Rb, t=None):
438
    """
439
    Find rotation Rd such that Rd*Rb is closest to Ra in chordal metric.
440
    
441
    Args:
442
        Ra (quaternion array): Target quaternions
443
        Rb (quaternion array): Base quaternions to align
444
        t (array_like, optional): Time weights for integral
445
        
446
    Returns:
447
        quaternion: Optimal alignment rotation
448
        
449
    Notes:
450
        Minimizes ∫|Rd*Rb - Ra|² dt by computing mean of Ra/Rb.
451
        Equivalent to mean_rotor_in_chordal_metric(Ra / Rb, t).
452
    """
453

454
def optimal_alignment_in_Euclidean_metric(avec, bvec, t=None):
455
    """
456
    Solve Wahba's problem: find R such that R*b⃗*R̄ ≈ a⃗ in Euclidean metric.
457
    
458
    Args:
459
        avec (array_like): Target 3D vectors with shape (..., 3)
460
        bvec (array_like): Base 3D vectors with shape (..., 3)
461
        t (array_like, optional): Time weights for integral
462
        
463
    Returns:
464
        quaternion: Optimal rotation using Davenport's method
465
        
466
    Notes:
467
        Uses Davenport's eigenvector method for optimal vector alignment.
468
        Constructs matrix from vector cross-correlations and finds dominant eigenvector.
469
    """
470
```
471

472
## Types
473

474
```python { .api }
475
# Core quaternion type (implemented in C extension)
476
class quaternion:
477
    """
478
    Quaternion number with w, x, y, z components representing q = w + xi + yj + zk.
479
    
480
    Attributes:
481
        w (float): Scalar component (real part)
482
        x (float): i component (first imaginary part)
483
        y (float): j component (second imaginary part)
484
        z (float): k component (third imaginary part)
485
        
486
    Notes:
487
        Components are stored as double-precision floating point.
488
        Supports all standard arithmetic operations and mathematical functions.
489
        Can be used in NumPy arrays with dtype=quaternion.
490
    """
491
    def __init__(self, w=0.0, x=0.0, y=0.0, z=0.0):
492
        """Initialize quaternion with given components."""
493
        
494
    def conjugate(self):
495
        """Return quaternion conjugate: q* = w - xi - yj - zk."""
496
        
497
    def norm(self):
498
        """Return quaternion norm (magnitude): |q| = sqrt(w² + x² + y² + z²)."""
499
        
500
    def normalized(self):
501
        """Return unit quaternion: q/|q|. Raises error if norm is zero."""
502
        
503
    def __add__(self, other):
504
        """Add quaternions or scalars component-wise."""
505
        
506
    def __mul__(self, other):
507
        """Multiply quaternions using Hamilton product or scale by scalar."""
508
        
509
    def __truediv__(self, other):
510
        """Divide by quaternion (multiply by inverse) or scalar."""
511
        
512
    def __pow__(self, other):
513
        """Raise quaternion to scalar power using q^p = exp(p * log(q))."""
514
        
515
    def exp(self):
516
        """Exponential: exp(q) = exp(w) * (cos(|v|) + sin(|v|)*v/|v|) where v=(x,y,z)."""
517
        
518
    def log(self):
519
        """Natural logarithm: log(q) = ln(|q|) + arccos(w/|q|) * v/|v| where v=(x,y,z)."""
520
```