Skip to content

Rotations, Quaternions and Geometry

Conversions between Euler angles, rotation matrices (real and complex), and quaternions (real and complex), quaternion arithmetic, and 2D/3D geometry primitives.

Euler / matrix / quaternion conversions

v_roteu2qr

V_ROTEU2QR - convert Euler angles to real unit quaternion.

v_roteu2qr 

v_roteu2qr(m, e=None) -> ndarray

Convert Euler angles to real unit quaternion.

Parameters:

Name Type Description Default
m str

Rotation code string.

 required
e array_like

Euler angles. Shape (K,) or (K, N) where K is the number of xyz rotations in m.

 None

Returns:

Name Type Description
q (ndarray, shape(4) or (4, N))

Quaternion(s) [w, x, y, z].
Source code in pyvoicebox/v_roteu2qr.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def v_roteu2qr(m, e=None) -> np.ndarray:
    """Convert Euler angles to real unit quaternion.

    Parameters
    ----------
    m : str
        Rotation code string.
    e : array_like, optional
        Euler angles. Shape (K,) or (K, N) where K is the number of
        xyz rotations in m.

    Returns
    -------
    q : ndarray, shape (4,) or (4, N)
        Quaternion(s) [w, x, y, z].
    """
    mv = v_roteucode(m)
    ne = int(mv[1, -1])  # number of euler angles needed

    if ne == 0:
        q = np.zeros((4, 1))
        q[0, :] = 1.0
        r = q.copy()
    else:
        e = np.asarray(e, dtype=float)
        original_shape = e.shape
        if e.ndim == 1:
            e = e.reshape(-1, 1)
        elif e.ndim == 2 and original_shape[0] == 1 and e.size == ne:
            e = e.reshape(-1, 1)
        else:
            e = e.reshape(original_shape[0], -1)
        q = np.zeros((4, e.shape[1]))
        q[0, :] = 1.0
        r = q.copy()

    ef = 0.5 * mv[3, -1]  # signed euler angle scale factor (includes 0.5)
    nm = mv.shape[1] - 1  # number of rotations

    for i in range(nm):
        mvi = mv[:, i]
        mi = int(mvi[0])  # 1-indexed rotation code
        x = _Y[mi - 1, :] - 1  # convert to 0-indexed

        if mi <= 3:
            ei = int(mvi[1]) - 1  # 0-indexed euler angle index
            b = ef * e[ei, :]  # semi-rotation angle in radians
            c = np.cos(b)
            s = np.sin(b)
        else:
            c = _CB[mi - 1]
            s = _SB[mi - 1]

        r[x[0:2], :] = q[x[2:4], :]
        r[x[4:6], :] = -q[x[6:8], :]
        if np.isscalar(c):
            q = c * q + s * r
        else:
            q = c[np.newaxis, :] * q + s[np.newaxis, :] * r

    q[0, :] *= mv[6, -1]  # invert rotation if necessary

    q = _qr_normalize(q)

    # Reshape output
    if ne == 0:
        return q.ravel()
    elif e.shape[1] == 1:
        return q.ravel()
    else:
        return q

v_rotqr2eu

V_ROTQR2EU - convert real quaternion to Euler angles.

v_rotqr2eu 

v_rotqr2eu(m, q) -> ndarray

Convert real quaternion to Euler angles.

Parameters:

Name Type Description Default
m str

Rotation code string.

 required
q (array_like, shape(4) or (4, N))

Quaternion(s) [w, x, y, z].

 required

Returns:

Name Type Description
e (ndarray, shape(K) or (K, N))

Euler angles.
Source code in pyvoicebox/v_rotqr2eu.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def v_rotqr2eu(m, q) -> np.ndarray:
    """Convert real quaternion to Euler angles.

    Parameters
    ----------
    m : str
        Rotation code string.
    q : array_like, shape (4,) or (4, N)
        Quaternion(s) [w, x, y, z].

    Returns
    -------
    e : ndarray, shape (K,) or (K, N)
        Euler angles.
    """
    return v_rotro2eu(m, v_rotqr2ro(q))

v_roteu2ro

V_ROTEU2RO - convert Euler angles to rotation matrix.

v_roteu2ro 

v_roteu2ro(m, e=None) -> ndarray

Convert Euler angles to rotation matrix.

Parameters:

Name Type Description Default
m str

Rotation code string.

 required
e array_like

Euler angles.

 None

Returns:

Name Type Description
r (ndarray, shape(3, 3) or (3, 3, N))

Rotation matrix/matrices.
Source code in pyvoicebox/v_roteu2ro.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def v_roteu2ro(m, e=None) -> np.ndarray:
    """Convert Euler angles to rotation matrix.

    Parameters
    ----------
    m : str
        Rotation code string.
    e : array_like, optional
        Euler angles.

    Returns
    -------
    r : ndarray, shape (3, 3) or (3, 3, N)
        Rotation matrix/matrices.
    """
    if e is None:
        return v_rotqr2ro(v_roteu2qr(m))
    else:
        return v_rotqr2ro(v_roteu2qr(m, e))

v_rotro2eu

V_ROTRO2EU - convert rotation matrix to Euler angles.

v_rotro2eu 

v_rotro2eu(m, r) -> ndarray

Convert rotation matrix to Euler angles.

Parameters:

Name Type Description Default
m str

Rotation code string.

 required
r (array_like, shape(3, 3) or (3, 3, N))

Rotation matrix/matrices.

 required

Returns:

Name Type Description
e (ndarray, shape(K) or (K, N))

Euler angles.
Source code in pyvoicebox/v_rotro2eu.py 
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def v_rotro2eu(m, r) -> np.ndarray:
    """Convert rotation matrix to Euler angles.

    Parameters
    ----------
    m : str
        Rotation code string.
    r : array_like, shape (3, 3) or (3, 3, N)
        Rotation matrix/matrices.

    Returns
    -------
    e : ndarray, shape (K,) or (K, N)
        Euler angles.
    """
    mv = v_roteucode(m)
    nm = mv.shape[1] - 1  # number of rotation codes

    r = np.asarray(r, dtype=float)
    original_shape = r.shape
    squeeze = (r.ndim == 2)

    # Vectorize rotation matrices to (9, N), column-major (copy to avoid modifying input)
    if r.ndim == 2:
        rv = r.ravel(order='F').copy().reshape(9, 1)
    else:
        n = int(np.prod(original_shape[2:]))
        rv = np.zeros((9, n))
        for i in range(n):
            rv[:, i] = r[:, :, i].ravel(order='F')

    nr = rv.shape[1]

    # Transpose rotation matrix if needed
    if mv[6, -1] < 0:
        rv = rv[_RTR, :]

    ne = int(mv[1, -1])  # number of euler angles
    e = np.zeros((ne, nr))
    ef = mv[3, -1]  # scale factor

    # Process rotations in reverse order
    for i in range(nm - 1, -1, -1):
        mvi = mv[:, i]
        mi = int(mvi[0])  # rotation code (1-indexed)

        if mi <= 3:  # rotation around x, y, or z
            if mvi[5] != 0:  # skip if redundant
                # mvi[3] and mvi[4] are 1-indexed into vectorized matrix
                idx4 = int(mvi[3]) - 1  # 0-indexed
                idx5 = int(mvi[4]) - 1  # 0-indexed
                sign7 = mvi[6]

                si, ci, ri, ti = v_atan2sc(sign7 * rv[idx4, :], sign7 * rv[idx5, :])

                ei = int(mv[1, i]) - 1  # 0-indexed euler angle index
                e[ei, :] = ti * mvi[5] / ef

                # si_arr: shape (2, nr), row 0 = +mvi[5]*si, row 1 = -mvi[5]*si
                si_arr = mvi[5] * _PMW[:, np.newaxis] * si[np.newaxis, :]

                # Apply reverse rotation
                ax = mi - 1  # 0-indexed axis
                ci_rep = np.tile(ci, (6, 1))
                rv[_RTCI[:, ax], :] = ci_rep * rv[_RTCI[:, ax], :] - si_arr[_X6, :] * rv[_RTSI[:, ax], :]
        else:
            # Fixed rotation (code 4-12)
            ai = int(_SCAI[3, mi - 1]) - 1  # 0-indexed axis

            # scai values for this rotation
            cos_vals = np.tile(_SCAI[2, mi - 1], (6, 1))  # cos values
            sin_vals = _SCAI[_X6, mi - 1][:, np.newaxis]  # sin values (indexed by x6)

            rv[_RTCI[:, ai], :] = cos_vals * rv[_RTCI[:, ai], :] - sin_vals * rv[_RTSI[:, ai], :]

    if squeeze:
        return e.ravel()
    else:
        return e

v_rotro2qr

V_ROTRO2QR - convert 3x3 rotation matrix to real quaternion.

v_rotro2qr 

v_rotro2qr(r) -> ndarray

Convert 3x3 rotation matrix to real quaternion.

Parameters:

Name Type Description Default
r (array_like, shape(3, 3) or (3, 3, N))

Rotation matrix/matrices.

 required

Returns:

Name Type Description
q (ndarray, shape(4) or (4, N))

Quaternion(s) [w, x, y, z].
Source code in pyvoicebox/v_rotro2qr.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def v_rotro2qr(r) -> np.ndarray:
    """Convert 3x3 rotation matrix to real quaternion.

    Parameters
    ----------
    r : array_like, shape (3, 3) or (3, 3, N)
        Rotation matrix/matrices.

    Returns
    -------
    q : ndarray, shape (4,) or (4, N)
        Quaternion(s) [w, x, y, z].
    """
    r = np.asarray(r, dtype=float)
    original_shape = r.shape
    squeeze = (r.ndim == 2)

    # Reshape to (9, N) - MATLAB column-major vectorization of 3x3
    if r.ndim == 2:
        # Single 3x3 matrix: vectorize column-major
        rv = r.ravel(order='F').reshape(9, 1)
    else:
        n = int(np.prod(original_shape[2:]))
        rv = np.zeros((9, n))
        for i in range(n):
            rv[:, i] = r[:, :, i].ravel(order='F')

    nq = rv.shape[1]
    q = np.zeros((4, nq))

    # d = 1 + r(1,1) + r(2,2) + r(3,3) = 4*cos^2(t/2)
    # In column-major vectorized form: r[0]=r(1,1), r[4]=r(2,2), r[8]=r(3,3)
    d = 1.0 + rv[0, :] + rv[4, :] + rv[8, :]

    mm = d > 1  # rotation angle < 120 degrees
    if np.any(mm):
        s = np.sqrt(d[mm]) * 2.0
        # r(2,3)-r(3,2) = rv[5]-rv[7], r(3,1)-r(1,3) = rv[6]-rv[2], r(1,2)-r(2,1) = rv[1]-rv[3]
        q[1, mm] = (rv[5, mm] - rv[7, mm]) / s
        q[2, mm] = (rv[6, mm] - rv[2, mm]) / s
        q[3, mm] = (rv[1, mm] - rv[3, mm]) / s
        q[0, mm] = 0.25 * s

    if np.any(~mm):
        # r(1,1) > r(2,2) and r(1,1) > r(3,3)
        m = (rv[0, :] > rv[4, :]) & (rv[0, :] > rv[8, :]) & ~mm
        if np.any(m):
            s = np.sqrt(1.0 + rv[0, m] - rv[4, m] - rv[8, m]) * 2.0
            q[1, m] = 0.25 * s
            q[2, m] = (rv[1, m] + rv[3, m]) / s
            q[3, m] = (rv[6, m] + rv[2, m]) / s
            q[0, m] = (rv[5, m] - rv[7, m]) / s
            mm = mm | m

        # r(2,2) > r(3,3)
        m = (rv[4, :] > rv[8, :]) & ~mm
        if np.any(m):
            s = np.sqrt(1.0 + rv[4, m] - rv[0, m] - rv[8, m]) * 2.0
            q[1, m] = (rv[1, m] + rv[3, m]) / s
            q[2, m] = 0.25 * s
            q[3, m] = (rv[5, m] + rv[7, m]) / s
            q[0, m] = (rv[6, m] - rv[2, m]) / s
            mm = mm | m

        # Remaining
        m = ~mm
        if np.any(m):
            s = np.sqrt(1.0 + rv[8, m] - rv[0, m] - rv[4, m]) * 2.0
            q[1, m] = (rv[6, m] + rv[2, m]) / s
            q[2, m] = (rv[5, m] + rv[7, m]) / s
            q[3, m] = 0.25 * s
            q[0, m] = (rv[1, m] - rv[3, m]) / s

    # Check normalization
    norm_err = np.max(np.abs(np.sum(q ** 2, axis=0) - 1.0))
    if norm_err > 1e-8:
        raise ValueError('Input to v_rotro2qr must be a rotation matrix with det(r)=+1')

    q = _qr_normalize(q)

    if squeeze:
        return q.ravel()
    else:
        return q

