Несложно изменить код Marching Cubes так, чтобы он подсчитывал только количество вершин и количество треугольников. Я думаю, что тот вариант, который вы предлагаете, может быть полезен. Вам нужно только определить, в какой части кода хранятся вершины (а также нормали и цвета) и индексы треугольников. В случае кода Marching Cubes 33 для языка C, который мы запрограммировали (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%.
Важное примечание: из соображений размера поста я удалил часть кода, переключатель (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);
}
person
David Vega
schedule
18.01.2021
