"""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 ---
[docs]def homogenize(vector):
"""Homogenize a coordinate vector by inserting a 1.
Parameters
----------
vector : numpy.ndarray of shape (n,)
The vector to be homogenized.
Returns
-------
numpy.ndarray of shape (n+1,)
"""
return np.append(vector, 1)
[docs]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,)
"""
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 ---
[docs]def in_general_position(points):
"""
Checks whether a list of k+1 points in n-dimensional projective
are in general position.
Parameters
----------
points : list
List 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.
"""
if points:
k = len(points) - 1
n = points[0].shape[0] - 1
if k <= n:
return True if rank(np.array(points)) == k + 1 else False
else:
subset_points = combinations(points, n + 1)
for points in list(subset_points):
if np.linalg.det(np.array(points)) == 0:
return False
return True
else:
raise ValueError("No points are given to test for general position.")
[docs]def is_projective_frame(points):
"""
Checks whether a list of n+2 points in n-dimensional projective
space form a projective frame.
Parameters
----------
points : list
List 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.
"""
if points:
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."
)
else:
return in_general_position(points)
else:
raise ValueError("No points are given to test for projective frame.")
[docs]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 : list
List of planes given in dual coordinates,
i.e. by a vector of the form (a,-b) for a.T * x - b = 0
Returns
-------
intersection : array
Intersection of the given planes
"""
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 ---
[docs]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.
"""
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
# --- translations ---
[docs]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`.
"""
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)
[docs]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)
"""
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
[docs]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)
"""
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