Source code for ddg.math.projective

"""Utility functions for the projective module."""
from itertools import combinations

import numpy as np

from ddg import nonexact
from ddg.math.linalg import rank


# --- homogenization ---
def homogenize(vector):
    """Homogenize a coordinate vector by inserting a 1.

    Parameters
    ----------
    vector : array_like of shape (n,)
        The vector to be homogenized.

    Returns
    -------
    numpy.ndarray of shape (n+1,)

    Examples
    --------
    >>> from ddg.math.projective import homogenize
    >>> homogenize((3, 4, 5))
    array([3, 4, 5, 1])
    """
    return np.append(vector, 1)


def dehomogenize(vector, atol=None, rtol=None):
    """Pick an affine representative for a homogeneous coordinate vector.

    Parameters
    ----------
    vector : numpy.ndarray of shape (n+1,)
        The vector to be dehomogenized.
    atol, rtol : float (default=None)
        This function uses the global tolerance defaults if `atol` or `rtol` are
        set to None. See :py:mod:`ddg.nonexact` for details.

    Returns
    -------
    numpy.ndarray of shape (n,)

    Raises
    ------
    ValueError
        If last value of vector is 0.

    Notes
    -----
    This function uses the global tolerance defaults if `atol` or `rtol` are
    set to None. See ddg.abc.NonExact for details.

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import dehomogenize
    >>> dehomogenize(np.array((1, 0, 1)))
    array([1., 0.])
    """
    vector = np.asanyarray(vector)
    if nonexact.isclose(vector[-1], 0, atol=atol, rtol=rtol):
        raise ValueError("Given point is at infinity.")
    return vector[0:-1] / vector[-1]


# --- projective utils ---
def in_general_position(points, atol=None, rtol=None):
    """
    Checks whether a list of k+1 points in n-dimensional projective
    are in general position.

    Parameters
    ----------
    points : array_like
        Array_like of numpy.ndarray of shape (n+1,).

    Returns
    -------
    bool
        True if the points are in general position, False otherwise.

    Raises
    ------
    ValueError
        If no points are given.

    Notes
    -----
    If k <= n, we test whether their lifts in (n+1)-dimensional
    space are linearly independent.
    If k > n, we test whether any (n+1)-subset of them is contained
    in a hyperplane of n-dimensional projective space.

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import in_general_position
    >>> points_2d = [
    ...     np.array([1.0, 0.0, 1.0]),
    ...     np.array([0.0, 1.0, 1.0]),
    ...     np.array([1.0, 1.0, 1.0]),
    ...     np.array([0.5, 0.5, 1.0]),
    ... ]
    >>> in_general_position(points_2d)
    False
    >>> points_3d = [
    ...     np.array([1.0, 0.0, 0.0, 1.0]),
    ...     np.array([0.0, 1.0, 0.0, 1.0]),
    ...     np.array([0.0, 0.0, 1.0, 1.0]),
    ...     np.array([1.0, 1.0, 1.0, 1.0]),
    ...     np.array([0.5, 0.5, 0.5, 1.0]),
    ... ]
    >>> in_general_position(points_3d)
    True
    """
    if not points:
        raise ValueError("No points are given to test for general position.")

    k = len(points) - 1
    n = points[0].shape[0] - 1
    if k <= n:
        return rank(np.array(points)) == k + 1

    subset_points = combinations(points, n + 1)
    return all(
        not nonexact.isclose(np.linalg.det(np.array(points)), 0, atol=atol, rtol=rtol)
        for points in list(subset_points)
    )


def is_projective_frame(points):
    """
    Checks whether a list of n+2 points in n-dimensional projective
    space form a projective frame.

    Parameters
    ----------
    points : array_like
        Array_like of numpy.ndarray of shape (n+1,)

    Returns
    -------
    bool
        True if the points are a projective frame, False otherwise.

    Raises
    ------
    ValueError
        If no points are given.

    Notes
    -----
    A list of n+1 fundamental points and one unit point in projective
    space are called a projective frame if they are in general position.

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import is_projective_frame
    >>> points_2d = [
    ...     np.array([1.0, 0.0, 1.0]),
    ...     np.array([0.0, 1.0, 1.0]),
    ...     np.array([1.0, 1.0, 1.0]),
    ...     np.array([0.5, 0.5, 1.0]),
    ... ]
    >>> is_projective_frame(points_2d)
    False
    >>> points_3d = [
    ...     np.array([1.0, 0.0, 0.0, 1.0]),
    ...     np.array([0.0, 1.0, 0.0, 1.0]),
    ...     np.array([0.0, 0.0, 1.0, 1.0]),
    ...     np.array([1.0, 1.0, 1.0, 1.0]),
    ...     np.array([0.5, 0.5, 0.5, 1.0]),
    ... ]
    >>> is_projective_frame(points_3d)
    True
    """
    if not points:
        raise ValueError("No points are given to test for projective frame.")

    n = len(points) - 2
    if np.shape(points) != (n + 2, n + 1):
        raise ValueError(
            "Only n+2 points can form a projective frame "
            "for n-dimensional projective space."
        )
    return in_general_position(points)


