如何在更高维度的超球面上均匀分布点? [英] How to distribute points evenly on the surface of hyperspheres in higher dimensions?
问题描述
我感兴趣的是将N个点均匀分布在3维及更高尺寸的球体表面上.
I am interested in evenly distributing N points on the surface of spheres in dimensions 3 and higher.
更具体地说:
- 给出点N和维数D(其中D> 1，N> 1)
- 每个点到原点的距离必须为1
- 两点之间的最小距离应尽可能大
- 每个点到它最近的邻居的距离不必对于每个点都相同(实际上，除非点的数量形成柏拉图式实体的顶点，或者如果N < = D).
我对以下内容不感兴趣:
I am not interested in:
- 在超球面上创建均匀的随机分布，因为我希望任意两点之间的最小距离尽可能大而不是随机分布.
- 粒子排斥模拟类型的方法，因为它们难以实现并且需要花费很长的时间才能运行较大的N(理想情况下，该方法应是确定性的，且应为O(n)).
一种满足这些条件的方法称为斐波那契晶格，但我只能在2d和3d中找到针对该条件的代码实现.
One method that satisfies these criteria is called the fibonacci lattice, but I have only been able to find code implementations for that in 2d and 3d.
斐波纳契晶格(也称为斐波纳契螺旋)背后的方法是生成围绕球体表面螺旋的一维线，以使该线所覆盖的表面积在每一圈都大致相同.然后，您可以下降N个点，它们均匀分布在螺旋上，并且它们将大致均匀地分布在球体的表面上.
The method behind the fibonacci lattice (also known as the fibonacci spiral) is to generate a 1d line that spirals around the surface of the sphere such that the surface area covered by the line is roughly the same at every turn. You can then drop N points equally distributed on the spiral and they will roughly be evenly distributed on the surface of the sphere.
在此答案中，有一个针对3个维度的python实现，会生成以下内容:
In this answer there is a python implementation for 3 dimensions that generates the following:
我想知道斐波那契螺旋是否可以扩展到大于3的尺寸，并在数学堆栈交换中发布了一个问题.令我惊讶的是，我收到了两个令人惊讶的答案，据我所知(因为我不完全理解显示的数学方法)表明确实有可能将此方法扩展到N个维度.
I wanted to know whether the fibonacci spiral could be extended to dimensions higher than 3 and posted a question on the maths stack exchange. To my surprise I received two amazing answers which as far as I can tell (because I don't fully understand the maths shown) shows it's indeed possible to extend this method to N dimensions.
不幸的是，我对所显示的数学知识还不够了解，无法将任何一个答案都转换为(伪)代码.我是一位经验丰富的计算机程序员，但我的数学背景仅此而已.
Unfortunately I don't understand enough of the maths shown to be able to turn either answer into (pseudo)code. I am an experienced computer programmer, but my maths background only goes so far.
我将复制我认为是以下答案之一最重要的部分(不幸的是，SO不支持mathjax，因此我必须复制为图像)
I will copy in what I believe to be the most important part of one of the answers below (unfortunately SO doesn't support mathjax so I had to copy as an image)
我遇到的上述困难:
- 如何解析用于ψn的反函数?
- 给出的示例是d = 3的.如何生成任意d的公式?
在座的任何人都可以理解所涉及的数学方法，从而能够在实现对链接的斐波那契晶格问题的任一答案的伪代码实现方面取得进展吗?我知道完整的实施可能会非常困难，所以我对部分实施感到满意，该实施使我能够足够自己完成其余的工作.
Would anyone here who understands the maths involved be able to make progress towards a pseudo code implementation of either answer to the linked fibonacci lattice question? I understand a full implementation may be quite difficult so I'd be happy with a part implementation that leads me far enough to be able to complete the rest myself.
为了简化操作，我已经编写了一个函数，该函数将N个维度的球面坐标转换为笛卡尔坐标，因此该实现可以输出任意一个，因为我可以轻松地进行转换.
To make it easier, I've already coded a function that takes spherical coordinates in N dimensions and turns them into cartesian coordinates, so the implementation can output either one as I can easily convert.
此外，我看到一个答案为每个附加维使用下一个质数.我可以轻松地编写一个输出每个连续素数的函数，因此可以假定已经实现了.
Additionally I see that one answer uses the next prime number for each additional dimension. I can easily code a function that outputs each successive prime, so you can assume that's already implemented.
未能在N维上实现斐波那契晶格，我很乐意接受满足上述约束的另一种方法.
Failing an implementation of the fibonacci lattice in N dimensions, I'd be happy to accept a different method that satisfies the above constraints.
推荐答案
非常有趣的问题.我想将其实现为我的4D渲染引擎，因为我很好奇它的外观，但我又懒又无能从数学角度处理ND超验问题.
Very interesting question. I wanted to implement this into mine 4D rendering engine as I was curious how would it look like but I was too lazy and incompetent to handle ND transcendent problems from the math side.
相反，我针对此问题提出了不同的解决方案.它不是 Fibonaci Latice ！！！ 相反，我扩展了将超球面"或"n球面" 转换为超螺旋线，然后仅拟合螺旋参数，以使这些点大致等距.
Instead I come up with different solution to this problem. Its not a Fibonaci Latice !!! Instead I expand the parametrical equation of a hypersphere or n-sphere into hyperspiral and then just fit the spiral parameters so the points are more or less equidistant.
我知道这听起来很可怕，但并不是那么难，并且在解决了一些愚蠢的错别字复制/粘贴错误之后，结果对我来说(最终:)正确)
It sounds horrible I know but its not that hard and the results look correct to me (finally :) after solving some silly typos copy/paste bugs)
主要思想是对超球面使用n维参数方程式，以从角度和半径计算其表面点.在这里实现:
The main idea is to use the n-dimensional parametrical equations for hypersphere to compute its surface points from angles and radius. Here implementation:
请参见 [edit2] .现在问题可以归结为两个主要问题:
see the [edit2]. Now the problem boils down to 2 main problems:
-
计算螺钉数量
因此，如果我们希望我们的点是等距的，则它们必须等距地躺在螺旋路径上(请参见项目符号#2 )，但螺钉本身之间的距离也应相同.为此，我们可以利用超球面的几何特性.让我们从2D开始:
so if we want that our points are equidistant so they must lye on the spiral path in equidistances (see bullet #2) but also the screws itself should have the same distance between each other. For that we can exploit geometrical properties of the hypersphere. Let start with 2D:
就是screws = r/d.点数也可以推断为points = area/d^2 = PI*r^2/d^2.
so simply screws = r/d. The number of points can be also inferred as points = area/d^2 = PI*r^2/d^2.
所以我们可以简单地将2D螺旋写为:
so we can simply write 2D spiral as:
t = <0.0,1.0>
a = 2.0*M_PI*screws*t;
x = r*t*cos(a);
y = r*t*sin(a);
为简单起见，我们可以假设r=1.0，所以假设d=d/r(稍后再缩放点).然后展开(每个维度仅添加angle参数)如下所示:
To be more simple we can assume r=1.0 so d=d/r (and just scale the points later). Then the expansions (each dimension just adds angle parameter) look like this:
2D:
screws=1.0/d; // radius/d
points=M_PI/(d*d); // surface_area/d^2
a = 2.0*M_PI*t*screws;
x = t*cos(a);
y = t*sin(a);
3D:
screws=M_PI/d; // half_circumference/d
points=4.0*M_PI/(d*d); // surface_area/d^2
a= M_PI*t;
b=2.0*M_PI*t*screws;
x=cos(a) ;
y=sin(a)*cos(b);
z=sin(a)*sin(b);
4D:
screws = M_PI/d;
points = 3.0*M_PI*M_PI*M_PI/(4.0*d*d*d);
a= M_PI*t;
b= M_PI*t*screws;
c=2.0*M_PI*t*screws*screws;
x=cos(a) ;
y=sin(a)*cos(b) ;
z=sin(a)*sin(b)*cos(c);
w=sin(a)*sin(b)*sin(c);
现在要提防4D的问题只是我的假设.我凭经验发现它们与constant/d^3有关，但不完全相关.每个角度的螺丝都不同.我的假设是，除了screws^i之外没有其他尺度，但可能需要不断调整(因为对我来说结果不错，因此不需要对所得的点云进行分析)
Now beware points for 4D are just my assumption. I empirically found out that they relate to constant/d^3 but not exactly. The screws are different for each angle. Mine assumption is that there is no other scale than screws^i but it might need some constant tweaking (did not do analysis of the resulting point-cloud as the result look ok to me)
现在我们可以从单个参数t=<0.0,1.0>生成螺旋上的任何点.
Now we can generate any point on spiral from single parameter t=<0.0,1.0>.
请注意，如果您逆转方程式，则d=f(points)您可以将点作为输入值，但要注意其近似的点数不准确！
Note if you reverse the equation so d=f(points) you can have points as input value but beware its just approximate number of points not exact !!!
在螺旋上生成台阶，使点等距
这是我跳过代数混乱而使用拟合代替的部分.我只是对二进制搜索增量t进行二进制搜索，因此所得的点与上一个点相距d.因此，只需生成点t=0，然后在估计位置附近生成二进制搜索t，直到d远离起点.然后重复此操作，直到t<=1.0 ...
This is the part I skip the algebraic mess and use fitting instead. I simply binary search delta t so the resulting point is d distant to previous point. So simply generate point t=0 and then binary search t near estimated position until is d distant to the start point. Then repeat this until t<=1.0 ...
您可以使用二进制搜索或其他方式.我知道它的速度不及O(1)代数方法，但不需要为每个维度派生东西.看起来10次迭代就可以拟合，因此也不慢.
You can use binary search or what ever. I know its not as fast as O(1) algebraic approach but no need to derive the stuff for each dimension... Looks 10 iterations are enough for fitting so its not that slow either.
这里是我的4D引擎 C ++/GL/VCL 的实现:
Here implementation from my 4D engine C++/GL/VCL:
void ND_mesh::set_HyperSpiral(int N,double r,double d)
{
int i,j;
reset(N);
d/=r; // unit hyper-sphere
double dd=d*d; // d^2
if (n==2)
{
// r=1,d=!,screws=?
// S = PI*r^2
// screws = r/d
// points = S/d^2
int i0,i;
double a,da,t,dt,dtt;
double x,y,x0,y0;
double screws=1.0/d;
double points=M_PI/(d*d);
dbg=points;
da=2.0*M_PI*screws;
x0=0.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
x=(t*cos(a))-x0; x*=x;
y=(t*sin(a))-y0; y*=y;
if ((!j)&&(x+y<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
x0=t*cos(a); pnt.add(x0);
y0=t*sin(a); pnt.add(y0);
as2(i0,i);
}
}
if (n==3)
{
// r=1,d=!,screws=?
// S = 4*PI*r^2
// screws = 2*PI*r/(2*d)
// points = S/d^2
int i0,i;
double a,b,da,db,t,dt,dtt;
double x,y,z,x0,y0,z0;
double screws=M_PI/d;
double points=4.0*M_PI/(d*d);
dbg=points;
da= M_PI;
db=2.0*M_PI*screws;
x0=1.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
z0=0.0; pnt.add(z0);
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
x=cos(a) -x0; x*=x;
y=sin(a)*cos(b)-y0; y*=y;
z=sin(a)*sin(b)-z0; z*=z;
if ((!j)&&(x+y+z<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
b=db*t;
x0=cos(a) ; pnt.add(x0);
y0=sin(a)*cos(b); pnt.add(y0);
z0=sin(a)*sin(b); pnt.add(z0);
as2(i0,i);
}
}
if (n==4)
{
// r=1,d=!,screws=?
// S = 2*PI^2*r^3
// screws = 2*PI*r/(2*d)
// points = 3*PI^3/(4*d^3);
int i0,i;
double a,b,c,da,db,dc,t,dt,dtt;
double x,y,z,w,x0,y0,z0,w0;
double screws = M_PI/d;
double points=3.0*M_PI*M_PI*M_PI/(4.0*d*d*d);
dbg=points;
da= M_PI;
db= M_PI*screws;
dc=2.0*M_PI*screws*screws;
x0=1.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
z0=0.0; pnt.add(z0);
w0=0.0; pnt.add(w0);
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
c=dc*t;
x=cos(a) -x0; x*=x;
y=sin(a)*cos(b) -y0; y*=y;
z=sin(a)*sin(b)*cos(c)-z0; z*=z;
w=sin(a)*sin(b)*sin(c)-w0; w*=w;
if ((!j)&&(x+y+z+w<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z+w>dd) t-=dtt;
} dt=dtt;
if (t>1.0) break;
a=da*t;
b=db*t;
c=dc*t;
x0=cos(a) ; pnt.add(x0);
y0=sin(a)*cos(b) ; pnt.add(y0);
z0=sin(a)*sin(b)*cos(c); pnt.add(z0);
w0=sin(a)*sin(b)*sin(c); pnt.add(w0);
as2(i0,i);
}
}
for (i=0;i<pnt.num;i++) pnt.dat[i]*=r;
for (i=0;i<s1.num;i++) s1.dat[i]*=n;
for (i=0;i<s2.num;i++) s2.dat[i]*=n;
for (i=0;i<s3.num;i++) s3.dat[i]*=n;
for (i=0;i<s4.num;i++) s4.dat[i]*=n;
}
其中n=N是设置的尺寸，r是半径，d是点之间的期望距离.我正在使用很多未在此处声明的内容，但重要的是pnt[]列出了对象的点列表，并且as2(i0,i1)将索引为i0,i1的点的行添加到网格中.
Where n=N are set dimensionality, r is radius and d is desired distance between points. I am using a lot of stuff not declared here but what isimportant is just that pnt[] list the list of points of the object and as2(i0,i1) adds line from points at indexes i0,i1 to the mesh.
这里有几个屏幕截图...
Here few screenshots...
3D透视:
4D视角:
具有超平面w=0.0的4D横截面:
4D cross-section with hyperplane w=0.0:
并且具有更多点和更大半径的情况:
and the same with more points and bigger radius:
形状随着动画的旋转而变化...
the shape changes with rotations in which its animated ...
[Edit1]更多代码/信息
这是我的引擎网格类的样子:
This is how my engine mesh class look like:
//---------------------------------------------------------------------------
//--- ND Mesh: ver 1.001 ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _ND_mesh_h
#define _ND_mesh_h
//---------------------------------------------------------------------------
#include "list.h" // my dynamic list you can use std::vector<> instead
#include "nd_reper.h" // this is just 5x5 transform matrix
//---------------------------------------------------------------------------
enum _render_enum
{
_render_Wireframe=0,
_render_Polygon,
_render_enums
};
const AnsiString _render_txt[]=
{
"Wireframe",
"Polygon"
};
enum _view_enum
{
_view_Orthographic=0,
_view_Perspective,
_view_CrossSection,
_view_enums
};
const AnsiString _view_txt[]=
{
"Orthographic",
"Perspective",
"Cross section"
};
struct dim_reduction
{
int view; // _view_enum
double coordinate; // cross section hyperplane coordinate or camera focal point looking in W+ direction
double focal_length;
dim_reduction() { view=_view_Perspective; coordinate=-3.5; focal_length=2.0; }
dim_reduction(dim_reduction& a) { *this=a; }
~dim_reduction() {}
dim_reduction* operator = (const dim_reduction *a) { *this=*a; return this; }
//dim_reduction* operator = (const dim_reduction &a) { ...copy... return this; }
};
//---------------------------------------------------------------------------
class ND_mesh
{
public:
int n; // dimensions
List<double> pnt; // ND points (x0,x1,x2,x3,...x(n-1))
List<int> s1; // ND points (i0)
List<int> s2; // ND wireframe (i0,i1)
List<int> s3; // ND triangles (i0,i1,i2,)
List<int> s4; // ND tetrahedrons (i0,i1,i2,i3)
DWORD col; // object color 0x00BBGGRR
int dbg; // debug/test variable
ND_mesh() { reset(0); }
ND_mesh(ND_mesh& a) { *this=a; }
~ND_mesh() {}
ND_mesh* operator = (const ND_mesh *a) { *this=*a; return this; }
//ND_mesh* operator = (const ND_mesh &a) { ...copy... return this; }
// add simplex
void as1(int a0) { s1.add(a0); }
void as2(int a0,int a1) { s2.add(a0); s2.add(a1); }
void as3(int a0,int a1,int a2) { s3.add(a0); s3.add(a1); s3.add(a2); }
void as4(int a0,int a1,int a2,int a3){ s4.add(a0); s4.add(a1); s4.add(a2); s4.add(a3); }
// init ND mesh
void reset(int N);
void set_HyperTetrahedron(int N,double a); // dimensions, side
void set_HyperCube (int N,double a); // dimensions, side
void set_HyperSphere (int N,double r,int points); // dimensions, radius, points per axis
void set_HyperSpiral (int N,double r,double d); // dimensions, radius, distance between points
// render
void glDraw(ND_reper &rep,dim_reduction *cfg,int render); // render mesh
};
//---------------------------------------------------------------------------
#define _cube(a0,a1,a2,a3,a4,a5,a6,a7) { as4(a1,a2,a4,a7); as4(a0,a1,a2,a4); as4(a2,a4,a6,a7); as4(a1,a2,a3,a7); as4(a1,a4,a5,a7); }
//---------------------------------------------------------------------------
void ND_mesh::reset(int N)
{
dbg=0;
if (N>=0) n=N;
pnt.num=0;
s1.num=0;
s2.num=0;
s3.num=0;
s4.num=0;
col=0x00AAAAAA;
}
//---------------------------------------------------------------------------
void ND_mesh::set_HyperSpiral(int N,double r,double d)
{
int i,j;
reset(N);
d/=r; // unit hyper-sphere
double dd=d*d; // d^2
if (n==2)
{
// r=1,d=!,screws=?
// S = PI*r^2
// screws = r/d
// points = S/d^2
int i0,i;
double a,da,t,dt,dtt;
double x,y,x0,y0;
double screws=1.0/d;
double points=M_PI/(d*d);
dbg=points;
da=2.0*M_PI*screws;
x0=0.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
x=(t*cos(a))-x0; x*=x;
y=(t*sin(a))-y0; y*=y;
if ((!j)&&(x+y<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
x0=t*cos(a); pnt.add(x0);
y0=t*sin(a); pnt.add(y0);
as2(i0,i);
}
}
if (n==3)
{
// r=1,d=!,screws=?
// S = 4*PI*r^2
// screws = 2*PI*r/(2*d)
// points = S/d^2
int i0,i;
double a,b,da,db,t,dt,dtt;
double x,y,z,x0,y0,z0;
double screws=M_PI/d;
double points=4.0*M_PI/(d*d);
dbg=points;
da= M_PI;
db=2.0*M_PI*screws;
x0=1.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
z0=0.0; pnt.add(z0);
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
x=cos(a) -x0; x*=x;
y=sin(a)*cos(b)-y0; y*=y;
z=sin(a)*sin(b)-z0; z*=z;
if ((!j)&&(x+y+z<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
b=db*t;
x0=cos(a) ; pnt.add(x0);
y0=sin(a)*cos(b); pnt.add(y0);
z0=sin(a)*sin(b); pnt.add(z0);
as2(i0,i);
}
}
if (n==4)
{
// r=1,d=!,screws=?
// S = 2*PI^2*r^3
// screws = 2*PI*r/(2*d)
// points = 3*PI^3/(4*d^3);
int i0,i;
double a,b,c,da,db,dc,t,dt,dtt;
double x,y,z,w,x0,y0,z0,w0;
double screws = M_PI/d;
double points=3.0*M_PI*M_PI*M_PI/(4.0*d*d*d);
dbg=points;
da= M_PI;
db= M_PI*screws;
dc=2.0*M_PI*screws*screws;
x0=1.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
z0=0.0; pnt.add(z0);
w0=0.0; pnt.add(w0);
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
c=dc*t;
x=cos(a) -x0; x*=x;
y=sin(a)*cos(b) -y0; y*=y;
z=sin(a)*sin(b)*cos(c)-z0; z*=z;
w=sin(a)*sin(b)*sin(c)-w0; w*=w;
if ((!j)&&(x+y+z+w<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z+w>dd) t-=dtt;
} dt=dtt;
if (t>1.0) break;
a=da*t;
b=db*t;
c=dc*t;
x0=cos(a) ; pnt.add(x0);
y0=sin(a)*cos(b) ; pnt.add(y0);
z0=sin(a)*sin(b)*cos(c); pnt.add(z0);
w0=sin(a)*sin(b)*sin(c); pnt.add(w0);
as2(i0,i);
}
}
for (i=0;i<pnt.num;i++) pnt.dat[i]*=r;
for (i=0;i<s1.num;i++) s1.dat[i]*=n;
for (i=0;i<s2.num;i++) s2.dat[i]*=n;
for (i=0;i<s3.num;i++) s3.dat[i]*=n;
for (i=0;i<s4.num;i++) s4.dat[i]*=n;
}
//---------------------------------------------------------------------------
void ND_mesh::glDraw(ND_reper &rep,dim_reduction *cfg,int render)
{
int N,i,j,i0,i1,i2,i3;
const int n0=0,n1=n,n2=n+n,n3=n2+n,n4=n3+n;
double a,b,w,F,*p0,*p1,*p2,*p3,_zero=1e-6;
vector<4> v;
List<double> tmp,t0; // temp
List<double> S1,S2,S3,S4; // reduced simplexes
#define _swap(aa,bb) { double *p=aa.dat; aa.dat=bb.dat; bb.dat=p; int q=aa.siz; aa.siz=bb.siz; bb.siz=q; q=aa.num; aa.num=bb.num; bb.num=q; }
// apply transform matrix pnt -> tmp
tmp.allocate(pnt.num); tmp.num=pnt.num;
for (i=0;i<pnt.num;i+=n)
{
v.ld(0.0,0.0,0.0,0.0);
for (j=0;j<n;j++) v.a[j]=pnt.dat[i+j];
rep.l2g(v,v);
for (j=0;j<n;j++) tmp.dat[i+j]=v.a[j];
}
// copy simplexes and convert point indexes to points (only due to cross section)
S1.allocate(s1.num*n); S1.num=0; for (i=0;i<s1.num;i++) for (j=0;j<n;j++) S1.add(tmp.dat[s1.dat[i]+j]);
S2.allocate(s2.num*n); S2.num=0; for (i=0;i<s2.num;i++) for (j=0;j<n;j++) S2.add(tmp.dat[s2.dat[i]+j]);
S3.allocate(s3.num*n); S3.num=0; for (i=0;i<s3.num;i++) for (j=0;j<n;j++) S3.add(tmp.dat[s3.dat[i]+j]);
S4.allocate(s4.num*n); S4.num=0; for (i=0;i<s4.num;i++) for (j=0;j<n;j++) S4.add(tmp.dat[s4.dat[i]+j]);
// reduce dimensions
for (N=n;N>2;)
{
N--;
if (cfg[N].view==_view_Orthographic){} // no change
if (cfg[N].view==_view_Perspective)
{
w=cfg[N].coordinate;
F=cfg[N].focal_length;
for (i=0;i<S1.num;i+=n)
{
a=S1.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S1.dat[i+j]*=a;
}
for (i=0;i<S2.num;i+=n)
{
a=S2.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S2.dat[i+j]*=a;
}
for (i=0;i<S3.num;i+=n)
{
a=S3.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S3.dat[i+j]*=a;
}
for (i=0;i<S4.num;i+=n)
{
a=S4.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S4.dat[i+j]*=a;
}
}
if (cfg[N].view==_view_CrossSection)
{
w=cfg[N].coordinate;
_swap(S1,tmp); for (S1.num=0,i=0;i<tmp.num;i+=n1) // points
{
p0=tmp.dat+i+n0;
if (fabs(p0[N]-w)<=_zero)
{
for (j=0;j<n;j++) S1.add(p0[j]);
}
}
_swap(S2,tmp); for (S2.num=0,i=0;i<tmp.num;i+=n2) // lines
{
p0=tmp.dat+i+n0; a=p0[N]; b=p0[N];// a=min,b=max
p1=tmp.dat+i+n1; if (a>p1[N]) a=p1[N]; if (b<p1[N]) b=p1[N];
if (fabs(a-w)+fabs(b-w)<=_zero) // fully inside
{
for (j=0;j<n;j++) S2.add(p0[j]);
for (j=0;j<n;j++) S2.add(p1[j]);
continue;
}
if ((a<=w)&&(b>=w)) // intersection -> points
{
a=(w-p0[N])/(p1[N]-p0[N]);
for (j=0;j<n;j++) S1.add(p0[j]+a*(p1[j]-p0[j]));
}
}
_swap(S3,tmp); for (S3.num=0,i=0;i<tmp.num;i+=n3) // triangles
{
p0=tmp.dat+i+n0; a=p0[N]; b=p0[N];// a=min,b=max
p1=tmp.dat+i+n1; if (a>p1[N]) a=p1[N]; if (b<p1[N]) b=p1[N];
p2=tmp.dat+i+n2; if (a>p2[N]) a=p2[N]; if (b<p2[N]) b=p2[N];
if (fabs(a-w)+fabs(b-w)<=_zero) // fully inside
{
for (j=0;j<n;j++) S3.add(p0[j]);
for (j=0;j<n;j++) S3.add(p1[j]);
for (j=0;j<n;j++) S3.add(p2[j]);
continue;
}
if ((a<=w)&&(b>=w)) // cross section -> t0
{
t0.num=0;
i0=0; if (p0[N]<w-_zero) i0=1; if (p0[N]>w+_zero) i0=2;
i1=0; if (p1[N]<w-_zero) i1=1; if (p1[N]>w+_zero) i1=2;
i2=0; if (p2[N]<w-_zero) i2=1; if (p2[N]>w+_zero) i2=2;
if (i0+i1==3){ a=(w-p0[N])/(p1[N]-p0[N]); for (j=0;j<n;j++) t0.add(p0[j]+a*(p1[j]-p0[j])); }
if (i1+i2==3){ a=(w-p1[N])/(p2[N]-p1[N]); for (j=0;j<n;j++) t0.add(p1[j]+a*(p2[j]-p1[j])); }
if (i2+i0==3){ a=(w-p2[N])/(p0[N]-p2[N]); for (j=0;j<n;j++) t0.add(p2[j]+a*(p0[j]-p2[j])); }
if (!i0) for (j=0;j<n;j++) t0.add(p0[j]);
if (!i1) for (j=0;j<n;j++) t0.add(p1[j]);
if (!i2) for (j=0;j<n;j++) t0.add(p2[j]);
if (t0.num==n1) for (j=0;j<t0.num;j++) S1.add(t0.dat[j]);// copy t0 to target simplex based on points count
if (t0.num==n2) for (j=0;j<t0.num;j++) S2.add(t0.dat[j]);
if (t0.num==n3) for (j=0;j<t0.num;j++) S3.add(t0.dat[j]);
}
}
_swap(S4,tmp); for (S4.num=0,i=0;i<tmp.num;i+=n4) // tetrahedrons
{
p0=tmp.dat+i+n0; a=p0[N]; b=p0[N];// a=min,b=max
p1=tmp.dat+i+n1; if (a>p1[N]) a=p1[N]; if (b<p1[N]) b=p1[N];
p2=tmp.dat+i+n2; if (a>p2[N]) a=p2[N]; if (b<p2[N]) b=p2[N];
p3=tmp.dat+i+n3; if (a>p3[N]) a=p3[N]; if (b<p3[N]) b=p3[N];
if (fabs(a-w)+fabs(b-w)<=_zero) // fully inside
{
for (j=0;j<n;j++) S4.add(p0[j]);
for (j=0;j<n;j++) S4.add(p1[j]);
for (j=0;j<n;j++) S4.add(p2[j]);
for (j=0;j<n;j++) S4.add(p3[j]);
continue;
}
if ((a<=w)&&(b>=w)) // cross section -> t0
{
t0.num=0;
i0=0; if (p0[N]<w-_zero) i0=1; if (p0[N]>w+_zero) i0=2;
i1=0; if (p1[N]<w-_zero) i1=1; if (p1[N]>w+_zero) i1=2;
i2=0; if (p2[N]<w-_zero) i2=1; if (p2[N]>w+_zero) i2=2;
i3=0; if (p3[N]<w-_zero) i3=1; if (p3[N]>w+_zero) i3=2;
if (i0+i1==3){ a=(w-p0[N])/(p1[N]-p0[N]); for (j=0;j<n;j++) t0.add(p0[j]+a*(p1[j]-p0[j])); }
if (i1+i2==3){ a=(w-p1[N])/(p2[N]-p1[N]); for (j=0;j<n;j++) t0.add(p1[j]+a*(p2[j]-p1[j])); }
if (i2+i0==3){ a=(w-p2[N])/(p0[N]-p2[N]); for (j=0;j<n;j++) t0.add(p2[j]+a*(p0[j]-p2[j])); }
if (i0+i3==3){ a=(w-p0[N])/(p3[N]-p0[N]); for (j=0;j<n;j++) t0.add(p0[j]+a*(p3[j]-p0[j])); }
if (i1+i3==3){ a=(w-p1[N])/(p3[N]-p1[N]); for (j=0;j<n;j++) t0.add(p1[j]+a*(p3[j]-p1[j])); }
if (i2+i3==3){ a=(w-p2[N])/(p3[N]-p2[N]); for (j=0;j<n;j++) t0.add(p2[j]+a*(p3[j]-p2[j])); }
if (!i0) for (j=0;j<n;j++) t0.add(p0[j]);
if (!i1) for (j=0;j<n;j++) t0.add(p1[j]);
if (!i2) for (j=0;j<n;j++) t0.add(p2[j]);
if (!i3) for (j=0;j<n;j++) t0.add(p3[j]);
if (t0.num==n1) for (j=0;j<t0.num;j++) S1.add(t0.dat[j]);// copy t0 to target simplex based on points count
if (t0.num==n2) for (j=0;j<t0.num;j++) S2.add(t0.dat[j]);
if (t0.num==n3) for (j=0;j<t0.num;j++) S3.add(t0.dat[j]);
if (t0.num==n4) for (j=0;j<t0.num;j++) S4.add(t0.dat[j]);
}
}
}
}
glColor4ubv((BYTE*)(&col));
if (render==_render_Wireframe)
{
// add points from higher primitives
for (i=0;i<S2.num;i++) S1.add(S2.dat[i]);
for (i=0;i<S3.num;i++) S1.add(S3.dat[i]);
for (i=0;i<S4.num;i++) S1.add(S4.dat[i]);
glPointSize(5.0);
glBegin(GL_POINTS);
glNormal3d(0.0,0.0,1.0);
if (n==2) for (i=0;i<S1.num;i+=n1) glVertex2dv(S1.dat+i);
if (n>=3) for (i=0;i<S1.num;i+=n1) glVertex3dv(S1.dat+i);
glEnd();
glPointSize(1.0);
glBegin(GL_LINES);
glNormal3d(0.0,0.0,1.0);
if (n==2)
{
for (i=0;i<S2.num;i+=n1) glVertex2dv(S2.dat+i);
for (i=0;i<S3.num;i+=n3)
{
glVertex2dv(S3.dat+i+n0); glVertex2dv(S3.dat+i+n1);
glVertex2dv(S3.dat+i+n1); glVertex2dv(S3.dat+i+n2);
glVertex2dv(S3.dat+i+n2); glVertex2dv(S3.dat+i+n0);
}
for (i=0;i<S4.num;i+=n4)
{
glVertex2dv(S4.dat+i+n0); glVertex2dv(S4.dat+i+n1);
glVertex2dv(S4.dat+i+n1); glVertex2dv(S4.dat+i+n2);
glVertex2dv(S4.dat+i+n2); glVertex2dv(S4.dat+i+n0);
glVertex2dv(S4.dat+i+n0); glVertex2dv(S4.dat+i+n3);
glVertex2dv(S4.dat+i+n1); glVertex2dv(S4.dat+i+n3);
glVertex2dv(S4.dat+i+n2); glVertex2dv(S4.dat+i+n3);
}
}
if (n>=3)
{
for (i=0;i<S2.num;i+=n1) glVertex3dv(S2.dat+i);
for (i=0;i<S3.num;i+=n3)
{
glVertex3dv(S3.dat+i+n0); glVertex3dv(S3.dat+i+n1);
glVertex3dv(S3.dat+i+n1); glVertex3dv(S3.dat+i+n2);
glVertex3dv(S3.dat+i+n2); glVertex3dv(S3.dat+i+n0);
}
for (i=0;i<S4.num;i+=n4)
{
glVertex3dv(S4.dat+i+n0); glVertex3dv(S4.dat+i+n1);
glVertex3dv(S4.dat+i+n1); glVertex3dv(S4.dat+i+n2);
glVertex3dv(S4.dat+i+n2); glVertex3dv(S4.dat+i+n0);
glVertex3dv(S4.dat+i+n0); glVertex3dv(S4.dat+i+n3);
glVertex3dv(S4.dat+i+n1); glVertex3dv(S4.dat+i+n3);
glVertex3dv(S4.dat+i+n2); glVertex3dv(S4.dat+i+n3);
}
}
glEnd();
}
if (render==_render_Polygon)
{
double nor[3],a[3],b[3],q;
#define _triangle2(ss,p0,p1,p2) \
{ \
glVertex2dv(ss.dat+i+p0); \
glVertex2dv(ss.dat+i+p1); \
glVertex2dv(ss.dat+i+p2); \
}
#define _triangle3(ss,p0,p1,p2) \
{ \
for(j=0;(j<3)&&(j<n);j++) \
{ \
a[j]=ss.dat[i+p1+j]-ss.dat[i+p0+j]; \
b[j]=ss.dat[i+p2+j]-ss.dat[i+p1+j]; \
} \
for(;j<3;j++){ a[j]=0.0; b[j]=0.0; } \
nor[0]=(a[1]*b[2])-(a[2]*b[1]); \
nor[1]=(a[2]*b[0])-(a[0]*b[2]); \
nor[2]=(a[0]*b[1])-(a[1]*b[0]); \
q=sqrt((nor[0]*nor[0])+(nor[1]*nor[1])+(nor[2]*nor[2])); \
if (q>1e-10) q=1.0/q; else q-0.0; \
for (j=0;j<3;j++) nor[j]*=q; \
glNormal3dv(nor); \
glVertex3dv(ss.dat+i+p0); \
glVertex3dv(ss.dat+i+p1); \
glVertex3dv(ss.dat+i+p2); \
}
#define _triangle3b(ss,p0,p1,p2) \
{ \
glNormal3dv(nor3.dat+(i/n)); \
glVertex3dv(ss.dat+i+p0); \
glVertex3dv(ss.dat+i+p1); \
glVertex3dv(ss.dat+i+p2); \
}
glBegin(GL_TRIANGLES);
if (n==2)
{
glNormal3d(0.0,0.0,1.0);
for (i=0;i<S3.num;i+=n3) _triangle2(S3,n0,n1,n2);
for (i=0;i<S4.num;i+=n4)
{
_triangle2(S4,n0,n1,n2);
_triangle2(S4,n3,n0,n1);
_triangle2(S4,n3,n1,n2);
_triangle2(S4,n3,n2,n0);
}
}
if (n>=3)
{
for (i=0;i<S3.num;i+=n3) _triangle3 (S3,n0,n1,n2);
for (i=0;i<S4.num;i+=n4)
{
_triangle3(S4,n0,n1,n2);
_triangle3(S4,n3,n0,n1);
_triangle3(S4,n3,n1,n2);
_triangle3(S4,n3,n2,n0);
}
glNormal3d(0.0,0.0,1.0);
}
glEnd();
#undef _triangle2
#undef _triangle3
}
#undef _swap
}
//---------------------------------------------------------------------------
#undef _cube
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
我使用我的动态列表模板是这样的:
I use mine dynamic list template so:
List<double> xxx;与double xxx[];相同
xxx.add(5);将5添加到列表的末尾
xxx[7]访问数组元素(安全)
xxx.dat[7]访问数组元素(不安全但快速的直接访问)
xxx.num是数组的实际使用大小
xxx.reset()清除数组并设置xxx.num=0
xxx.allocate(100)为100个项目预分配空间
List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items
，因此您需要将其移植到您可以使用的任何列表中(例如std:vector<>).我还使用5x5变换矩阵，其中
so you need to port it to any list you have at disposal (like std:vector<>). I also use 5x5 transform matrix where
void ND_reper::g2l (vector<4> &l,vector<4> &g); // global xyzw -> local xyzw
void ND_reper::l2g (vector<4> &g,vector<4> &l); // global xyzw <- local xyzw
将点转换为全局或局部坐标(通过将正矩阵或逆矩阵乘以点).您可以忽略它，因为它仅在渲染中使用过一次，您可以复制这些点(不旋转)...在同一标头中还包含一些常量:
convert point either to global or local coordinates (by multiplying direct or inverse matrix by point). You can ignore it as its used just once in the rendering and you can copy the points instead (no rotation)... In the same header are also some constants:
const double pi = M_PI;
const double pi2 =2.0*M_PI;
const double pipol=0.5*M_PI;
const double deg=M_PI/180.0;
const double rad=180.0/M_PI;
我还将矢量和矩阵数学模板集成在转换矩阵头中，因此vector<n>是n维向量，matrix<n>是n*n方阵，但是它仅用于渲染，因此您可以再次忽略它.如果您对这里感兴趣，那么所有链接都来自其中:
I got also vector and matrix math template integrated in the transform matrix header so vector<n> is n dimensional vector and matrix<n> is n*n square matrix but its used only for rendering so again you can ignore it. If youre interested here few links from whic all this was derived:
- How to use 4d rotors
- ND vector and square matrix math template
枚举和降维仅用于渲染. cfg包含如何将每个尺寸缩小到2D.
The enums and dimension reductions are used only for rendering. The cfg holds how should be each dimension reduced down to 2D.
AnsiString是来自 VCL 的自重定位字符串，因此可以使用char*或环境中的字符串类. DWORD只是无符号的32位int.希望我没有忘记什么...
AnsiString is a self relocating string from VCL so either use char* or string class you got in your environment. DWORD is just unsigned 32 bit int. Hope I did not forget something ...
这篇关于如何在更高维度的超球面上均匀分布点?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
