四维空间球极投影与霍普夫纤维丛具体算法实现

四维空间球极投影与霍普夫纤维丛具体算法实现

作者Wxy

本文是一篇根据一部数学科普CG电影《维度：数学漫步》[1]对四维空间与霍普夫纤维丛的讲解基础上，讲解如何在支持3D软件编程下绘制出霍普夫纤维丛的算法，实现绘制《维度：数学漫步》影片上所展示的图形。

球极投影(Stereographic Projection)

首先，我们来看一下三维空间中球极投影的公式：设一个球心在原点的半径=R的球。设投影极点为（0，0，1），投影平面：z=-1

我们不难得出：给一球面上的点（x0,y0,z0），它在投影平面：z=-1上投影坐标为： x = -((2*x0)/(-2+R*R+R*z0));

y = -((2*y0)/(-2+R*R+R*z0));

这是一段伪代码：

function stgpro(x0,y0,z0) {

x = -((2*x0)/(-2+R*R+R*z0));

y = -((2*y0)/(-2+R*R+R*z0));

return [x, y, -1];

}

如果要画一个正多面体的球极投影，得先把它“膨胀”到球心在原点的球面上：

function proSphere(x0,y0,z0) {

l = Math.sqrt(x0*x0+y0*y0+z0*z0);

return [R*x0/l, R*y0/l, R*z0/l];

}

而通过类比的思想，我们可以得出四维空间球极投影公式：

x = -((2*x0)/(-2+R*R+R*t0));

y = -((2*y0)/(-2+R*R+R*t0));

z = -((2*z0)/(-2+R*R+R*t0));

function stgpro4D(x0,y0,z0,t0) {

x = -((2*x0)/(-2+R*R+R*t0));

y = -((2*y0)/(-2+R*R+R*t0));

z = -((2*z0)/(-2+R*R+R*t0));

return [x, y, z, -1];

}

有了球极投影公式，再加上坐标旋转变换，

及多胞体顶点坐标数据[2]，我们就可以把

多胞体在四维空间中滚动做出来了。此图是我用3d max的Maxscript实现的正120胞体球极投影模型。

霍普夫纤维丛中的圆(Hopf Circle)

定义一个欧式复二维空间：C2，这个空间中的任意一点都能用两个有序复数对（Z1,Z2）表示，即用四个有序实数对表示的空间——实四维空间。

在《维度：数学漫步》第七八章的详细说明[3]中说，我们用Z2=k*Z1来表示一个四维空间

里的二维平面。（但注意用Z2=k*Z1并不能表示所有过原点的二维平面，比如xOz坐标面就不能写成这样的形式，能写成Z2=k*Z1的过原点平面是很特殊的，所以我们不能用此方法直接研究一般情形的平面，但我们用这些特殊的平面来截S3已足够了，因为它们所截得的圆周两两不相交）

所以只要给定一个k，就给定了一个二维平面，而这个二维平面与S3相交出来的圆周就确定了。现在我们要解决的是给一个k值，我们该如何让程序找到这个交出来的圆，并球极投影到我们的三维空间当中，把这个圆周描出来。

首先，设k=a+bi Z1=x+yi Z2=z+ti

则我们得到：

z=ax-by

t=ay+bx

将这两个式子与S3联立：x2+y2+z2+t2=1

我们即可以得出这个圆周在四维空间中的方程：

做一下整理，可得：

x2+y2=1/(1+a2+b2)

从这个式子中我们可以发现，给定一个复数k=a+bi，我们即定义了一个四维空间当中的过圆心的圆周，而且已知圆周上的一点的一个坐标分量（如x），根据方程我们就可以解出其余的坐标（y,z,t坐标值）而我们正是想让计算机来把圆周上的点描出来。怎样来选取圆周上的点呢？

可以考虑转换成参数方程:引入一个角θ（0≤θ≤2π）

令x = cosθ/sqrt(1+a2+b2)，则y = sinθ/sqrt(1+a2+b2)，又因为z=ax-by，t=ay+bx，我们就可以把这个圆在四维空间中描出来了。

function drawcircles(k) {

a = k[0];//get the value of Re k

b = k[1];//get the value of Im k

var x, y, z, t;

for(var i=0;i<=Math.PI*2;i+= Math.PI/20){

r = Math.sqrt(1+a*a+b*b);

x = Math.cos(i)/r;

y = Math.sin(i)/r;

z=a*x-b*y;

t=a*y+b*x;

P=stgpro4D(x,y,z,t); //得到此点在三维空间的球极投影

lineTo(P);//描点连线

}

} Array霍普夫圆环（Hopf T orus）

一个霍普夫圆环面是由无数个霍普夫圆两

两嵌套但互不相交组成。而我们把k的值绕