def point_of_intersection(planes, homogeneous_coords=True):
    """
    Computes the point of intersection of given hyperplanes (in dual coordinates).
    If more hyperplanes are given than the dimension of the space,
    a least square solution is computed.

    E.g. in RP^3 it computes the intersection of (at least) three given planes.

    Parameters
    ----------
    planes : numpy.ndarray of shape (n, k)
        Amount of planes `n` and dimension of each plane `k`.
    homogeneous: bool, (default=True)
        Determines whether a `1.0` will be added to the output vector.

    Returns
    -------
    intersection : numpy.ndarray of shape (k + 1,) or (k,)
        Intersection of the given planes

    Examples
    --------
    >>> from ddg.math.projective import point_of_intersection
    >>> plane1 = np.array([1.0, 0.0, 0.0, -1.0])
    >>> plane2 = np.array([0.0, 1.0, 0.0, -2.0])
    >>> plane3 = np.array([0.0, 0.0, 1.0, -3.0])
    >>> planes = np.array((plane1, plane2, plane3))
    >>> point_of_intersection(planes)
    array([1., 2., 3., 1.])
    """
    if np.linalg.matrix_rank(planes[:, 0:-1]) < len(planes[0]) - 1:
        raise ValueError("Intersection can only be computed if rank(planes) >= 3")
    A = planes[:, 0:-1]
    b = -planes[:, -1]
    x = np.linalg.lstsq(A, b, rcond=None)[0]
    if homogeneous_coords:
        x = np.hstack((x, [1.0]))
    return x


# --- affine transformations ---
def affine_transformation(A, b):
    """Assemble affine transformation.

    Given a matrix `A` and a vector `b`, assemble the projective
    transformation ::

             A | b
        F =  --|---
             0 | 1

    It corresponds in affine coordinates to the
    affine transformation ``x -> Ax + b``.

    Parameters
    ----------
    A : numpy.ndarray of shape (n, m)
    b : numpy.ndarray of shape (n,)

    Returns
    -------
    F : numpy.ndarray of shape (n+1, m+1)

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import affine_transformation
    >>> A = np.arange(9).reshape((3, 3))
    >>> A
    array([[0, 1, 2],
           [3, 4, 5],
           [6, 7, 8]])
    >>> b = np.ones(3) * 15
    >>> b
    array([15., 15., 15.])
    >>> affine_transformation(A, b)
    array([[ 0.,  1.,  2., 15.],
           [ 3.,  4.,  5., 15.],
           [ 6.,  7.,  8., 15.],
           [ 0.,  0.,  0.,  1.]])
    """
    n, m = np.shape(A)
    F = np.zeros((n + 1, m + 1))
    F[:-1, :-1] = A
    F[:-1, -1] = b
    F[-1, -1] = 1
    return F


def decompose_affine_transformation(F, atol=None, rtol=None):
    """Decompose affine transformation.

    Given a projective transformation ::

             A | b
        F =  --|---
             0 | s

    disassemble it into a matrix `A/s` and a vector `b/s`.

    Parameters
    ----------
    F : numpy.ndarray of shape (n+1, m+1)
    atol, rtol : float (default=None)
        This function uses the global tolerance defaults if `atol` or `rtol` are
        set to None. See :py:mod:`ddg.nonexact` for details.

    Returns
    -------
    A : numpy.ndarray of shape (n, m)
    b : numpy.ndarray of shape (n,)

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import decompose_affine_transformation
    >>> F = np.array([[0, 1, 2, 15], [3, 4, 5, 15], [6, 7, 8, 15], [0, 0, 0, 1]])
    >>> F
    array([[ 0,  1,  2, 15],
           [ 3,  4,  5, 15],
           [ 6,  7,  8, 15],
           [ 0,  0,  0,  1]])
    >>> A, b = decompose_affine_transformation(F)
    >>> A
    array([[0., 1., 2.],
           [3., 4., 5.],
           [6., 7., 8.]])
    >>> b
    array([15., 15., 15.])
    """
    if nonexact.isclose(F[-1, -1], 0, atol=atol, rtol=rtol) or not nonexact.allclose(
        F[-1, :-1], 0, atol=atol, rtol=rtol
    ):
        raise ValueError("Not an affine transformation")
    F_dehomogenized = F / F[-1, -1]
    A = F_dehomogenized[:-1, :-1]
    b = F_dehomogenized[:-1, -1]
    return A, b