v_rotqr2ro

V_ROTQR2RO - convert real quaternion to 3x3 rotation matrix.

v_rotqr2ro 

v_rotqr2ro(q) -> ndarray

Convert real quaternion to 3x3 rotation matrix.

Parameters:

Name Type Description Default
q (array_like, shape(4) or (4, N))

Real quaternion(s) [w, x, y, z].

 required

Returns:

Name Type Description
r (ndarray, shape(3, 3) or (3, 3, N))

Rotation matrix/matrices.
Source code in pyvoicebox/v_rotqr2ro.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def v_rotqr2ro(q) -> np.ndarray:
    """Convert real quaternion to 3x3 rotation matrix.

    Parameters
    ----------
    q : array_like, shape (4,) or (4, N)
        Real quaternion(s) [w, x, y, z].

    Returns
    -------
    r : ndarray, shape (3, 3) or (3, 3, N)
        Rotation matrix/matrices.
    """
    q = np.asarray(q, dtype=float)
    if q.ndim == 1:
        q = q.reshape(4, 1)
        squeeze = True
    else:
        squeeze = False

    nq = q.shape[1]

    # Compute quadratic terms: p = 2*q[h]*q[m] / sum(q^2)
    qnorm = np.sum(q ** 2, axis=0, keepdims=True)  # (1, nq)
    p = 2.0 * q[_H, :] * q[_M, :] / np.tile(qnorm, (16, 1))  # (16, nq)

    r = np.zeros((9, nq))
    r[_A, :] = 1.0 - p[_B, :] - p[_C, :]
    r[_D, :] = p[_E, :] - p[_F, :]
    r[_G, :] = p[_E, :] + p[_F, :]

    if squeeze:
        return r[:, 0].reshape(3, 3, order='F')
    else:
        result = np.zeros((3, 3, nq))
        for i in range(nq):
            result[:, :, i] = r[:, i].reshape(3, 3, order='F')
        return result

v_rotmr2qr

V_ROTMR2QR - convert real quaternion matrices to quaternion vectors.

v_rotmr2qr 

v_rotmr2qr(mr) -> ndarray

Convert real quaternion matrix form to real quaternion vector form.

Parameters:

Name Type Description Default
mr array_like, shape (4m, 4n, ...)

Real quaternion matrices.

 required

Returns:

Name Type Description
qr ndarray, shape (4m,) or (4m, n, ...)

Real quaternion vectors.
Source code in pyvoicebox/v_rotmr2qr.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def v_rotmr2qr(mr) -> np.ndarray:
    """Convert real quaternion matrix form to real quaternion vector form.

    Parameters
    ----------
    mr : array_like, shape (4m, 4n, ...)
        Real quaternion matrices.

    Returns
    -------
    qr : ndarray, shape (4m,) or (4m, n, ...)
        Real quaternion vectors.
    """
    mr = np.asarray(mr, dtype=float)
    s = list(mr.shape)
    n_out = s[1] // 4
    mr2 = mr.reshape(s[0], -1)
    qr = mr2[:, 0::4]
    if n_out == 1:
        return qr.ravel()
    else:
        out_s = list(s)
        out_s[1] = n_out
        return qr.reshape(out_s)

v_rotqr2mr

V_ROTQR2MR - convert real quaternion vectors to quaternion matrices.

v_rotqr2mr 

v_rotqr2mr(qr) -> ndarray

Convert real quaternion vector form to real quaternion matrix form.

Each quaternion [w, x, y, z] becomes a 4x4 matrix: [[ w, -x, -y, -z], [ x, w, -z, y], [ y, z, w, -x], [ z, -y, x, w]]

Parameters:

Name Type Description Default
qr array_like, shape (4m,) or (4m, n, ...)

Real quaternion vectors.

 required

Returns:

Name Type Description
mr ndarray, shape (4m, 4) or (4m, 4n, ...)

Real quaternion matrices.
Source code in pyvoicebox/v_rotqr2mr.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def v_rotqr2mr(qr) -> np.ndarray:
    """Convert real quaternion vector form to real quaternion matrix form.

    Each quaternion [w, x, y, z] becomes a 4x4 matrix:
    [[ w, -x, -y, -z],
     [ x,  w, -z,  y],
     [ y,  z,  w, -x],
     [ z, -y,  x,  w]]

    Parameters
    ----------
    qr : array_like, shape (4m,) or (4m, n, ...)
        Real quaternion vectors.

    Returns
    -------
    mr : ndarray, shape (4m, 4) or (4m, 4n, ...)
        Real quaternion matrices.
    """
    qr = np.asarray(qr, dtype=float)
    squeeze = (qr.ndim == 1)
    if squeeze:
        qr = qr.reshape(-1, 1)

    s = list(qr.shape)
    m = s[0]
    qa = qr.reshape(m, -1)
    n_total = qa.shape[1]
    n_quats = m // 4

    mr = np.zeros((m, 4 * n_total))

    for qi in range(n_quats):
        for col_idx in range(n_total):
            w = qa[4 * qi + 0, col_idx]
            x = qa[4 * qi + 1, col_idx]
            y = qa[4 * qi + 2, col_idx]
            z = qa[4 * qi + 3, col_idx]

            row_base = 4 * qi
            col_base = 4 * col_idx

            mr[row_base:row_base + 4, col_base:col_base + 4] = np.array([
                [w, -x, -y, -z],
                [x, w, -z, y],
                [y, z, w, -x],
                [z, -y, x, w],
            ])

    if squeeze:
        return mr
    else:
        out_s = list(s)
        out_s[1] = 4 * s[1]
        return mr.reshape(out_s)

v_rotmc2qc

V_ROTMC2QC - convert complex quaternion matrices to complex quaternion vectors.

v_rotmc2qc 

v_rotmc2qc(mc) -> ndarray

Convert complex quaternion matrix form to complex quaternion vector form.

Parameters:

Name Type Description Default
mc array_like, shape (2m, 2n, ...)

Complex quaternion matrices.

 required

Returns:

Name Type Description
qc ndarray, shape (2m, n, ...) or (2m,) if n=1

Complex quaternion vectors.
Source code in pyvoicebox/v_rotmc2qc.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def v_rotmc2qc(mc) -> np.ndarray:
    """Convert complex quaternion matrix form to complex quaternion vector form.

    Parameters
    ----------
    mc : array_like, shape (2m, 2n, ...)
        Complex quaternion matrices.

    Returns
    -------
    qc : ndarray, shape (2m, n, ...) or (2m,) if n=1
        Complex quaternion vectors.
    """
    mc = np.asarray(mc, dtype=complex)
    s = list(mc.shape)
    s[1] = s[1] // 2
    mc2 = mc.reshape(s[0], -1)
    qc = mc2[:, 0::2]
    if s[1] == 1:
        return qc.ravel()
    return qc.reshape(s)

v_rotqc2mc

V_ROTQC2MC - convert complex quaternion vectors to complex quaternion matrices.

v_rotqc2mc 

v_rotqc2mc(qc) -> ndarray

Convert complex quaternion vector form to complex quaternion matrix form.

Each quaternion r+ai+bj+ck is stored as a 2x1 complex vector [r+bi, a+ci]. The matrix form is [[r+bi, -a+ci], [a+ci, r-bi]].

Parameters:

Name Type Description Default
qc array_like, shape (2m,) or (2m, n, ...)

Complex quaternion vectors.

 required

Returns:

Name Type Description
mc ndarray, shape (2m, 2) or (2m, 2n, ...)

Complex quaternion matrices.
Source code in pyvoicebox/v_rotqc2mc.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def v_rotqc2mc(qc) -> np.ndarray:
    """Convert complex quaternion vector form to complex quaternion matrix form.

    Each quaternion r+ai+bj+ck is stored as a 2x1 complex vector [r+bi, a+ci].
    The matrix form is [[r+bi, -a+ci], [a+ci, r-bi]].

    Parameters
    ----------
    qc : array_like, shape (2m,) or (2m, n, ...)
        Complex quaternion vectors.

    Returns
    -------
    mc : ndarray, shape (2m, 2) or (2m, 2n, ...)
        Complex quaternion matrices.
    """
    qc = np.asarray(qc, dtype=complex)
    squeeze = (qc.ndim == 1)
    if squeeze:
        qc = qc.reshape(-1, 1)

    s = list(qc.shape)
    m = s[0]
    qa = qc.reshape(m, -1)
    n_cols = qa.shape[1]
    mc = np.zeros((m, 2 * n_cols), dtype=complex)

    ix = np.arange(0, m, 2)  # even rows (0-indexed)

    # Odd columns (0-indexed: 0, 2, 4, ...) = qa
    mc[:, 0::2] = qa
    # Even columns (1, 3, 5, ...)
    mc[ix, 1::2] = -np.conj(qa[ix + 1, :])
    mc[ix + 1, 1::2] = np.conj(qa[ix, :])

    if not squeeze:
        s[1] = 2 * s[1]
        mc = mc.reshape(s)
    return mc

v_rotqc2qr

V_ROTQC2QR - convert complex quaternion to real quaternion.

v_rotqc2qr 

v_rotqc2qr(qc) -> ndarray

Convert complex quaternion [r+jb, a+jc] to real [r, a, b, c].

Parameters:

Name Type Description Default
qc array_like, shape (2m, ...)

Complex quaternion vectors.

 required

Returns:

Name Type Description
qr ndarray, shape (4m, ...)

Real quaternion vectors.
Source code in pyvoicebox/v_rotqc2qr.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def v_rotqc2qr(qc) -> np.ndarray:
    """Convert complex quaternion [r+j*b, a+j*c] to real [r, a, b, c].

    Parameters
    ----------
    qc : array_like, shape (2m, ...)
        Complex quaternion vectors.

    Returns
    -------
    qr : ndarray, shape (4m, ...)
        Real quaternion vectors.
    """
    qc = np.asarray(qc, dtype=complex)
    s = list(qc.shape)
    s[0] = 2 * s[0]

    # Interleave real and imaginary parts: [real, imag] for each element
    flat = qc.ravel()
    pairs = np.array([flat.real, flat.imag])  # (2, total)
    qr = pairs.reshape(4, -1)  # (4, total/2)
    # Reorder from [r, b, a, c] to [r, a, b, c]
    a = np.array([0, 2, 1, 3])
    qr = qr[a, :]
    qr = qr.reshape(s)
    return qr

