BUG: spatial: Rotation.align_vectors() incorrect for anti-parallel vectors

[kalekundert]

Describe your issue.

Rotation.align_vectors(a, b) is supposed to calculate the "best estimate of the rotation that transforms b to a". However, with a = np.array([-1, 0, 0]) and b = np.array([1, 0, 0]), it incorrectly returns the identity transformation.

Reproducing Code Example

from scipy.spatial.transform import Rotation

a = np.array([-1, 0, 0])
b = np.array([ 1, 0, 0])

rot, _ = Rotation.align_vectors(a, b)

rot.as_matrix()
# array([[ 1., -0.,  0.],
#        [ 0.,  1., -0.],
#        [ 0.,  0.,  1.]])

Error message

N/A

SciPy/NumPy/Python version and system information

1.13.0 1.26.4 sys.version_info(major=3, minor=12, micro=2, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-giifv1os/overlay/lib/python3.12/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp312-cp312/bin/python
  version: '3.12'

[kalekundert]

I looked into this more closely, and I found that I can get the right answer by copying the exact cython code from scipy/spatial/transform/_rotation.pyx (minus some irrelevant code paths, for clarity) into a normal python script:

import numpy as np
from scipy.spatial.transform import Rotation
from math import sqrt, atan2

def py_align_vectors(a, b):
    a_original = np.array(a, dtype=float)
    b_original = np.array(b, dtype=float)
    a = np.atleast_2d(a_original)
    b = np.atleast_2d(b_original)

    N = len(a)

    a_pri = a
    b_pri = b
    a_pri_norm = _norm3(a_pri[0])
    b_pri_norm = _norm3(b_pri[0])
    if a_pri_norm == 0 or b_pri_norm == 0:
        raise ValueError("Cannot align zero length primary vectors")
    a_pri /= a_pri_norm
    b_pri /= b_pri_norm

    # We first find the minimum angle rotation between the primary
    # vectors.
    cross = np.cross(b_pri[0], a_pri[0])
    cross_norm = _norm3(cross)
    theta = atan2(cross_norm, _dot3(a_pri[0], b_pri[0]))
    tolerance = 1e-3  # tolerance for small angle approximation (rad)
    R_flip = Rotation.identity()
    if (np.pi - theta) < tolerance:
        # Near pi radians, the Taylor series appoximation of x/sin(x)
        # diverges, so for numerical stability we flip pi and then
        # rotate back by the small angle pi - theta
        if cross_norm == 0:
            # For antiparallel vectors, cross = [0, 0, 0] so we need to
            # manually set an arbitrary orthogonal axis of rotation
            i = np.argmin(np.abs(a_pri[0]))
            r = np.zeros(3)
            r[i - 1], r[i - 2] = a_pri[0][i - 2], -a_pri[0][i - 1]
        else:
            r = cross  # Shortest angle orthogonal axis of rotation
        R_flip = Rotation.from_rotvec(r / np.linalg.norm(r) * np.pi)
        theta = np.pi - theta
        cross = -cross
    if abs(theta) < tolerance:
        # Small angle Taylor series approximation for numerical stability
        theta2 = theta * theta
        r = cross * (1 + theta2 / 6 + theta2 * theta2 * 7 / 360)
    else:
        r = cross * theta / np.sin(theta)
    R_pri = Rotation.from_rotvec(r) * R_flip

    # No secondary vectors, so we are done
    R_opt = R_pri

    return R_opt, None

def _empty1(n):
    return np.array(shape=(n,))

def _cross3(a, b):
    result = _empty1(3)
    result[0] = a[1]*b[2] - a[2]*b[1]
    result[1] = a[2]*b[0] - a[0]*b[2]
    result[2] = a[0]*b[1] - a[1]*b[0]
    return result

def _dot3(a, b):
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

def _norm3(elems):
    return sqrt(_dot3(elems, elems))

a = np.array([-1, 0, 0])
b = np.array([ 1, 0, 0])

R1, _ = py_align_vectors(a, b)
# array([[-1.0000000e+00, -1.2246468e-16,  0.0000000e+00],
#        [ 1.2246468e-16, -1.0000000e+00,  0.0000000e+00],
#        [ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00]])

R2, _ = Rotation.align_vectors(a, b)
# array([[ 1., -0.,  0.],
#        [ 0.,  1., -0.],
#        [ 0.,  0.,  1.]])

So it seems that the issue is something relating to cython, but that's pretty much the limit of my debugging ability.

[scottshambaugh]

This was fixed with #20573, should be out in 1.14.0, and I believe will be backported to 1.13.1. But good catch, this was a tricky edge case!

I believe your pasting into a separate script working was just a floating point thing.

[kalekundert]

Thanks for the fix, and sorry I missed the already-open issue!