有没有一种方法可以快速估计行进立方体算法将产生多少几何体?就像对num vert+num个面或STL文件大小的大致估计一样?
我正在尝试实现一个预先三角测量检查,告诉用户降低分辨率(即使用较低数量的样本,跨越较大的空间区域(,或在无意中尝试生成隐式曲面的多GB-大STL文件网格之前降低域大小,这是任何切片软件都无法处理的。
在最坏的情况下,我可能只生成一大堆具有不同设置的STL来找出这些值。
修改Marching Cubes代码,使其只计算顶点数和三角形数并不困难。我认为你提出的这个选择是有用的。您只需要确定三角形的顶点(以及法线和颜色(和索引存储在代码的哪个部分。在我们编程的C语言的Marching Cubes 33代码(Marching_Cubes_33_C_library_v4,MC33_libraries(的情况下,顶点存储在MC33_storePointNormal函数中。顶点计数器不是调用MC33_storePointNormal,而是简单地增加1。在这种特定情况下,存储三角形索引的代码被删除,只剩下三角形数量的计数器增量。
我做了一个函数,计算顶点和三角形的数量,并计算曲面占用的内存量:
// modified from calculate_isosurface function
unsigned long long memory_of_isosurface (MC33 * M, float iso, unsigned int * nV, unsigned int * nT);
此函数调用MC33_findCase_count函数,而不是调用MC33_findCase函数。您只需要将这两个函数的代码复制到source/marching_ubes_33.c文件的末尾,并将memory_of_isosurface函数的声明放在include/marching_cubes_33.h文件中。
我对只计数的函数进行了基准测试,执行时间减少了12%到32%。
重要提示:由于帖子的大小原因,我删除了部分代码switch(c>>8(。这部分代码必须从MC33_findCase函数中复制。
// modified from MC33_findCase
void MC33_findCase_count(MC33 *M, unsigned int x, unsigned int y, unsigned int z, unsigned int i, float *v)
{
unsigned int p[13] = {FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF};
int f[6];
unsigned int ti[3];
unsigned int c, m, k;
const unsigned short int *pcase = MC33_all_tables;
c = pcase[i&0x80? i^0xFF: i];
m = (i^c)&0x80;
k = c&0x7F;
switch (c>>8)
{
/* copy here the missing code part from the MC33_findCase function */
}
while (i)
{
i = *(++pcase);
for (k = 3; k;)
{
c = i&0x0F;
if (p[c] == FF)
{
switch (c)
{
case 0:
if (z || x)
p[0] = M->Oy[y][x];
else
{
if (v[0] == 0)
{
if (p[3] != FF)
p[0] = p[3];
else if (p[8] != FF)
p[0] = p[8];
else if (y && signbf(v[3]))
p[0] = M->Lz[y][0];
else if (y && signbf(v[4]))
p[0] = M->Ox[y][0];
else if (y? sbf(M->iso - M->F[0][y - 1][0]): 0)
p[0] = M->Oy[y - 1][0];
else
p[0] = M->nV++;
}
else if (v[1] == 0)
{
if (p[1] != FF)
p[0] = p[1];
else if (p[9] != FF)
p[0] = p[9];
else
p[0] = M->nV++;
}
else
p[0] = M->nV++;
M->Oy[y][0] = p[0];
}
break;
case 1:
if (x)
p[1] = M->Lz[y + 1][x];
else
{
if (v[1] == 0)
{
if (p[0] != FF)
p[1] = p[0];
else if (p[9] != FF)
p[1] = p[9];
else if (z && signbf(v[0]))
p[1] = M->Oy[y][0];
else if (z && signbf(v[5]))
p[1] = M->Ox[y + 1][0];
else if (z && y + 1 < M->ny? sbf(M->iso - M->F[z][y + 2][0]): 0)
p[1] = M->Oy[y + 1][0];
else if (z? sbf(M->iso - M->F[z - 1][y + 1][0]): 0)
p[1] = M->Lz[y + 1][0];
else
p[1] = M->nV++;
}
else if (v[2] == 0)
{
if (p[2] != FF)
p[1] = p[2];
else if (p[10] != FF)
p[1] = p[10];
else
p[1] = M->nV++;
}
else
p[1] = M->nV++;
M->Lz[y + 1][0] = p[1];
}
break;
case 2:
if (x)
p[2] = M->Ny[y][x];
else
{
if (v[3] == 0)
{
if (p[3] != FF)
p[2] = p[3];
else if (p[11] != FF)
p[2] = p[11];
else if (y && signbf(v[0]))
p[2] = M->Lz[y][0];
else if (y && signbf(v[7]))
p[2] = M->Nx[y][0];
else if (y? sbf(M->iso - M->F[z + 1][y - 1][0]): 0)
p[2] = M->Ny[y - 1][0];
else
p[2] = M->nV++;
}
else if (v[2] == 0)
{
if (p[1] != FF)
p[2] = p[1];
else if (p[10] != FF)
p[2] = p[10];
else
p[2] = M->nV++;
}
else
p[2] = M->nV++;
M->Ny[y][0] = p[2];
}
break;
case 3:
if (y || x)
p[3] = M->Lz[y][x];
else
{
if (v[0] == 0)
{
if (p[0] != FF)
p[3] = p[0];
else if (p[8] != FF)
p[3] = p[8];
else if (z && signbf(v[1]))
p[3] = M->Oy[0][0];
else if (z && signbf(v[4]))
p[3] = M->Ox[0][0];
else if (z? sbf(M->iso - M->F[z - 1][0][0]): 0)
p[3] = M->Lz[0][0];
else
p[3] = M->nV++;
}
else if (v[3] == 0)
{
if (p[2] != FF)
p[3] = p[2];
else if (p[11] != FF)
p[3] = p[11];
else
p[3] = M->nV++;
}
else
p[3] = M->nV++;
M->Lz[0][0] = p[3];
}
break;
case 4:
if (z)
p[4] = M->Oy[y][x + 1];
else
{
if (v[4] == 0)
{
if (p[8] != FF)
p[4] = p[8];
else if (p[7] != FF)
p[4] = p[7];
else if (y && signbf(v[0]))
p[4] = M->Ox[y][x];
else if (y && signbf(v[7]))
p[4] = M->Lz[y][x + 1];
else if (y? sbf(M->iso - M->F[0][y - 1][x + 1]): 0)
p[4] = M->Oy[y - 1][x + 1];
else if (y && x + 1 < M->nx? sbf(M->iso - M->F[0][y][x + 2]): 0)
p[4] = M->Ox[y][x + 1];
else
p[4] = M->nV++;
}
else if (v[5] == 0)
{
if (p[9] != FF)
p[4] = p[9];
else if (p[5] != FF)
p[4] = p[5];
else
p[4] = M->nV++;
}
else
p[4] = M->nV++;
M->Oy[y][x + 1] = p[4];
}
break;
case 5:
if (v[5] == 0)
{
if (signbf(v[4]))
{
if (z)
{
p[5] = p[4] = M->Oy[y][x + 1];
if (signbf(v[1]))
p[9] = p[5];
}
else
{
if (p[4] != FF)
p[5] = p[4];
else if (p[9] != FF)
p[5] = p[9];
else
p[5] = M->nV++;
}
}
else if (signbf(v[1]))
{
if (z)
p[5] = p[9] = M->Ox[y + 1][x];
else
{
if (p[9] != FF)
p[5] = p[9];
else
p[5] = M->Ox[y + 1][x] = p[9] = M->nV++;
}
}
else
{
if (z && x + 1 < M->nx? sbf(M->iso - M->F[z][y + 1][x + 2]): 0)
p[5] = M->Ox[y + 1][x + 1];
else if (z && y + 1 < M->ny? sbf(M->iso - M->F[z][y + 2][x + 1]): 0)
p[5] = M->Oy[y + 1][x + 1];
else if (z? sbf(M->iso - M->F[z - 1][y + 1][x + 1]): 0)
p[5] = M->Lz[y + 1][x + 1];
else
p[5] = M->nV++;
}
}
else if (v[6] == 0)
{
if (p[6] == FF)
{
if (p[10] == FF)
{
p[5] = M->nV++;
if (signbf(v[2])) M->Nx[y + 1][x] = p[10] = p[5];
if (signbf(v[7])) M->Ny[y][x + 1] = p[6] = p[5];
}
else
{
p[5] = p[10];
if (signbf(v[7])) M->Ny[y][x + 1] = p[6] = p[5];
}
}
else
{
p[5] = p[6];
if (signbf(v[2])) M->Nx[y + 1][x] = p[10] = p[5];
}
}
else
p[5] = M->nV++;
M->Lz[y + 1][x + 1] = p[5];
break;
case 6:
if (v[7] == 0)
{
if (signbf(v[3]))
{
if (y)
{
p[6] = p[11] = M->Nx[y][x];
if (signbf(v[4]))
p[7] = p[6];
}
else
{
if (p[11] != FF)
p[6] = p[11];
else if (p[7] != FF)
p[6] = p[7];
else
p[6] = M->nV++;
}
}
else if (signbf(v[4]))
{
if (y)
p[6] = p[7] = M->Lz[y][x + 1];
else if (p[7] != FF)
p[6] = p[7];
else
p[6] = M->Lz[y][x + 1] = p[7] = M->nV++;
}
else
{
if (y? sbf(M->iso - M->F[z + 1][y - 1][x + 1]): 0)
p[6] = M->Ny[y - 1][x + 1];
else if (y && x + 1 < M->nx? sbf(M->iso - M->F[z + 1][y][x + 2]): 0)
p[6] = M->Nx[y][x + 1];
else
p[6] = M->nV++;
}
}
else if (v[6] == 0)
{
if (p[5] == FF)
{
if (p[10] == FF)
{
p[6] = M->nV++;
if (signbf(v[5])) M->Lz[y + 1][x + 1] = p[5] = p[6];
if (signbf(v[2])) M->Nx[y + 1][x] = p[10] = p[6];
}
else
{
p[6] = p[10];
if (signbf(v[5])) M->Lz[y + 1][x + 1] = p[5] = p[6];
}
}
else
{
p[6] = p[5];
if (signbf(v[2])) M->Nx[y + 1][x] = p[10] = p[6];
}
}
else
p[6] = M->nV++;
M->Ny[y][x + 1] = p[6];
break;
case 7:
if (y)
p[7] = M->Lz[y][x + 1];
else
{
if (v[4] == 0)
{
if (p[8] != FF)
p[7] = p[8];
else if (p[4] != FF)
p[7] = p[4];
else if (z && signbf(v[0]))
p[7] = M->Ox[0][x];
else if (z && x + 1 < M->nx? sbf(M->iso - M->F[z][0][x + 2]): 0)
p[7] = M->Ox[0][x + 1];
else if (z? sbf(M->iso - M->F[z - 1][0][x + 1]): 0)
p[7] = M->Lz[0][x + 1];
else
p[7] = M->nV++;
}
else if (v[7] == 0)
{
if (p[6] != FF)
p[7] = p[6];
else if (p[11] != FF)
p[7] = p[11];
else
p[7] = M->nV++;
}
else
p[7] = M->nV++;
M->Lz[0][x + 1] = p[7];
}
break;
case 8:
if (z || y)
p[8] = M->Ox[y][x];
else
{
if (v[0] == 0)
{
if (p[3] != FF)
p[8] = p[3];
else if (p[0] != FF)
p[8] = p[0];
else if (x && signbf(v[3]))
p[8] = M->Lz[0][x];
else if (x && signbf(v[1]))
p[8] = M->Oy[0][x];
else if (x? sbf(M->iso - M->F[0][0][x - 1]): 0)
p[8] = M->Ox[0][x - 1];
else
p[8] = M->nV++;
}
else if (v[4] == 0)
{
if (p[4] != FF)
p[8] = p[4];
else if (p[7] != FF)
p[8] = p[7];
else
p[8] = M->nV++;
}
else
p[8] = M->nV++;
M->Ox[0][x] = p[8];
}
break;
case 9:
if (z)
p[9] = M->Ox[y + 1][x];
else
{
if (v[1] == 0)
{
if (p[0] != FF)
p[9] = p[0];
else if (p[1] != FF)
p[9] = p[1];
else if (x && signbf(v[0]))
p[9] = M->Oy[y][x];
else if (x && signbf(v[2]))
p[9] = M->Lz[y + 1][x];
else if (x? sbf(M->iso - M->F[0][y + 1][x - 1]): 0)
p[9] = M->Ox[y + 1][x - 1];
else
p[9] = M->nV++;
}
else if (v[5] == 0)
{
if (p[5] != FF)
p[9] = p[5];
else if (p[4] != FF)
p[9] = p[4];
else
p[9] = M->nV++;
}
else
p[9] = M->nV++;
M->Ox[y + 1][x] = p[9];
}
break;
case 10:
if (v[2] == 0)
{
if (signbf(v[1]))
{
if (x)
{
p[10] = p[1] = M->Lz[y + 1][x];
if (signbf(v[3])) p[2] = p[10];
}
else
{
if (p[1] != FF)
p[10] = p[1];
else if (p[2] != FF)
p[10] = p[2];
else
p[10] = M->nV++;
}
}
else if (signbf(v[3]))
{
if (x)
p[10] = p[2] = M->Ny[y][x];
else
p[10] = M->nV++;
}
else
{
if (x? sbf(M->iso - M->F[z + 1][y + 1][x - 1]): 0)
p[10] = M->Nx[y + 1][x - 1];
else
p[10] = M->nV++;
}
}
else if (v[6] == 0)
{
if (p[6] == FF)
{
if (p[5] == FF)
{
p[10] = M->nV++;
if (signbf(v[5])) M->Lz[y + 1][x + 1] = p[5] = p[10];
if (signbf(v[7])) M->Ny[y][x + 1] = p[6] = p[10];
}
else
{
p[10] = p[5];
if (signbf(v[7])) M->Ny[y][x + 1] = p[6] = p[10];
}
}
else
{
p[10] = p[6];
if (signbf(v[5])) M->Lz[y + 1][x + 1] = p[5] = p[10];
}
}
else
p[10] = M->nV++;
M->Nx[y + 1][x] = p[10];
break;
case 11:
if (y)
p[11] = M->Nx[y][x];
else
{
if (v[3] == 0)
{
if (p[3] != FF)
p[11] = p[3];
else if (p[2] != FF)
p[11] = p[2];
else if (x && signbf(v[0]))
p[11] = M->Lz[0][x];
else if (x? sbf(M->iso - M->F[z + 1][0][x - 1]): 0)
p[11] = M->Nx[0][x - 1];
else
p[11] = M->nV++;
}
else if (v[7] == 0)
{
if (p[6] != FF)
p[11] = p[6];
else if (p[7] != FF)
p[11] = p[7];
else
p[11] = M->nV++;
}
else
p[11] = M->nV++;
M->Nx[0][x] = p[11];
}
break;
default:
p[12] = M->nV++;
}
}
ti[--k] = p[c];
i >>= 4;
}
if (ti[0] != ti[1] && ti[0] != ti[2] && ti[1] != ti[2])
M->nT++;
}
}
// modified from calculate_isosurface function
unsigned long long memory_of_isosurface(MC33 *M, float iso, unsigned int *nV, unsigned int *nT)
{
unsigned int x, y, z, nx = M->nx, ny = M->ny, nz = M->nz;
GRD_data_type ***F = M->F, **F0, **F1, *V00, *V01, *V11, *V10;
float V[12];
float *v1 = V, *v2 = V + 4;
M->nT = M->nV = 0;
M->iso = iso;
for (z = 0; z != nz; ++z)
{
F0 = *F;
F1 = *(++F);
for (y = 0; y != ny; ++y)
{
V00 = *F0;
V01 = *(++F0);
V10 = *F1;
V11 = *(++F1);
v2[0] = iso - *V00;
v2[1] = iso - *V01;
v2[2] = iso - *V11;
v2[3] = iso - *V10;
unsigned int i = signbf(v2[3]) != 0;
if (signbf(v2[2])) i |= 2;
if (signbf(v2[1])) i |= 4;
if (signbf(v2[0])) i |= 8;
for (x = 0; x != nx; ++x)
{
{float *P = v1; v1 = v2; v2 = P;}
v2[0] = iso - *(++V00);
v2[1] = iso - *(++V01);
v2[2] = iso - *(++V11);
v2[3] = iso - *(++V10);
i = ((i&0x0F)<<4)|(signbf(v2[3]) != 0);
if (signbf(v2[2])) i |= 2;
if (signbf(v2[1])) i |= 4;
if (signbf(v2[0])) i |= 8;
if (i && i^0xFF)
{
if (v1 > v2) {float *t = v2; float *s = t + 8; *s = *t; *(++s) = *(++t); *(++s) = *(++t); *(++s) = *(++t);}
MC33_findCase_count(M,x,y,z,i,v1);
}
}
}
{unsigned int** P = M->Ox; M->Ox = M->Nx; M->Nx = P;}
{unsigned int** P = M->Oy; M->Oy = M->Ny; M->Ny = P;}
}
*nV = M->nV;
*nT = M->nT;
// number of vertices * (size of vertex and normal + size of color ) + number of triangle * size of triangle + size of struct surface
return (*nV) * (6 * sizeof(float) + sizeof(int)) + (*nT) * (3 * sizeof(int)) + sizeof(surface);
}