v_rotqr2qc

V_ROTQR2QC - convert real quaternion to complex quaternion.

v_rotqr2qc 

v_rotqr2qc(qr) -> ndarray

Convert real quaternion [r, a, b, c] to complex [r+jb, a+jc].

Parameters:

Name Type Description Default
qr array_like, shape (4m, ...)

Real quaternion vectors.

 required

Returns:

Name Type Description
qc ndarray, shape (2m, ...)

Complex quaternion vectors.
Source code in pyvoicebox/v_rotqr2qc.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def v_rotqr2qc(qr) -> np.ndarray:
    """Convert real quaternion [r, a, b, c] to complex [r+j*b, a+j*c].

    Parameters
    ----------
    qr : array_like, shape (4m, ...)
        Real quaternion vectors.

    Returns
    -------
    qc : ndarray, shape (2m, ...)
        Complex quaternion vectors.
    """
    qr = np.asarray(qr, dtype=float)
    s = list(qr.shape)
    qq = qr.reshape(4, -1)  # (4, n)
    # Reorder: [r, a, b, c] -> [r, b, a, c]
    a = np.array([0, 2, 1, 3])
    qq = qq[a, :]
    # Now [r, b, a, c]: pair (r, b) -> r + j*b, (a, c) -> a + j*c
    qc = qq[0::2, :] + 1j * qq[1::2, :]  # (2, n)
    s[0] = s[0] // 2
    qc = qc.reshape(s)
    return qc

v_rotax2qr

V_ROTAX2QR - convert rotation axis and angle to quaternion.

v_rotax2qr 

v_rotax2qr(a, t) -> ndarray

Convert rotation axis and angle to quaternion.

Parameters:

Name Type Description Default
a (array_like, shape(3))

Rotation axis vector (need not be unit length).

 required
t float

Rotation angle in radians.

 required

Returns:

Name Type Description
q (ndarray, shape(4))

Quaternion [w, x, y, z].
Source code in pyvoicebox/v_rotax2qr.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def v_rotax2qr(a, t) -> np.ndarray:
    """Convert rotation axis and angle to quaternion.

    Parameters
    ----------
    a : array_like, shape (3,)
        Rotation axis vector (need not be unit length).
    t : float
        Rotation angle in radians.

    Returns
    -------
    q : ndarray, shape (4,)
        Quaternion [w, x, y, z].
    """
    a = np.asarray(a, dtype=float).ravel()
    if np.all(a == 0):
        a = np.array([1.0, 0.0, 0.0])
    m = np.sqrt(a @ a)
    q = np.zeros(4)
    q[0] = np.cos(0.5 * t)
    q[1:] = np.sin(0.5 * t) * a / m
    return q

v_rotqr2ax

V_ROTQR2AX - convert quaternion to rotation axis and angle.

v_rotqr2ax 

v_rotqr2ax(q) -> tuple[ndarray, float]

Convert quaternion to rotation axis and angle.

Parameters:

Name Type Description Default
q (array_like, shape(4))

Quaternion [w, x, y, z].

 required

Returns:

Name Type Description
a (ndarray, shape(3))

Unit rotation axis vector.
t float

Rotation angle in radians (0 to 2*pi).
Source code in pyvoicebox/v_rotqr2ax.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def v_rotqr2ax(q) -> tuple[np.ndarray, float]:
    """Convert quaternion to rotation axis and angle.

    Parameters
    ----------
    q : array_like, shape (4,)
        Quaternion [w, x, y, z].

    Returns
    -------
    a : ndarray, shape (3,)
        Unit rotation axis vector.
    t : float
        Rotation angle in radians (0 to 2*pi).
    """
    q = np.asarray(q, dtype=float).ravel()
    a = q[1:4].copy()
    m = np.sqrt(a @ a)
    t = 2.0 * np.arctan2(m, q[0])
    if m > 0:
        a = a / m
    else:
        a = np.array([0.0, 0.0, 1.0])
    return a, t

v_rotpl2ro

V_ROTPL2RO - find rotation matrix from plane vectors.

v_rotpl2ro 

v_rotpl2ro(u, v, t=None) -> ndarray

Find rotation matrix to rotate in the plane containing u and v.

Parameters:

Name Type Description Default
u (array_like, shape(n))

First vector defining the rotation plane.

 required
v (array_like, shape(n))

Second vector defining the rotation plane.

 required
t float

Rotation angle in radians. If omitted, defaults to the angle between u and v.

 None

Returns:

Name Type Description
r (ndarray, shape(n, n))

Rotation matrix.
Source code in pyvoicebox/v_rotpl2ro.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def v_rotpl2ro(u, v, t=None) -> np.ndarray:
    """Find rotation matrix to rotate in the plane containing u and v.

    Parameters
    ----------
    u : array_like, shape (n,)
        First vector defining the rotation plane.
    v : array_like, shape (n,)
        Second vector defining the rotation plane.
    t : float, optional
        Rotation angle in radians. If omitted, defaults to the angle
        between u and v.

    Returns
    -------
    r : ndarray, shape (n, n)
        Rotation matrix.
    """
    u = np.asarray(u, dtype=float).ravel()
    n = len(u)
    v = np.asarray(v, dtype=float).ravel()

    l = np.sqrt(u @ u)
    if l == 0:
        raise ValueError('input u is a zero vector')
    u = u / l

    q = v - (v @ u) * u  # q is orthogonal to u
    l = np.sqrt(q @ q)
    if l == 0:
        # u and v are colinear or v=zero
        _, i = np.max(np.abs(u)), np.argmax(np.abs(u))
        q = np.zeros(n)
        q[(i + 1) % n] = 1.0
        q = q - (q @ u) * u
        l = np.sqrt(q @ q)
    q = q / l

    if t is None:
        s, c, _, _ = v_atan2sc(v @ q, v @ u)
        r = (np.eye(n) + (c - 1) * (np.outer(u, u) + np.outer(q, q))
             + s * (np.outer(q, u) - np.outer(u, q)))
    else:
        r = (np.eye(n) + (np.cos(t) - 1) * (np.outer(u, u) + np.outer(q, q))
             + np.sin(t) * (np.outer(q, u) - np.outer(u, q)))
    return r

v_rotro2pl

V_ROTRO2PL - find plane and rotation angle of a rotation matrix.

v_rotro2pl 

v_rotro2pl(r) -> tuple[ndarray, ndarray, float]

Find the plane and rotation angle of a rotation matrix.

Parameters:

Name Type Description Default
r (array_like, shape(n, n))

Rotation matrix.

 required

Returns:

Name Type Description
u (ndarray, shape(n))

First orthonormal vector defining the rotation plane.
v (ndarray, shape(n))

Second orthonormal vector or rotated u.
t float

Rotation angle in radians (0 <= t <= pi).
Source code in pyvoicebox/v_rotro2pl.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
def v_rotro2pl(r) -> tuple[np.ndarray, np.ndarray, float]:
    """Find the plane and rotation angle of a rotation matrix.

    Parameters
    ----------
    r : array_like, shape (n, n)
        Rotation matrix.

    Returns
    -------
    u : ndarray, shape (n,)
        First orthonormal vector defining the rotation plane.
    v : ndarray, shape (n,)
        Second orthonormal vector or rotated u.
    t : float
        Rotation angle in radians (0 <= t <= pi).
    """
    r = np.asarray(r, dtype=float)
    n = r.shape[0]

    # scipy.linalg.schur returns (T, Z) where A = Z @ T @ Z^H
    # MATLAB's [q, e] = schur(r) returns Z then T (i.e., q=Z, e=T)
    e, q = schur(r, output='real')

    # Find the 2x2 block with the largest subdiagonal element
    # e(2:n+1:n^2) in MATLAB = subdiagonal elements
    sub_diag = np.abs(e[np.arange(1, n), np.arange(n - 1)])
    i = np.argmax(sub_diag)  # 0-indexed position of largest subdiag

    z = int(e[i + 1, i] < 0)

    if z == 0:
        uv = q[:, i:i + 2]
    else:
        uv = q[:, i + 1:i - 1:-1] if i > 0 else q[:, i + 1::-1]

    u = uv[:, 0]
    v_vec = uv[:, 1]
    t_val = np.arctan2(np.abs(e[i + 1, i]), e[i, i])

    return u, v_vec, t_val

v_rotlu2ro

V_ROTLU2RO - convert look and up directions to rotation matrix.

v_rotlu2ro 

v_rotlu2ro(l, u=None) -> ndarray

Convert look and up directions to a rotation matrix.

The rotation maps the look direction to the negative z-axis and the up direction to lie in the y-z plane with a positive y component.

Parameters:

Name Type Description Default
l (array_like, shape(3) or (3, N))

Look direction vector(s).

 required
u (array_like, shape(3) or (3, N))

Up direction vector(s). Default is [0, 0, 1] unless l is a multiple of this, in which case [0, 1, 0].

 None

Returns:

Name Type Description
r (ndarray, shape(3, 3) or (3, 3, N))

Rotation matrix/matrices.
Source code in pyvoicebox/v_rotlu2ro.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def v_rotlu2ro(l, u=None) -> np.ndarray:
    """Convert look and up directions to a rotation matrix.

    The rotation maps the look direction to the negative z-axis and the
    up direction to lie in the y-z plane with a positive y component.

    Parameters
    ----------
    l : array_like, shape (3,) or (3, N)
        Look direction vector(s).
    u : array_like, shape (3,) or (3, N), optional
        Up direction vector(s). Default is [0, 0, 1] unless l is a
        multiple of this, in which case [0, 1, 0].

    Returns
    -------
    r : ndarray, shape (3, 3) or (3, 3, N)
        Rotation matrix/matrices.
    """
    l = np.asarray(l, dtype=float)
    if l.ndim == 1:
        l = l.reshape(3, 1)
        squeeze = True
    else:
        squeeze = False

    n = l.shape[1]

    # mk: sign correction masks, shape (3, 3, 4)
    mk = np.zeros((3, 3, 4))
    mk[:, :, 0] = np.tile(np.array([-1, 1, -1]).reshape(3, 1), (1, 3))
    mk[:, :, 1] = np.tile(np.array([1, -1, -1]).reshape(3, 1), (1, 3))
    mk[:, :, 2] = np.tile(np.array([1, 1, 1]).reshape(3, 1), (1, 3))
    mk[:, :, 3] = np.tile(np.array([-1, -1, 1]).reshape(3, 1), (1, 3))

    if u is None:
        u = np.tile(np.array([0, 0, 1], dtype=float).reshape(3, 1), (1, n))
        # If l is along z-axis, use [0, 1, 0] instead
        along_z = (l[0, :] == 0) & (l[1, :] == 0)
        u[1, along_z] += 1.0

    u = np.asarray(u, dtype=float)
    if u.ndim == 1:
        u = u.reshape(3, 1)

    if squeeze:
        r = np.zeros((3, 3))
    else:
        r = np.zeros((3, 3, n))

    for i in range(n):
        li = l[:, i:i + 1]
        ui = u[:, i:i + 1]
        # QR decomposition
        q_mat, t_mat = np.linalg.qr(np.hstack([li, ui]))

        # rx = [cross(q2, q1), q2, q1]^T
        q1 = q_mat[:, 0]
        q2 = q_mat[:, 1]
        rx = np.vstack([np.cross(q2, q1), q2, q1])

        # Select sign correction based on orientation
        cond1 = int(rx[2, :] @ li[:, 0] < 0)  # rx(3,:)*l < 0
        cond2 = int(rx[1, :] @ ui[:, 0] < 0)  # rx(2,:)*u < 0
        mk_idx = 2 * cond1 + cond2

        ri = rx * mk[:, :, mk_idx]
        if squeeze:
            r = ri
        else:
            r[:, :, i] = ri

    return r

