"""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
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.spheres import SphereLike
from ddg.geometry.subspaces import Point, Subspace, coordinate_hyperplane
[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 "
"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)} "
"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)} "
"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 "
"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 = 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
)
[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 = 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.
"""
raise NotImplementedError(
f"Inverse stereographic projection of objects of type {type(object_)} "
"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 = 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)
[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 = 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)
[docs]@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 : 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}].
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.
"""
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 : CayleyKleinSphere
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."
)
c = sphere.center.point
b = quadric.inner_product
q = projection_point.point
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)."
)
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
[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)