def decompose_similarity(F, atol=None, rtol=None):
    """Decompose a similarity transformation into scaling, orthogonal matrix
    and translation.

    Parameters
    ----------
    F : numpy.ndarray of shape (n+1, m+1)
    atol, rtol : float (default=None)

    Returns
    -------
    s : float
        Scaling factor
    O : numpy.ndarray of shape (n, m)
        Orthogonal matrix
    b : numpy.ndarray of shape (n,)
        Translation vector
    atol, rtol : float (default=None)

    Raises
    ------
    ValueError
        If `F` does not represent a similarity embedding.

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import decompose_similarity
    >>> F = np.array(
    ...     [
    ...         [2.0, 0.0, 0.0, 0.0],
    ...         [0.0, 2.0, 0.0, 0.0],
    ...         [0.0, 0.0, 2.0, 0.0],
    ...         [0.0, 0.0, 0.0, 1.0],
    ...     ]
    ... )
    >>> s, O, b = decompose_similarity(F)
    >>> s
    2.0
    >>> O
    array([[1., 0., 0.],
           [0., 1., 0.],
           [0., 0., 1.]])
    >>> b
    array([0., 0., 0.])
    """
    A, b = decompose_affine_transformation(F, atol=atol, rtol=rtol)
    D = A.T @ A
    d = D[0, 0]
    if not nonexact.allclose(D / d, np.eye(D.shape[0]), atol=atol, rtol=rtol):
        raise ValueError("Not a similarity")
    s = np.sqrt(d)
    O = A / s
    return s, O, b


def affine_component_transformation(n, before, after=-1):
    """Transformation transforming between affine pictures.

    If an object dehomogenized wth affine component `before` produces a certain
    affine picture, the object transformed with `F` will produce the same
    picture when dehomogenized with affine component `after`.

    Parameters
    ----------
    n : int
        Size of returned matrix
    before : int
    after : int (default=-1)

    Returns
    -------
    F : numpy.ndarray of shape (n, n)

    Examples
    --------
    >>> from ddg.math.projective import affine_component_transformation
    >>> affine_component_transformation(4, 2)
    array([[1., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., 0., 1.],
           [0., 0., 1., 0.]])
    >>> affine_component_transformation(5, 4, 2)
    array([[1., 0., 0., 0., 0.],
           [0., 1., 0., 0., 0.],
           [0., 0., 0., 0., 1.],
           [0., 0., 1., 0., 0.],
           [0., 0., 0., 1., 0.]])
    """
    F = np.eye(n)
    ac_entries = F[before]
    F = np.delete(F, before, axis=0)
    # Convert to non-negative index because numpy inserts BEFORE the index,
    # which we don't always want.
    if after < 0:
        after += n
    F = np.insert(F, after, ac_entries, axis=0)
    return F


# --- translations ---
def translation(normal, direction, atol=None, rtol=None):
    """Projective translation.

    Returns the matrix ::

        I + direction @ normal.T

    This transformation fixes all vectors v in the hyperplane orthogonal to
    `normal`. It can also be thought of as mapping ::

        v -> v + dot(normal, v) * direction

    Parameters
    ----------
    normal, direction : numpy.ndarray of shape (n,)
        `direction` must lie in the hyperplane orthogonal to `normal`.
    atol, rtol : float (default=None)
        This function uses the global tolerance defaults if `atol` or `rtol` are
        set to None. See :py:mod:`ddg.nonexact` for details.

    Returns
    -------
    numpy.ndarray of shape (n, n)

    Raises
    ------
    ValueError
        * If the dimensions of `normal` and `direction` do not match.
        * If `direction` is not orthogonal to `normal`.

    Examples
    --------
    >>> from ddg.math.projective import translation
    >>> translation(np.array((1, 0, 1)), np.array((-1, 0, 1)))
    array([[ 0.,  0., -1.],
           [ 0.,  1.,  0.],
           [ 1.,  0.,  2.]])
    >>> translation(np.array((1, 0, 1, 0)), np.array((-1, 0, 1, 0)))
    array([[ 0.,  0., -1.,  0.],
           [ 0.,  1.,  0.,  0.],
           [ 1.,  0.,  2.,  0.],
           [ 0.,  0.,  0.,  1.]])
    """
    if np.size(normal) != np.size(direction):
        raise ValueError("normal and direction must have equal dimensions.")
    if not nonexact.isclose(np.dot(normal, direction), 0.0, atol=atol, rtol=rtol):
        raise ValueError("direction must be orthogonal to normal.")
    return np.eye(np.size(normal)) + np.outer(direction, normal)