v_rotro2lu

V_ROTRO2LU - convert rotation matrix to look and up directions.

v_rotro2lu 

v_rotro2lu(r) -> tuple[ndarray, ndarray]

Convert rotation matrix to look and up directions.

The rotation maps the look direction to the negative z-axis and the up direction to lie in the y-z plane with a positive y component.

Parameters:

Name Type Description Default
r (array_like, shape(3, 3) or (3, 3, N))

Rotation matrix/matrices.

 required

Returns:

Name Type Description
l (ndarray, shape(3) or (3, N))

Look direction vector(s).
u (ndarray, shape(3) or (3, N))

Up direction vector(s).
Source code in pyvoicebox/v_rotro2lu.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def v_rotro2lu(r) -> tuple[np.ndarray, np.ndarray]:
    """Convert rotation matrix to look and up directions.

    The rotation maps the look direction to the negative z-axis and the
    up direction to lie in the y-z plane with a positive y component.

    Parameters
    ----------
    r : array_like, shape (3, 3) or (3, 3, N)
        Rotation matrix/matrices.

    Returns
    -------
    l : ndarray, shape (3,) or (3, N)
        Look direction vector(s).
    u : ndarray, shape (3,) or (3, N)
        Up direction vector(s).
    """
    r = np.asarray(r, dtype=float)
    original_shape = r.shape

    if r.ndim == 2:
        # Single rotation matrix
        rv = r.ravel(order='F')  # column-major vectorization
        l = -rv[[2, 5, 8]]
        u = rv[[1, 4, 7]]
    else:
        n = int(np.prod(original_shape[2:]))
        rv = np.zeros((9, n))
        for i in range(n):
            rv[:, i] = r[:, :, i].ravel(order='F')
        out_shape = [3] + list(original_shape[2:])
        l = -rv[[2, 5, 8], :].reshape(out_shape)
        u = rv[[1, 4, 7], :].reshape(out_shape)

    return l, u

v_roteucode

V_ROTEUCODE - decode Euler angle rotation code string.

v_roteucode 

v_roteucode(m) -> ndarray

Decode a string specifying a rotation axis sequence.

Parameters:

Name Type Description Default
m str

Rotation code string.

 required

Returns:

Name Type Description
mv (ndarray, shape(7, k))

Decoded rotation parameters where k-1 is the number of rotations. mv[0, j] = rotation code (1-indexed: 1,2,3 for x,y,z; 4-12 for fixed) mv[1, j] = index into euler angle array (1-indexed) mv[2, j] = rotation class (pattern ID, 1-indexed) mv[3, j] = index into vectorized matrix (1-indexed) mv[4, j] = index of other changing element mv[5, j] = sign of sine term mv[6, j] = sign of entry before rotation / inversion flag (last col) mv[3, -1] = scale factor for euler angles mv[6, -1] = -1 to invert rotation (transpose matrix) or +1
Source code in pyvoicebox/v_roteucode.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
def v_roteucode(m) -> np.ndarray:
    """Decode a string specifying a rotation axis sequence.

    Parameters
    ----------
    m : str
        Rotation code string.

    Returns
    -------
    mv : ndarray, shape (7, k)
        Decoded rotation parameters where k-1 is the number of rotations.
        mv[0, j] = rotation code (1-indexed: 1,2,3 for x,y,z; 4-12 for fixed)
        mv[1, j] = index into euler angle array (1-indexed)
        mv[2, j] = rotation class (pattern ID, 1-indexed)
        mv[3, j] = index into vectorized matrix (1-indexed)
        mv[4, j] = index of other changing element
        mv[5, j] = sign of sine term
        mv[6, j] = sign of entry before rotation / inversion flag (last col)
        mv[3, -1] = scale factor for euler angles
        mv[6, -1] = -1 to invert rotation (transpose matrix) or +1
    """
    # Convert characters to integers (m - 'w')
    mm_raw = np.array([ord(c) - ord('w') for c in m], dtype=int)

    # Find characters XYZ (uppercase) and convert to xyz
    mi_upper = (mm_raw >= -31) & (mm_raw <= -29)
    mm_raw[mi_upper] += 32

    # Find digits 1-9 and convert to codes 4-12
    mi_digit = (mm_raw >= -70) & (mm_raw <= -62)
    mm_raw[mi_digit] += 74

    # Separate control characters (<=0) from rotations (>0)
    mi_ctrl = mm_raw <= 0
    mc = mm_raw[mi_ctrl]
    mm = mm_raw[~mi_ctrl]

    # Process control characters
    ef = 1.0  # angle scale factor
    es = 1    # angle sign
    fl = 1    # no rotation matrix transposing

    for c in mc:
        if c == -5:     # 'r' = radians
            pass
        elif c == -19:  # 'd' = degrees
            ef = np.pi / 180.0
        elif c == -37:  # 'R' = negated radians
            ef = -1.0
        elif c == -51:  # 'D' = negated degrees
            ef = -np.pi / 180.0
        elif c == -8:   # 'o' = object-extrinsic
            pass
        elif c == -40:  # 'O' = object-intrinsic
            fl = -1
            es = -1
        elif c == -22:  # 'a' = axes-extrinsic
            fl = -1
        elif c == -54:  # 'A' = axes-intrinsic
            es = -1
        else:
            raise ValueError(f'Invalid character: {chr(c + ord("w"))}')

    ef = ef * es  # change sign of scale factor if necessary
    if es < 0:
        mm = MES[mm - 1]  # sign-reverse: interchange 4,5,6 with 10,11,12

    nm = len(mm)
    mv = np.zeros((7, nm + 1), dtype=float)
    mv[0, :nm] = mm
    mv[0, nm] = 0

    # Cumulative sum for euler angle indices
    euler_mask = np.concatenate([(mm <= 3).astype(int), [0]])
    mv[1, :] = np.cumsum(euler_mask)

    mv[2, 0] = 1  # initial pattern is identity matrix (pattern 1)
    for i in range(nm):
        mmi = int(mm[i])  # rotation code (1-indexed)
        mv[2, i + 1] = _TRMAP[int(mv[2, i]) - 1, mmi - 1]  # pattern ID after rotation
        if mmi <= 3:
            # zel is (4, 3, 52) with 1-indexed pattern in dim 2
            mv[3:7, i] = ZEL[:, mmi - 1, int(mv[2, i]) - 1]

    mv[6, nm] = fl
    mv[3, nm] = ef

    return mv

v_rotation

V_ROTATION - encode and decode rotation matrices.

v_rotation 

v_rotation(x, y=None, z=None) -> ndarray

Create or decompose a rotation matrix.

(1) r = v_rotation(x, y, angle): rotation in the plane of x and y (2) axis, angle = v_rotation®: decompose 3x3 rotation matrix (3) r = v_rotation(axis, angle): rotation about axis (3D only) (4) r = v_rotation(axis_times_angle): rotation about axis (3D only)

Parameters:

Name Type Description Default
x array_like

First vector, rotation matrix, or axis*angle vector.

 required
y array_like or float

Second vector or angle.

 None
z float

Rotation angle.

 None

Returns:

Type Description
Depends on usage mode; see above.
Source code in pyvoicebox/v_rotation.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def v_rotation(x, y=None, z=None) -> np.ndarray:
    """Create or decompose a rotation matrix.

    (1) r = v_rotation(x, y, angle): rotation in the plane of x and y
    (2) axis, angle = v_rotation(r): decompose 3x3 rotation matrix
    (3) r = v_rotation(axis, angle): rotation about axis (3D only)
    (4) r = v_rotation(axis_times_angle): rotation about axis (3D only)

    Parameters
    ----------
    x : array_like
        First vector, rotation matrix, or axis*angle vector.
    y : array_like or float, optional
        Second vector or angle.
    z : float, optional
        Rotation angle.

    Returns
    -------
    Depends on usage mode; see above.
    """
    x = np.asarray(x, dtype=float)

    if z is not None:
        # Mode (1): r = v_rotation(x, y, angle)
        x_vec = x.ravel()
        l = len(x_vec)
        x_vec = x_vec / np.sqrt(x_vec @ x_vec)
        y_vec = np.asarray(y, dtype=float).ravel()
        y_vec = y_vec - (y_vec @ x_vec) * x_vec
        y_vec = y_vec / np.sqrt(y_vec @ y_vec)
        angle = float(z)
        r = (np.eye(l) + (np.cos(angle) - 1) * (np.outer(x_vec, x_vec) + np.outer(y_vec, y_vec))
             + np.sin(angle) * (np.outer(y_vec, x_vec) - np.outer(x_vec, y_vec)))
        return r

    elif x.ndim == 1 or (x.ndim == 2 and min(x.shape) == 1):
        # x is a vector
        x_vec = x.ravel()
        l = len(x_vec)
        if l == 3:
            # Mode (3) or (4): rotation about axis
            lx = np.sqrt(x_vec @ x_vec)
            if y is None:
                angle = lx
            else:
                angle = float(y)
            x_vec = x_vec / lx
            xx = np.outer(x_vec, x_vec)
            # Build skew-symmetric matrix
            # MATLAB: s([6 7 2])=x means s(row2,col1)=x(0), s(row0,col2)=x(1), s(row1,col0)=x(2)
            # in 0-based column-major flat indexing: 5->r(2,1), 6->r(0,2), 1->r(1,0)
            s = np.zeros((3, 3))
            s[2, 1] = x_vec[0]
            s[0, 2] = x_vec[1]
            s[1, 0] = x_vec[2]
            s[1, 2] = -x_vec[0]
            s[2, 0] = -x_vec[1]
            s[0, 1] = -x_vec[2]
            r = xx + np.cos(angle) * (np.eye(3) - xx) + np.sin(angle) * s
            return r
        else:
            raise ValueError("For non-3D vectors, use the 3-argument form")
    else:
        # x is a matrix: decomposition mode
        r_mat = x
        n = r_mat.shape[0]

        # Use eigendecomposition
        from scipy.linalg import schur
        # scipy.linalg.schur returns (T, Z) where A = Z @ T @ Z^H
        e, q = schur(r_mat, output='complex')
        d = np.diag(e)
        j = np.argsort(np.real(d))
        j1 = j[0]
        ang = np.angle(d[j1])

        sq = np.sqrt(2)
        r_vec = np.imag(q[:, j1]) * sq
        if np.all(np.abs(r_vec) < 1e-10):
            p_vec = q[:, j1].real
            r_vec = q[:, j[1]].real
        else:
            p_vec = np.real(q[:, j1]) * sq

        if n == 3:
            axis = np.cross(r_vec, p_vec)
            return axis, ang
        else:
            return r_vec, p_vec, ang

Quaternion arithmetic and mean

v_rotqrmean

V_ROTQRMEAN - calculate mean rotation of quaternion array.

v_rotqrmean 

v_rotqrmean(q) -> tuple[ndarray, ndarray, float]

Calculate the mean rotation of a quaternion array.

Parameters:

Name Type Description Default
q (array_like, shape(4, n))

Normalized real quaternion array.

 required

Returns:

Name Type Description
y (ndarray, shape(4))

Normalized mean quaternion.
s (ndarray, shape(n))

Sign vector such that y = normalize(q @ s).
v float

