Assume you have a sphere of radius R1 centered at the origin of the standard coordinate system O e1 e2 e3
. Then the sphere is given by all points x = [x[0], x[1], x[2]]
in 3D that satisfy the equation x[0]^2 + x[1]^2 + x[2]^2 = R1^2
. On this sphere you have a patch and the patch has a center c = [c[0], c[1], c[2]]
.
First, rotate the patch so that the center c
goes to the north pole, then project it onto a plane, using an area preserving map for the sphere of radius R1
, then map it back using the analogous area preserving map but for radius R2
sphere and finally rotate back the north pole to the scaled position of the center.
Functions you may need to define:
Function 1: Define spherical coordinates
x = sc(u, v, R):
return
x[0] = R*sin(u)*sin(v)
x[1] = R*sin(u)*cos(v)
x[2] = R*cos(u)
where
0 <= u <= pi and 0 <= v < 2*pi
Function 2: Define inverse spherical coordinates:
[u, v] = inv_sc(x, R):
return
u = arccos( x[2] / R )
if x[1] > 0
v = arccot(x[0] / x[1]) if x[1] > 0
else if x[1] < 0
v = 2*pi - arccot(x[1] / x[0])
else if x[1] = 0 and x[0] > 0
v = 0
else if x[1] = 0 and x[0] < 0
v = pi
where x[0]^2 + x[1]^2 + x[2]^2 = R^2
Function 3: Rotation matrix that rotates the center c
to the north pole:
Assume the center c
is given in spherical coordinates [uc, vc]
. Then apply function 1
c = [c[0], c[2], c[3]] = sc(uc, vc, R1)
Then, find for which index i
we have c[i] = min( abs(c[0]), abs(c[1]), abs(c[2]))
. Say i=2
and take the coordinate vector e2 = [0, 1, 0]
.
Calculate the cross-product vectors cross(c, e2)
and cross(cross(c, e2), c)
, think of them as row-vectors, and form the 3 by 3
rotation matrix
A3 = c / norm(c)
A2 = cross(c, e2) / norm(cross(c, e2))
A1 = cross(A2, A3)
A = [ A1,
A2,
A3 ]
Functions 4:
[w,z] = area_pres(u,v,R1,R2):
return
w = arccos( 1 - (R1/R2)^2 * (1 - cos(u)) )
z = v
Now if you re-scale the sphere from radius R1
to radius R2
then any point x
from the patch on the sphere with radius R1
gets transformed to the point y
on the sphere of radius R2
by the following chain of transformations:
If x is given in spherical coordinates `[ux, vx]`, first apply
x = [x[0], x[1], x[2]] = sc(ux, vx, R1)
Then rotate with the matrix A:
x = matrix_times_vector(A, x)
Then apply the chain of transformations:
[u,v] = inv_sc(x, R1)
[w,z] = area_pres(u,v,R1,R2)
y = sc(w,z,R2)
Now y is on the R2 sphere.
Finally,
y = matrix_times_vector(transpose(A), y)
As a result all of these points y
fill-in the corresponding transformed patch on the sphere of radius R2
and the patch-area on R2 equals the patch-area of the original patch on sphere R1
. Plus the center point c
gets just scaled up or down along a ray emanating from the center of the sphere.
The general idea behind this appriach is that, basically, the area element of the R1 sphere is R1^2*sin(u) du dv
and we can look for a transformation of the latitude-longitude coordinates [u,v]
of the R1 sphere into latitude-longitude coordinates [w,z]
of the R2 sphere where we have the functions w = w(u,v)
and z = z(u,v)
such that
R2^2*sin(w) dw dz = R1^2*sin(u) du dv
When you expand the derivatives of [w,z]
with respect to [u,v]
, you get
dw = dw/du(u,v) du + dw/dv(u,v) dv
dz = dz/du(u,v) du + dz/dv(u,v) dv
Plug them in the first formula, and you get
R2^2*sin(w) dw dz = R2^2*sin(w) * ( dw/du(u,v) du + dw/dv(u,v) dv ) wedge ( dz/du(u,v) du + dz/dv(u,v) dv )
= R1^2*sin(u) du dv
which simplifies to the equation
R2^2*sin(w) * ( dw/du(u,v) dz/dv(u,v) - dw/dv(u,v) dz/du(u,v) ) du dv = R^2*sin(u) du dv
So the general differential equation that guarantees the area preserving property of the transformation between the spherical patch on R1 and its image on R2 is
R2^2*sin(w) * ( dw/du(u,v) dz/dv(u,v) - dw/dv(u,v) dz/du(u,v) ) = R^2*sin(u)
Now, recall that the center of the patch has been rotated to the north pole of the R1 sphere, so you can think the center of the patch is the north pole. If you want a nice transformation of the patch so that it is somewhat homogeneous and isotropic from the patch's center, i.e. when standing at the center c
of the patch (c = north pole
) you see the patch deformed so that longitudes (great circles passing through c
) are preserved (i.e. all points from a longitude get mapped to points of the same longitude), you get the restriction that the longitude coordinate v
of point [u, v]
gets transformed to a new point [w, z]
which should be on the same longitude, i.e. z = v
. Therefore such longitude preserving transformation should look like this:
w = w(u,v)
z = v
Consequently, the area-preserving equation simplifies to the following partial differential equation
R2^2*sin(w) * dw/du(u,v) = R1^2*sin(u)
because dz/dv = 1
and dz/du = 0
.
To solve it, first fix the variable v
, and you get the ordinary differential equation
R2^2*sin(w) * dw = R1^2*sin(u) du
whose solution is
R2^2*(1 - cos(w)) = R1^2*(1 - cos(u)) + const
Therefore, when you let v
vary, the general solution for the partial differential equation
R2^2*sin(w) * dw/du(u,v) = R^2*sin(u)
in implicit form (equation that links the variables w, u, v
) should look like
R2^2*(1 - cos(w)) = R1^2*(1 - cos(u)) + f(v)
for any function f(v)
However, let us not forget that the north pole stays fixed during this transformation, i.e. we have the restriction that w= 0
whenever u = 0
. Plug this condition into the equation above and you get the restriction for the function f(v)
R2^2*(1 - cos(0)) = R1^2*(1 - cos(0)) + f(v)
R2^2*(1 - 1) = R1^2*(1 - 1) + f(v)
0 = f(v)
for every longitude v
Therefore, as soon as you impose longitudes to be transformed to the same longitudes and the north pole to be preserved, the only option you are left with is the equation
R2^2*(1 - cos(w)) = R1^2*(1 - cos(u))
which means that when you solve for w
you get
w = arccos( 1 - (R1/R2)^2 * (1 - cos(u)) )
and thus, the corresponding area preserving transformation between the patch on sphere R1 and the patch on sphere R2 with the same area, fixed center and a uniform deformation at the center so that longitudes are transformed to the same longitudes, is
w = arccos( 1 - (R1/R2)^2 * (1 - cos(u)) )
z = v
Here I implemented some of these functions in Python and ran a simple simulation:
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
def trig(uv):
return np.cos(uv), np.sin(uv)
def sc_trig(cos_uv, sin_uv, R):
n, dim = cos_uv.shape
x = np.empty((n,3), dtype=float)
x[:,0] = sin_uv[:,0]*cos_uv[:,1] #cos_u*sin_v
x[:,1] = sin_uv[:,0]*sin_uv[:,1] #cos_u*cos_v
x[:,2] = cos_uv[:,0] #sin_u
return R*x
def sc(uv,R):
cos_uv, sin_uv = trig(uv)
return sc_trig(cos_uv, sin_uv, R)
def inv_sc_trig(x):
n, dim = x.shape
cos_uv = np.empty((n,2), dtype=float)
sin_uv = np.empty((n,2), dtype=float)
Rad = np.sqrt(x[:,0]**2 + x[:,1]**2 + x[:,2]**2)
r_xy = np.sqrt(x[:,0]**2 + x[:,1]**2)
cos_uv[:,0] = x[:,2]/Rad #cos_u = x[:,2]/R
sin_uv[:,0] = r_xy/Rad #sin_v = x[:,1]/R
cos_uv[:,1] = x[:,0]/r_xy
sin_uv[:,1] = x[:,1]/r_xy
return cos_uv, sin_uv
def center_x(x,R):
n, dim = x.shape
c = np.sum(x, axis=0)/n
return R*c/math.sqrt(c.dot(c))
def center_uv(uv,R):
x = sc(uv,R)
return center_x(x,R)
def center_trig(cos_uv, sin_uv, R):
x = sc_trig(cos_uv, sin_uv, R)
return center_x(x,R)
def rot_mtrx(c):
i = np.where(c == min(c))[0][0]
e_i = np.zeros(3)
e_i[i] = 1
A = np.empty((3,3), dtype=float)
A[2,:] = c/math.sqrt(c.dot(c))
A[1,:] = np.cross(A[2,:], e_i)
A[1,:] = A[1,:]/math.sqrt(A[1,:].dot(A[1,:]))
A[0,:] = np.cross(A[1,:], A[2,:])
return A.T # ready to apply to a n x 2 matrix of points from the right
def area_pres(cos_uv, sin_uv, R1, R2):
cos_wz = np.empty(cos_uv.shape, dtype=float)
sin_wz = np.empty(sin_uv.shape, dtype=float)
cos_wz[:,0] = 1 - (R1/R2)**2 * (1 - cos_uv[:,0])
cos_wz[:,1] = cos_uv[:,1]
sin_wz[:,0] = np.sqrt(1 - cos_wz[:,0]**2)
sin_wz[:,1] = sin_uv[:,1]
return cos_wz, sin_wz
def sym_patch_0(n,m):
u = math.pi/2 + np.linspace(-math.pi/3, math.pi/3, num=n)
v = math.pi/2 + np.linspace(-math.pi/3, math.pi/3, num=m)
uv = np.empty((n, m, 2), dtype=float)
uv[:,:,0] = u[:, np.newaxis]
uv[:,:,1] = v[np.newaxis,:]
uv = np.reshape(uv, (n*m, 2), order='F')
return uv, u, v
uv, u, v = sym_patch_0(18,18)
r1 = 1
r2 = 2/3
r3 = 2
limits = max(r1,r2,r3)
p = math.pi
x = sc(uv,r1)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x[:,0], x[:,1], x[:,2])
ax.set_xlim(-limits, limits)
ax.set_ylim(-limits, limits)
ax.set_zlim(-limits, limits)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
B = rot_mtrx(center_x(x,r1))
x = x.dot(B)
cs, sn = inv_sc_trig(x)
cs1, sn1 = area_pres(cs, sn, r1, r2)
y = sc_trig(cs1, sn1, r2)
y = y.dot(B.T)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(y[:,0], y[:,1], y[:,2])
ax.set_xlim(-limits, limits)
ax.set_ylim(-limits, limits)
ax.set_zlim(-limits, limits)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
cs1, sn1 = area_pres(cs, sn, r1, r3)
y = sc_trig(cs1, sn1, r3)
y = y.dot(B.T)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(y[:,0], y[:,1], y[:,2])
ax.set_xlim(-limits, limits)
ax.set_ylim(-limits, limits)
ax.set_zlim(-limits, limits)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
plt.show()
One can see three figures of how a patch gets deformed when the radius of the sphere changes from radius 2/3, through radius 1 and finally to radius 2. The patch's area doesn't change and the transformation of the patch is homogeneous in all direction with no excessive deformation.