def translations_from_quad_3d(x00, x10, x11, x01, atol=None, rtol=None):
    """Translations that map along a quadrilateral in RP2.

    Returns two projective translations `T1`, `T2` of RP2 such that ::

               T1
          x01 ---> x11
           ^        ^
        T2 |        | T2
           |        |
          x00 ---> x10
               T1

    Parameters
    ----------
    x00, x10, x11, x01 : numpy.ndarray of shape (3,)
    atol, rtol : float (default=None)
        This function uses the global tolerance defaults if `atol` or `rtol` are
        set to None. See :py:mod:`ddg.nonexact` for details.

    Returns
    -------
    T1, T2 : numpy.ndarray of shape (4, 4)

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import translations_from_quad_3d
    >>> x00 = np.array([0.0, 0.0, 1.0])
    >>> x10 = np.array([1.0, 0.0, 1.0])
    >>> x11 = np.array([1.0, 1.0, 1.0])
    >>> x01 = np.array([0.0, 1.0, 1.0])
    >>> T1, T2 = translations_from_quad_3d(x00, x10, x11, x01)
    >>> np.array(np.round(T1), dtype=int)
    array([[1, 0, 1],
           [0, 1, 0],
           [0, 0, 1]])
    >>> np.array(np.round(T2), dtype=int)
    array([[1, 0, 0],
           [0, 1, 1],
           [0, 0, 1]])
    """
    quad = np.array([x00, x10, x11, x01])
    U, _, _ = np.linalg.svd(quad)
    dependency = U[:, 3]
    y1 = dependency[0] * x00 + dependency[1] * x10
    y2 = dependency[0] * x00 + dependency[3] * x01
    plane = np.cross(y1, y2)
    plane = -1.0 / dependency[0] / plane.dot(x00) * plane
    T1 = translation(plane, y1, atol=atol, rtol=rtol)
    T2 = translation(plane, y2, atol=atol, rtol=rtol)
    return T1, T2


def translations_from_quad_4d(x00, x10, x11, x01, n, atol=None, rtol=None):
    """Translations that map along a quadrilateral in a plane in RP3.

    Returns two projective translations `T1`, `T2` of RP3 such that ::

               T1
          x01 ---> x11
           ^        ^
        T2 |        | T2
           |        |
          x00 ---> x10
               T1

    Parameters
    ----------
    x00, x10, x11, x01, n : numpy.ndarray of shape (4,)
        The corners of the quadrilateral must lie in the plane orthogonal to n.
    atol, rtol : float (default=None)
        This function uses the global tolerance defaults if `atol` or `rtol` are
        set to None. See :py:mod:`ddg.nonexact` for details.

    Returns
    -------
    T1, T2 : numpy.ndarray of shape (4, 4)

    Examples
    --------
    >>> import numpy as np
    >>> from ddg.math.projective import translations_from_quad_4d
    >>> x00 = np.array([0.0, 0.0, 0.0, 1.0])
    >>> x10 = np.array([1.0, 0.0, 0.0, 1.0])
    >>> x11 = np.array([1.0, 1.0, 0.0, 1.0])
    >>> x01 = np.array([0.0, 1.0, 0.0, 1.0])
    >>> n = np.array([0.0, 0.0, 1.0, 0.0])
    >>> T1, T2 = translations_from_quad_4d(x00, x10, x11, x01, n)
    >>> np.array(np.round(T1), dtype=int)
    array([[1, 0, 0, 1],
           [0, 1, 0, 0],
           [0, 0, 1, 0],
           [0, 0, 0, 1]])
    >>> np.array(np.round(T2), dtype=int)
    array([[1, 0, 0, 0],
           [0, 1, 0, 1],
           [0, 0, 1, 0],
           [0, 0, 0, 1]])
    """
    quad = np.array([x00, x10, x11, x01])
    U, _, _ = np.linalg.svd(quad)
    dependency = U[:, 3]
    y1 = dependency[0] * x00 + dependency[1] * x10
    y2 = dependency[0] * x00 + dependency[3] * x01
    _, _, V2 = np.linalg.svd([y1, y2, n])
    plane = V2[3, :]
    plane = -1.0 / dependency[0] / plane.dot(x00) * plane
    T1 = translation(plane, y1, atol=atol, rtol=rtol)
    T2 = translation(plane, y2, atol=atol, rtol=rtol)
    return T1, T2