"""Euclidean geometry module.
Contains model classes and functions for conversion between the models.
"""
import numpy as np
import ddg.math.euclidean as em
import ddg.math.linalg as la
import ddg.math.projective as pmath
import ddg.math.symmetric_matrices as sm
from ddg import nonexact
from ddg.geometry import moebius_models
from ddg.geometry.geometry_model_templates import CayleyKleinGeometry, MetricGeometry
from ddg.geometry.quadrics import Quadric, intersect_quadric_subspace
from ddg.geometry.spheres import QuadricSphere, _QuadricSphereGeometryBridge
from ddg.geometry.subspaces import (
Point,
coordinate_hyperplane,
subspace_from_affine_points,
subspace_from_affine_points_and_directions,
)
__all__ = [
"ProjectiveModel",
"DualProjectiveModel",
"MoebiusModel",
"projective_to_moebius",
"moebius_to_projective",
]
[docs]class DualProjectiveModel(CayleyKleinGeometry):
"""Dual projective model of Euclidean geometry.
.. rubric:: Model space
The model space is RP^n without [0,...,0, 1]. The points in this space can
be thought of as hyperplanes in R^n.
.. rubric:: Representation of objects
Objects are just projective points, subspaces, quadrics and Cayley-Klein
spheres that do not contain the origin.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
"""
@property
def absolute(self):
"""The absolute quadric with matrix ``diag([1,...,1, 0])``.
Returns
-------
Quadric
"""
matrix = np.eye(self.dimension + 1)
matrix[-1, -1] = 0
return Quadric(matrix)
[docs] def angle(self, subspace1, subspace2):
"""Euclidean angle.
Computes the angle between the duals of the subspaces.
Parameters
----------
subspace1, subspace2 : Subspace
Returns
-------
float
Raises
------
NotImplementedError
If the subspaces have unsupported dimensions.
"""
if subspace1.dimension == 0 and subspace1.dimension == 0:
return np.arccos(np.sqrt(self.cayley_klein_distance(subspace1, subspace2)))
if subspace1.codimension == 1 and subspace2.codimension == 1:
return ProjectiveModel(self.dimension).angle(
subspace1.dualize(), subspace2.dualize()
)
raise NotImplementedError(
"Angle can currently only be computed between two vectors or "
"between two hyperplanes."
)
def __contains__(self, point):
return self.inner_product(point, point) > 0
[docs]class ProjectiveModel(CayleyKleinGeometry, MetricGeometry):
"""Euclidean geometry.
.. rubric:: Model space
The model space is RP^n without the points at infinity.
.. rubric:: Representation of objects
Objects are just projective points, subspaces and quadrics not at infinity
and Euclidean spheres.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
"""
@property
def absolute(self):
"""The absolute quadric.
This is the dual quadric of the quadric with matrix ``diag([1,...,1, 0])``.
As the dual of a degenerate quadric, it is contained in the hyperplane at
infinity.
Returns
-------
Quadric
"""
return DualProjectiveModel(self.dimension).absolute.dualize()
[docs] @staticmethod
def d(x, y):
"""Euclidean distance.
Dehomogenizes and computes the norm of the difference of `x` and `y`.
Parameters
----------
x, y : Point or array_like of shape (dimension+1,)
Returns
-------
float
"""
if not isinstance(x, Point):
x = Point(x)
if not isinstance(y, Point):
y = Point(y)
return np.linalg.norm(x.affine_point - y.affine_point)
[docs] def angle(self, subspace1, subspace2):
"""Euclidean angle.
If both subspaces are points, the angle between them as vectors will be
returned, up to pi. Otherwise, the smaller angle between the subspaces
will be returned, up to pi/2. Currently only hyperplanes are supported
in the latter case.
Parameters
----------
subspace1, subspace2 : Subspace
Returns
-------
float
Raises
------
NotImplementedError
If the subspaces have unsupported dimensions.
"""
if subspace1.dimension == 0 and subspace2.dimension == 0:
p1 = subspace1.affine_point
p2 = subspace2.affine_point
p1 /= np.linalg.norm(p1)
p2 /= np.linalg.norm(p2)
return np.arccos(np.dot(p1, p2))
if subspace1.codimension == 1 and subspace2.codimension == 1:
return DualProjectiveModel(self.dimension).angle(
subspace1.dualize(), subspace2.dualize()
)
raise NotImplementedError(
"Angle can currently only be computed between two vectors or "
"between two hyperplanes."
)
[docs] def ellipse_from_foci(self, f1, f2, a):
"""Get ellipse in 2D from foci and distance parameter.
The returned ellipse will be the set::
{x : |x - f1| + |x - f2| = 2a}
Parameters
----------
f1, f2 : Point
a : float
2 * a must ge greater than the distance between the foci.
Returns
-------
Quadric
Raises
------
ValueError
If 2 * a is not greater than the distance between the foci.
"""
if self.dimension != 2:
raise NotImplementedError("Function ellipse_from_foci only works in 2D.")
f1 = f1.affine_point
f2 = f2.affine_point
c = np.linalg.norm(f1 - f2) / 2
if a <= c:
raise ValueError(
"2 * a must be greater than the distance between the foci."
)
# We first create the standard ellipse centered at the origin with the
# x-axis being the major axis. Then we rotate and translate it
# appropriately.
Q = Quadric(np.diag([1 / (a**2), 1 / (a**2 - c**2), -1]))
# It doesn't matter which way the difference goes because of symmetries of the
# ellipse.
angle = np.angle(complex(*(f2 - f1)))
A = em.rotation_2d(angle)
F = pmath.affine_transformation(A, (f1 + f2) / 2)
Q = Q.transform(F)
return Q
[docs] @classmethod
def sphere(cls, center, radius, subspace=None, atol=None, rtol=None):
if not isinstance(center, Point):
center = Point(center)
dim = center.ambient_dimension
geometry_bridge = _EuclideanQSGB(cls(dim))
return QuadricSphere(
center,
radius,
subspace=subspace,
geometry_bridge=geometry_bridge,
atol=atol,
rtol=rtol,
)
[docs] @classmethod
def sphere_from_affine_point_and_normals(
cls, point, radius, normals=None, atol=None, rtol=None
):
"""Create a Euclidean sphere from an affine point and normal directions.
Parameters
----------
point : array_like of shape (n,)
Affine coordinate vector in n-dimensional space.
radius : float
normals : array_like of shape (k, n) (default=None)
The rows/elements of `normals` are vectors in affine coordinates.
`sphere.subspace` will be the affine subspace ::
orthogonal_complement(span(normals)) + point.
The default is an empty array, i.e. `sphere.subspace` will be the whole
space.
atol, rtol : float (default=None)
Tolerances for the sphere.
Returns
-------
sphere : Sphere
Let k'=dim(span(normals)). Returns a Euclidean (n-k'-1)-sphere contained
in the (n-k')-dimensional projective subspace described above.
"""
center = subspace_from_affine_points(point)
if normals is None:
normals = np.array([])
else:
normals = np.array(normals)
if normals.size == 0:
return cls.sphere(center, radius, atol=atol, rtol=rtol)
directions = list(la.nullspace(normals).T)
subspace = subspace_from_affine_points_and_directions(point, directions)
return cls.sphere(center, radius, subspace=subspace, atol=atol, rtol=rtol)
[docs] @staticmethod
def sphere_to_quadric(sphere):
"""Convert Euclidean sphere to quadric
Parameters
----------
sphere : Sphere
Returns
-------
Quadric
"""
# This is necessary because quadric.transform actually transforms
# the quadric's subspace and can only handle invertible matrices. This
# should be fixed once transform returns a new object instead of
# mutating.
atol = sphere.atol
rtol = sphere.rtol
if nonexact.isclose(sphere.radius, 0, atol=atol, rtol=rtol):
return Quadric(np.array([[0]]), subspace=sphere.center)
r = sphere.radius
c = sphere.center.affine_point
k = sphere.ambient_dimension
Q = np.eye(k + 1)
Q[-1, -1] = -1
Q = Quadric(Q)
F = pmath.affine_transformation(r * np.eye(k), c)
Q = Q.transform(F)
if sphere.subspace.codimension > 0:
Q = intersect_quadric_subspace(Q, sphere.subspace)
Q.atol = atol
Q.rtol = rtol
return Q
[docs] @classmethod
def quadric_to_sphere(cls, quadric, atol=None, rtol=None):
"""Convert quadric to a Euclidean sphere, if it is one.
Parameters
----------
quadric : Quadric
atol, rtol : float (default=None)
If None is given, the global defaults are used. See ddg.nonexact for
details.
Returns
-------
Sphere
Raises
------
ValueError
If quadric is not a sphere.
"""
# Check signature
sgn = quadric.signature(affine=True)
k = quadric.subspace.dimension
if sgn != sm.AffineSignature(k, 1, last_entry=-1):
raise ValueError(
f"Signature of quadric is {sgn}, but has to be (k, 1) with -1 "
f"at the last position for it to be a sphere."
)
Q = quadric.normalize(affine=True)
B = Q.subspace.matrix
try:
r, _, c = pmath.decompose_similarity(B, atol=atol, rtol=rtol)
except ValueError as exc:
raise ValueError("Quadric is not a Euclidean sphere.") from exc
c = subspace_from_affine_points(c)
return cls.sphere(c, r, subspace=Q.subspace)
def __contains__(self, point):
return not point.at_infinity()
class _EuclideanQSGB(_QuadricSphereGeometryBridge):
"""Auxiliary object for spheres. See the spheres developer's guide for details."""
def __init__(self, geometry):
self.geometry = geometry
def validate_sphere(self, sphere):
if sphere.center.at_infinity():
raise ValueError("Center of Euclidean spheres can not be at infinity.")
@property
def model_dimension(self):
return self.geometry.dimension
def sphere_to_quadric(self, sphere):
return self.geometry.sphere_to_quadric(sphere)
def embed(self, sphere, subspace=None):
if subspace is None:
subspace = coordinate_hyperplane(sphere.ambient_dimension + 1)
try:
# The last basis vector must be a point not at infinity and the others
# must be a scaled orthonormal basis of the subspace. Determine the
# scaling s.
s, _, _ = pmath.decompose_similarity(
subspace.matrix, atol=sphere.atol, rtol=sphere.rtol
)
except ValueError as exc:
raise ValueError(
"This embedding would not result in a Euclidean sphere. If "
"this is intentional, convert the sphere to a Quadric before "
"embedding."
) from exc
return self.geometry.sphere(
sphere.center.embed(subspace),
s * sphere.radius,
sphere.subspace.embed(subspace),
)
def unembed(self, sphere, subspace=None):
if subspace is None:
subspace = coordinate_hyperplane(sphere.ambient_dimension)
if not sphere.subspace <= subspace:
raise ValueError("Given subspace must contain sphere.")
try:
# The last basis vector must be a point not at infinity and the others
# must be a scaled orthonormal basis of the subspace. Determine the
# scaling s.
s, _, _ = pmath.decompose_similarity(
subspace.matrix, atol=sphere.atol, rtol=sphere.rtol
)
except ValueError as exc:
raise ValueError(
"This coordinate computation would not result in a Euclidean "
"sphere. If this is intentional, convert the sphere to a "
" Quadric and compute its coordinates."
) from exc
return self.geometry.sphere(
sphere.center.unembed(subspace),
sphere.radius / s,
sphere.subspace.unembed(subspace),
)
def __repr__(self):
return f"_EuclideanQSGB({self.geometry})"
[docs]class MoebiusModel(moebius_models._ProjectiveSubgeometry):
"""Euclidean geometry as a subgeometry of Möbius geometry.
.. rubric:: Model space
The model space is the standard Möbius quadric with matrix
`diag([1,...,1, -1])` without `fixed_point`.
.. rubric:: Representation of objects
Subspaces and spheres are represented as Möbius spheres, i.e. Quadric
objects contained in the Möbius quadric.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
Notes
-----
This class supports comparison with `==`. Two geometries are equal
if and only if they have the same type and `dimension`.
This class also supports the `in` operator, which checks whether a point is in the
model space of this geometry, i.e. the Möbius quadric.
"""
@property
def fixed_point(self):
"""Projection point on the Möbius quadric. Used to compute metric.
This is the north pole with representative [0,...,0, 0.5, 0.5].
Returns
-------
Point
"""
fixed_point = np.zeros(self.ambient_dimension + 1)
fixed_point[[-2, -1]] = 0.5
return Point(fixed_point)
[docs] def d(self, v, w):
"""Euclidean distance in the Möbius model.
Given by ::
d([x], [y])^2 = (-1/2) * <x, y> / (<x, p>, <y, p>)
where `p` is `self.fixed_point` and `x` and `y` are points in the
Möbius quadric not equal to `p`.
Parameters
----------
x, y : Point or array_like of shape (dimension+2,)
Returns
-------
float
"""
if isinstance(v, Point):
v = v.point
if isinstance(w, Point):
w = w.point
p = self.fixed_point.point
b = self.inner_product
return np.sqrt(-0.5 * b(v, w) / (b(v, p) * b(w, p)))
def __contains__(self, point):
return point in self.absolute and point != self.fixed_point
[docs]def projective_to_moebius(object_, embedded=False):
"""Alias for :py:func:`.moebius_models.euclidean_to_projective`"""
return moebius_models.euclidean_to_projective(object_, embedded=embedded)
[docs]def moebius_to_projective(object_, embedded=False):
"""Alias for :py:func:`.moebius_models.projective_to_euclidean`"""
return moebius_models.projective_to_euclidean(object_, embedded=embedded)