Average squared deviation from the mean quaternion.
Source code in pyvoicebox/v_rotqrmean.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def v_rotqrmean(q) -> tuple[np.ndarray, np.ndarray, float]:
    """Calculate the mean rotation of a quaternion array.

    Parameters
    ----------
    q : array_like, shape (4, n)
        Normalized real quaternion array.

    Returns
    -------
    y : ndarray, shape (4,)
        Normalized mean quaternion.
    s : ndarray, shape (n,)
        Sign vector such that y = normalize(q @ s).
    v : float
        Average squared deviation from the mean quaternion.
    """
    q = np.asarray(q, dtype=float)
    if q.ndim == 1:
        q = q.reshape(4, 1)

    mmax = 10  # number of n-best hypotheses to keep
    nq = q.shape[1]

    mkx = np.zeros((nq, mmax), dtype=int)
    msum = np.zeros((4, 2 * mmax))
    msum[:, 0] = q[:, 0]

    ix = np.arange(mmax)
    jx = np.arange(mmax, 2 * mmax)

    for i in range(1, nq):
        qi = q[:, i:i + 1]  # (4, 1)
        msum[:, jx] = msum[:, ix] - qi
        msum[:, ix] = msum[:, ix] + qi

        # Sort by squared norm (descending)
        norms = np.sum(msum ** 2, axis=0)
        kx = np.argsort(-norms)

        mkx[i, :] = kx[:mmax]
        msum[:, ix] = msum[:, kx[:mmax]]

    y = msum[:, 0]
    y = y / np.sqrt(y @ y)

    # Traceback to find signs
    s = np.zeros(nq)
    k = 0
    for i in range(nq - 1, 0, -1):
        neg = int(mkx[i, k] >= mmax)
        s[i] = neg
        k = mkx[i, k] - mmax * neg

    s = 1.0 - 2.0 * s

    v = float(np.sum(np.mean((q - np.outer(y, s)) ** 2, axis=1)))

    return y, s, v

v_rotqrvec

V_ROTQRVEC - rotate vectors by quaternion.

v_rotqrvec 

v_rotqrvec(q, x) -> ndarray

Apply a quaternion rotation to a vector array.

Parameters:

Name Type Description Default
q (array_like, shape(4))

Quaternion [w, x, y, z] (possibly unnormalized).

 required
x array_like

Array of 3D column vectors, shape (3n, ...).

 required

Returns:

Name Type Description
y ndarray

Rotated vectors, same shape as x.
Source code in pyvoicebox/v_rotqrvec.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def v_rotqrvec(q, x) -> np.ndarray:
    """Apply a quaternion rotation to a vector array.

    Parameters
    ----------
    q : array_like, shape (4,)
        Quaternion [w, x, y, z] (possibly unnormalized).
    x : array_like
        Array of 3D column vectors, shape (3n, ...).

    Returns
    -------
    y : ndarray
        Rotated vectors, same shape as x.
    """
    q = np.asarray(q, dtype=float)
    x = np.asarray(x, dtype=float)
    original_shape = x.shape

    r = v_rotqr2ro(q)  # (3, 3)
    y = r @ x.reshape(3, -1)
    return y.reshape(original_shape)

v_qrmult

V_QRMULT - Multiply two real quaternion matrices.

v_qrmult 

v_qrmult(q1, q2) -> ndarray

Multiply two real quaternion matrices.

Parameters:

Name Type Description Default
q1 array_like, shape (4m, n)

First quaternion matrix.

 required
q2 array_like, shape (4n, r)

Second quaternion matrix.

 required

Returns:

Name Type Description
q ndarray, shape (4m, r)

