"""Module containing projections, e.g. central and stereographic.
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)).
"""
from functools import singledispatch
import numpy as np
import ddg.geometry.spheres
from ddg import nonexact
from ddg.geometry import moebius_models
from ddg.geometry._intersection import join, meet
from ddg.geometry._quadrics import Quadric, touching_cone
from ddg.geometry._subspaces import Point, Subspace, coordinate_hyperplane
from ddg.geometry.spheres import SphereLike
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])
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))
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
-------
Quadric
"""
return meet(quadric, join(subspace, center))
[docs]def central_project(object_, target, center):
"""
Project the `object_` onto the `target`
with respect to given `center` of projection.
Currently supported types of objects and targets for central projection:
* Subspace onto Subspace
* Quadric onto Subspace
* Subspace onto Quadric
| The projection is obtained as the intersection
of `target` and `join(object_, center)`.
| Thus, in general case, a point of the `object`
is mapped to a subspace or to a quadric.
| Such projection maps points to points when
`center` and `target` are disjoint Subspaces
with the sum `center.dimension + target.dimension + 1`
equal to the dimension of their common ambient space.
| If, in addition, the `object_` is a Subspace,
the function `central_project` returns the Subspace
spanned by the images of `object_.points`.
Parameters
----------
object_: ddg.geometry.Subspace or ddg.geometry.Quadric
Object to be projected.
target: ddg.geometry.Subspace or ddg.geometry.Quadric
Subspace or Quadric to be projected onto.
center: ddg.geometry.Subspace
| The center of projection.
| The most common example of `center` is a Point
outside the `target` Subspace.
However, the dimension of `center` may also be positive.
Returns
-------
ddg.geometry.Subspace or ddg.geometry.Quadric
Raises
------
NotImplementedError
If the pair of given objects' types is currently not supported.
If a quadric is given as `object_` and is not contained in a proper subspace.
"""
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 "
"is not implemented (yet)."
)
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))
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))
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):
"""
| Take the contour of image of `object_` under projection onto the `target`
with respect to given `center` of projection.
| In other words, the function `central_project_contour` returns the boundary
of the shadow (obtained via `central_project` with the same input)
instead of the shadow itself.
Currently supported types of objects and targets
for central projection in contour mode:
* Subspace onto Subspace
* Quadric onto Subspace
* Subspace onto Quadric
The contour of a proper subspace is considered to be the proper subspace itself.
The contour of the whole space is the empty subspace.
| When projecting the contour of a quadric contained in the whole space
(i.e. not contained in any proper subspace), the center must be a point,
but not as a subspace of larger dimension.
| In this case the function `central_project_contour` returns the intersection
of `target` and `touching_cone` to given quadric.
Parameters
----------
object_: ddg.geometry.Subspace or ddg.geometry.Quadric
Object to be projected.
target: ddg.geometry.Subspace or ddg.geometry.Quadric
Subspace or Quadric to be projected onto.
center: ddg.geometry.Subspace
| The center of projection.
| If `object_` is Quadric whose `subspace` attribute is the whole space,
then `center` must be a Point (i.e. a Subspace of zero dimension).
Returns
-------
ddg.geometry.Subspace or ddg.geometry.Quadric
If the image of `object_` has no boundary
under central projection onto `target`,
the function `central_project_contour` returns the empty object of
expected type (i.e. Subspace or Quadric) in expected subspace.
Raises
------
NotImplementedError
If the pair of given objects' types is currently not supported.
If `object_` given as a quadric not contained
in a proper subspace and `center` is not a point.
"""
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)} "
"objects is not implemented (yet)."
)
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):
"""
| The complex mode of `central_project` function.
| Every real projective subspace or quadric can be treated
as a slice of complex subspace or quadric.
The `central_project_complex` implementation
is equivalent the following procedure:
#. consider complexified `object_`, `target` and `center`;
#. take the image of complexified `object_` under central projection
with respect to complexified `center` onto complexified `target`;
#. bring the image back to real space by forgetting its complex dimensions.
Currently supported types of objects and targets
for central projection in complex mode:
* Quadric onto Subspace
While projecting a quadric onto a subspace, the result will always be a subspace:
either the whole `target` subspace
or its intersection with the subspace of `object_` quadric.
Parameters
----------
object_ : ddg.geometry.Quadric
Quadric to project.
target : ddg.geometry.Subspace
Subspace to be projected onto.
center : ddg.geometry.Subspace
The center of projection.
Returns
-------
ddg.geometry.Subspace
Raises
------
NotImplementedError
If the pair of given objects' types is currently not supported.
"""
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)} "
"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.
Currently supported types of objects to be projected:
* Quadric
* Subspace
Parameters
----------
object_ : `ddg.geometry.Quadric` or `ddg.geometry.Subspace`
| Object (Quadric or Subspace) to be projected.
| Note that the object does not have to be contained by the Quadric
from which the stereographic projection is implemented,
since the image defined just by `projection_point` and
target `hyperplane`.
| If given the Subspace coinciding with `projection_point`,
then `np.inf` is returned.
hyperplane: `ddg.geometry.Subspace` (default=None)
Hyperplane not containing `projection_point` to project onto.
Default is the plane spanned by the homogeneous vectors
e_1,...,e_n-1, e_{n+1}.
projection_point: `ddg.geometry.Point` (default=None)
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
`ddg.geometry.Quadric` (if Quadric is given)
or `ddg.geometry.Subspace` or `float` (if Subspace is given)
Notes
-----
In order to obtain `ddg.geometry.spheres.SphereLike` from a Moebius sphere
use `ddg.geometry.moebius_models.projective_to_euclidean`
or `ddg.geometry.euclidean_models.moebius_to_projective`.
Raises
------
NotImplementedError
If given type of object is currently not supported.
"""
raise NotImplementedError(
f"Stereographic projection of objects of type {type(object_)} is not "
"implemented (yet)."
)
@stereographic_project.register(Subspace)
def stereographic_project_subspace(subspace, hyperplane=None, projection_point=None):
"""
Parameters
----------
subspace : ddg.geometry.Subspace
Subspace to be projected.
If this is the projection point itself, returns `inf`.
hyperplane : ddg.geometry.Subspace (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 : ddg.geometry.Point or None (default=None)
Point to project from. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
ddg.geometry.Subspace or float
"""
n = subspace.ambient_dimension
if projection_point is None:
projection_point = north_pole(n)
if hyperplane is None:
hyperplane = coordinate_hyperplane(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
)
@stereographic_project.register(Quadric)
def stereographic_project_quadric(quadric, hyperplane=None, projection_point=None):
"""
Parameters
----------
quadric : ddg.geometry.Quadric
Quadric to be projected.
hyperplane : ddg.geometry.Subspace (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 : ddg.geometry.Point (default=None)
Point to project from. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
ddg.geometry.Quadric
"""
n = quadric.ambient_dimension
if projection_point is None:
projection_point = north_pole(n)
if hyperplane is None:
hyperplane = coordinate_hyperplane(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.
Currently supported types of objects to be projected:
* Point
* Subspace
* Sphere
Parameters
----------
object_ : `ddg.geometry.Subspace` or `ddg.geometry.spheres.SphereLike`
| Object (Point, Subspace or SphereLike) to be projected.
| The cases Point and Subspace are differentiated because of
different types of output (Point and Quadric respectively).
quadric: `ddg.geometry.Quadric` (default=None)
Quadric to project onto. Default is the Moebius quadric.
projection_point: `ddg.geometry.Point` (default=None)
Center of projection. Must lie on `quadric`.
Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
ddg.geometry.Point (if Point is given) or ddg.geometry.Quadric (otherwise)
Raises
------
NotImplementedError
If given type of object is currently not supported.
If `SphereLike` is given and
* the intersection of `sphere.subspace` and the tangent plane at
`projection_point` is not at infinity.
* the induced metric is distorted, i.e. planar sections of the
`quadric` would not map to spheres in `sphere.subspace`.
* 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`.
"""
raise NotImplementedError(
f"Inverse stereographic projection of objects of type {type(object_)} "
"is not implemented (yet)."
)
@inverse_stereographic_project.register(Point)
def inverse_stereographic_project_point(point, quadric=None, projection_point=None):
"""
Parameters
----------
point : ddg.geometry.Point
Point to project.
quadric : ddg.geometry.Quadric (default=None)
Quadric to project onto. Default is the Möbius quadric.
projection_point : ddg.geometry.Point
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
ddg.geometry.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 = moebius_models.ProjectiveModel(n - 1).absolute
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)
@inverse_stereographic_project.register(Subspace)
def inverse_stereographic_project_subspace(
subspace, quadric=None, projection_point=None
):
"""
Parameters
----------
subspace : ddg.geometry.Subspace
Subspace to project.
quadric : ddg.geometry.Quadric or None (default=None)
Quadric to project onto. Default is the Möbius quadric.
projection_point : ddg.geometry.Point
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
ddg.geometry.Quadric
Raises
------
ValueError
If `projection_point` does not lie on `quadric`.
"""
n = subspace.ambient_dimension
if quadric is None:
quadric = moebius_models.ProjectiveModel(n - 1).absolute
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)
@inverse_stereographic_project.register(SphereLike)
def inverse_stereographic_project_sphere(sphere, quadric=None, projection_point=None):
"""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 : Euclidean sphere
quadric : ddg.geometry.Quadric(default=None)
projection_point : ddg.geometry.Point (default=None)
Center of projection. Default is the "north pole" [e_n+e_{n+1}].
Returns
-------
ddg.geometry.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.
"""
n = sphere.ambient_dimension
if quadric is None:
quadric = moebius_models.ProjectiveModel(n - 1).absolute
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():
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, rtol = nonexact.combine_tols(sphere, quadric)
if not nonexact.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():
raise NotImplementedError("Origin is at infinity.")
O = O.dehomogenize().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 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 : ddg.geometry.spheres.CayleyKleinSphere
quadric : ddg.geometry.Quadric
Quadric to lift to.
projection_point : ddg.geometry.Point
A point not on the quadric.
Returns
-------
ddg.geometry.Quadric, ddg.geometry.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."
)
c = sphere.center.point
b = quadric.inner_product
q = projection_point.point
if isinstance(sphere, ddg.geometry.spheres.GeneralizedCayleyKleinSphere):
# This is for lifting horospheres
mu = sphere.generalized_radius()
else:
mu = sphere.cayley_klein_radius()
# For Cayley-Klein spheres, if the radius is < 0, the sphere and its center
# are on different sides of the absolute. So if the center is outside the
# projection image and the radius is non-negative, the sphere also lies
# outside the image and the lift is empty.
# For mu==0, the "sphere" is just the polar hyperplane and the lift still
# works.
preimage = central_project(sphere.center, quadric, projection_point)
if preimage.signature().is_definite and mu > 0:
# We could also return some empty quadric.
raise ValueError(
"Sphere must lie in the central projection image of the quadric "
"via projection_point onto quadric.polarize(projection_point)."
)
if isinstance(sphere, ddg.geometry.spheres.GeneralizedCayleyKleinSphere):
polar1 = Point(c + np.sqrt(-mu / b(q, q)) * q)
polar2 = Point(c - np.sqrt(-mu / b(q, q)) * q)
else:
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)
subspace = join(sphere.subspace, projection_point)
Q1 = meet(quadric.polarize(polar1), quadric, subspace)
Q2 = meet(quadric.polarize(polar2), quadric, subspace)
return Q1, Q2
def north_pole(n):
""" "North pole" in n-dimensional space.
ddg.geometry.Point(en + en+1)
"""
p = np.zeros(n + 1)
p[[-1, -2]] = 1
return Point(p)