"""Utility functions for the projective module."""
import numpy as np
from itertools import combinations
from ddg.math.linalg import e, rank
from ddg.abc import NonExact
# --- homogenization ---
[docs]def homogenize(vector, affine_component=-1):
"""Homogenize a coordinate vector by inserting a 1.
Parameters
----------
vector : numpy.ndarray of shape (n,)
The vector to be homogenized.
affine_component : int (default=-1)
Index of inserted 1 in the output vector.
Returns
-------
numpy.ndarray of shape (n+1,)
"""
i = affine_component
if -np.size(vector) <= i <= -1:
i += np.size(vector) + 1
return np.insert(vector, i, 1)
[docs]@NonExact.nonexact_function
def dehomogenize(vector, affine_component=-1, atol=None, rtol=None):
"""Pick an affine representative for a homogenous coordinate vector.
Parameters
----------
vector : numpy.ndarray of shape (n+1,)
The vector to be dehomogenized.
affine_component : int (default=-1)
Index of affine component normalized to 1
in the output vector.
Returns
-------
numpy.ndarray of shape (n,)
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
i = affine_component
if np.isclose(vector[i], 0, atol=atol, rtol=rtol):
raise ValueError("Given point is at infinity.")
if i == -1:
return vector[0:-1] / vector[-1]
else:
return np.concatenate((vector[:i], vector[i + 1:])) / vector[i]
# --- 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]@NonExact.nonexact_function
def is_similarity(F, affine_component=-1, atol=None, rtol=None):
"""Check whether an affine transformation F is a similarity transformation.
Parameters
----------
F : numpy.ndarray of shape (n+1, n+1)
The affine transformation to check.
affine_component : int, default=-1
The diagonal index of affine component.
atol, rtol : float
Returns
-------
bool
True if F is a similarity transformation.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
A, b = decompose_affine_transformation(
F,
affine_component=affine_component
)
I = A.T @ A
I /= I[0, 0]
if np.allclose(I, np.eye(I.shape[0]), atol=atol, rtol=rtol):
return True
else:
return False
# --- translations ---
[docs]@NonExact.nonexact_function
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
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`.
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
if np.size(normal) != np.size(direction):
raise ValueError("normal and direction must have equal dimensions.")
if not np.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]@NonExact.nonexact_function
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,)
Returns
-------
T1, T2 : numpy.ndarray of shape (4, 4)
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
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]@NonExact.nonexact_function
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.
Returns
-------
T1, T2 : numpy.ndarray of shape (4, 4)
Notes
-----
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
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