"""Conversion module for geometries.
This module contains functions for conversion between different models of a
Geometry as well as functions that deal with geometric objects of varying type,
e.g. projection functions.
A word on projection functions: Central projection is usually defined as a map
between two k-planes S1, S2 in RP\\ :sup:`n` via another, (n-k-1)-dim. subspace C
disjoint from S1, S2. The projection maps x in S1 to meet(S2, join(x, S1)) and
this is a well-defined bijective map and a projective transformation. In this
module, we take a more general approach: If objects Y and C are given, we map
any object X to meet(Y, join(X, C)).
"""
import numpy as np
from functools import singledispatch
import ddg
from ddg.geometry.geometries import _KNOWN_GEOMETRIES, get_geometry, Geometry
from ddg.geometry.subspaces import Subspace, Point, normal_with_level,\
subspace_from_columns
from ddg.geometry.intersection import join, meet
from ddg.geometry.quadrics import Quadric, intersect_quadric_subspace,\
cayley_klein_sphere, touching_cone
from ddg.geometry.spheres import Sphere
import ddg.math.linalg as la
import ddg.math.projective as pmath
import ddg.math.symmetric_matrices as sm
from ddg.abc import NonExact
_error_msg = 'Could not determine geometry. Neither source geometry was given,\
nor does the object have any associated geometry.'
[docs]def convert(obj, source_geometry, target_geometry):
"""General convert function for geometries.
Parameters
----------
obj : ddg Object
Object to convert to target geometry
target_geometry : ddg.Geometry or string
Target geometry. When give as a string, the dimension has to be
included. (This might change at some point)
source_geometry : ddg.Geometry (optional, default=None)
When given, the object will be interpreted as being of this geometry.
Returns
-------
ddg Object
converted object
"""
#We have to make sure that we get our standard names for the geometries
if source_geometry is None:
if not hasattr(obj, 'geometry') or obj.geometry is None:
raise Exception(_error_msg)
source_geometry = obj.geometry
elif not isinstance(source_geometry, Geometry):
if isinstance(source_geometry, str):
source_geometry = get_geometry(source_geometry)
else:
source_geometry = get_geometry(*source_geometry)
if not isinstance(target_geometry, Geometry):
if isinstance(target_geometry, str):
target_geometry = get_geometry(target_geometry)
else:
target_geometry = get_geometry(*target_geometry)
source_geometry = source_geometry.name
target_geometry = target_geometry.name
#Remove dual if it occurs. Conversion helper function have to handle this
#anyway
if 'dual ' in target_geometry:
target_geometry = target_geometry.replace('dual ', '')
elif 'dual' in target_geometry:
target_geometry = target_geometry.replace('dual', '')
if 'dual ' in source_geometry:
source_geometry = source_geometry.replace('dual ', '')
elif 'dual' in source_geometry:
source_geometry = source_geometry.replace('dual', '')
#Create expected function name
fct_name = '_' + ''.join(source_geometry.split())
fct_name += '_to_'
fct_name += ''.join(target_geometry.split())
if hasattr(ddg.geometry.conversion, fct_name):
fct = getattr(ddg.geometry.conversion, fct_name)
else:
#Function was not found. Find alternative name and try again
alt1 = None
alt2 = None
for key in _KNOWN_GEOMETRIES:
if key in source_geometry:
alt1 = key
if key in target_geometry:
alt2 = key
if (alt1 != None) and (alt2 != None):
break
fct_name = '_' + alt1 + '_to_' + alt2
if hasattr(ddg.geometry.conversion, fct_name):
fct = getattr(ddg.geometry.conversion, fct_name)
else:
raise NotImplementedError('Could not find conversion function.')
return fct(obj)
def _euclidean_to_laguerre(obj):
if obj.codimension == 1:
return Point(np.array([*normal_with_level(obj), 1]))
raise NotImplementedError
def _laguerre_to_euclidean(obj):
if isinstance(obj, Point):
p = obj.point[:-1]
return Point(p).dualize()
raise NotImplementedError
def _euclidean_to_moebius(object_, subspace=None, projection_point=None,
affine_component=-1):
"""Convert from Euclidean to Möbius geometry.
This function will embed `object_` in `subspace` and then stereographically
project onto the Möbius quadric via `projection_point`.
Parameters
----------
object_ : Subspace or Sphere
Objects in n-dimensional ambient space. All homogeneous coordinate
vectors are interpreted as coordinates w.r.t. the basis
`subspace.points`.
subspace : Subspace or None (default=None)
Subspace of dimension at most n in (n+1)-dimensional ambient space.
Default is the equatorial plane e_n+1 = 0.
projection_point : Point or None (default=None)
Point on the Möbius quadric acting as center of projection. Default is
the "north pole" e_n+1 + e_n+2.
affine_component : int (default=-1)
Affine component used in
:py:func:`.inverse_stereographic_project_sphere`.
Returns
-------
Quadric
An intersection of the Möbius quadric with a subspace.
Raises
------
ValueError
If `subspace.ambient_dimension != object_.ambient_dimension + 1`
TypeError
If object_ type is not supported.
"""
n = object_.ambient_dimension
# projection_point is handled by inverse_stereographic_project_sphere
# subspace would also be handled by this function, except that we need to
# know it to embed into it.
quadric = get_geometry('moebius', n+1).quadric
if subspace is None:
subspace = equatorial_plane(n+1)
elif subspace.ambient_dimension != n+1:
raise ValueError("Dimension mismatch: subspace.ambient dimension must "
"be equal to object_.ambient_dimension + 1.")
B = subspace.matrix
if isinstance(object_, Subspace):
object_embedded = subspace_from_columns(B @ object_.matrix)
return inverse_stereographic_project(object_embedded, quadric,
projection_point)
if isinstance(object_, Sphere):
subspace_embedded = subspace_from_columns(B @ object_.subspace.matrix)
center_embedded = Point(B @ object_.center.point)
object_embedded = Sphere(center_embedded, object_.radius,
subspace_embedded)
return inverse_stereographic_project(
object_embedded, quadric, projection_point,
affine_component=affine_component
)
raise TypeError(f"Object of type {type(object_)} can't be converted "
f"from Euclidean to Möbius goemetry.")
def _moebius_to_euclidean(object_, subspace=None, projection_point=None):
"""Convert a Möbius sphere to a Euclidean sphere or subspace.
Spheres will be represented as quadrics since they can look like ellipses
depending on the choice of `projection_point` and `subspace`.
This function stereographically projects and then returns the object
expressed in the given basis of `subspace`. More precisely, the embedded
version of a returned subspace `proj_obj` would be given by ::
subspace_from_columns(subspace.matrix @ proj_obj.matrix)
and the embedded version of a returned quadric `proj_obj` would be given by
::
Quadric(
proj_obj.matrix,
subspace_from_columns(
subspace.matrix @ proj_obj.subspace.matrix
)
)
Parameters
----------
object_ : Quadric
Intersection of Möbius quadric with a subspace, all in (n+1)-dim.
ambient space.
subspace : Subspace
Subspace to project onto. The returned object will be expressed in the
coordinates of this subspace basis.
projection_point : Point
Point on the quadric acting as center of projection. Default is the
"north pole" e_n+1 + e_n+2.
Returns
-------
Quadric or Subspace
Object in n-dim. ambient space.
"""
n = object_.ambient_dimension
if subspace is None:
subspace = equatorial_plane(n)
img_embedded = stereographic_project(object_, subspace,
projection_point)
if isinstance(img_embedded, Subspace):
return subspace_from_columns(
la.coordinates(img_embedded.matrix, subspace.matrix)
)
if isinstance(img_embedded, Quadric):
A = la.coordinates(img_embedded.subspace.matrix, subspace.matrix)
return Quadric(img_embedded.matrix, subspace=A.T)
raise Exception("Something went wrong (stereographic projection returned "
"neither a quadric nor a subspace).")
######################## Projections and lifts ################################
[docs]def central_project_subspace_onto_subspace(subspace, target_subspace, center):
"""
Parameters
----------
subspace : Subspace
If possible, the basis of the resulting subspace will consist of the
projections of the basis points of `subspace`. For this to work,
target_subspace.dimension+center.dimension must be equal to
ambient_dimension-1 and they must be disjoint. Otherwise, this function
simply returns meet(target_subspace, join(subspace, center)).
target_subspace : Subspace
center : Subspace
Returns
-------
Subspace
"""
# Dimension requirement for standard point-by-point central projection is
# not fulfilled, so we can't project basis points individually. But we can
# still project the subspace as a whole.
if center.codimension != target_subspace.dimension + 1:
return meet(target_subspace, join(subspace, center))
image_points = []
for p in subspace.points:
p = Point(p)
image_point = meet(target_subspace, join(p, center))
image_points.append(image_point)
return Subspace(*[p.point for p in image_points])
[docs]def central_project_quadric_onto_subspace(quadric, subspace, center):
"""
Parameters
----------
quadric : Quadric
target_subspace : Subspace
center : Subspace
Returns
-------
Quadric
Raises
------
NotImplementedError
If quadric is not contained in a subspace.
"""
if quadric.subspace.codimension == 0:
raise NotImplementedError(
"Can't project a quadric in the whole space onto a subspace. Use "
"central_project_contour to do this."
)
return meet(subspace, join(center, quadric))
[docs]def central_project_subspace_onto_quadric(subspace, quadric, center):
"""
Parameters
----------
subspace : Subspace
Subspace to be projected.
quadric : Quadric
Quadric to project onto.
center : Subspace
Projection center.
Returns
-------
Subspace
"""
return meet(quadric, join(subspace, center))
[docs]def central_project(object_, target, center):
"""Central projection dispatcher.
"""
if isinstance(object_, Subspace) and isinstance(target, Subspace):
return central_project_subspace_onto_subspace(object_, target, center)
if isinstance(object_, Quadric) and isinstance(target, Subspace):
return central_project_quadric_onto_subspace(object_, target, center)
if isinstance(object_, Subspace) and isinstance(target, Quadric):
return central_project_subspace_onto_quadric(object_, target, center)
raise NotImplementedError(
f"Projection of {type(object_)} objects onto {type(target)} objects "
f"is not implemented (yet)."
)
[docs]def central_project_subspace_onto_subspace_contour(subspace, target_subspace,
center):
"""Central project subspace onto subspace (contour mode).
If image is a proper subspace of `target_subspace`, returns the regular
central projection. If image is all of target_subspace, returns the empty
subspace.
Parameters
----------
subspace, target_subspace, center : Subspace
Returns
-------
Subspace
"""
image = central_project_subspace_onto_subspace(subspace, target_subspace,
center)
if image.dimension < target_subspace.dimension:
return image
return Subspace(np.zeros(subspace.ambient_dimension + 1))
[docs]def central_project_quadric_onto_subspace_contour(quadric, subspace, center):
"""Central project quadric onto subspace (contour version).
This function also supports quadrics in the whole space.
Parameters
----------
quadric : Quadric
subspace : Subspace
center : Point
Returns
-------
Quadric
Raises
------
NotImplementedError
If `quadric` is in the whole space and `center` is not a point.
"""
if quadric.subspace.codimension > 0:
return central_project_quadric_onto_subspace(quadric, subspace, center)
if not isinstance(center, Point):
raise NotImplementedError(
"Quadrics in the whole space must be projected from a point "
"(because touching_cone only supports points)."
)
return meet(subspace, touching_cone(center, quadric))
[docs]def central_project_subspace_onto_quadric_contour(subspace, quadric, center):
"""Central project subspace onto quadric (contour mode).
If image is a proper subset of `quadric`, returns image. If image is all of
`quadric`, returns an empty quadric.
Parameters
----------
subspace : Subspace
quadric : Quadric
center : Subspace
Returns
-------
Quadric
"""
image = central_project_subspace_onto_quadric(subspace, quadric, center)
if image.dimension() < quadric.dimension():
return image
return Quadric(np.eye(quadric.subspace.dimension + 1),
subspace=quadric.subspace)
[docs]def central_project_contour(object_, target, center):
"""Central projection dispatcher (contour mode).
This returns the boundary of the shadow instead of the shadow itself, so to
speak.
"""
if isinstance(object_, Subspace) and isinstance(target, Subspace):
return central_project_subspace_onto_subspace_contour(object_, target,
center)
if isinstance(object_, Quadric) and isinstance(target, Subspace):
return central_project_quadric_onto_subspace_contour(object_, target,
center)
if isinstance(object_, Subspace) and isinstance(target, Quadric):
return central_project_subspace_onto_quadric_contour(object_, target,
center)
raise NotImplementedError(
f"Contour projection of {type(object_)} objects onto {type(target)} "
f"objects is not implemented (yet)."
)
[docs]def central_project_quadric_onto_subspace_complex(quadric, subspace, center):
"""Central project quadric onto subspace (complex version).
If we also consider complex intersections of lines through `center` with
`quadric`, the projection image will always be all of subspace (or the
intersection of `subspace` with `quadric.subspace`).
Parameters
----------
quadric : Quadric
subspace : Subspace
center : Subspace
Returns
-------
Subspace
"""
if quadric.subspace.codimension == 0:
return subspace
return meet(subspace, quadric.subspace)
[docs]def central_project_complex(object_, target, center):
"""Central projection dispatcher (complex mode).
"""
if isinstance(object_, Quadric) and isinstance(target, Subspace):
return central_project_quadric_onto_subspace_complex(object_, target,
center)
raise NotImplementedError(
f"Complex projection of {type(object_)} objects onto {type(target)} "
f"objects is not implemented (yet)."
)
[docs]@singledispatch
def stereographic_project(object_):
"""Stereographic projection (Quadric to subspace) dispatcher.
Accepts different keyword arguments based on type. See the type-specific
functions for details.
"""
raise NotImplementedError(
f"Stereographic projection of objects of type {type(object_)} is not "
f"implemented (yet)."
)
[docs]@stereographic_project.register(Subspace)
def stereographic_project_subspace(subspace, hyperplane=None,
projection_point=None):
"""
Parameters
----------
subspace : Subspace
Subspace to be projected.
If this is the projection point itself, returns `inf`.
hyperplane : Subspace or None (default=None)
Hyperplane not containing `projection_point`.
Default is the plane spanned by the homogeneous vectors
e_1,...,e_n-1, e_{n+1}.
projection_point : Point or None (default=None)
Point to project from. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
Subspace or float
"""
n = subspace.ambient_dimension
if projection_point is None:
projection_point = north_pole(n)
if hyperplane is None:
hyperplane = equatorial_plane(n)
# TODO Is this a good choice? I could also imagine it being the whole line
# at infinity for example.
if subspace == projection_point:
return np.inf
return central_project_subspace_onto_subspace(
subspace,
target_subspace=hyperplane,
center=projection_point
)
[docs]@stereographic_project.register(Quadric)
def stereographic_project_quadric(quadric, hyperplane=None,
projection_point=None):
"""
Parameters
----------
quadric : Quadric
Quadric to be projected.
hyperplane : Subspace or None (default=None)
Hyperplane not containing `projection_point`.
Default is the plane spanned by the homogeneous vectors
e_1,...,e_n-1, e_{n+1}.
projection_point : Point or None (default=None)
Point to project from. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
Quadric
"""
n = quadric.ambient_dimension
if projection_point is None:
projection_point = north_pole(n)
if hyperplane is None:
hyperplane = equatorial_plane(n)
return central_project_quadric_onto_subspace(
quadric,
subspace=hyperplane,
center=projection_point
)
[docs]@singledispatch
def inverse_stereographic_project(object_):
"""Inverse stereographic projection (Subspace to Quadric) dispatcher.
Accepts different keyword arguments based on type. See the type-specific
functions for details.
"""
raise NotImplementedError(
f"Inverse stereographic projection of objects of type {type(object_)} "
f"is not implemented (yet)."
)
[docs]@inverse_stereographic_project.register(Point)
def inverse_stereographic_project_point(point, quadric=None,
projection_point=None):
"""
Parameters
----------
point : Point
Point to project.
quadric : Quadric or None (default=None)
Quadric to project onto. Default is the Möbius quadric.
projection_point : Point
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
Point
Notes
-----
This is separate from inverse_stereographic_project_subspace because in
this case there is an explicit formula and because of some special cases.
Raises
------
ValueError
If `projection_point` does not lie on `quadric`.
"""
n = point.ambient_dimension
if quadric is None:
quadric = get_geometry('moebius', n).quadric
if projection_point is None:
projection_point = north_pole(n)
if projection_point not in quadric:
raise ValueError("projection_point must lie on quadric.")
if quadric.conjugate(projection_point, point):
if point in quadric:
# The whole line is contained in the quadric.
return join(point, projection_point)
# point is "at infinity", no regular inverse image.
return projection_point
b = quadric.inner_product
c = projection_point.point
p = point.point
return Point(-b(p, p) / (2 * b(c, p)) * c + p)
[docs]@inverse_stereographic_project.register(Subspace)
def inverse_stereographic_project_subspace(subspace, quadric=None,
projection_point=None):
"""
Parameters
----------
subspace : Subspace
Subspace to project.
quadric : Quadric or None (default=None)
Quadric to project onto. Default is the Möbius quadric.
projection_point : Point
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
Quadric
Raises
------
ValueError
If `projection_point` does not lie on `quadric`.
"""
n = subspace.ambient_dimension
if quadric is None:
quadric = get_geometry('moebius', n).quadric
if projection_point is None:
projection_point = north_pole(n)
if projection_point not in quadric:
raise ValueError("projection_point must lie on quadric.")
return central_project_subspace_onto_quadric(subspace, quadric,
projection_point)
[docs]@inverse_stereographic_project.register(Sphere)
def inverse_stereographic_project_sphere(sphere, quadric=None,
projection_point=None,
affine_component=-1):
"""Stereographically project a sphere contained in a subspace to a quadric.
Works in all cases in which planar sections of `quadric` are
stereographically projected to spheres (in the sense of the Euclidean
metric in affine coordinates) in sphere.subspace.
Parameters
----------
sphere : Sphere
quadric : Quadric or None (default=None)
projection_point : Point or None (default=None)
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
affine_component : int
Returns
-------
Quadric
Raises
------
NotImplementedError
* If the intersection of sphere.subspace and the tangent plane at
projection_point is not at infinity.
* If the induced metric is distorted, i.e. planar sections of the
quadric would not map to spheres in sphere.subspace.
* If the intersection of the line through projection_point and
polarize(sphere.subspace) with sphere.subspace is at infinity. This
point would normally act as a sort of origin.
ValueError
If `projection_point` does not lie on `quadric`.
Notes
-----
This function uses the euclidean distance on a quadric of signature
(n+1, 1) given by::
d([x], [y])^2 = -1/2 * <x,y>/(<x,p><y,p>)
for a point [p] on the quadric with a fixed representative vector. We use
that formula to find all points a constant distance from a point. This set
will be a planar section of the quadric.
This metric corresponds to the Euclidean metric in affine coordinates after
stereographic projection, but only to certain subspaces via certain points.
We have to determine whether we have such a case and then find the right
representative vector `p` to get the right scaling.
See :ref:`supplemental material <euclidean_moebius>` for details.
"""
i = affine_component
n = sphere.ambient_dimension
if quadric is None:
quadric = get_geometry('moebius', n).quadric
if projection_point is None:
projection_point = north_pole(n)
if projection_point not in quadric:
raise ValueError("Projection point must be on the quadric.")
l_inf = meet(quadric.polarize(projection_point), sphere.subspace)
# It's possible that this is not necessary: If l_inf is not at infinity,
# the metric will probably also be distorted, which we check below anyway.
# But let's be safe.
if not l_inf.at_infinity(i):
raise NotImplementedError(
"The intersection of sphere.subspace and the tangent plane at "
"projection_point is not at infinity."
)
Q_restricted = meet(l_inf, quadric)
# The columns of this matrix are an orthonormal basis of l_inf w.r.t.
# Q.inner_product. We want to check if they are a scaled orthonormal basis
# of l_inf w.r.t. the standard Euclidean inner product as well.
A = Q_restricted.normalize().subspace.matrix
A = A.T @ A
l = A[0, 0]
atol = max([sphere.atol, quadric.atol])
rtol = max([sphere.rtol, quadric.rtol])
if not np.allclose(A / l, np.eye(A.shape[0]), atol=atol, rtol=rtol):
raise NotImplementedError(
"The image of this sphere on the quadric would not be a planar "
"section."
)
# If B is any hyperplane containing sphere.subspace, b is its pole and we
# centrally project b to B via projection_point, we always get the same
# point, independently of the hyperplane. Namely, we get this O.
O = central_project(quadric.polarize(sphere.subspace), sphere.subspace,
projection_point)
if O.at_infinity(i):
raise NotImplementedError("Origin is at infinity.")
O = O.dehomogenize(i).point
e_inf = projection_point.point
b = quadric.inner_product
# This makes it so that b(O, e_inf) = +/- 1/2.
e_inf /= 2 * b(O, e_inf)
# Currently, "vectors" that have Euclidean length l have Q-length 1, so to
# match, we have to scale the Q-metric by l.
# Scaling e_inf by a scales the metric by (1/a)^2.
e_inf /= np.sqrt(l)
c = inverse_stereographic_project(sphere.center, quadric, projection_point)
c = c.point
r = sphere.radius
pole = Point(b(c, e_inf) * r * r * e_inf + c / 2)
polar = quadric.polarize(pole)
proj_sphere = meet(polar, quadric)
if sphere.subspace.codimension > 1:
proj_sphere = meet(
proj_sphere,
join(sphere.subspace, projection_point)
)
return proj_sphere
[docs]def north_pole(n):
""" "North pole" in n-dimensional space.
Point(en + en+1)
"""
p = np.zeros(n + 1)
p[[-1, -2]] = 1
return Point(p)
[docs]def equatorial_plane(n):
""" "Equatorial plane" in n-dimensional space.
Subspace(e1,...,en-1, en+1)
"""
S = np.eye(n + 1)
S[[-1, -2]] = S[[-2, -1]]
return Subspace(*S[:-1])
[docs]def lift_sphere_to_quadric(sphere, quadric, projection_point):
"""Find the lift of a sphere to a quadric via a point.
This is the inverse of the following map: Map a planar section of `quadric`
to `quadric.polarize(projection_point)` via `projection_point`, giving a
Cayley-Klein sphere.
Because the inverse map (the projection) is a double cover, this will
always return two quadrics.
Parameters
----------
sphere : Sphere
Sphere. Must have a geometry that allows computation of Cayley-Klein
"distance" from metric distance. Its center must lie in the polar plane
of `projection_point`.
quadric : Quadric
Quadric to lift to.
projection_point : Point
A point not on the quadric.
Returns
-------
Quadric, Quadric
Raises
------
ValueError
* If projection_point is on the quadric
* If center of sphere does not lie in
`quadric.polarize(projection_point)`.
* If no preimage exists. This is the case if the center of the sphere
does not lie in the projection image of the quadric.
Notes
-----
Taken from `Non-Euclidean Laguerre Geometry and Incircular Nets`, page 43.
"""
if projection_point in quadric:
raise ValueError("projection_point must not lie in quadric.")
if not quadric.conjugate(sphere.center, projection_point):
raise ValueError(
"Center of sphere must lie in polar plane of projection_point."
)
# For Cayley-Klein spheres, if the radius is < 0, the sphere and its center
# are on different sides of the absolute. Because our spheres are given in
# terms of the metric, this can not happen.
preimage = central_project(sphere.center, quadric, projection_point)
if preimage.signature().is_definite:
# We could also return some empty quadric.
raise ValueError(
"Center of sphere must lie in the central projection image of the "
"quadric via projection_point onto "
"quadric.polarize(projection_point)."
)
c = sphere.center.point
b = quadric.inner_product
q = projection_point.point
mu = sphere.cayley_klein_radius()
polar1 = Point(c + np.sqrt(-mu * b(c, c) / b(q, q)) * q)
polar2 = Point(c - np.sqrt(-mu * b(c, c) / b(q, q)) * q)
Q1 = meet(quadric.polarize(polar1), quadric)
Q2 = meet(quadric.polarize(polar2), quadric)
return Q1, Q2
### Type conversions
[docs]@NonExact.nonexact_function
def quadric_to_euclidean_sphere(quadric, affine_component=-1, atol=None,
rtol=None):
"""Convert quadric to a Euclidean sphere, if it is one.
Parameters
----------
quadric : Quadric
affine_component : int (default=-1)
atol, rtol : float (default=None)
If None is given, the global defaults are used. See ddg.abc.NonExact
for details.
Returns
-------
Sphere
Raises
------
ValueError
If quadric is not a sphere.
"""
ac = affine_component
# Check signature
sgn = quadric.signature(affine=True, affine_component=ac)
k = quadric.subspace.dimension
if sgn != sm.AffineSignature(k, 1, affine_component_entry=-1):
raise ValueError(
f"Signature of quadric is {sgn}, but has to be (k, 1) with -1 at "
f"position affine_component for it to be a sphere."
)
Q = quadric.normalize(affine=True, affine_component=ac)
B = Q.subspace.matrix
# The unit sphere in a standard subspace is mapped to the quadric by an
# affine transformation. It is a sphere iff the transformation is a
# similarity.
B /= B[ac, -1]
A = np.delete(B, ac, axis=0)
A = np.delete(A, -1, axis=1)
I = A.T @ A
tmp = I / I[0, 0]
if not np.allclose(tmp, np.eye(tmp.shape[0]), atol=atol, rtol=rtol):
raise ValueError("Quadric is not a Euclidean sphere.")
c = B[:, -1]
r = np.sqrt(I[0, 0])
return Sphere(c, r, subspace=Q.subspace)
[docs]def euclidean_sphere_to_quadric(sphere, affine_component=-1):
"""Convert Euclidean sphere to quadric
Parameters
----------
sphere : Sphere
affine_component : int
Returns
-------
Quadric
"""
ac = affine_component
r = sphere.radius
c = sphere.center.affine_point(ac)
k = sphere.ambient_dimension
Q = np.eye(k + 1)
Q[ac, ac] = -1
Q = Quadric(Q)
F = pmath.affine_transformation(r * np.eye(k), c, affine_component=ac)
Q.transform(F)
if sphere.subspace.codimension > 0:
Q = intersect_quadric_subspace(Q, sphere.subspace)
return Q
[docs]def cayley_klein_sphere_to_quadric(sphere):
"""Convert Sphere in a Cayley-Klein geometry to quadric.
For this to work, the sphere's geometry must implement
metric_to_cayley_klein_distance.
Parameters
----------
sphere : Sphere
Returns
-------
quadric
"""
return cayley_klein_sphere(sphere.center, sphere.cayley_klein_radius(),
absolute=sphere.geometry.quadric)