Remove dependency on transforms3d (#1392)

* Add reimplemented transformations

* Add temporary, proof-of-concept tests

* Remove transforms3d implementation and sanity check

* Remove transforms3d dependency and mentions

* Explain math in transformation docstrings

* Add parameters and examples to transformations

* Add dedicated transformations tests

* Use exact check for the fast path

environment.yml

@@ -8,7 +8,6 @@ dependencies:
   - scooby>=0.5.1
   - meshio>=4.0.3, <5.0
   - matplotlib
-  - transforms3d==0.3.1
   - appdirs
   - pyqt
   - imageio>=2.5.0

pyvista/utilities/transformations.py

@@ -1,20 +1,214 @@
-"""Transformations leveraging transforms3d library."""
+"""Module implementing point transformations and their matrices."""
 import numpy as np
 
 
-
 def axis_angle_rotation(axis, angle, point=None, deg=True):
-    """Return a 4x4 matrix for rotation about an axis by given angle, optionally about a given point."""
+    r"""Return a 4x4 matrix for rotation about any axis by given angle.
+
+    Rotations around an axis that contains the origin can easily be
+    computed using Rodrigues' rotation formula. The key quantity is
+    the ``K`` cross product matrix for the unit vector ``n`` defining
+    the axis of the rotation:
+
+             /   0  -nz   ny \
+        K =  |  nz    0  -nx |
+             \ -ny   nx    0 /
+
+    For a rotation angle ``phi`` around the vector ``n`` the rotation
+    matrix is given by
+
+        R = I + sin(phi) K  + (1 - cos(phi)) K^2
+
+    where ``I`` is the 3-by-3 unit matrix and ``K^2`` denotes the matrix
+    square of ``K``.
+
+    If the rotation axis doesn't contain the origin, we have to first
+    shift real space to transform the axis' ``p0`` reference point into
+    the origin, then shift the points back after rotation:
+
+        p' = R @ (p - p0) + p0 = R @ p + (p0 - R @ p0)
+
+    This means that the rotation in general consists of a 3-by-3
+    rotation matrix ``R``, and a translation given by
+    ``b = p0 - R @ p0``. These can be encoded in a 4-by-4 transformation
+    matrix by filling the 3-by-3 leading principal submatrix with ``R``,
+    and filling the top 3 values in the last column with ``b``.
+
+    Parameters
+    ----------
+    axis : 3-length sequence
+        The direction vector of the rotation axis. It need not be a
+        unit vector, but it must not be a zero vector.
+
+    angle : float
+        Angle of rotation around the axis. The angle is defined as a
+        counterclockwise rotation when facing the normal vector of the
+        rotation axis. Passed either in degrees or radians depending on
+        the value of ``deg``.
+
+    point : 3-length sequence, optional
+        The origin of the rotation (a reference point through which the
+        rotation axis passes). By default the rotation axis contains the
+        origin.
+
+    deg : bool, optional
+        Whether the angle is specified in degrees. ``False`` implies
+        radians.
+
+    Examples
+    --------
+    Generate a transformation matrix for rotation around a cube's body
+    diagonal by 120 degrees.
+
+    >>> import numpy as np
+    >>> from pyvista import transformations
+    >>> trans = transformations.axis_angle_rotation([1, 1, 1], 120)
+
+    Check that the transformation cycles the cube's three corners.
+
+    >>> corners = np.array([
+    ...     [1, 0, 0],
+    ...     [0, 1, 0],
+    ...     [0, 0, 1],
+    ... ])
+    >>> rotated = transformations.apply_transformation_to_points(trans, corners)
+    >>> np.allclose(rotated, corners[[1, 2, 0], :])
+    True
+
+    """
     if deg:
+        # convert to radians
         angle *= np.pi / 180
-    import transforms3d as tf3d
-    return tf3d.axangles.axangle2aff(axis, angle, point=point)
+
+    # return early for no rotation; play it safe and check only exact equality
+    if angle % (2 * np.pi) == 0:
+        return np.eye(4)
+
+    axis = np.asarray(axis, dtype='float64')
+    if axis.shape != (3,):
+        raise ValueError('Axis must be a 3-length array-like.')
+    if point is not None:
+        point = np.asarray(point)
+        if point.shape != (3,):
+            raise ValueError('Rotation center must be a 3-length array-like.')
+
+    # check and normalize
+    axis_norm = np.linalg.norm(axis)
+    if np.isclose(axis_norm, 0):
+        raise ValueError('Cannot rotate around zero vector axis.')
+    if not np.isclose(axis_norm, 1):
+        axis = axis / axis_norm
+
+    # build Rodrigues' rotation matrix
+    K = np.zeros((3, 3))
+    K[[2, 0, 1], [1, 2, 0]] = axis
+    K += -K.T
+    R = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K
+    augmented = np.eye(4)
+    augmented[:-1, :-1] = R
+
+    if point is not None:
+        # rotation of point p would be R @ (p - point) + point
+        # which is R @ p + (point - R @ point)
+        augmented[:-1, -1] = point - R @ point
+
+    return augmented
 
 
 def reflection(normal, point=None):
-    """Return a 4x4 matrix for reflection across a normal about a point."""
-    import transforms3d as tf3d
-    return tf3d.reflections.rfnorm2aff(normal, point=point)
+    """Return a 4x4 matrix for reflection across a normal about a point.
+
+    Projection to a unit vector ``n`` can be computed using the dyadic
+    product (or outer product) ``P`` of ``n`` with itself, which is a
+    3-by-3 symmetric matrix.
+
+    Reflection across a plane that contains the origin amounts to
+    reversing the components of real space points that are perpendicular
+    to the reflection plane. This gives us the transformation ``R``
+    acting on a point ``p`` as
+
+        p' = R @ p = p - 2 P @ p = (I - 2 P) @ p
+
+    so the reflection's transformation matrix is the unit matrix minus
+    twice the dyadic product ``P``.
+
+    If additionally we want to compute a reflection to a plane that does
+    not contain the origin, we can we can first shift every point in
+    real space by ``-p0`` (if ``p0`` is a point that lies on the plane)
+
+        p' = R @ (p - p0) + p0 = R @ p + (p0 - R @ p0)
+
+    This means that the reflection in general consists of a 3-by-3
+    reflection matrix ``R``, and a translation given by
+    ``b = p0 - R @ p0``. These can be encoded in a 4-by-4 transformation
+    matrix by filling the 3-by-3 leading principal submatrix with ``R``,
+    and filling the top 3 values in the last column with ``b``.
+
+    Parameters
+    ----------
+    normal : 3-length sequence
+        The normal vector of the reflection plane. It need not be a unit
+        vector, but it must not be a zero vector.
+
+    point : 3-length sequence, optional
+        The origin of the reflection (a reference point through which
+        the reflection plane passes). By default the reflection plane
+        contains the origin.
+
+    Examples
+    --------
+    Generate a transformation matrix for reflection over the XZ plane.
+
+    >>> import numpy as np
+    >>> from pyvista import transformations
+    >>> trans = transformations.reflection([0, 1, 0])
+
+    Check that the reflection transforms corners of a cube among one
+    another.
+
+    >>> verts = np.array([
+    ...     [ 1, -1,  1],
+    ...     [-1, -1,  1],
+    ...     [-1, -1, -1],
+    ...     [-1, -1,  1],
+    ...     [ 1,  1,  1],
+    ...     [-1,  1,  1],
+    ...     [-1,  1, -1],
+    ...     [-1,  1,  1],
+    ... ])
+    >>> mirrored = transformations.apply_transformation_to_points(trans, verts)
+    >>> np.allclose(mirrored, verts[[np.r_[4:8, 0:4]], :])
+    True
+
+    """
+    normal = np.asarray(normal, dtype='float64')
+    if normal.shape != (3,):
+        raise ValueError('Normal must be a 3-length array-like.')
+    if point is not None:
+        point = np.asarray(point)
+        if point.shape != (3,):
+            raise ValueError('Plane reference point must '
+                             'be a 3-length array-like.')
+
+    # check and normalize
+    normal_norm = np.linalg.norm(normal)
+    if np.isclose(normal_norm, 0):
+        raise ValueError('Plane normal cannot be zero.')
+    if not np.isclose(normal_norm, 1):
+        normal = normal / normal_norm
+
+    # build reflection matrix
+    projection = np.outer(normal, normal)
+    R = np.eye(3) - 2 * projection
+    augmented = np.eye(4)
+    augmented[:-1, :-1] = R
+
+    if point is not None:
+        # reflection of point p would be R @ (p - point) + point
+        # which is R @ p + (point - R @ point)
+        augmented[:-1, -1] = point - R @ point
+
+    return augmented
 
 
 def apply_transformation_to_points(transformation, points, inplace=False):
@@ -50,6 +244,7 @@ def apply_transformation_to_points(transformation, points, inplace=False):
     >>> tf[3, 3,] = 1
     >>> pyvista.transformations.apply_transformation_to_points(tf, points, inplace=True)
     >>> assert np.all(np.isclose(points, scale_factor * points_orig))
+
     """
     transformation_shape = transformation.shape
     if transformation_shape not in ((3, 3), (4, 4)):

setup.py

@@ -22,7 +22,6 @@
                     'scooby>=0.5.1',
                     'meshio>=4.0.3, <5.0',
                     'vtk',
-                    'transforms3d==0.3.1'
                     ]
 
 readme_file = os.path.join(filepath, 'README.rst')

tests/test_utilities.py

@@ -515,3 +515,74 @@ def test_vtk_error_catcher():
     error_catcher = pyvista.utilities.errors.VtkErrorCatcher(raise_errors=True)
     with error_catcher:
         pass
+
+
+def test_axis_angle_rotation():
+    # rotate cube corners around body diagonal
+    points = np.array([
+        [1, 0, 0],
+        [0, 1, 0],
+        [0, 0, 1],
+    ])
+    axis = [1, 1, 1]
+
+    # no-op case
+    angle = 360
+    trans = transformations.axis_angle_rotation(axis, angle)
+    actual = transformations.apply_transformation_to_points(trans, points)
+    assert np.array_equal(actual, points)
+
+    # default origin
+    angle = np.radians(120)
+    expected = points[[1, 2, 0], :]
+    trans = transformations.axis_angle_rotation(axis, angle, deg=False)
+    actual = transformations.apply_transformation_to_points(trans, points)
+    assert np.allclose(actual, expected)
+
+    # non-default origin
+    p0 = [-2, -3, 4]
+    points += p0
+    expected += p0
+    trans = transformations.axis_angle_rotation(axis, angle, point=p0, deg=False)
+    actual = transformations.apply_transformation_to_points(trans, points)
+    assert np.allclose(actual, expected)
+
+    # invalid cases
+    with pytest.raises(ValueError):
+        transformations.axis_angle_rotation([1, 0, 0, 0], angle)
+    with pytest.raises(ValueError):
+        transformations.axis_angle_rotation(axis, angle, point=[1, 0, 0, 0])
+    with pytest.raises(ValueError):
+        transformations.axis_angle_rotation([0, 0, 0], angle)
+
+
+def test_reflection():
+    # reflect points of a square across a diagonal
+    points = np.array([
+        [ 1,  1, 0],
+        [-1,  1, 0],
+        [-1, -1, 0],
+        [ 1, -1, 0],
+    ])
+    normal = [1, 1, 0]
+
+    # default origin
+    expected = points[[2, 1, 0, 3], :]
+    trans = transformations.reflection(normal)
+    actual = transformations.apply_transformation_to_points(trans, points)
+    assert np.allclose(actual, expected)
+
+    # non-default origin
+    p0 = [1, 1, 0]
+    expected += 2 * np.array(p0)
+    trans = transformations.reflection(normal, point=p0)
+    actual = transformations.apply_transformation_to_points(trans, points)
+    assert np.allclose(actual, expected)
+
+    # invalid cases
+    with pytest.raises(ValueError):
+        transformations.reflection([1, 0, 0, 0])
+    with pytest.raises(ValueError):
+        transformations.reflection(normal, point=[1, 0, 0, 0])
+    with pytest.raises(ValueError):
+        transformations.reflection([0, 0, 0])