如何从 2D 面创建 3D X、Y、Z 数组，以便保留点之间的连续性
How to create a 3D X,Y,Z array from 2D faces so that contiguity between points is preserved
我正在尝试实现一种算法，该算法将 X、Y、Z 点的六个 2-D 矩阵存储在 X、Y、Z 点的 3D 矩阵中，从而保留它们的连通性。一个例子将阐明问题。
我要用三个 3-D 矩阵表示
( M3D_X(:,:,:), M3D_Y(:,:,:),M_3D_Z(:,:,:) )
space 的六面体区域。在那个六面体区域中，只有它的面是已知的，存储在三个二维矩阵中
( [face1_x(:,:,:),face1_y(:,:,:),face1_z(:,:,:)], [face2_x.....]....,],...,[...,face6_Z(:,:,:)] )
关于面孔之间连通性的唯一已知信息是 face_i 与 face_j 有共同的一面。没有给出哪一侧是共同的。
那么，如何将 2-D 矩阵存储在 3-D 矩阵中，以便保留面之间的连续性？
澄清一下，六面体区域的内部将使用 3-D 超限插值创建。
示例（只有两个面而不是六个）：
! []-->matrix ,--> column separator ;--> row separator
! face_i_x= [p11_x, p12_x,..,p1n_x;
....................
pm1_x,..........,pmn_x]
face1_x = [0.0,0.5,1.0;
0.0,0.5,1.0;
0.0,0.5,1.0]
face1_y = [0.0,0.0,0.0;
0.5,0.5,0.5;
1.0,1.0,1.0]
face1_z = [0.0,0.0,0.0;
0.0,0.0,0.0;
0.0,0.0,0.0]
face2_x = [0.0,0.5,1.0;
0.0,0.5,1.0;
0.0,0.5,1.0]
face2_y = [1.0,1.0,1.0;
1.0,1.0,1.0;
1.0,1.0,1.0]
face2_z = [0.0,0.0,0.0;
0.5,0.5,0.5;
1.0,1.0,1.0]
cube3D_x = (:,:,:)
cube3D_y = (:,:,:)
cube3D_z = (:,:,:)
face1表示单位立方体z=0处的面。 face2 表示单位立方体在 y=1.0 处的面。为了保持这个例子的简单性，我不会考虑 face3,..,face6。现在我需要创建一个 3D 矩阵来表示单位立方体。
我可以这样开始
cube3D_x(1,:,:)=face1_x
cube3D_y(1,:,:)=face1_y
cube3D_z(1,:,:)=face1_z
但是第二个（以及随后的插入必须特别小心处理，因为
face1 与 face2 有共同的边（y=1.0，z=0.0 处的边）。
如何在立方体中插入第二个面？
为了更清晰，这里有一个 python 脚本，可以更好地可视化问题：
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
face1x=np.array( [ [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] ])
face1y=np.array( [ [0.0, 0.0, 0.0, 0.0], [0.3,0.3,0.3,0.3], [0.7,0.7,0.7,0.7],[1.0,1.0,1.0,1.0]])
face1z=np.array( [ [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0]] )
face2x=np.array( [ [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] ])
face2y=np.array( [ [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0]] )
face2z=np.array( [ [0.0, 0.0, 0.0, 0.0], [0.3,0.3,0.3,0.3],[0.5,0.5,0.5,0.5] ,[0.7,0.7,0.7,0.7],[1.0,1.0,1.0,1.0]])
# because face2 is stored with arbitrary direction of indices, just simulate this randomness with some flipping and rotation
np.flip(face2x)
np.flip(face2y)
np.flip(face2z)
np.rot90(face2x)
np.rot90(face2y)
np.rot90(face2z)
# now plot the two faces
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(face1x,face1y,face1z,color='b')
ax.plot_wireframe(face2x,face2y,face2z,color='r')
plt.show()
还有：https://en.wikipedia.org/wiki/Regular_grid
我对性能不感兴趣。
非常感谢。
----编辑---
同样的问题可以用线而不是面和面而不是 3D 体积提出。
line1 = [1,2,3]
line2 = [5,4,3]
line3 = [5,6,7]
line4 = [9,8,7]
face= matrix(3,3)
可以从
开始
face(1,:)= line1
他获得了
face = [1,2,3 ;
*,*,* ;
*,*,* ]
现在如果line2是这样插入的
face(:,3)=line2
结果会是
face = [1,2,5 ;
*,*,4 ;
*,*,3 ]
当正确的解决方案是
face(:,3)= reverse(line2)
（因为第3点与line1和line2有共同点）得到
face = [1,2,3 ;
*,*,4 ;
*,*,5 ]
希望这能让问题更清楚。
所以如果我做对了，你得到了 6 个 3D 点的 2D 矩阵，代表 6 个具有某种形状（连接在一起）的表面，处于彼此之间的任何 flip/mirror/order 状态，并且想要构建单个 3D 矩阵空的内部（例如点 (0,0,0) 及其表面将包含重新排序的 6 个表面，重新定向以使它们的边缘匹配。
由于您总是有 6 个具有 2D 拓扑的 3D 表面，您的 3D 拓扑可以固定，这大大简化了事情。
定义 3D 拓扑
可以是任何我决定使用这个：
-------------
/ /
/ 5 /
y / /------
------------- |
x | |
| 3 |
/| | | /|
/ | | | / |
/ | y ------------- / |
| | x | |
| 0 | | 1 |
| / | /
| / ------------- | /
y|/x | | y|/x
| |
| 2 |
| |------------
| | /
y ------------- 4 /
x y / /
-------------
x
数字 0..5 代表您的面部 (surfaces/faces f0,f1,...f5) 索引和 x,y 每个二维矩阵中的索引原点和方向。
转换你的 6 个表面以匹配拓扑结构
您将需要 2 个包含 6 个值的数组。一个告诉输入矩阵是否已在拓扑中使用 (us[6])，另一个保持输入表面的索引映射到拓扑 f0...f5 (id[6])。在开始时全部设置为 -1。如果需要，us 也可以保持 id 的反向。
第一个面 f0 可以硬编码为第一个输入面
id[0]=0;
us[0]=0;
现在只需遍历 f0 的所有边并在所有未使用的表面中找到匹配的边。所以每个表面有 4 条边，但每条边有 2 个可能的方向。所以这是每个面 8 个边缘测试。
找到匹配边后，我们需要调整匹配面的方向，使其与我们选择的拓扑相匹配。所以我们需要镜像 x，镜像 y 和/或交换 x/y.
这将获得f2,f3,f4,f5.
现在只需选取任何新发现的面孔，然后通过匹配其共享边来找到 f0。再次定位 f0 以匹配拓扑。
将曲面复制到 3D 矩阵中
因此分配 3D 矩阵大小（从 f0(x,y),f2(x) 的分辨率获得，现在只需将所有 f0..f5 复制到它们在矩阵中的位置...您可以用 [=15 填充内部点=] 或在曲面之间进行插值...
这里小C++/VCL/OpenGL例子：
arrays.h 1D/2D/3D 容器
//---------------------------------------------------------------------------
#ifndef _arrays_h
#define _arrays_h
/*---------------------------------------------------------------------------
// inline safety constructors/destructors add this to any class/struct used as T
T() {}
T(T& a) { *this=a; }
~T() {}
T* operator = (const T *a) { *this=*a; return this; }
//T* operator = (const T &a) { ...copy... return this; }
//-------------------------------------------------------------------------*/
template <class T> class array1D
{
public:
T *dat; // items array
int nx; // allocated size
array1D() { dat=NULL; free(); }
array1D(array1D& a) { *this=a; }
~array1D() { free(); }
array1D* operator = (const array1D *a) { *this=*a; return this; }
array1D* operator = (const array1D &a)
{
int x;
allocate(a.nx);
for (x=0;x<nx;x++) dat[x]=a.dat[x];
return this;
}
T& operator[] (int x){ return dat[x]; }
void free()
{
if (dat!=NULL) delete[] dat;
nx=0; dat=NULL;
}
int allocate(int _nx)
{
free();
dat=new T[_nx];
nx=_nx;
}
};
//---------------------------------------------------------------------------
template <class T> class array2D
{
public:
T **dat; // items array
int nx,ny; // allocated size
array2D() { dat=NULL; free(); }
array2D(array2D& a) { *this=a; }
~array2D() { free(); }
array2D* operator = (const array2D *a) { *this=*a; return this; }
array2D* operator = (const array2D &a)
{
int x,y;
allocate(a.nx,a.ny);
for (x=0;x<nx;x++)
for (y=0;y<ny;y++)
dat[x][y]=a.dat[x][y];
return this;
}
T* operator[] (int x){ return dat[x]; }
void free()
{
if (dat!=NULL)
{
if (dat[0]!=NULL) delete[] dat[0];
delete[] dat;
}
nx=0; ny=0; dat=NULL;
}
int allocate(int _nx,int _ny)
{
free();
dat=new T*[_nx];
dat[0]=new T[_nx*_ny];
nx=_nx; ny=_ny;
for (int x=1;x<nx;x++) dat[x]=dat[x-1]+ny;
}
};
//---------------------------------------------------------------------------
template <class T> class array3D
{
public:
T ***dat; // items array
int nx,ny,nz; // allocated size
array3D() { dat=NULL; free(); }
array3D(array3D& a) { *this=a; }
~array3D() { free(); }
array3D* operator = (const array3D *a) { *this=*a; return this; }
array3D* operator = (const array3D &a)
{
int x,y,z;
allocate(a.nx,a.ny,a.nz);
for (x=0;x<nx;x++)
for (y=0;y<ny;y++)
for (z=0;z<nz;z++)
dat[x][y][z]=a.dat[x][y][z];
return this;
}
T* operator[] (int x){ return dat[x]; }
void free()
{
if (dat!=NULL)
{
if (dat[0]!=NULL)
{
if (dat[0][0]!=NULL) delete[] dat[0][0];
delete[] dat[0];
}
delete[] dat;
}
nx=0; ny=0; nz=0; dat=NULL;
}
int allocate(int _nx,int _ny,int _nz)
{
int x,y;
free();
dat=new T**[_nx];
dat[0]=new T*[_nx*_ny];
dat[0][0]=new T[_nx*_ny*_nz];
nx=_nx; ny=_ny; nz=_nz;
for (x=1;x<nx;x++) dat[x]=dat[x-1]+ny;
for (x=0;x<nx;x++)
for (y=0;y<ny;y++)
dat[x][y]=dat[0][0]+(x*ny*nz)+(y*nz);
}
};
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
只是简单的 OpenGL
parametric_grid.h这是你需要消化的主要来源
//---------------------------------------------------------------------------
#ifndef _parametric_grid_h
#define _parametric_grid_h
//---------------------------------------------------------------------------
#include "arrays.h"
//---------------------------------------------------------------------------
struct point
{
float x,y,z;
point() {}
point(point& a) { *this=a; }
~point() {}
point* operator = (const point *a) { *this=*a; return this; }
//point* operator = (const point &a) { ...copy... return this; }
bool point::operator==(const point &a)
{
float d;
d =(x-a.x)*(x-a.x);
d+=(y-a.y)*(y-a.y);
d+=(z-a.z)*(z-a.z);
return (d<=0.001);
}
bool point::operator!=(const point &a)
{
float d;
d =(x-a.x)*(x-a.x);
d+=(y-a.y)*(y-a.y);
d+=(z-a.z)*(z-a.z);
return (d>0.001);
}
};
typedef array1D<float> param[3]; // parametrers <0,1>
typedef array2D<point> surface[6]; // surface
typedef array3D<point> volume; // volume
//---------------------------------------------------------------------------
void swap_xy(array2D<point> &f)
{
int ix,iy;
array2D<point> f0;
f0=f;
f.allocate(f0.ny,f0.nx);
for (ix=0;ix<f.nx;ix++)
for (iy=0;iy<f.ny;iy++)
f.dat[ix][iy]=f0.dat[iy][ix];
}
//---------------------------------------------------------------------------
void mirror_x(array2D<point> &f)
{
int ix0,ix1,iy;
point p;
for (ix0=0,ix1=f.nx-1;ix0<ix1;ix0++,ix1--)
for (iy=0;iy<f.ny;iy++)
{ p=f.dat[ix0][iy]; f.dat[ix0][iy]=f.dat[ix1][iy]; f.dat[ix1][iy]=p; }
}
//---------------------------------------------------------------------------
void mirror_y(array2D<point> &f)
{
int ix,iy0,iy1;
point p;
for (ix=0;ix<f.nx;ix++)
for (iy0=0,iy1=f.ny-1;iy0<iy1;iy0++,iy1--)
{ p=f.dat[ix][iy0]; f.dat[ix][iy0]=f.dat[ix][iy1]; f.dat[ix][iy1]=p; }
}
//---------------------------------------------------------------------------
void surface2volume(volume &v,surface &s) // normalize surface and convert it to volume
{
point p;
float *pp=(float*)&p;
int ix,iy,iz,nx,ny,nz,mx,my,i,j,k,id[6],us[6];
// find which surface belongs to which face f0..f5 and store it to id[]
// also mirror/rotate to match local x,y orientations
id[0]=0; us[0]=0; // first face hard coded as 0
for (i=1;i<6;i++) id[i]=-1; // all the rest not known yet
for (i=1;i<6;i++) us[i]=-1; // unused
nx=s[0].nx;
ny=s[0].ny;
for (iy=0,j=id[0],k=4,i=1;i<6;i++) // find f4 so it share f0(iy=0) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
}
} if (id[k]<0) return; // fail
nz=s[id[k]].nx; // 3th resolution
v.allocate(nx,ny,nz); // allocate volume container
for (iy=ny-1,j=id[0],k=5,i=1;i<6;i++)// find f5 so it share f0(iy=ny-1) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
}
} if (id[k]<0) return; // fail
for (ix=0,j=id[0],k=2,i=1;i<6;i++) // find f2 so it share f0(ix=0) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
}
} if (id[k]<0) return; // fail
for (ix=nx-1,j=id[0],k=3,i=1;i<6;i++)// find f3 so it share f0(ix=nx-1) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
}
} if (id[k]<0) return; // fail
for (ix=nz-1,j=id[2],k=1,i=1;i<6;i++)// find f1 so it share f2(ix=nz-1) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
}
} if (id[k]<0) return; // fail
i=0;
// copy surfaces into volume matrix
for (ix=0;ix<nx;ix++)
for (iy=0;iy<ny;iy++)
for (iz=0;iz<nz;iz++)
{
// interior points
p.x=0.0;
p.y=0.0;
p.z=0.0;
// surface points
if (iz== 0) p=s[id[0]].dat[ix][iy];
if (iz==nz-1) p=s[id[1]].dat[ix][iy];
if (ix== 0) p=s[id[2]].dat[iz][iy];
if (ix==nx-1) p=s[id[3]].dat[iz][iy];
if (iy== 0) p=s[id[4]].dat[iz][ix];
if (iy==ny-1) p=s[id[5]].dat[iz][ix];
v.dat[ix][iy][iz]=p;
}
}
//---------------------------------------------------------------------------
void draw(const surface &s) // render surface
{
int ix,iy,i;
DWORD col[6]= // colors per face
{
0x00202060,
0x00202030,
0x00206020,
0x00203020,
0x00602020,
0x00302020
};
glBegin(GL_LINES);
for (i=0;i<6;i++)
{
glColor4ubv((BYTE*)(&col[i]));
for (ix=0;ix<s[i].nx;ix++)
for (iy=1;iy<s[i].ny;iy++)
{
glVertex3fv((float*)&(s[i].dat[ix][iy ]));
glVertex3fv((float*)&(s[i].dat[ix][iy-1]));
}
for (ix=1;ix<s[i].nx;ix++)
for (iy=0;iy<s[i].ny;iy++)
{
glVertex3fv((float*)&(s[i].dat[ix ][iy]));
glVertex3fv((float*)&(s[i].dat[ix-1][iy]));
}
}
glEnd();
}
//---------------------------------------------------------------------------
void draw(const volume &v)
{
int ix,iy,iz,i;
// all points
glBegin(GL_POINTS);
for (ix=0;ix<v.nx;ix++)
for (iy=0;iy<v.ny;iy++)
for (iz=0;iz<v.nz;iz++)
glVertex3fv((float*)&(v.dat[ix][iy][iz]));
glEnd();
// surface edges for visual check if points are ordered as should
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz= 0;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy=v.ny-1,iz= 0;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy=v.ny-1,iz=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz= 0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy= 0,iz= 0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy= 0,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy=v.nz-1,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy= 0,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy=v.nz-1,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
}
//---------------------------------------------------------------------------
void cube(surface &s,const param &t) // cube -> surface
{
point p;
float *pp=(float*)&p;
int iu,iv,nx,ny,nz,i,u,v;
// count of points per axis
nx=t[0].nx;
ny=t[1].nx;
nz=t[2].nx;
if ((!nx)||(!ny)||(!nz)) return;
// create 6 faces grid
for (i=0;i<6;i++)
{
if (i==0){ p.z=t[2].dat[ 0]; u=0; v=1; } // z = first coordinate
if (i==1){ p.z=t[2].dat[nz-1]; u=0; v=1; } // z = last coordinate
if (i==2){ p.y=t[0].dat[ 0]; u=0; v=2; } // y = first coordinate
if (i==3){ p.y=t[0].dat[ny-1]; u=0; v=2; } // y = last coordinate
if (i==4){ p.x=t[1].dat[ 0]; u=2; v=1; } // x = first coordinate
if (i==5){ p.x=t[2].dat[nx-1]; u=2; v=1; } // x = last coordinate
s[i].allocate(t[u].nx,t[v].nx); // allocate face resolution
for (iu=0;iu<s[i].nx;iu++) // process all cells
for (iv=0;iv<s[i].ny;iv++)
{
pp[u]=t[u].dat[iu];
pp[v]=t[v].dat[iv];
s[i].dat[iu][iv]=p;
}
}
}
//---------------------------------------------------------------------------
void cubic(surface &s,const param &t) // hardcoded cubic curved object -> surface
{
point p;
float *pp=(float*)&p;
int ix,iy,iz,nx,ny,nz,i,j;
// 3 hard coded cubic control points
float cp[3][4][3]= // [parameter][point][coordinate]
{
{{ 0.00,0.00,0.00 },
{ 0.50,0.50,0.00 },
{ 1.00,0.50,0.00 },
{ 1.50,0.00,0.00 }},
{{ 0.00,0.00,0.00 },
{ 0.10,0.15,0.00 },
{ 0.15,0.30,0.00 },
{ 0.16,0.45,0.00 }},
{{ 0.00,0.00,0.00 },
{ 0.00,0.10,0.15 },
{ 0.00,0.10,0.30 },
{ 0.00,0.00,0.45 }}
};
float ce[3][4][3]; // cubic coeficients
float tx[4],ty[4],tz[4],d1,d2;
for (ix=0;ix<3;ix++)
for (iy=0;iy<3;iy++)
{
d1=0.5*(cp[ix][2][iy]-cp[ix][0][iy]);
d2=0.5*(cp[ix][3][iy]-cp[ix][1][iy]);
ce[ix][0][iy]=cp[ix][1][iy];
ce[ix][1][iy]=d1;
ce[ix][2][iy]=(3.0*(cp[ix][2][iy]-cp[ix][1][iy]))-(2.0*d1)-d2;
ce[ix][3][iy]=d1+d2+(2.0*(-cp[ix][2][iy]+cp[ix][1][iy]));
}
// count of points per axis
nx=t[0].nx;
ny=t[1].nx;
nz=t[2].nx;
if ((!nx)||(!ny)||(!nz)) return;
// allocate faces
s[0].allocate(nx,ny);
s[1].allocate(nx,ny);
s[2].allocate(ny,nz);
s[3].allocate(ny,nz);
s[4].allocate(nz,nx);
s[5].allocate(nz,nx);
// process all points
for (ix=0;ix<nx;ix++)
for (iy=0;iy<ny;iy++)
for (iz=0;iz<nz;iz++)
{
// compute t,t^2,t^3 per each parameter
for (tx[0]=1.0,i=1;i<4;i++) tx[i]=tx[i-1]*t[0].dat[ix];
for (ty[0]=1.0,i=1;i<4;i++) ty[i]=ty[i-1]*t[1].dat[iy];
for (tz[0]=1.0,i=1;i<4;i++) tz[i]=tz[i-1]*t[2].dat[iz];
// compute position as superposition of 3 cubics
for (i=0;i<3;i++)
{
pp[i]=0.0;
for (j=0;j<4;j++) pp[i]+=ce[0][j][i]*tx[j];
for (j=0;j<4;j++) pp[i]+=ce[1][j][i]*ty[j];
for (j=0;j<4;j++) pp[i]+=ce[2][j][i]*tz[j];
}
// copy posiiton to corresponding surfaces
if (iz== 0) s[0].dat[ix][iy]=p;
if (iz==nz-1) s[1].dat[ix][iy]=p;
if (ix== 0) s[2].dat[iy][iz]=p;
if (ix==nx-1) s[3].dat[iy][iz]=p;
if (iy== 0) s[4].dat[iz][ix]=p;
if (iy==ny-1) s[5].dat[iz][ix]=p;
}
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
用法：没有任何组件的简单单窗体 VCL 应用程序
//---------------------------------------------------------------------------
#include <vcl.h> // VCL stuff (ignore)
#pragma hdrstop // VCL stuff (ignore)
#include "Unit1.h" // VCL stuff (header of this window)
#include "gl_simple.h" // my GL init (source included)
#include "parametric_grid.h"
//---------------------------------------------------------------------------
#pragma package(smart_init) // VCL stuff (ignore)
#pragma resource "*.dfm" // VCL stuff (ignore)
TForm1 *Form1; // VCL stuff (this window)
surface sur; // 6 2D matrices (surfaces in any order/orientation before normalization)
volume vol; // 3D matrix (volume)
//---------------------------------------------------------------------------
void gl_draw()
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDisable(GL_CULL_FACE);
glEnable(GL_DEPTH_TEST);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(-0.9,-0.75,-2.0);
glRotatef(45.0,0.5,0.5,0.0);
glPointSize(4.0);
glColor3f(0.1,0.5,0.7); draw(vol);
glPointSize(1.0);
// draw(sur);
glFlush();
SwapBuffers(hdc);
}
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
// this is called on window startup
gl_init(Handle); // init OpenGL 1.0
int i;
param t; // source coordinates
/*
t[0].allocate(5); i=0;
t[0][i]=0.000; i++;
t[0][i]=0.125; i++;
t[0][i]=0.250; i++;
t[0][i]=0.500; i++;
t[0][i]=1.000; i++;
t[1].allocate(5); i=0;
t[1][i]=0.000; i++;
t[1][i]=0.250; i++;
t[1][i]=0.500; i++;
t[1][i]=0.750; i++;
t[1][i]=1.000; i++;
t[2].allocate(6); i=0;
t[2][i]=0.0; i++;
t[2][i]=0.2; i++;
t[2][i]=0.4; i++;
t[2][i]=0.6; i++;
t[2][i]=0.8; i++;
t[2][i]=1.0; i++;
*/
t[0].allocate(6); i=0;
t[0][i]=0.0; i++;
t[0][i]=0.2; i++;
t[0][i]=0.4; i++;
t[0][i]=0.6; i++;
t[0][i]=0.8; i++;
t[0][i]=1.0; i++;
t[1].allocate(6); i=0;
t[1][i]=0.0; i++;
t[1][i]=0.2; i++;
t[1][i]=0.4; i++;
t[1][i]=0.6; i++;
t[1][i]=0.8; i++;
t[1][i]=1.0; i++;
t[2].allocate(6); i=0;
t[2][i]=0.0; i++;
t[2][i]=0.2; i++;
t[2][i]=0.4; i++;
t[2][i]=0.6; i++;
t[2][i]=0.8; i++;
t[2][i]=1.0; i++;
// cube(sur,t);
cubic(sur,t);
surface2volume(vol,sur);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
{
// this is called before window exits
gl_exit(); // exit OpenGL
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
{
// this is called on each window resize (and also after startup)
gl_resize(ClientWidth,ClientHeight);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
{
// this is called whnewer app needs repaint
gl_draw();
}
//---------------------------------------------------------------------------
这是上面代码的预览：
我通过为每个表面尝试 mirror/swap 的所有组合进行调试（首先不按原样使用）并渲染共享边缘（如果它们完全重叠）并调整找到的面的重新排序直到它们完成所以代码应该适用于任何输入（连接）
我正在尝试实现一种算法，该算法将 X、Y、Z 点的六个 2-D 矩阵存储在 X、Y、Z 点的 3D 矩阵中，从而保留它们的连通性。一个例子将阐明问题。
我要用三个 3-D 矩阵表示
( M3D_X(:,:,:), M3D_Y(:,:,:),M_3D_Z(:,:,:) )
space 的六面体区域。在那个六面体区域中，只有它的面是已知的，存储在三个二维矩阵中
( [face1_x(:,:,:),face1_y(:,:,:),face1_z(:,:,:)], [face2_x.....]....,],...,[...,face6_Z(:,:,:)] )
关于面孔之间连通性的唯一已知信息是 face_i 与 face_j 有共同的一面。没有给出哪一侧是共同的。 那么，如何将 2-D 矩阵存储在 3-D 矩阵中，以便保留面之间的连续性？ 澄清一下，六面体区域的内部将使用 3-D 超限插值创建。
示例（只有两个面而不是六个）：
! []-->matrix ,--> column separator ;--> row separator
! face_i_x= [p11_x, p12_x,..,p1n_x;
....................
pm1_x,..........,pmn_x]
face1_x = [0.0,0.5,1.0;
0.0,0.5,1.0;
0.0,0.5,1.0]
face1_y = [0.0,0.0,0.0;
0.5,0.5,0.5;
1.0,1.0,1.0]
face1_z = [0.0,0.0,0.0;
0.0,0.0,0.0;
0.0,0.0,0.0]
face2_x = [0.0,0.5,1.0;
0.0,0.5,1.0;
0.0,0.5,1.0]
face2_y = [1.0,1.0,1.0;
1.0,1.0,1.0;
1.0,1.0,1.0]
face2_z = [0.0,0.0,0.0;
0.5,0.5,0.5;
1.0,1.0,1.0]
cube3D_x = (:,:,:)
cube3D_y = (:,:,:)
cube3D_z = (:,:,:)
face1表示单位立方体z=0处的面。 face2 表示单位立方体在 y=1.0 处的面。为了保持这个例子的简单性，我不会考虑 face3,..,face6。现在我需要创建一个 3D 矩阵来表示单位立方体。 我可以这样开始
cube3D_x(1,:,:)=face1_x
cube3D_y(1,:,:)=face1_y
cube3D_z(1,:,:)=face1_z
但是第二个（以及随后的插入必须特别小心处理，因为
face1 与 face2 有共同的边（y=1.0，z=0.0 处的边）。
如何在立方体中插入第二个面？
为了更清晰，这里有一个 python 脚本，可以更好地可视化问题：
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
face1x=np.array( [ [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] ])
face1y=np.array( [ [0.0, 0.0, 0.0, 0.0], [0.3,0.3,0.3,0.3], [0.7,0.7,0.7,0.7],[1.0,1.0,1.0,1.0]])
face1z=np.array( [ [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0]] )
face2x=np.array( [ [0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] , [0.0,0.4,0.6,1.0],[0.0,0.4,0.6,1.0] ])
face2y=np.array( [ [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0]] )
face2z=np.array( [ [0.0, 0.0, 0.0, 0.0], [0.3,0.3,0.3,0.3],[0.5,0.5,0.5,0.5] ,[0.7,0.7,0.7,0.7],[1.0,1.0,1.0,1.0]])
# because face2 is stored with arbitrary direction of indices, just simulate this randomness with some flipping and rotation
np.flip(face2x)
np.flip(face2y)
np.flip(face2z)
np.rot90(face2x)
np.rot90(face2y)
np.rot90(face2z)
# now plot the two faces
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(face1x,face1y,face1z,color='b')
ax.plot_wireframe(face2x,face2y,face2z,color='r')
plt.show()
还有：https://en.wikipedia.org/wiki/Regular_grid 我对性能不感兴趣。 非常感谢。
----编辑--- 同样的问题可以用线而不是面和面而不是 3D 体积提出。
line1 = [1,2,3]
line2 = [5,4,3]
line3 = [5,6,7]
line4 = [9,8,7]
face= matrix(3,3)
可以从
开始face(1,:)= line1
他获得了
face = [1,2,3 ;
*,*,* ;
*,*,* ]
现在如果line2是这样插入的
face(:,3)=line2
结果会是
face = [1,2,5 ;
*,*,4 ;
*,*,3 ]
当正确的解决方案是
face(:,3)= reverse(line2)
（因为第3点与line1和line2有共同点）得到
face = [1,2,3 ;
*,*,4 ;
*,*,5 ]
希望这能让问题更清楚。
所以如果我做对了，你得到了 6 个 3D 点的 2D 矩阵，代表 6 个具有某种形状（连接在一起）的表面，处于彼此之间的任何 flip/mirror/order 状态，并且想要构建单个 3D 矩阵空的内部（例如点 (0,0,0) 及其表面将包含重新排序的 6 个表面，重新定向以使它们的边缘匹配。
由于您总是有 6 个具有 2D 拓扑的 3D 表面，您的 3D 拓扑可以固定，这大大简化了事情。
定义 3D 拓扑
可以是任何我决定使用这个：
------------- / / / 5 / y / /------ ------------- | x | | | 3 | /| | | /| / | | | / | / | y ------------- / | | | x | | | 0 | | 1 | | / | / | / ------------- | / y|/x | | y|/x | | | 2 | | |------------ | | / y ------------- 4 / x y / / ------------- x数字
0..5代表您的面部 (surfaces/faces f0,f1,...f5) 索引和x,y每个二维矩阵中的索引原点和方向。转换你的 6 个表面以匹配拓扑结构
您将需要 2 个包含 6 个值的数组。一个告诉输入矩阵是否已在拓扑中使用 (
us[6])，另一个保持输入表面的索引映射到拓扑f0...f5(id[6])。在开始时全部设置为-1。如果需要，us也可以保持id的反向。第一个面
f0可以硬编码为第一个输入面id[0]=0; us[0]=0;现在只需遍历
f0的所有边并在所有未使用的表面中找到匹配的边。所以每个表面有 4 条边，但每条边有 2 个可能的方向。所以这是每个面 8 个边缘测试。找到匹配边后，我们需要调整匹配面的方向，使其与我们选择的拓扑相匹配。所以我们需要镜像 x，镜像 y 和/或交换 x/y.
这将获得
f2,f3,f4,f5.现在只需选取任何新发现的面孔，然后通过匹配其共享边来找到
f0。再次定位f0以匹配拓扑。将曲面复制到 3D 矩阵中
因此分配 3D 矩阵大小（从
f0(x,y),f2(x)的分辨率获得，现在只需将所有f0..f5复制到它们在矩阵中的位置...您可以用 [=15 填充内部点=] 或在曲面之间进行插值...
这里小C++/VCL/OpenGL例子：
arrays.h 1D/2D/3D 容器
//---------------------------------------------------------------------------
#ifndef _arrays_h
#define _arrays_h
/*---------------------------------------------------------------------------
// inline safety constructors/destructors add this to any class/struct used as T
T() {}
T(T& a) { *this=a; }
~T() {}
T* operator = (const T *a) { *this=*a; return this; }
//T* operator = (const T &a) { ...copy... return this; }
//-------------------------------------------------------------------------*/
template <class T> class array1D
{
public:
T *dat; // items array
int nx; // allocated size
array1D() { dat=NULL; free(); }
array1D(array1D& a) { *this=a; }
~array1D() { free(); }
array1D* operator = (const array1D *a) { *this=*a; return this; }
array1D* operator = (const array1D &a)
{
int x;
allocate(a.nx);
for (x=0;x<nx;x++) dat[x]=a.dat[x];
return this;
}
T& operator[] (int x){ return dat[x]; }
void free()
{
if (dat!=NULL) delete[] dat;
nx=0; dat=NULL;
}
int allocate(int _nx)
{
free();
dat=new T[_nx];
nx=_nx;
}
};
//---------------------------------------------------------------------------
template <class T> class array2D
{
public:
T **dat; // items array
int nx,ny; // allocated size
array2D() { dat=NULL; free(); }
array2D(array2D& a) { *this=a; }
~array2D() { free(); }
array2D* operator = (const array2D *a) { *this=*a; return this; }
array2D* operator = (const array2D &a)
{
int x,y;
allocate(a.nx,a.ny);
for (x=0;x<nx;x++)
for (y=0;y<ny;y++)
dat[x][y]=a.dat[x][y];
return this;
}
T* operator[] (int x){ return dat[x]; }
void free()
{
if (dat!=NULL)
{
if (dat[0]!=NULL) delete[] dat[0];
delete[] dat;
}
nx=0; ny=0; dat=NULL;
}
int allocate(int _nx,int _ny)
{
free();
dat=new T*[_nx];
dat[0]=new T[_nx*_ny];
nx=_nx; ny=_ny;
for (int x=1;x<nx;x++) dat[x]=dat[x-1]+ny;
}
};
//---------------------------------------------------------------------------
template <class T> class array3D
{
public:
T ***dat; // items array
int nx,ny,nz; // allocated size
array3D() { dat=NULL; free(); }
array3D(array3D& a) { *this=a; }
~array3D() { free(); }
array3D* operator = (const array3D *a) { *this=*a; return this; }
array3D* operator = (const array3D &a)
{
int x,y,z;
allocate(a.nx,a.ny,a.nz);
for (x=0;x<nx;x++)
for (y=0;y<ny;y++)
for (z=0;z<nz;z++)
dat[x][y][z]=a.dat[x][y][z];
return this;
}
T* operator[] (int x){ return dat[x]; }
void free()
{
if (dat!=NULL)
{
if (dat[0]!=NULL)
{
if (dat[0][0]!=NULL) delete[] dat[0][0];
delete[] dat[0];
}
delete[] dat;
}
nx=0; ny=0; nz=0; dat=NULL;
}
int allocate(int _nx,int _ny,int _nz)
{
int x,y;
free();
dat=new T**[_nx];
dat[0]=new T*[_nx*_ny];
dat[0][0]=new T[_nx*_ny*_nz];
nx=_nx; ny=_ny; nz=_nz;
for (x=1;x<nx;x++) dat[x]=dat[x-1]+ny;
for (x=0;x<nx;x++)
for (y=0;y<ny;y++)
dat[x][y]=dat[0][0]+(x*ny*nz)+(y*nz);
}
};
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
parametric_grid.h这是你需要消化的主要来源
//---------------------------------------------------------------------------
#ifndef _parametric_grid_h
#define _parametric_grid_h
//---------------------------------------------------------------------------
#include "arrays.h"
//---------------------------------------------------------------------------
struct point
{
float x,y,z;
point() {}
point(point& a) { *this=a; }
~point() {}
point* operator = (const point *a) { *this=*a; return this; }
//point* operator = (const point &a) { ...copy... return this; }
bool point::operator==(const point &a)
{
float d;
d =(x-a.x)*(x-a.x);
d+=(y-a.y)*(y-a.y);
d+=(z-a.z)*(z-a.z);
return (d<=0.001);
}
bool point::operator!=(const point &a)
{
float d;
d =(x-a.x)*(x-a.x);
d+=(y-a.y)*(y-a.y);
d+=(z-a.z)*(z-a.z);
return (d>0.001);
}
};
typedef array1D<float> param[3]; // parametrers <0,1>
typedef array2D<point> surface[6]; // surface
typedef array3D<point> volume; // volume
//---------------------------------------------------------------------------
void swap_xy(array2D<point> &f)
{
int ix,iy;
array2D<point> f0;
f0=f;
f.allocate(f0.ny,f0.nx);
for (ix=0;ix<f.nx;ix++)
for (iy=0;iy<f.ny;iy++)
f.dat[ix][iy]=f0.dat[iy][ix];
}
//---------------------------------------------------------------------------
void mirror_x(array2D<point> &f)
{
int ix0,ix1,iy;
point p;
for (ix0=0,ix1=f.nx-1;ix0<ix1;ix0++,ix1--)
for (iy=0;iy<f.ny;iy++)
{ p=f.dat[ix0][iy]; f.dat[ix0][iy]=f.dat[ix1][iy]; f.dat[ix1][iy]=p; }
}
//---------------------------------------------------------------------------
void mirror_y(array2D<point> &f)
{
int ix,iy0,iy1;
point p;
for (ix=0;ix<f.nx;ix++)
for (iy0=0,iy1=f.ny-1;iy0<iy1;iy0++,iy1--)
{ p=f.dat[ix][iy0]; f.dat[ix][iy0]=f.dat[ix][iy1]; f.dat[ix][iy1]=p; }
}
//---------------------------------------------------------------------------
void surface2volume(volume &v,surface &s) // normalize surface and convert it to volume
{
point p;
float *pp=(float*)&p;
int ix,iy,iz,nx,ny,nz,mx,my,i,j,k,id[6],us[6];
// find which surface belongs to which face f0..f5 and store it to id[]
// also mirror/rotate to match local x,y orientations
id[0]=0; us[0]=0; // first face hard coded as 0
for (i=1;i<6;i++) id[i]=-1; // all the rest not known yet
for (i=1;i<6;i++) us[i]=-1; // unused
nx=s[0].nx;
ny=s[0].ny;
for (iy=0,j=id[0],k=4,i=1;i<6;i++) // find f4 so it share f0(iy=0) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
}
} if (id[k]<0) return; // fail
nz=s[id[k]].nx; // 3th resolution
v.allocate(nx,ny,nz); // allocate volume container
for (iy=ny-1,j=id[0],k=5,i=1;i<6;i++)// find f5 so it share f0(iy=ny-1) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; break; }
for (ix=0;ix<s[j].nx;ix++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ ix=-1; break; } if (ix>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
}
} if (id[k]<0) return; // fail
for (ix=0,j=id[0],k=2,i=1;i<6;i++) // find f2 so it share f0(ix=0) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
}
} if (id[k]<0) return; // fail
for (ix=nx-1,j=id[0],k=3,i=1;i<6;i++)// find f3 so it share f0(ix=nx-1) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
}
} if (id[k]<0) return; // fail
for (ix=nz-1,j=id[2],k=1,i=1;i<6;i++)// find f1 so it share f2(ix=nz-1) edge
if (us[i]<0) // only in unknonw yet surfaces
{
mx=s[i].nx-1;
my=s[i].ny-1;
if ((s[j].nx==s[i].nx)&&(s[j].ny==s[i].ny)) // try x,y matched
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][ iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-ix][my-iy]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); break; }
}
if ((s[j].nx==s[i].ny)&&(s[j].ny==s[i].nx)) // try x,y swapped
{
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[ iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][ ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); mirror_y(s[i]); swap_xy(s[i]); break; }
for (iy=0;iy<s[j].ny;iy++) if (s[j].dat[ix][iy]!=s[i].dat[mx-iy][my-ix]){ iy=-1; break; } if (iy>=0){ id[k]=i; us[i]=k; mirror_x(s[i]); swap_xy(s[i]); break; }
}
} if (id[k]<0) return; // fail
i=0;
// copy surfaces into volume matrix
for (ix=0;ix<nx;ix++)
for (iy=0;iy<ny;iy++)
for (iz=0;iz<nz;iz++)
{
// interior points
p.x=0.0;
p.y=0.0;
p.z=0.0;
// surface points
if (iz== 0) p=s[id[0]].dat[ix][iy];
if (iz==nz-1) p=s[id[1]].dat[ix][iy];
if (ix== 0) p=s[id[2]].dat[iz][iy];
if (ix==nx-1) p=s[id[3]].dat[iz][iy];
if (iy== 0) p=s[id[4]].dat[iz][ix];
if (iy==ny-1) p=s[id[5]].dat[iz][ix];
v.dat[ix][iy][iz]=p;
}
}
//---------------------------------------------------------------------------
void draw(const surface &s) // render surface
{
int ix,iy,i;
DWORD col[6]= // colors per face
{
0x00202060,
0x00202030,
0x00206020,
0x00203020,
0x00602020,
0x00302020
};
glBegin(GL_LINES);
for (i=0;i<6;i++)
{
glColor4ubv((BYTE*)(&col[i]));
for (ix=0;ix<s[i].nx;ix++)
for (iy=1;iy<s[i].ny;iy++)
{
glVertex3fv((float*)&(s[i].dat[ix][iy ]));
glVertex3fv((float*)&(s[i].dat[ix][iy-1]));
}
for (ix=1;ix<s[i].nx;ix++)
for (iy=0;iy<s[i].ny;iy++)
{
glVertex3fv((float*)&(s[i].dat[ix ][iy]));
glVertex3fv((float*)&(s[i].dat[ix-1][iy]));
}
}
glEnd();
}
//---------------------------------------------------------------------------
void draw(const volume &v)
{
int ix,iy,iz,i;
// all points
glBegin(GL_POINTS);
for (ix=0;ix<v.nx;ix++)
for (iy=0;iy<v.ny;iy++)
for (iz=0;iz<v.nz;iz++)
glVertex3fv((float*)&(v.dat[ix][iy][iz]));
glEnd();
// surface edges for visual check if points are ordered as should
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz= 0;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy=v.ny-1,iz= 0;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy=v.ny-1,iz=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz= 0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy= 0,iz= 0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy= 0,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy= 0,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iy=v.nz-1,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy= 0,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iy=v.nz-1,iz= 0;iz<v.nz;iz++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
}
//---------------------------------------------------------------------------
void cube(surface &s,const param &t) // cube -> surface
{
point p;
float *pp=(float*)&p;
int iu,iv,nx,ny,nz,i,u,v;
// count of points per axis
nx=t[0].nx;
ny=t[1].nx;
nz=t[2].nx;
if ((!nx)||(!ny)||(!nz)) return;
// create 6 faces grid
for (i=0;i<6;i++)
{
if (i==0){ p.z=t[2].dat[ 0]; u=0; v=1; } // z = first coordinate
if (i==1){ p.z=t[2].dat[nz-1]; u=0; v=1; } // z = last coordinate
if (i==2){ p.y=t[0].dat[ 0]; u=0; v=2; } // y = first coordinate
if (i==3){ p.y=t[0].dat[ny-1]; u=0; v=2; } // y = last coordinate
if (i==4){ p.x=t[1].dat[ 0]; u=2; v=1; } // x = first coordinate
if (i==5){ p.x=t[2].dat[nx-1]; u=2; v=1; } // x = last coordinate
s[i].allocate(t[u].nx,t[v].nx); // allocate face resolution
for (iu=0;iu<s[i].nx;iu++) // process all cells
for (iv=0;iv<s[i].ny;iv++)
{
pp[u]=t[u].dat[iu];
pp[v]=t[v].dat[iv];
s[i].dat[iu][iv]=p;
}
}
}
//---------------------------------------------------------------------------
void cubic(surface &s,const param &t) // hardcoded cubic curved object -> surface
{
point p;
float *pp=(float*)&p;
int ix,iy,iz,nx,ny,nz,i,j;
// 3 hard coded cubic control points
float cp[3][4][3]= // [parameter][point][coordinate]
{
{{ 0.00,0.00,0.00 },
{ 0.50,0.50,0.00 },
{ 1.00,0.50,0.00 },
{ 1.50,0.00,0.00 }},
{{ 0.00,0.00,0.00 },
{ 0.10,0.15,0.00 },
{ 0.15,0.30,0.00 },
{ 0.16,0.45,0.00 }},
{{ 0.00,0.00,0.00 },
{ 0.00,0.10,0.15 },
{ 0.00,0.10,0.30 },
{ 0.00,0.00,0.45 }}
};
float ce[3][4][3]; // cubic coeficients
float tx[4],ty[4],tz[4],d1,d2;
for (ix=0;ix<3;ix++)
for (iy=0;iy<3;iy++)
{
d1=0.5*(cp[ix][2][iy]-cp[ix][0][iy]);
d2=0.5*(cp[ix][3][iy]-cp[ix][1][iy]);
ce[ix][0][iy]=cp[ix][1][iy];
ce[ix][1][iy]=d1;
ce[ix][2][iy]=(3.0*(cp[ix][2][iy]-cp[ix][1][iy]))-(2.0*d1)-d2;
ce[ix][3][iy]=d1+d2+(2.0*(-cp[ix][2][iy]+cp[ix][1][iy]));
}
// count of points per axis
nx=t[0].nx;
ny=t[1].nx;
nz=t[2].nx;
if ((!nx)||(!ny)||(!nz)) return;
// allocate faces
s[0].allocate(nx,ny);
s[1].allocate(nx,ny);
s[2].allocate(ny,nz);
s[3].allocate(ny,nz);
s[4].allocate(nz,nx);
s[5].allocate(nz,nx);
// process all points
for (ix=0;ix<nx;ix++)
for (iy=0;iy<ny;iy++)
for (iz=0;iz<nz;iz++)
{
// compute t,t^2,t^3 per each parameter
for (tx[0]=1.0,i=1;i<4;i++) tx[i]=tx[i-1]*t[0].dat[ix];
for (ty[0]=1.0,i=1;i<4;i++) ty[i]=ty[i-1]*t[1].dat[iy];
for (tz[0]=1.0,i=1;i<4;i++) tz[i]=tz[i-1]*t[2].dat[iz];
// compute position as superposition of 3 cubics
for (i=0;i<3;i++)
{
pp[i]=0.0;
for (j=0;j<4;j++) pp[i]+=ce[0][j][i]*tx[j];
for (j=0;j<4;j++) pp[i]+=ce[1][j][i]*ty[j];
for (j=0;j<4;j++) pp[i]+=ce[2][j][i]*tz[j];
}
// copy posiiton to corresponding surfaces
if (iz== 0) s[0].dat[ix][iy]=p;
if (iz==nz-1) s[1].dat[ix][iy]=p;
if (ix== 0) s[2].dat[iy][iz]=p;
if (ix==nx-1) s[3].dat[iy][iz]=p;
if (iy== 0) s[4].dat[iz][ix]=p;
if (iy==ny-1) s[5].dat[iz][ix]=p;
}
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
用法：没有任何组件的简单单窗体 VCL 应用程序
//---------------------------------------------------------------------------
#include <vcl.h> // VCL stuff (ignore)
#pragma hdrstop // VCL stuff (ignore)
#include "Unit1.h" // VCL stuff (header of this window)
#include "gl_simple.h" // my GL init (source included)
#include "parametric_grid.h"
//---------------------------------------------------------------------------
#pragma package(smart_init) // VCL stuff (ignore)
#pragma resource "*.dfm" // VCL stuff (ignore)
TForm1 *Form1; // VCL stuff (this window)
surface sur; // 6 2D matrices (surfaces in any order/orientation before normalization)
volume vol; // 3D matrix (volume)
//---------------------------------------------------------------------------
void gl_draw()
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDisable(GL_CULL_FACE);
glEnable(GL_DEPTH_TEST);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glTranslatef(-0.9,-0.75,-2.0);
glRotatef(45.0,0.5,0.5,0.0);
glPointSize(4.0);
glColor3f(0.1,0.5,0.7); draw(vol);
glPointSize(1.0);
// draw(sur);
glFlush();
SwapBuffers(hdc);
}
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
{
// this is called on window startup
gl_init(Handle); // init OpenGL 1.0
int i;
param t; // source coordinates
/*
t[0].allocate(5); i=0;
t[0][i]=0.000; i++;
t[0][i]=0.125; i++;
t[0][i]=0.250; i++;
t[0][i]=0.500; i++;
t[0][i]=1.000; i++;
t[1].allocate(5); i=0;
t[1][i]=0.000; i++;
t[1][i]=0.250; i++;
t[1][i]=0.500; i++;
t[1][i]=0.750; i++;
t[1][i]=1.000; i++;
t[2].allocate(6); i=0;
t[2][i]=0.0; i++;
t[2][i]=0.2; i++;
t[2][i]=0.4; i++;
t[2][i]=0.6; i++;
t[2][i]=0.8; i++;
t[2][i]=1.0; i++;
*/
t[0].allocate(6); i=0;
t[0][i]=0.0; i++;
t[0][i]=0.2; i++;
t[0][i]=0.4; i++;
t[0][i]=0.6; i++;
t[0][i]=0.8; i++;
t[0][i]=1.0; i++;
t[1].allocate(6); i=0;
t[1][i]=0.0; i++;
t[1][i]=0.2; i++;
t[1][i]=0.4; i++;
t[1][i]=0.6; i++;
t[1][i]=0.8; i++;
t[1][i]=1.0; i++;
t[2].allocate(6); i=0;
t[2][i]=0.0; i++;
t[2][i]=0.2; i++;
t[2][i]=0.4; i++;
t[2][i]=0.6; i++;
t[2][i]=0.8; i++;
t[2][i]=1.0; i++;
// cube(sur,t);
cubic(sur,t);
surface2volume(vol,sur);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
{
// this is called before window exits
gl_exit(); // exit OpenGL
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
{
// this is called on each window resize (and also after startup)
gl_resize(ClientWidth,ClientHeight);
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
{
// this is called whnewer app needs repaint
gl_draw();
}
//---------------------------------------------------------------------------
这是上面代码的预览：
我通过为每个表面尝试 mirror/swap 的所有组合进行调试（首先不按原样使用）并渲染共享边缘（如果它们完全重叠）并调整找到的面的重新排序直到它们完成所以代码应该适用于任何输入（连接）