估计行进立方体输出几何的大小

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);

【讨论】:

以上是关于估计行进立方体输出几何的大小的主要内容,如果未能解决你的问题,请参考以下文章

行进立方体等值

行进立方体 - 获得倾斜的表面

计算行进立方体网格的法线、索引和 UV 坐标

几何的三大问题

学术-几何-维:超级立方体 (五维超级超立方体)

OCAF中Reference-key model概念