THE STEREOGRAPHIC PROJECTION IN SAGE
JEFF BEORSE, KEVIN LINDEMAN & SIMON SPICER
WED JUNE 2, 2010
The stereographoc projection is a 1-1 mapping from the extended real plane $\mathbb{R}^2 \cup \{\infty\}$ to the unit sphere $S = \{(x,y,z) \in \mathbb{R}^3$ s.t. $x^2+y^2+z^2=1\}$.
If $(X,Y) \in \mathbb{R}^2 \cup \{\infty\}$ and $(x,y,z) \in S$, then the bijection is given by
$\begin{eqnarray*}
(x,y,z) &=& \left(\frac{2X}{X^2+Y^2+1}, \;\frac{2Y}{X^2+Y^2+1}, \;\frac{X^2+Y^2-1}{X^2+Y^2+1}\right) \\
(X,Y) &=& \left(\frac{x}{1-z}, \;\frac{y}{1-z}\right)
\end{eqnarray*}$
This allows us to 'wrap' any two-dimensional graph or image onto the unit sphere, or vice versa:
{{{id=10|
def symbolic_project_up(A):
V = 1 + A[0]^2 + A[1]^2
return vector([2*A[0]/V, 2*A[1]/V, (V-2)/V])
def symbolic_project_down(B):
return vector([B[0]/(1-B[2]), B[1]/(1-B[2])])
///
}}}
{{{id=15|
v = vector((3/2,3/2))
show(symbolic_project_up(v))
///
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{6}{11},\frac{6}{11},\frac{7}{11}\right)
}}}
{{{id=14|
w = vector((6/11,6/11,7/11))
show(symbolic_project_down(w))
///
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{3}{2},\frac{3}{2}\right)
}}}
{{{id=1|
%hide
P = sphere(opacity=0.4,color='lightblue')
Q = polygon3d([[-2,-2,0], [-2,2,0], [2,2,0], [2,-2,0]],opacity=0.2)
R = line([(3/2,3/2,0),(0,0,1)],thickness=2,color='black')
S = point3d((3/2,3/2,0),size=10,color='red') + point3d((0,0,1),size=10,color='green')
S += point3d((6/11,6/11,7/11),size=10,color='yellow')
A = P+Q+R+S
A.show(aspect_ratio=1)
///
}}}
- Points on the plane inside the unit disc get mapped to the southern hemisphere
- Points on the unit circle stay put
- Everything else goes to the northern hemisphere
- The point at infinity goes to the north pole
The stereographic projection is an example of a conformal map: angles are preserved:
{{{id=6|
%hide
import numpy as np
P = sphere(opacity=0.4,color='lightblue')
Q = polygon3d([[-2,-2,0], [-2,2,0], [2,2,0], [2,-2,0]],opacity=0.2)
D = np.arange(0,1,0.01)
O = np.zeros(len(D))
X1 = D + 1
Y1 = O - 3/2
X2 = O + 3/2
Y2 = D - 2
R = line(zip(X1,Y1),color='red',thickness=4) + line(zip(X2,Y2),color='red',thickness=4)
v1 = 1 + X1^2 + Y1^2
v2 = 1 + X2^2 + Y2^2
x1,y1,z1 = 2*X1/v1, 2*Y1/v1, (v1-2)/v1
x2,y2,z2 = 2*X2/v2, 2*Y2/v2, (v2-2)/v2
S = line(zip(x1,y1,z1),color='black',thickness=3) + line(zip(x2,y2,z2),color='black',thickness=3)
A = P+Q+R+S
T = line([(1,-1.5,0),(0,0,1)],color='yellow',thickness=1)
T += line([(1.5,-1,0),(0,0,1)],color='yellow',thickness=1)
T += line([(2,-1.5,0),(0,0,1)],color='yellow',thickness=1)
T += line([(1.5,-2,0),(0,0,1)],color='yellow',thickness=1)
T += point3d((0,0,1),size=10,color='green')
A += T
A.show(aspect_ratio=1)
///
}}}
If one identifies $\mathbb{R}^2$ with the complex plane (i.e. $Z = X + iY$), then the stereographic transformations become
$\begin{eqnarray*}
P(Z) &=& \left(\frac{Z+\bar{Z}}{Z \bar{Z}+1}, \;\frac{-i(Z-\bar{Z})}{Z \bar{Z}+1}, \;\frac{Z \bar{Z}-1}{Z \bar{Z}+1}\right) \\
P^{-1}((x,y,z)) &=& \left(\frac{x}{1-z}\right) + i \left(\frac{y}{1-z}\right),
\end{eqnarray*}$
What happens when you project up, roll the sphere around, then project back down?
You get something called a linear fractional transformation from the plane back onto itself. These are conformal maps of the form
$\begin{equation*}
f(Z) = \frac{aZ + b}{cZ + d}
\end{equation*}$
for some complex constants $a,b,c$ and $d$.
And here's what they look like:
Some directions for the following interact:
- Type in any explicit function $f(x)$ or any elliptic curve (e.g. $y^2 == x^3 - 4x + 1$ or `389a' in quote marks) to see its pre- and post-transformation plots
- Upload a .png image to the worksheet and type in its name in quote marks to see it being transformed. 'lion.png' has already been uploaded, so try that one out first.
{{{id=20|
%hide
# Rotates a given numpy array of points on the unit sphere about one of the three principal axes
def rotate_about_axis(B,theta,axis=2):
# Create the rotation matrix. Rotation about the x-axis corresponds to axis==0, y-axis to axis==1 and
# z-axis to axis==2
M = np.array([])
if axis==0:
M = np.array([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]],dtype='float128')
elif axis==1:
M = np.array([[cos(theta),0,-sin(theta)],[0,1,0],[sin(theta),0,cos(theta)]],dtype='float128')
elif axis==2:
M = np.array([[cos(theta),-sin(theta),0],[sin(theta),cos(theta),0],[0,0,1]],dtype='float128')
# Numpy makes large floating point matrix manipulations easy
return np.dot(M,B)
# Projects a numpy array of plane points up onto the sphere
def project_up(A):
V = 1 + A[0]^2 + A[1]^2
return np.array([2*A[0]/V, 2*A[1]/V, (V-2)/V])
# Projects a numpy array of sphere points down onto the plane
def project_down(B):
return np.array([B[0]/(1-B[2]), B[1]/(1-B[2])])
# Projects the numpy array of plate points up onto the sphere, rotates the sphere,
# and projects them back down to the plane
def project(A, x_axis_rot, y_axis_rot, z_axis_rot):
H = project_up(A)
H = rotate_about_axis(H,x_axis_rot,0)
H = rotate_about_axis(H,y_axis_rot,1)
H = rotate_about_axis(H,z_axis_rot,2)
return project_down(H), H
# Extracts the x,y data from the image file and puts it into one dimensional
# arrays that can be used with the projection functions
def get_image_data(name):
image = pylab.imread(DATA + name)
y_max, x_max, pixSize = image.shape
# Create the empty one dimensional arrays with length enough to hold all the pixels
x_data = np.zeros(x_max*y_max)
y_data = np.zeros(x_max*y_max)
# Fill in the indices of the pixels, row by row
for i in range(y_max):
for j in range(x_max):
x_data[i*x_max+j] = j
y_data[i*x_max+j] = i
# Scale the indices down to 0 to 1 range to make the transformations look better
x_data /= x_max
y_data /= y_max
return [x_data, y_data, image]
# Wraps the projected x,y data back up into the pixel matrix so it can be displayed
# as a graphics object
def build_image(x_data, y_data, name):
image = pylab.imread(DATA + name)
old_y_max, old_x_max, pixSize = image.shape
# Scale the indices back up to their original sizes
x_data *= old_x_max
y_data *= old_y_max
# Find the maximum and minimum transforme indices so that they can be shifted
# back towards the origin
new_x_max = int(max(x_data))+1
new_x_min = int(min(x_data))
new_y_max = int(max(y_data))+1
new_y_min = int(min(y_data))
width = new_x_max - new_x_min
height = new_y_max - new_y_min
new_image = np.zeros((height, width, pixSize))
# Iterate through all the indices, filling in the original pixels at the
# corresponding transformed indices
for i in range(0, old_y_max):
for j in range(0, old_x_max):
xInd = int(x_data[i*old_x_max+j]) - new_x_min
yInd = int(y_data[i*old_x_max+j]) - new_y_min
new_image[yInd][xInd] = image[i][j]
return new_image
///
}}}
{{{id=12|
%hide
import numpy as np
import pylab
var('x y')
@interact
def plotSterographicProjection(curve=input_box(x*cos(x), 'Function'), \
zoom=slider(1, 100, 1, 30, 'Zoom'), \
plot_resolution=slider(1, 10000, 1, 1000, 'Plot Resolution'), \
x_axis_rot=slider(-pi, pi, pi/100, 0, 'Rotate Sphere About x-axis'), \
y_axis_rot=slider(-pi, pi, pi/100, 0, 'Rotate Sphere About y-axis'), \
z_axis_rot=slider(-pi, pi, pi/100, 0, 'Rotate Sphere About z-axis'), \
show_3d_plot=checkbox(default=False, label="Show 3d Projection")):
"""
This interact plots a stereographic projection of an explicit function, elliptic curve, or image provided by
the user.
Inputs
Function
A string that specifies the explicit function, elliptic curve, or image.
Elliptic curves and images must be wrapped in quotes
Images must be uploaded to the worksheet before they can be used
Zoom
A slider that specifies the zoom level of the plot
This input has no function for images
Plot Resolution
A slider that specifies the number of points to plot
Increasing this will make the projection calculate slower, but improve the accuracy of the plot
This input has no function for images
Rotate Sphere About x-axis, y-axis, z-axis
Sliders that specify the number of radians to rotate the unit sphere about the specified axis
Show 3-d projection
A check box that specifies whether or not to show the projection of the plot onto the unit sphere
This input has no function for images
Output
Graphics showing the original function and its stereographic projection. Optional 3-d graphics showing the
original function projected onto the unit sphere
EXAMPLES
Explicit function: x*cos(x)
Elliptic curve: '389a'
Image: 'lion.png'
"""
# Figure out what object we're dealing with
try:
X, Y, image = get_image_data(curve);
P = None
except:
try:
E = EllipticCurve(curve)
P = plot(E, plot_points=plot_resolution,xmin=-zoom,xmax=zoom)
except:
try:
P = plot(curve, (x,-zoom,zoom),plot_points=plot_resolution)
except:
print "Please provide one of the following:\n\tElliptic Curve\n\tExplicit function in terms of x \
\n\tImage file (note that zoom and plot resolution do not affect images)"
return
# R, S and T are empty graphics objects that we will fill with the line data for the original,
# transformed and 3d graphs respectively
R,S,T = Graphics(), Graphics(), Graphics()
if P == None:
# Pull the X and Y data from the image, project it
A = np.array([X, Y])
B, H = project(A, x_axis_rot, y_axis_rot, z_axis_rot)
# Reconstruct the image matrix
new_image = build_image(B[0], B[1], curve)
# Plot
g1 = graphics_array([matrix_plot(image)])
g2 = graphics_array([matrix_plot(new_image)])
show(g1, axes=False, figsize=(8,3))
show(g2, axes=False, figsize=(8,3))
else:
# Iterate over the plot objects in P (elliptic curves can have two branches)
for p in P:
# Pull the X and Y data from P, project it
A = np.array([p.xdata,p.ydata])
B, H = project(A, x_axis_rot, y_axis_rot, z_axis_rot)
# Create line objects from the 3 sets of data
R += line(zip(A[0],A[1]),rgbcolor=(0,1/2,0),thickness=2)
S += line(zip(B[0],B[1]),rgbcolor=(1/2,0,0),thickness=2)
if show_3d_plot: T += line(zip(H[0],H[1],H[2]),rgbcolor=(1,0,0),thickness=3)
#Plot
G = graphics_array([R,S])
show(G, figsize=(9,3))
if show_3d_plot:
T += sphere(center=(0,0,0),size=1,opacity=0.4)
T.show(aspect_ratio=[1,1,1])
///
}}}
A majority of the time was spent iterating from extremely simple code (plotting a single explicit function onto a unit sphere) to being able to do the full 2D-3D-2D mapping on elliptic curves.
Had to find a few tricks to do that.
- An elliptic curve plot object is made up of 1 or more line objects (for each chunk).
- If you have a plot p, you can get the individual line objects with p[0], p[1] and so on.
- Once you have that, you can use .xdata and .ydata to get lists of the x,y values of the line
The original code was also extremely slow.
- We were able to optimize by only having one code path
- Using numpy for doing the transformations on the matrices of x,y points greatly sped it up
Lots of room for future work.
- Work with implicit functions
- Hoped implicit plots would help - but the data structure for x,y data is really confusing!
- We asked in the irc channel and the sage mailing list but no one replied :(
- Performing transformations by manipulating the sphere itself instead of interacts.
- More dynamic behavior
- showing the result of the transformation as you move the slider instead of destroying and rebuilding the plots.
- Make it even faster (especially for large x,y arrays)
- Compute the end LFT and apply instead of doing a series of transformations in succession.
{{{id=23|
///
}}}
{{{id=25|
///
}}}
{{{id=26|
///
}}}