估计行进立方体输出几何的大小
Posted
技术标签:
【中文标题】估计行进立方体输出几何的大小【英文标题】：Estimating size of marching cubes output geometry 【发布时间】：2020-11-29 23:13:48 【问题描述】：有没有一种方法可以快速估计行进立方体算法将产生多少几何图形？就像对 num verts + num faces 或 STL 文件大小的粗略估计一样？
我正在尝试实施预三角检查，告诉用户降低分辨率（即使用更少的样本数，跨越更大的空间区域），或者在无意中尝试生成多区域之前减小域大小GB-large STL 文件的隐式曲面网格，任何切片软件都无法处理。
在最坏的情况下，我可能会生成一大堆具有不同设置的 STL 来找出这些值。
【问题讨论】：
【参考方案1】：修改 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_cubes_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);
【讨论】：
以上是关于估计行进立方体输出几何的大小的主要内容，如果未能解决你的问题，请参考以下文章