Precision in Programming Möbius transformations of the Riemann Sphere

Question

Apologies if this is better asked on StackOverflow, but I thought the community here would be more likely to have an answer, even though it's related to coding.

I'm coding up a little library for experimenting with the inversive geometry of the sphere. A point on the sphere is represented by $(x, y, z)$ where $\sqrt{x^2 + y^2 + z^2} = 1$. I want to apply Möbius transformations to sets of points on the sphere. Obviously, I could accomplish this by first applying a stereographic projection of each point $(x, y, z)$ to a point $z$ in the complex plane, applying a Möbius transformation $z \mapsto \frac{az + b}{cz + d}$ to each point, and then projecting back onto the sphere.

The problem with this approach, is that because I'm using double precision floating point numbers for the coordinates, points that are near the South pole of the sphere will have large precision errors on projection (since they are "near" infinity), whereas points near the equator will experience a lot less of a precision issue.

So I'm wondering if there is a more natural method for applying Möbius transformations directly to the points $(x, y, z)$ rather than composing with stereographic projections?

Answer 1

If you interpret a point of the unit sphere (a.k.a., the Riemann sphere, or complex projective line) as the complex ordered pair $\left[\begin{array}{c}u\\ v\end{array}\right]$, then a Möbius transformation is induced by the complex linear transformation $$ \left[\begin{array}{c} u \\ v \\ \end{array}\right] \mapsto \left[\begin{array}{cc} a & b \\ c & d \\ \end{array}\right] \left[\begin{array}{c} u \\ v \\ \end{array}\right] = \left[\begin{array}{c} au + bv \\ cu + dv \\ \end{array}\right]. \tag{*} $$ The north pole corresponds to $u = 0$ and the south pole to $v = 0$.

If you're transforming a point $(x, y, z)$ with $0 \leq z \leq 1$, then $v$ is "far from $0$". Dividing (*) by $v/v$ expresses the Möbius transformation in terms of the numerically stable complex parameter $u/v$, for which your stereographic projection scheme works well.

If you're transforming a point $(x, y, z)$ with $-1 \leq z < 0$, then $u$ is "far from $0$". Dividing (*) by $u/u$ expresses the Möbius transformation in terms of the numerically stable complex parameter $v/u$, for which the complex conjugate of stereographic projection from the south pole works well.

Answer 2

My suggestion is to implement your first approach and test whether or not it exhibits the behavior you're afraid of! It will probably be fine.

points that are near the South pole of the sphere will have large precision errors on projection (since they are "near" infinity)

The double-precision format gives you full precision on large numbers up to $10^{308}$. You might eke out a few more orders of magnitude by playing with an alternate coordinate chart and exploiting subnormal numbers, but I doubt you'll see a qualitative difference.

Answer 3

I write the stereographic projection $\sigma:\ S^2\to\bar{\mathbb C}$ in the form $$\sigma:\quad (x,y,z)\mapsto\zeta:={x+iy\over 1-z}={1+z\over x-iy}\ .\tag{1}$$ At each of the poles $N$ and $S$ one of the expressions on the right hand side is ${0\over0}$, hence undefined.

In order to facilitate matters we introduce homogeneous coordinates $(\zeta_1,\zeta_2)$ in $\bar{\mathbb C}$, whereby $\zeta={\zeta_1\over\zeta_2}$. The map $\sigma$ then appears as $$\sigma:\quad(x,y,z)\mapsto(\zeta_1,\zeta_2):=\left\{\eqalign{&\ (x+iy,\ 1-z)\qquad(z<1)\cr &\ (1+z,\ x-iy)\qquad(z>-1)\cr}\right.\quad.\tag{2}$$ This means that away from $N$ you may use the upper formula, and away from $S$ you may use the lower formula.

In terms of homogeneous coordinates the Moebius transformation $T(\zeta):={a\zeta+b\over c\zeta+d}$ appears as linear map $$T:\quad(\zeta_1,\zeta_2)\mapsto\left\{\eqalign{\zeta_1'&:=a\zeta_1+b\zeta_2 \cr \zeta_2'&:=c\zeta_1+d\zeta_2 \cr}\right.\quad.$$ The inverse of $(1)$ is given by $$\sigma^{-1}:\quad\zeta\mapsto\left\{\eqalign{&x:={\zeta+\bar\zeta\over|\zeta^2|+1}\cr &y:={\zeta-\bar\zeta\over i(|\zeta|^2+1)}\cr &z:={|\zeta|^2-1\over|\zeta|^2+1}\cr}\right.\quad.$$ In terms of the homogeneous coordinates $(\zeta_1,\zeta_2)$ this becomes $$\sigma^{-1}:\quad\zeta\mapsto\left\{\eqalign{ &x:={\zeta_1\bar\zeta_2+\bar\zeta_1\zeta_2\over|\zeta_1|^2+|\zeta_2|^2}\cr &y:={\zeta_1\bar\zeta_2-\bar\zeta_1\zeta_2\over i(|\zeta_1|^2+|\zeta_2|^2)}\cr &z:= {|\zeta_1|^2-|\zeta_2|^2\over|\zeta_1|^2+|\zeta_2|^2}\cr}\right.\quad,$$ whereby it plays no rôle which of the formulas $(2)$ you have used in the first step.

Apart from the alternative in $(2)$ no extra measures have to be taken. As long as $a$, $b$, $c$, $d$ are constant there will be no problems of small denominators.

Answer 4

In addition to the answers already given to my question, here's another method that should work:

Select an even number of spheres, each of which is orthogonal to the unit sphere, and perform a sphere inversion on the points one after another. Such an inversion fixes the unit sphere, and is thus a circle inversion of the points of the unit sphere through the circle of intersection between the inversion sphere and the unit sphere. Thus, the composition of an even number of such inversions is a Möbius transformation of the sphere to itself.

Answer 5

The question asks whether it is possible to compute a Möbius transformation directly on the unit sphere $\Sigma$, without recourse to stereoscopic projection. But so far the upvoted and accepted answers resort to stereoscopic projection.

Let $\sigma$ be the stereoscopic projection $\Sigma\to \mathbb C \cup \infty$. A Möbius transformations $M:z \mapsto \frac{az + b}{cz + d}$ on the complex plane corresponds to a 3D homograpy (aka projective transformation) $T$ that leaves the unit sphere invariant. What is $T(p)$ for a given plane $p$ that intersects $\Sigma$?

$C=p\cap\Sigma$ is a circle. Since $M,\sigma$ and $\sigma^{-1}$ all map circles to circles, $C'=\sigma^{-1}M\sigma(C)$ is a circle on $\Sigma$. Then $T(p)$ is the plane $p'$ of $C'$.

The above skips some details, but granting that $T$ is a homography we can compute $T$ by the same methods that computer vision folks compute 2D homographies from 4 pairs of points. A very clear explanation can be found on this site at Finding the Transform matrix from 4 projected points (with Javascript). In 3D, of course, we need 5 projected points, and the matrices are 4x4.

$T$ now lets you map points $(x,y,z)\mapsto(x',y',z')$ without doing a computation on the plane.