问题描述
我正在尝试实现一种算法，该算法将 X、Y、Z 点的六个二维矩阵存储在 X、Y、Z 点的 3D 矩阵中，以保持它们的连接性。一个例子将说明问题。
我必须用三个 3-D 矩阵来表示
( M3D_X(:,:,:),M3D_Y(:,M_3D_Z(:,:) )
空间的六面体区域。该六面体区域将只知道它的面，存储在三个二维矩阵中
( [face1_x(:,face1_y(:,face1_z(:,:)],[face2_x.....]....,],...,[...,face6_Z(:,:)] )
关于人脸之间连通性的唯一已知信息是 face_i 与 face_j 有共同的一面。没有给出哪一方是共同的。 那么，如何在 3-D 矩阵中存储 2-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,1.0]
face1_y = [0.0,0.0,0.0;
0.5,0.5;
1.0,1.0,1.0]
face1_z = [0.0,0.0;
0.0,0.0]
face2_x = [0.0,1.0]
face2_y = [1.0,1.0;
1.0,1.0]
face2_z = [0.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,1.0] ])
face1y=np.array( [ [0.0,0.0],[0.3,0.3,0.3],[0.7,0.7,0.7],[1.0,1.0]])
face1z=np.array( [ [0.0,0.0]] )
face2x=np.array( [ [0.0,1.0] ])
face2y=np.array( [ [0.0,0.0]] )
face2z=np.array( [ [0.0,[0.5,0.5],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,3 ;
*,*,* ;
*,* ]
现在如果 line2 像这样插入
face(:,3)=line2
结果是
face = [1,5 ;
*,4 ;
*,3 ]
当正确的解决方案是
face(:,3)= reverse(line2)
（因为点 3 与 line1 和 line2 相同）得到
face = [1,5 ]
希望这能让问题更清楚。
解决方法
所以，如果我做对了，你会得到 6 个 3D 点的 2D 矩阵，它们代表 6 个形状（连接在一起）的表面，处于彼此之间的任何翻转/镜像/有序状态，并且想要构建具有空内部的单个 3D 矩阵（例如，点 (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代表每个二维矩阵内的人脸（表面/面 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复制到它们在矩阵中的位置...您可以用(0,0)填充内部点或在曲面之间进行插值...
这里是 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
//---------------------------------------------------------------------------
gl_simple.h 只是简单的 OpenGL
parametric_grid.h这是您需要消化的主要来源
//---------------------------------------------------------------------------
#ifndef _parametric_grid_h
#define _parametric_grid_h
//---------------------------------------------------------------------------
#include "arrays.h"
//---------------------------------------------------------------------------
struct point
{
float x,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,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,nz); // allocate volume container
for (iy=ny-1,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 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,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,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 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,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,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,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=v.nz-1;ix<v.nx;ix++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 0,iz= 0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 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,iz= 0;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix=v.nx-1,iz=v.nz-1;iy<v.ny;iy++) glVertex3fv((float*)&(v.dat[ix][iy][iz])); glEnd();
glBegin(GL_LINE_STRIP); for (ix= 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,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,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,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,{ 1.00,{ 1.50,0.00 }},{{ 0.00,{ 0.10,0.15,{ 0.15,0.30,{ 0.16,0.45,{ 0.00,0.10,0.15 },0.30 },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.0);
glPointSize(4.0);
glColor3f(0.1,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();
}
//---------------------------------------------------------------------------
这是上面代码的预览：
我通过尝试每个表面的镜像/交换的所有组合进行调试（首先不计为按原样使用）并渲染共享边缘（如果它们完全重叠）并调整找到的面的重新排序，直到它们这样做，代码应该为任何输入（连接）工作