Matrix product.
Source code in pyvoicebox/v_qrmult.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def v_qrmult(q1, q2) -> np.ndarray:
    """Multiply two real quaternion matrices.

    Parameters
    ----------
    q1 : array_like, shape (4m, n)
        First quaternion matrix.
    q2 : array_like, shape (4n, r)
        Second quaternion matrix.

    Returns
    -------
    q : ndarray, shape (4m, r)
        Matrix product.
    """
    q1 = np.asarray(q1, dtype=float)
    q2 = np.asarray(q2, dtype=float)
    s1 = q1.shape
    s2 = q2.shape

    if s1 == (4,) or (len(s1) == 2 and s1 == (4, 1)):
        # q1 is a scalar quaternion
        q1c = q1.reshape(4, 1) if q1.ndim == 1 else q1
        q2c = q2 if q2.ndim == 2 else q2.reshape(4, 1)
        return v_qrdotmult(np.tile(q1c, (q2c.shape[0] // 4, q2c.shape[1])), q2c)
    elif s2 == (4,) or (len(s2) == 2 and s2 == (4, 1)):
        # q2 is a scalar quaternion
        q2c = q2.reshape(4, 1) if q2.ndim == 1 else q2
        q1c = q1 if q1.ndim == 2 else q1.reshape(4, 1)
        return v_qrdotmult(q1c, np.tile(q2c, (q1c.shape[0] // 4, q1c.shape[1])))
    else:
        return v_rotqr2mr(q1) @ q2

v_qrdivide

V_QRDIVIDE - Divide two real quaternions.

v_qrdivide 

v_qrdivide(q1, q2=None) -> ndarray

Divide two real quaternions: q = q1/q2 such that q1 = q*q2.

Parameters:

Name Type Description Default
q1 (array_like, shape(4))

First quaternion [r, i, j, k].

 required
q2 (array_like, shape(4))

Second quaternion. If omitted, returns inverse of q1.

 None

Returns:

Name Type Description
q (ndarray, shape(4))

Quotient quaternion.
Source code in pyvoicebox/v_qrdivide.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def v_qrdivide(q1, q2=None) -> np.ndarray:
    """Divide two real quaternions: q = q1/q2 such that q1 = q*q2.

    Parameters
    ----------
    q1 : array_like, shape (4,)
        First quaternion [r, i, j, k].
    q2 : array_like, shape (4,), optional
        Second quaternion. If omitted, returns inverse of q1.

    Returns
    -------
    q : ndarray, shape (4,)
        Quotient quaternion.
    """
    q1 = np.asarray(q1, dtype=float).ravel()

    if q2 is None:
        # Just invert q1
        q = q1 / (q1 @ q1)
        q[1:4] = -q[1:4]
        return q

    q2 = np.asarray(q2, dtype=float).ravel()

    # Invert q2
    qi = q2 / (q2 @ q2)
    qi[1:4] = -qi[1:4]

    # Multiply q1 * qi
    t = np.outer(q1, qi)
    s = np.zeros((4, 4))

    # MATLAB indices (1-based): a=[5 8 9 10 15 13], b=[6 7 11 12 14 16]
    # These are linearized indices into a 4x4 matrix (column-major)
    # a: (0,1),(3,1),(0,2),(1,2),(2,3),(0,3) -> negate t at b positions
    # b: (1,1),(2,1),(2,2),(3,2),(1,3),(3,3)

    # c=[1 2 3 4 6 7 11 12 16 14] (1-based, col-major) -> positive terms
    # d=[1 2 3 4 5 8 9 10 13 15] (1-based, col-major)

    # Direct quaternion multiplication: q1 * qi
    r1, i1, j1, k1 = q1
    r2, i2, j2, k2 = qi

    q = np.zeros(4)
    q[0] = r1 * r2 - i1 * i2 - j1 * j2 - k1 * k2
    q[1] = r1 * i2 + i1 * r2 + j1 * k2 - k1 * j2
    q[2] = r1 * j2 - i1 * k2 + j1 * r2 + k1 * i2
    q[3] = r1 * k2 + i1 * j2 - j1 * i2 + k1 * r2
    return q

v_qrdotmult

V_QRDOTMULT - Element-wise quaternion multiplication.

v_qrdotmult 

v_qrdotmult(q1, q2) -> ndarray

Multiply two real quaternion arrays element-wise (Hadamard product).

Parameters:

Name Type Description Default
q1 array_like, shape (4n, ...)

First quaternion array.

 required
q2 array_like, shape (4n, ...)

Second quaternion array.

 required

Returns:

Name Type Description
q ndarray, shape (4n, ...)

Element-wise quaternion product.
Source code in pyvoicebox/v_qrdotmult.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def v_qrdotmult(q1, q2) -> np.ndarray:
    """Multiply two real quaternion arrays element-wise (Hadamard product).

    Parameters
    ----------
    q1 : array_like, shape (4n, ...)
        First quaternion array.
    q2 : array_like, shape (4n, ...)
        Second quaternion array.

    Returns
    -------
    q : ndarray, shape (4n, ...)
        Element-wise quaternion product.
    """
    q1 = np.asarray(q1, dtype=float)
    q2 = np.asarray(q2, dtype=float)
    s = q1.shape
    # Use Fortran order to match MATLAB column-major reshape
    qa = q1.reshape(4, -1, order='F')
    qb = q2.reshape(4, -1, order='F')

    # Index arrays (0-based)
    a = np.array([0, 1, 2, 3, 1, 0, 3, 2, 2, 3, 0, 1, 3, 2, 1, 0])
    b = np.array([0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3])
    # Sign negation indices (0-based from MATLAB c=[2 3 4 7 12 14] -> 1,2,3,6,11,13)
    c = np.array([1, 2, 3, 6, 11, 13])

    prod = qa[a, :] * qb[b, :]
    prod[c, :] = -prod[c, :]
    # MATLAB: reshape(sum(reshape(q,4,[]),1),s) -- all column-major
    q = prod.reshape(4, -1, order='F').sum(axis=0).reshape(s, order='F')
    return q

v_qrdotdiv

V_QRDOTDIV - Element-wise quaternion division.

v_qrdotdiv 

v_qrdotdiv(x, y=None) -> ndarray

Divide two real quaternion arrays element-wise.

Parameters:

Name Type Description Default
x array_like, shape (4n, ...)

First quaternion array.

 required
y array_like, shape (4n, ...)

Second quaternion array. If omitted, returns inverse of x.

 None

Returns:

Name Type Description
q ndarray, shape (4n, ...)

Element-wise quaternion quotient.
Source code in pyvoicebox/v_qrdotdiv.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def v_qrdotdiv(x, y=None) -> np.ndarray:
    """Divide two real quaternion arrays element-wise.

    Parameters
    ----------
    x : array_like, shape (4n, ...)
        First quaternion array.
    y : array_like, shape (4n, ...), optional
        Second quaternion array. If omitted, returns inverse of x.

    Returns
    -------
    q : ndarray, shape (4n, ...)
        Element-wise quaternion quotient.
    """
    x = np.asarray(x, dtype=float)
    s = x.shape
    # Use Fortran order to match MATLAB column-major reshape
    xr = x.reshape(4, -1, order='F')

    if y is None:
        # Invert x
        m = np.sum(xr ** 2, axis=0)
        result = xr / m[np.newaxis, :]
        result[1:, :] = -result[1:, :]
        return result.reshape(s, order='F')

    y = np.asarray(y, dtype=float)
    yr = y.reshape(4, -1, order='F')
    m = np.sum(yr ** 2, axis=0)  # shape (N,)

    # Index arrays (0-based)
    a = np.array([0, 1, 2, 3, 1, 0, 3, 2, 2, 3, 0, 1, 3, 2, 1, 0])
    b = np.array([0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3])
    # Sign negation indices (0-based from MATLAB c=[6 8 10 11 15 16] -> 5,7,9,10,14,15)
    c = np.array([5, 7, 9, 10, 14, 15])

    prod = xr[a, :] * yr[b, :]
    prod[c, :] = -prod[c, :]
    # MATLAB: reshape(sum(reshape(q,4,[]),1),s)./m(ones(4,1),:)
    # After sum, get (4N,) reshaped to s. Divide by m replicated to match.
    q_4xN = prod.reshape(4, -1, order='F').sum(axis=0).reshape(4, -1, order='F')
    q_4xN = q_4xN / m[np.newaxis, :]
    q = q_4xN.reshape(s, order='F')
    return q

v_qrabs

V_QRABS - Absolute value and normalization of real quaternions.

v_qrabs 

v_qrabs(q1) -> tuple[ndarray, ndarray]

Absolute value and normalization of real quaternions.

Parameters:

Name Type Description Default
q1 array_like, shape (4n, ...)

Real quaternion array.

 required

Returns:

Name Type Description
m (ndarray, shape(n, ...))

Quaternion magnitudes.
q ndarray, shape (4n, ...)

Normalized quaternions (unit magnitude).
Source code in pyvoicebox/v_qrabs.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def v_qrabs(q1) -> tuple[np.ndarray, np.ndarray]:
    """Absolute value and normalization of real quaternions.

    Parameters
    ----------
    q1 : array_like, shape (4n, ...)
        Real quaternion array.

    Returns
    -------
    m : ndarray, shape (n, ...)
        Quaternion magnitudes.
    q : ndarray, shape (4n, ...)
        Normalized quaternions (unit magnitude).
    """
    q1 = np.asarray(q1, dtype=float)
    s = list(q1.shape)
    q = q1.reshape(4, -1).copy()
    m = np.sqrt(np.sum(q ** 2, axis=0))

    # Normalize non-zero quaternions
    nz = m > 0
    q[:, nz] = q[:, nz] / m[np.newaxis, nz]
    # Zero quaternions become [1, 0, 0, 0]
    q[0, ~nz] = 1.0
    q[1:, ~nz] = 0.0

    q = q.reshape(s)
    s[0] = s[0] // 4
    m = m.reshape(s)
    return m, q

v_qrpermute

V_QRPERMUTE - Transpose or permute a quaternion array.

v_qrpermute 

v_qrpermute(x, p=None) -> ndarray

Transpose or permute a quaternion array.

Parameters:

Name Type Description Default
x array_like, shape (4m, ...)

Real quaternion array.

 required
p array_like

New order of dimensions (0-based). Default transposes first two dims.

 None

Returns:

Name Type Description
y ndarray

Permuted quaternion array.
Source code in pyvoicebox/v_qrpermute.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def v_qrpermute(x, p=None) -> np.ndarray:
    """Transpose or permute a quaternion array.

    Parameters
    ----------
    x : array_like, shape (4m, ...)
        Real quaternion array.
    p : array_like, optional
        New order of dimensions (0-based). Default transposes first two dims.

    Returns
    -------
    y : ndarray
        Permuted quaternion array.
    """
    x = np.asarray(x, dtype=float)
    s = list(x.shape)
    if p is None:
        # Default: transpose first two dimensions
        if len(s) == 1:
            return x.copy()
        p = list(range(len(s)))
        p[0], p[1] = p[1], p[0]

    p = list(p)
    s[0] = s[0] // 4
    t = [s[i] for i in p]
    t[0] = 4 * t[0]

    # Reshape to (4, s0, s1, ...), permute, then reshape back
    new_shape = [4] + s
    perm = [0] + [i + 1 for i in p]
    y = np.transpose(x.reshape(new_shape), perm).reshape(t)
    return y

Geometry

v_polygonarea

V_POLYGONAREA - calculate polygon area.

v_polygonarea 

v_polygonarea(p) -> float

Calculate the area of a polygon using the shoelace formula.

Parameters:

Name Type Description Default
p (array_like, shape(n, 2))

Polygon vertices.

 required

Returns:

Name Type Description
a float

Signed area (positive if vertices are counter-clockwise).
Source code in pyvoicebox/v_polygonarea.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def v_polygonarea(p) -> float:
    """Calculate the area of a polygon using the shoelace formula.

    Parameters
    ----------
    p : array_like, shape (n, 2)
        Polygon vertices.

    Returns
    -------
    a : float
        Signed area (positive if vertices are counter-clockwise).
    """
    p = np.asarray(p, dtype=float)
    # Append first point
    p_closed = np.vstack([p, p[0:1, :]])
    a = 0.5 * np.sum(
        (p_closed[:-1, 0] - p_closed[1:, 0]) * (p_closed[:-1, 1] + p_closed[1:, 1])
    )
    return float(a)

v_polygonwind

V_POLYGONWIND - test if points are inside a polygon.

v_polygonwind 

v_polygonwind(p, x) -> ndarray

Calculate the winding number for points with respect to a polygon.

Parameters:

Name Type Description Default
p (array_like, shape(n, 2))

Polygon vertices.

 required
x (array_like, shape(m, 2))

Points to test.

 required

Returns:

Name Type Description
w (ndarray, shape(m))

Winding number for each point. For a CCW polygon, 0=outside, 1=inside.
Source code in pyvoicebox/v_polygonwind.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def v_polygonwind(p, x) -> np.ndarray:
    """Calculate the winding number for points with respect to a polygon.

    Parameters
    ----------
    p : array_like, shape (n, 2)
        Polygon vertices.
    x : array_like, shape (m, 2)
        Points to test.

    Returns
    -------
    w : ndarray, shape (m,)
        Winding number for each point. For a CCW polygon, 0=outside, 1=inside.
    """
    p = np.asarray(p, dtype=float)
    x = np.asarray(x, dtype=float)
    n = p.shape[0]
    m = x.shape[0]

    # Close the polygon
    q = np.zeros((2, n + 1))
    q[:, :n] = p.T
    q[:, n] = q[:, 0]

    # Indices for edges: i -> j
    ii = np.arange(n)
    jj = np.arange(1, n + 1)

    # Cross product term for each (point, edge) pair
    # q(1,i)*q(2,j) - q(2,i)*q(1,j) + x(:,1)*(q(2,i)-q(2,j)) + x(:,2)*(q(1,j)-q(1,i))
    cross_const = q[0, ii] * q[1, jj] - q[1, ii] * q[0, jj]  # (n,)
    dy = q[1, ii] - q[1, jj]  # (n,)
    dx = q[0, jj] - q[0, ii]  # (n,)

    # cross_val(m, n)
    cross_val = (cross_const[np.newaxis, :]
                 + x[:, 0:1] * dy[np.newaxis, :]
                 + x[:, 1:2] * dx[np.newaxis, :])

    sign_term = 2 * (cross_val > 0).astype(float) - 1  # (m, n)

    # Edge crossing test: y coordinates
    above_i = (q[1, ii][np.newaxis, :] > x[:, 1:2])  # (m, n)
    above_j = (q[1, jj][np.newaxis, :] > x[:, 1:2])  # (m, n)

    crossing = np.abs(above_j.astype(float) - above_i.astype(float))  # (m, n)

    w = np.sum(sign_term * crossing, axis=1) / 2.0
    return w.astype(int) if np.all(w == np.round(w)) else w

v_polygonxline

V_POLYGONXLINE - find where a line crosses a polygon.

v_polygonxline 

v_polygonxline(
    p, l
) -> tuple[ndarray, ndarray, ndarray, ndarray]

Find where a line crosses a polygon.

Parameters:

Name Type Description Default
p (array_like, shape(n, 2))

Polygon vertices.

 required
l (array_like, shape(3))

Line in the form l @ [x, y, 1] = 0.

 required

Returns:

Name Type Description
xc (ndarray, shape(k, 2))

Crossing point coordinates.
ec (ndarray, shape(k))

Edge indices (1-indexed) where crossings occur.
tc (ndarray, shape(k))

Parametric positions along the line.
xy0 (ndarray, shape(2))

Starting point of the parametric line.
Source code in pyvoicebox/v_polygonxline.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def v_polygonxline(p, l) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Find where a line crosses a polygon.

    Parameters
    ----------
    p : array_like, shape (n, 2)
        Polygon vertices.
    l : array_like, shape (3,)
        Line in the form l @ [x, y, 1] = 0.

    Returns
    -------
    xc : ndarray, shape (k, 2)
        Crossing point coordinates.
    ec : ndarray, shape (k,)
        Edge indices (1-indexed) where crossings occur.
    tc : ndarray, shape (k,)
        Parametric positions along the line.
    xy0 : ndarray, shape (2,)
        Starting point of the parametric line.
    """
    p = np.asarray(p, dtype=float)
    l = np.asarray(l, dtype=float).ravel()
    n = p.shape[0]

    # Close polygon
    q = np.ones((n + 1, 3))
    q[:n, :2] = p
    q[n, :] = q[0, :]

    # Signed distance of each vertex from the line
    cdist = q @ l  # (n+1,)
    cside = cdist > 0
    cside[n] = cside[0]

    # Find edges that cross the line (where side changes)
    ec = np.where(cside[1:n + 1] != cside[:n])[0]  # 0-indexed edge indices
    nc = len(ec)

    if nc == 0:
        return np.empty((0, 2)), np.empty(0, dtype=int), np.empty(0), np.empty(0)

    l2 = l[:2]
    l2m = l2 @ l2
    l3 = l[2]
    a = np.array([-l[1], l[0]])
    xy0 = -l3 / l2m * l2

    # Parametric position along the line
    tn = (q[:, :2] - xy0[np.newaxis, :]) @ a / l2m  # (n+1,)
    tc = (cdist[ec + 1] * tn[ec] - cdist[ec] * tn[ec + 1]) / (cdist[ec + 1] - cdist[ec])

    # Sort crossings
    idx = np.argsort(tc)
    tc = tc[idx]
    ec = ec[idx]

    # Remove duplicate pairs
    if len(tc) > 1:
        tm = tc[1:] == tc[:-1]
        tm_full = np.concatenate([[False], tm]) | np.concatenate([tm, [False]])
        tc = tc[~tm_full]
        ec = ec[~tm_full]

    nc = len(ec)
    xc = xy0[np.newaxis, :] + tc[:, np.newaxis] * a[np.newaxis, :]

    # Convert to 1-indexed edge indices (MATLAB convention)
    ec = ec + 1

    return xc, ec, tc, xy0

v_minspane

V_MINSPANE - minimum spanning tree using Euclidean distance.

v_minspane 

v_minspane(x) -> ndarray

Calculate minimum spanning tree using Euclidean distance.

Uses Delaunay triangulation to find candidate edges, then applies Kruskal's algorithm.

Parameters:

Name Type Description Default
x (array_like, shape(n, d))

d-dimensional data points.

 required

Returns:

Name Type Description
p (ndarray, shape(n - 1))

Parent node indices (1-indexed). Node n is the root.
s (ndarray, shape(n - 1))

Edge indices sorted by ascending Euclidean distance (1-indexed).
Source code in pyvoicebox/v_minspane.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
def v_minspane(x) -> np.ndarray:
    """Calculate minimum spanning tree using Euclidean distance.

    Uses Delaunay triangulation to find candidate edges, then applies
    Kruskal's algorithm.

    Parameters
    ----------
    x : array_like, shape (n, d)
        d-dimensional data points.

    Returns
    -------
    p : ndarray, shape (n-1,)
        Parent node indices (1-indexed). Node n is the root.
    s : ndarray, shape (n-1,)
        Edge indices sorted by ascending Euclidean distance (1-indexed).
    """
    x = np.asarray(x, dtype=float)
    np_pts, nd = x.shape

    # Delaunay triangulation
    tri = Delaunay(x)
    simplices = tri.simplices  # (nt, nd+1)

    # Extract all edges from simplices
    from itertools import combinations
    edges = set()
    for simplex in simplices:
        for i, j in combinations(sorted(simplex), 2):
            edges.add((i, j))

    ee = np.array(sorted(edges))  # (ne, 2), 0-indexed
    ne = ee.shape[0]

    # Compute edge lengths squared
    sz = np.sum((x[ee[:, 0], :] - x[ee[:, 1], :]) ** 2, axis=1)
    mz = np.argsort(sz)
    ee = ee[mz, :]  # sort by ascending length

    # Kruskal's algorithm with union-find
    parent = np.full(np_pts, -1, dtype=int)  # -1 means root
    rank = np.zeros(np_pts, dtype=int)

    def find(i):
        while parent[i] >= 0:
            parent[i] = parent[parent[i]] if parent[parent[i]] >= 0 else parent[i]  # path compression
            i = parent[i]
        return i

    ei = []
    for i in range(ne):
        r1 = find(ee[i, 0])
        r2 = find(ee[i, 1])
        if r1 != r2:
            ei.append(i)
            # Union by rank
            if rank[r1] > rank[r2]:
                parent[r2] = r1
            elif rank[r1] < rank[r2]:
                parent[r1] = r2
            else:
                parent[r2] = r1
                rank[r1] += 1
            if len(ei) == np_pts - 1:
                break

    mst_edges = ee[ei, :]  # (n-1, 2), 0-indexed

    # Build tree with node (np_pts-1) as root (0-indexed), then convert to 1-indexed
    # Use BFS from root
    from collections import defaultdict, deque
    adj = defaultdict(list)
    edge_map = {}
    for idx, (u, v) in enumerate(mst_edges):
        adj[u].append(v)
        adj[v].append(u)
        edge_map[(min(u, v), max(u, v))] = idx

    root = np_pts - 1  # 0-indexed root
    p = np.zeros(np_pts - 1, dtype=int)  # parent array (1-indexed output)
    visited = np.zeros(np_pts, dtype=bool)
    visited[root] = True
    queue = deque([root])

    edge_lengths = np.sqrt(sz[np.array(ei)])

    # Build parent relationships and edge lengths for sorting
    node_edge_len = np.zeros(np_pts - 1)

    while queue:
        node = queue.popleft()
        for neighbor in adj[node]:
            if not visited[neighbor]:
                visited[neighbor] = True
                p[neighbor if neighbor < np_pts - 1 else neighbor] = node + 1  # 1-indexed parent
                if neighbor < np_pts - 1:
                    p[neighbor] = node + 1
                    key = (min(node, neighbor), max(node, neighbor))
                    edge_idx = edge_map[key]
                    node_edge_len[neighbor] = sz[ei[edge_idx]]
                queue.append(neighbor)

    # Sort edges by ascending length
    # s gives nodes sorted by edge length to parent
    lengths = np.sqrt(np.sum((x[:np_pts - 1, :] - x[p - 1, :]) ** 2, axis=1))
    s = np.argsort(lengths) + 1  # 1-indexed

    return p, s

v_imagehomog

V_IMAGEHOMOG - apply homography transformation to an image.

v_imagehomog 

v_imagehomog(
    im, h=None, m="", clip=None
) -> tuple[ndarray, ndarray, ndarray]

Apply a homography transformation to an image with bilinear interpolation.

Parameters:

Name Type Description Default
im (ndarray, shape(ny, nx) or (ny, nx, nc))

Input image (uint8).

 required
h (ndarray, shape(3, 3))

Homography matrix. Default is identity.

 None
m str

Mode string: 'c' central coordinates, 'k' clip to original, 'x' extend, 't' trim blank rows/cols.

 ''
clip array_like

Clipping specification.

 None

Returns:

Name Type Description
ih ndarray

Transformed image (uint8).
xa ndarray

X axis coordinates.
ya ndarray

Y axis coordinates.
Source code in pyvoicebox/v_imagehomog.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
def v_imagehomog(im, h=None, m='', clip=None) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Apply a homography transformation to an image with bilinear interpolation.

    Parameters
    ----------
    im : ndarray, shape (ny, nx) or (ny, nx, nc)
        Input image (uint8).
    h : ndarray, shape (3, 3), optional
        Homography matrix. Default is identity.
    m : str, optional
        Mode string: 'c' central coordinates, 'k' clip to original, 'x' extend,
        't' trim blank rows/cols.
    clip : array_like, optional
        Clipping specification.

    Returns
    -------
    ih : ndarray
        Transformed image (uint8).
    xa : ndarray
        X axis coordinates.
    ya : ndarray
        Y axis coordinates.
    """
    maxby = int(1e7)

    if im.ndim == 2:
        ny, nx = im.shape
        nc = 1
        im = im[:, :, np.newaxis]
    else:
        ny, nx, nc = im.shape

    if clip is None:
        clip = np.array([2.0])
    else:
        clip = np.asarray(clip, dtype=float).ravel()
    if h is None:
        h = np.eye(3)
    else:
        h = np.asarray(h, dtype=float)

    imr = im.reshape(nx * ny, nc).astype(np.float32)
    t = np.eye(3)

    if 'c' in m:
        t[0, 2] = 0.5 + nx / 2.0
        t[1, 2] = 0.5 + ny / 2.0
        h = t @ h @ np.linalg.inv(t)
        if len(clip) == 4:
            clip = clip + np.array([t[0, 2], t[0, 2], t[1, 2], t[1, 2]])

    if 'k' in m:
        clip = np.array([1.0, 1.0])
    elif len(clip) == 1:
        clip = np.array([clip[0], clip[0]])

    if len(clip) == 2:
        clip = np.array([
            1 * clip[0] - (clip[0] - 1) * (1 + nx) / 2,
            nx * clip[0] - (clip[0] - 1) * (1 + nx) / 2,
            1 * clip[1] - (clip[1] - 1) * (1 + ny) / 2,
            ny * clip[1] - (clip[1] - 1) * (1 + ny) / 2,
        ])

    clip = clip.ravel()
    clip[0] = np.floor(clip[0])
    clip[2] = np.floor(clip[2])
    clip[1] = np.ceil(clip[1])
    clip[3] = np.ceil(clip[3])

    # Determine image of source corners
    bi = np.array([[1, 1, nx, nx], [1, ny, ny, 1], [1, 1, 1, 1]], dtype=float)
    box = h @ bi
    b3 = box[2, :]

    if np.any(b3 <= 0):
        ib = np.where(b3 > 0)[0]
        nb = len(ib)
        if nb == 0:
            raise ValueError('image invisible')
        bb = np.ones((3, nb + 2))
        bb[:, :nb] = box[:, ib]
        px = np.array([3, 0, 1, 2])
        cross_idx = np.where(b3 * b3[px] <= 0)[0]
        ip = px[cross_idx]
        af = b3[ip]
        bf = b3[cross_idx]
        pof = np.array([[-1, 0, 1, 0], [0, 1, 0, -1], [0, 0, 0, 0]], dtype=float)
        for k, ci in enumerate(cross_idx):
            frac = bf[k] / (bf[k] - af[k])
            pt = bi[:, ip[k]] * (1 - frac) + bi[:, ci] * frac
            bb[:, nb + k] = h @ pt
        box = bb

    box_xy = box[:2, :] / box[2:3, :]
    box_bounds = np.array([
        np.min(box_xy[0, :]), np.max(box_xy[0, :]),
        np.min(box_xy[1, :]), np.max(box_xy[1, :])
    ])

    box_bounds[0] = np.floor(max(clip[0], box_bounds[0]))
    box_bounds[2] = np.floor(max(clip[2], box_bounds[2]))
    box_bounds[1] = np.ceil(min(clip[1], box_bounds[1]))
    box_bounds[3] = np.ceil(min(clip[3], box_bounds[3]))

    g = np.linalg.inv(h)
    mx = int(box_bounds[1] - box_bounds[0] + 1)
    my = int(box_bounds[3] - box_bounds[2] + 1)

    if mx <= 0 or my <= 0:
        ih = np.zeros((max(my, 0), max(mx, 0), nc), dtype=np.uint8)
        xa = np.arange(box_bounds[0], box_bounds[1] + 1) - t[0, 2]
        ya = np.arange(box_bounds[2], box_bounds[3] + 1) - t[1, 2]
        return ih, xa, ya

    ih = np.zeros((my * mx, nc), dtype=np.uint8)

    # Process all pixels at once (simpler than chunked approach)
    yy, xx = np.meshgrid(
        np.arange(box_bounds[2], box_bounds[3] + 1),
        np.arange(box_bounds[0], box_bounds[1] + 1),
        indexing='ij'
    )
    coords = np.vstack([xx.ravel(), yy.ravel(), np.ones(mx * my)])
    src = g @ coords
    gn = src[:2, :] / src[2:3, :]

    # Mask valid pixels
    mn = (gn[0, :] > -0.5) & (gn[1, :] > -0.5) & (gn[0, :] < nx + 0.5) & (gn[1, :] < ny + 0.5)

    if np.any(mn):
        fn1 = np.clip(np.floor(gn[0, mn]).astype(int), 1, nx - 1)
        fn2 = np.clip(np.floor(gn[1, mn]).astype(int), 1, ny - 1)
        dn1 = np.clip(gn[0, mn] - fn1, 0, 1)
        dn2 = np.clip(gn[1, mn] - fn2, 0, 1)
        dn1c = 1 - dn1
        dn2c = 1 - dn2

        # MATLAB uses 1-based indexing: imr(fn2 + ny*(fn1-1), :)
        # Python 0-based: imr[(fn2-1) + ny*(fn1-1), :]
        idx00 = (fn2 - 1) + ny * (fn1 - 1)
        idx01 = fn2 + ny * (fn1 - 1)
        idx10 = (fn2 - 1) + ny * fn1
        idx11 = fn2 + ny * fn1

        # Clip indices
        max_idx = nx * ny - 1
        idx00 = np.clip(idx00, 0, max_idx)
        idx01 = np.clip(idx01, 0, max_idx)
        idx10 = np.clip(idx10, 0, max_idx)
        idx11 = np.clip(idx11, 0, max_idx)

        val = (dn1c[:, np.newaxis] * (dn2c[:, np.newaxis] * imr[idx00, :]
                                       + dn2[:, np.newaxis] * imr[idx01, :])
               + dn1[:, np.newaxis] * (dn2c[:, np.newaxis] * imr[idx10, :]
                                        + dn2[:, np.newaxis] * imr[idx11, :]))
        jx = np.where(mn)[0]
        ih[jx, :] = np.clip(val, 0, 255).astype(np.uint8)

    ih = ih.reshape(my, mx, nc)

    if 'x' in m:
        b0 = int(box_bounds[0])
        b2 = int(box_bounds[2])
        c0 = int(clip[0])
        c1 = int(clip[1])
        c2 = int(clip[2])
        c3 = int(clip[3])
        new_ih = np.zeros((c3 - c2 + 1, c1 - c0 + 1, nc), dtype=np.uint8)
        row_start = b2 - c2
        col_start = b0 - c0
        new_ih[row_start:row_start + my, col_start:col_start + mx, :] = ih
        ih = new_ih
        xa = np.arange(c0, c1 + 1) - t[0, 2]
        ya = np.arange(c2, c3 + 1) - t[1, 2]
    else:
        xa = np.arange(box_bounds[0], box_bounds[1] + 1) - t[0, 2]
        ya = np.arange(box_bounds[2], box_bounds[3] + 1) - t[1, 2]

    if 't' in m:
        ihs = np.sum(ih, axis=2)
        azx = np.all(ihs == 0, axis=0)
        ix1 = np.argmax(~azx)
        if np.all(azx):
            return np.array([], dtype=np.uint8), np.array([]), np.array([])
        ix2 = len(azx) - 1 - np.argmax(~azx[::-1])
        azy = np.all(ihs == 0, axis=1)
        iy1 = np.argmax(~azy)
        iy2 = len(azy) - 1 - np.argmax(~azy[::-1])
        ih = ih[iy1:iy2 + 1, ix1:ix2 + 1, :]
        xa = xa[ix1:ix2 + 1]
        ya = ya[iy1:iy2 + 1]

    if nc == 1:
        ih = ih[:, :, 0]

    return ih, xa, ya

v_rectifyhomog

V_RECTIFYHOMOG - apply rectifying homographies to an image set.

v_rectifyhomog 

v_rectifyhomog(
    ims, roc=None, k0=None, mode=""
) -> tuple[ndarray, ndarray, ndarray]

Apply rectifying homographies to an image set.

Parameters:

Name Type Description Default
ims list of ndarray

List of input images.

 required
roc (ndarray, shape(3, 3) or (3, 3, nc))

Rotation matrices from world to camera coordinates.

 None
k0 float or ndarray

Camera matrix or focal length.

 None
mode str

Mode string.

 ''

Returns:

Name Type Description
imr list of ndarray

Rectified images.
xa list of ndarray

X axis for each image.
ya list of ndarray

Y axis for each image.
Source code in pyvoicebox/v_rectifyhomog.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def v_rectifyhomog(ims, roc=None, k0=None, mode='') -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Apply rectifying homographies to an image set.

    Parameters
    ----------
    ims : list of ndarray
        List of input images.
    roc : ndarray, shape (3, 3) or (3, 3, nc)
        Rotation matrices from world to camera coordinates.
    k0 : float or ndarray
        Camera matrix or focal length.
    mode : str
        Mode string.

    Returns
    -------
    imr : list of ndarray
        Rectified images.
    xa : list of ndarray
        X axis for each image.
    ya : list of ndarray
        Y axis for each image.
    """
    if not isinstance(ims, list):
        ims = [ims]
    nc = len(ims)

    if roc is None:
        roc = np.tile(np.eye(3)[:, :, np.newaxis], (1, 1, nc))
    roc = np.asarray(roc, dtype=float)
    if roc.ndim == 2:
        roc = roc[:, :, np.newaxis]

    if k0 is None:
        k0 = 0.8
    if mode is None:
        mode = ''

    vv = 'v' in mode

    # Determine mean camera orientation
    if 'a' in mode:
        qrc = np.zeros((4, nc))
        for i in range(nc):
            qrc[:, i] = v_rotro2qr(roc[:, :, i])
        rocmean = v_rotqr2ro(v_rotqrmean(qrc)[0])
    else:
        rocmean = np.eye(3)

    modeh = 'kxt' if 'k' in mode else 't'

    imr = []
    xa = []
    ya = []
    for i in range(nc):
        imsz = ims[i].shape
        if np.isscalar(k0) or (isinstance(k0, np.ndarray) and k0.size < 9):
            fe = float(k0) if np.isscalar(k0) else float(k0.ravel()[0])
            if fe < 0.1 * imsz[1]:
                fe = fe * imsz[1]
            xy0 = np.array([(imsz[1] + 1) / 2.0, (imsz[0] + 1) / 2.0])
            k_mat = np.eye(3)
            k_mat[0, 0] = fe
            k_mat[1, 1] = fe
            k_mat[0, 2] = xy0[0]
            k_mat[1, 2] = xy0[1]
        else:
            k_mat = np.asarray(k0, dtype=float).reshape(3, 3)

        rocall = rocmean @ roc[:, :, i].T
        h = k_mat @ rocall @ np.linalg.inv(k_mat)
        img_out, xa_out, ya_out = v_imagehomog(
            np.asarray(ims[i], dtype=np.uint8), h, modeh)
        imr.append(img_out)
        xa.append(xa_out)
        ya.append(ya_out)

    return imr, xa, ya

v_skew3d

V_SKEW3D - convert between vector and skew-symmetric matrix.

v_skew3d 

v_skew3d(x, m='') -> ndarray

Convert between a vector and the corresponding skew-symmetric matrix.

v_skew3d is its own inverse: v_skew3d(v_skew3d(x)) == x.

Parameters:

Name Type Description Default
x array_like

Input vector or matrix. Size must be (3,), (3, 3), (6,), or (4, 4).

 required
m str

Mode string: 'n' normalize, 'z' orthogonalize.

 ''

Returns:

Name Type Description
y ndarray

Output matrix or vector.
Source code in pyvoicebox/v_skew3d.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def v_skew3d(x, m='') -> np.ndarray:
    """Convert between a vector and the corresponding skew-symmetric matrix.

    v_skew3d is its own inverse: v_skew3d(v_skew3d(x)) == x.

    Parameters
    ----------
    x : array_like
        Input vector or matrix. Size must be (3,), (3, 3), (6,), or (4, 4).
    m : str, optional
        Mode string: 'n' normalize, 'z' orthogonalize.

    Returns
    -------
    y : ndarray
        Output matrix or vector.
    """
    x = np.asarray(x, dtype=float)
    mn = 'n' in m
    mz = 'z' in m

    j, k = x.shape if x.ndim == 2 else (x.shape[0], 1)

    if j == 3 and k == 1:
        # 3x1 vector -> 3x3 skew-symmetric matrix
        xv = x.ravel()
        if mn and xv @ xv > 0:
            xv = xv / np.sqrt(xv @ xv)
        y = np.zeros((3, 3))
        # MATLAB column-major flat indices: [6,7,2] -> (2,0), (0,1), (2,0) in row-major
        # Actually MATLAB flat indices (1-based, column-major): 6->row2,col1; 7->row0,col2; 2->row1,col0
        # 0-based column-major: index 5 -> (2,1), index 6 -> (0,2), index 1 -> (1,0)
        y[2, 1] = xv[0]
        y[0, 2] = xv[1]
        y[1, 0] = xv[2]
        y[1, 2] = -xv[0]
        y[2, 0] = -xv[1]
        y[0, 1] = -xv[2]
        return y
    elif j == 3 and k == 3:
        # 3x3 skew-symmetric matrix -> 3x1 vector
        # Extract from positions [6,7,2] (MATLAB 1-based column-major) = (2,1), (0,2), (1,0)
        y = np.array([x[2, 1], x[0, 2], x[1, 0]])
        if mn and y @ y > 0:
            y = y / np.sqrt(y @ y)
        return y
    elif j == 6 and k == 1:
        # 6x1 Plucker vector -> 4x4 Plucker matrix
        xv = x.ravel().copy()
        u = xv[:3]
        v = xv[5:2:-1]  # [x[5], x[4], x[3]]
        if mz and u @ u > 0 and v @ v > 0:
            v = v - (u @ v) / (2 * (u @ u)) * u
            xv = np.concatenate([u - (v @ u) / (v @ v) * v, v[::-1]])
        if mn and xv @ xv > 0:
            xv = xv / np.sqrt(xv @ xv)
        y = np.zeros((4, 4))
        # MATLAB 1-based column-major indices for positive entries: [5,9,13,10,8,15]
        # 0-based column-major: [4,8,12,9,7,14]
        # -> (0,1), (0,2), (0,3), (1,2), (3,1), (2,3)
        y[0, 1] = xv[0]
        y[0, 2] = xv[1]
        y[0, 3] = xv[2]
        y[1, 2] = xv[3]
        y[3, 1] = xv[4]
        y[2, 3] = xv[5]
        # Negative entries: [2,3,4,7,14,12] -> 0-based col-major: [1,2,3,6,13,11]
        # -> (1,0), (2,0), (3,0), (2,1), (1,3), (3,2)
        y[1, 0] = -xv[0]
        y[2, 0] = -xv[1]
        y[3, 0] = -xv[2]
        y[2, 1] = -xv[3]
        y[1, 3] = -xv[4]
        y[3, 2] = -xv[5]
        return y
    elif j == 4 and k == 4:
        # 4x4 Plucker matrix -> 6x1 Plucker vector
        # MATLAB 1-based column-major: u = x([5,9,13]) = (0,1), (0,2), (0,3)
        u = np.array([x[0, 1], x[0, 2], x[0, 3]])
        # v = x([15,8,10]) = (2,3), (3,1), (1,2)
        v = np.array([x[2, 3], x[3, 1], x[1, 2]])
        if mz and u @ u > 0 and v @ v > 0:
            v = v - (u @ v) / (2 * (u @ u)) * u
            y = np.concatenate([u - (v @ u) / (v @ v) * v, v[::-1]])
        else:
            y = np.concatenate([u, v[::-1]])
        if mn and y @ y > 0:
            y = y / np.sqrt(y @ y)
        return y
    else:
        raise ValueError('size(x) must be (3,), (3, 3), (6,), or (4, 4)')

v_upolyhedron

V_UPOLYHEDRON - calculate uniform polyhedron characteristics.

This is a simplified implementation that supports the most common polyhedra by name or index. The full MATLAB implementation supports all 75+ uniform polyhedra with Wythoff symbol computation.

v_upolyhedron 

v_upolyhedron(w, md='') -> ndarray

Calculate uniform polyhedron characteristics.

This is a simplified implementation supporting common polyhedra.

Parameters:

Name Type Description Default
w int or str

Polyhedron specification (index or name).

 required
md str

Mode string.

 ''

Returns:

Name Type Description
vlist ndarray

Vertex list with columns [x, y, z, d, n, e, t].
Source code in pyvoicebox/v_upolyhedron.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def v_upolyhedron(w, md='') -> np.ndarray:
    """Calculate uniform polyhedron characteristics.

    This is a simplified implementation supporting common polyhedra.

    Parameters
    ----------
    w : int or str
        Polyhedron specification (index or name).
    md : str, optional
        Mode string.

    Returns
    -------
    vlist : ndarray
        Vertex list with columns [x, y, z, d, n, e, t].
    """
    if isinstance(w, str):
        name = w.lower()
        if name not in _POLYHEDRA:
            raise ValueError(f'Unknown polyhedron: {w}')
        vertices = _POLYHEDRA[name]()
    elif isinstance(w, (int, float)):
        idx = int(w)
        if idx not in _INDEX_MAP:
            raise ValueError(f'Polyhedron index {idx} not supported in simplified implementation')
        vertices = _POLYHEDRA[_INDEX_MAP[idx]]()
    else:
        raise ValueError('w must be a string name or integer index')

    nv = vertices.shape[0]
    d = np.sqrt(np.sum(vertices ** 2, axis=1))
    vlist = np.column_stack([
        vertices,
        d,
        np.zeros(nv),  # valency placeholder
        np.zeros(nv),  # edge index placeholder
        np.ones(nv),   # type placeholder
    ])

    return vlist

v_sphrharm

V_SPHRHARM - Forward and inverse spherical harmonic transform (stub).

The full implementation is very complex. A basic stub is provided.

v_sphrharm 

v_sphrharm(*args, **kwargs) -> None

Forward and inverse spherical harmonic transform.

This is a complex function with many modes. A basic version is available through scipy.special.sph_harm for individual harmonics.

Raises:

Type Description
NotImplementedError

Full spherical harmonic transform implementation pending. Use scipy.special.sph_harm for individual harmonics.
Source code in pyvoicebox/v_sphrharm.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def v_sphrharm(*args, **kwargs) -> None:
    """Forward and inverse spherical harmonic transform.

    This is a complex function with many modes. A basic version is
    available through scipy.special.sph_harm for individual harmonics.

    Raises
    ------
    NotImplementedError
        Full spherical harmonic transform implementation pending.
        Use scipy.special.sph_harm for individual harmonics.
    """
    raise NotImplementedError(
        "v_sphrharm full implementation pending. "
        "Use scipy.special.sph_harm(m, n, theta, phi) for individual harmonics."
    )