原点取一周得到的一簇圆就会组成一个霍

普夫圆环面。而这个圆的半径r，决定了这

个霍普夫圆环面的大小。一般我们令r=tan

β，因为k所在的这个复平面可以认为是一

个S2球极投影的射影平面，而k取值的这

一簇簇的同心圆则是上的一圈圈纬线。引入的角β则可理解为纬度。（但注意，这样定义的话南极为0°，赤道45°，北极90°），当圆的半径为0时（即南极点，k=0），霍普夫圆环会缩小成一个圆。当圆的半径无限大时（即北极点，k=0），霍普夫圆环会无限膨胀，只留下一根过四维空间S3球面极点的“直线”。最后，我们把整个纤维丛画出来：

function show(k) {

for(var j=0;j<=Math.PI/2;j+= Math.PI/6){ //移动S2上的纬线圈

r = Math.tan(j); //计算k所在圆的半径

for(var k=0;k<=Math.PI*2;k+= Math.PI/20){

//旋转一周，把k点在复平面上取遍描出来

drawcircles([r*Math.cos(k), r*Math.sin(k)]);

//画一个Hopf圆

}

}

}

附：3d max Max script 脚本完整代码(直接粘贴到Maxscript上全部运行)：效果：

sphereR=1

gl4_x=0.0

gl4_y=0.0

gl4_z=0.0

--定义坐标旋转角度

fn Ls p1 p2=sqrt((p2[2]-p1[2])*(p2[2]-p1[2])+(p2[1]-p1[1])*(p2[1]-p1[1]))

fn gl4_P p = (--四维坐标旋转函数

_point4d=#()

_point4d = rotates gl4_y #(p[1], p[4]) #(0,0)

p[1] = _point4d[1]

p[4] = _point4d[2]

_point4d = rotates gl4_z #(p[2], p[4]) #(0,0)

p[2] = _point4d[1]

p[4] = _point4d[2]

_point4d = rotates gl4_x #(p[3], p[4]) #(0,0)

p[3] = _point4d[1]

p[4] = _point4d[2]

return proplane p

)

fn proplane Ponit4D =(--四维球极投影函数

x0 = Ponit4D[1]

y0 = Ponit4D[2]

z0 = Ponit4D[3]

t0 = Ponit4D[4]

l = sqrt(x0*x0+y0*y0+z0*z0+t0*t0)

x0 = Ponit4D[1]*sphereR/l

y0 = Ponit4D[2]*sphereR/l

z0 = Ponit4D[3]*sphereR/l

t0 = Ponit4D[4]*sphereR/l

x1 = -((2*x0)/(-2+sphereR*sphereR+sphereR*t0))

y1 = -((2*y0)/(-2+sphereR*sphereR+sphereR*t0))

z1 = -((2*z0)/(-2+sphereR*sphereR+sphereR*t0))

return #(x1,y1,z1)

)

fn rotates sita p1 p2 =(

_an_ = (Patan2 p2 p1)+sita

_r_ = Ls p1 p2

_NX_ = cos(_an_)*_r_+p2[1]

_NY_ = sin(_an_)*_r_+p2[2]

return #(_NX_, _NY_)

)

fn Patan2 p1 p2=atan2 (p1[2]-p2[2]) (p1[1]-p2[1])

fn L3 p1 p2 = ((p2[3]-p1[3])*(p2[3]-p1[3])+(p2[2]-p1[2])*(p2[2]-p1[2])+(p2[1]-p1[1]) *(p2[1]-p1[1]))

for man=0 to 90 by 15 do(

s = splineShape() --create a spline shape

snum=1

for ian=0 to 360 by 20 do(

mm = tan man*cos ian

mn = tan man*sin ian

addNewSpline s --add a spline to it

for jan=0 to 360 by 5 do(

yuanx = cos(jan)/sqrt(1+mm*mm+mn*mn)

yuany = sin(jan)/sqrt(1+mm*mm+mn*mn)

yuanz = yuanx*mm-yuany*mn

yuant = yuany*mm+yuanx*mn

p0=gl4_P #(yuanx, yuany, yuanz, yuant)

addKnot s snum #smooth #curve [p0[1],p0[2],p0[3]]

)

snum=snum+1

)

updateShape s --update the shape

)

参阅：

1.《维度：数学漫步》（Dimensions: A Walk Through Mathematics）官方网站:

https://www.360docs.net/doc/6311579986.html,/Dim_ZH_si.htm

2.顶点数据来自：四维之美https://www.360docs.net/doc/6311579986.html,/atyuwen/archive/2009/11/12/tesseract.html

3. 《维度：数学漫步》详细说明第七、八章：纤维丛

https://www.360docs.net/doc/6311579986.html,/Dim_CH7_ZH_si.htm