mirror of
https://github.com/bkaradzic/bgfx.git
synced 2026-10-29 12:15:30 +00:00
752 lines
25 KiB
C++
752 lines
25 KiB
C++
// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
|
|
#include "meshoptimizer.h"
|
|
|
|
#include <assert.h>
|
|
#include <float.h>
|
|
#include <math.h>
|
|
#include <string.h>
|
|
|
|
#ifndef TRACE
|
|
#define TRACE 0
|
|
#endif
|
|
|
|
#if TRACE
|
|
#include <stdio.h>
|
|
#endif
|
|
|
|
// Note: this is only exposed for development purposes; do *not* use
|
|
enum
|
|
{
|
|
meshopt_RemeshInternalDebug = 1 << 30
|
|
};
|
|
|
|
// This work is based on:
|
|
// William Lorensen, Harvey Cline. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. 1987
|
|
// Michael Garland and Paul S. Heckbert. Surface simplification using quadric error metrics. 1997
|
|
namespace meshopt
|
|
{
|
|
|
|
// triangles are encoded as three vertex codes in three nibbles
|
|
// vertex codes match the bit order XYZ (e.g. 3 maps to X=1 Y=1 Z=0)
|
|
struct Case
|
|
{
|
|
unsigned char code;
|
|
unsigned char altmode;
|
|
unsigned short triangles[5]; // up to 4 triangles and a null terminator
|
|
unsigned short alternate[5]; // alternate triangulation, restricted to have <= triangles vs primary
|
|
};
|
|
|
|
// the base cases establish a canonical triangulation; unlike classical Marching Cubes, we connect cube *corners* instead of edges
|
|
// each case may have an alternate triangulation, used based on deciders; crucially, open boundary edges must stay consistent as decider varies per cell
|
|
static const Case kBaseCases[] = {
|
|
{0x00, 0, {}, {}},
|
|
{0x01, 0, {}, {}},
|
|
{0x03, 0, {}, {}},
|
|
{0x06, 0, {}, {}},
|
|
{0x07, 0, {0x120}, {}},
|
|
{0x0f, 0, {0x130, 0x320}, {}},
|
|
{0x16, 1, {0x124, 0x142}, {}},
|
|
{0x17, 0, {0x124}, {}},
|
|
{0x18, 0, {}, {}},
|
|
{0x19, 1, {0x043, 0x403}, {}},
|
|
{0x1b, 1, {0x134, 0x430}, {0x410, 0x013}},
|
|
{0x1d, 1, {0x243, 0x340}, {0x320, 0x024}},
|
|
{0x1e, 1, {0x142, 0x432, 0x341}, {0x132}},
|
|
{0x1f, 2, {0x412, 0x213}, {0x134, 0x432}},
|
|
{0x3c, 1, {0x345, 0x324, 0x543, 0x423}, {}},
|
|
{0x3d, 1, {0x253, 0x524, 0x350}, {0x240, 0x504, 0x320}},
|
|
{0x3f, 0, {0x253, 0x452}, {}},
|
|
{0x69, 1, {0x035, 0x560, 0x306, 0x536}, {}},
|
|
{0x6b, 1, {0x063, 0x605, 0x365}, {0x305}},
|
|
{0x6f, 1, {0x365, 0x560}, {0x530, 0x036}},
|
|
{0x7e, 0, {0x142, 0x536}, {}},
|
|
{0x7f, 0, {0x365}, {}},
|
|
{0xf0, 0, {0x547, 0x746}, {}},
|
|
{0xff, 0, {}, {}},
|
|
};
|
|
|
|
// for each cube case, a triangulation is derived from the base cases via 90 degree rotations
|
|
static unsigned short kTriangleTable[256][2][5];
|
|
static unsigned char kTriangleCount[256];
|
|
static unsigned char kTriangleAlt[256];
|
|
|
|
static bool buildRemeshTables()
|
|
{
|
|
// corner remap for 90 degree rotations around X and Y
|
|
static const unsigned char rotations[2][8] = {{2, 3, 6, 7, 0, 1, 4, 5}, {4, 0, 6, 2, 5, 1, 7, 3}};
|
|
|
|
unsigned char filled[256] = {};
|
|
|
|
// copy base cases as is; we will generate the rest via rotations
|
|
for (size_t i = 0; i < sizeof(kBaseCases) / sizeof(kBaseCases[0]); ++i)
|
|
{
|
|
const Case& c = kBaseCases[i];
|
|
memcpy(kTriangleTable[c.code][0], c.triangles, sizeof(c.triangles));
|
|
memcpy(kTriangleTable[c.code][1], c.altmode ? c.alternate : c.triangles, sizeof(c.triangles));
|
|
kTriangleAlt[c.code] = c.altmode;
|
|
filled[kBaseCases[i].code] = 1;
|
|
}
|
|
|
|
// propagate base cases iteratively to fill all cube cases via rotations
|
|
// note: the worst case is 5 passes but in practice we only need 3 to fill all cases
|
|
for (int pass = 0; pass < 3; ++pass)
|
|
for (int code = 0; code < 256; ++code)
|
|
{
|
|
if (filled[code] != 1)
|
|
continue;
|
|
|
|
for (int r = 0; r < 2; ++r)
|
|
{
|
|
unsigned char rotated = 0;
|
|
for (int i = 0; i < 8; ++i)
|
|
rotated |= ((code >> i) & 1) << rotations[r][i];
|
|
|
|
if (filled[rotated])
|
|
continue;
|
|
|
|
for (int k = 0; k < 2; ++k)
|
|
for (int i = 0; kTriangleTable[code][k][i]; ++i)
|
|
{
|
|
unsigned short tri = kTriangleTable[code][k][i];
|
|
kTriangleTable[rotated][k][i] = (rotations[r][(tri >> 8) & 0xf] << 8) | (rotations[r][(tri >> 4) & 0xf] << 4) | rotations[r][tri & 0xf];
|
|
}
|
|
|
|
kTriangleAlt[rotated] = kTriangleAlt[code];
|
|
filled[rotated] = 1; // mark as pending
|
|
}
|
|
|
|
filled[code] = 2; // mark as processed
|
|
}
|
|
|
|
// finalize auxiliary tables
|
|
for (int code = 0; code < 256; ++code)
|
|
{
|
|
assert(filled[code]);
|
|
|
|
int count = 0;
|
|
while (kTriangleTable[code][0][count])
|
|
count++;
|
|
|
|
kTriangleCount[code] = (unsigned char)count;
|
|
|
|
// counting pass does not use alternate triangulations, so they need to stay at or below primary triangle count
|
|
assert(kTriangleTable[code][1][count] == 0);
|
|
|
|
// quadric decider expects a canonical quad encoding (0xabc 0xcbd)
|
|
if (kTriangleAlt[code] == 2)
|
|
assert(count == 2 && (kTriangleTable[code][0][0] & 0xf0) == (kTriangleTable[code][0][1] & 0xf0) && (kTriangleTable[code][0][0] & 0xf) == (kTriangleTable[code][0][1] >> 8));
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
static bool gRemeshTablesInitialized = buildRemeshTables();
|
|
|
|
struct Voxel
|
|
{
|
|
unsigned int coord;
|
|
unsigned char octants;
|
|
|
|
float px, py, pz;
|
|
float w;
|
|
|
|
// a00*x^2 + a11*y^2 + a22*z^2 + 2*a10*xy + 2*a20*xz + 2*a21*yz + 2*b0*x + 2*b1*y + 2*b2*z + c
|
|
float a00, a11, a22;
|
|
float a10, a20, a21;
|
|
float b0, b1, b2, c;
|
|
};
|
|
|
|
static void voxelAccumulate(Voxel& vox, float px, float py, float pz, float w)
|
|
{
|
|
vox.px += px * w;
|
|
vox.py += py * w;
|
|
vox.pz += pz * w;
|
|
vox.w += w;
|
|
}
|
|
|
|
static void voxelAccumulateQ(Voxel& vox, float a, float b, float c, float d, float w)
|
|
{
|
|
float aw = a * w;
|
|
float bw = b * w;
|
|
float cw = c * w;
|
|
float dw = d * w;
|
|
|
|
vox.a00 += a * aw;
|
|
vox.a11 += b * bw;
|
|
vox.a22 += c * cw;
|
|
vox.a10 += a * bw;
|
|
vox.a20 += a * cw;
|
|
vox.a21 += b * cw;
|
|
vox.b0 += a * dw;
|
|
vox.b1 += b * dw;
|
|
vox.b2 += c * dw;
|
|
vox.c += d * dw;
|
|
}
|
|
|
|
static float voxelError(const Voxel& vox, float x, float y, float z)
|
|
{
|
|
float rx = (vox.b0 + vox.a10 * y) * 2.f + vox.a00 * x;
|
|
float ry = (vox.b1 + vox.a21 * z) * 2.f + vox.a11 * y;
|
|
float rz = (vox.b2 + vox.a20 * x) * 2.f + vox.a22 * z;
|
|
|
|
return fabsf(vox.c + rx * x + ry * y + rz * z);
|
|
}
|
|
|
|
static bool voxelSolve(float& rx, float& ry, float& rz, const Voxel& vox, float lambda)
|
|
{
|
|
// solve A*p = -b where A is the quadric matrix and b is the linear term
|
|
float rw = lambda * vox.w;
|
|
float a00 = rw + vox.a00, a11 = rw + vox.a11, a22 = rw + vox.a22;
|
|
float a10 = vox.a10, a20 = vox.a20, a21 = vox.a21;
|
|
float x0 = lambda * vox.px - vox.b0, x1 = lambda * vox.py - vox.b1, x2 = lambda * vox.pz - vox.b2;
|
|
|
|
float eps = 1e-6f * vox.w;
|
|
|
|
// LDL decomposition: A = LDL^T
|
|
float d0 = a00;
|
|
float l10 = a10 / d0;
|
|
float l20 = a20 / d0;
|
|
|
|
float d1 = a11 - a10 * l10;
|
|
float dl21 = a21 - a20 * l10;
|
|
float l21 = dl21 / d1;
|
|
|
|
float d2 = a22 - a20 * l20 - dl21 * l21;
|
|
|
|
// solve L*y = x
|
|
float y0 = x0;
|
|
float y1 = x1 - l10 * y0;
|
|
float y2 = x2 - l20 * y0 - l21 * y1;
|
|
|
|
// solve D*z = y
|
|
float z0 = y0 / d0;
|
|
float z1 = y1 / d1;
|
|
float z2 = y2 / d2;
|
|
|
|
// substitute L^T*p = z
|
|
rz = z2;
|
|
ry = z1 - l21 * rz;
|
|
rx = z0 - l10 * ry - l20 * rz;
|
|
|
|
return fabsf(d0) > eps && fabsf(d1) > eps && fabsf(d2) > eps;
|
|
}
|
|
|
|
static float measureGrid(const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, int resolution, float* out_offset)
|
|
{
|
|
size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
|
|
|
|
float minv[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
|
|
float maxv[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
|
|
|
|
for (size_t i = 0; i < vertex_count; ++i)
|
|
{
|
|
const float* v = vertex_positions_data + i * vertex_stride_float;
|
|
|
|
for (int j = 0; j < 3; ++j)
|
|
{
|
|
float vj = v[j];
|
|
|
|
minv[j] = minv[j] > vj ? vj : minv[j];
|
|
maxv[j] = maxv[j] < vj ? vj : maxv[j];
|
|
}
|
|
}
|
|
|
|
float extent = 0.f;
|
|
|
|
extent = (maxv[0] - minv[0]) < extent ? extent : (maxv[0] - minv[0]);
|
|
extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]);
|
|
extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]);
|
|
|
|
// small subvoxel padding keeps vertices inside the last valid voxel to avoid clamping correctness issues
|
|
extent *= (resolution + 1e-2f) / float(resolution);
|
|
|
|
// rescale extents to [0..resolution - 2], because the first and last voxels are used as empty padding
|
|
float scale = extent == 0 ? 0.f : (resolution - 2) / extent;
|
|
|
|
// center mesh in the grid; this improves voxelization for symmetric inputs
|
|
for (int j = 0; j < 3; ++j)
|
|
out_offset[j] = minv[j] - (extent - (maxv[j] - minv[j])) * 0.5f;
|
|
|
|
return scale;
|
|
}
|
|
|
|
static void voxelize(unsigned char* grid, Voxel* voxels, const unsigned int* voxel_rows, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, float scale, const float offset[3], unsigned int options)
|
|
{
|
|
(void)vertex_count;
|
|
|
|
size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
|
|
|
|
for (size_t i = 0; i < index_count; i += 3)
|
|
{
|
|
unsigned int a = indices[i + 0], b = indices[i + 1], c = indices[i + 2];
|
|
assert(a < vertex_count && b < vertex_count && c < vertex_count);
|
|
|
|
const float* va = vertex_positions + a * vertex_stride_float;
|
|
const float* vb = vertex_positions + b * vertex_stride_float;
|
|
const float* vc = vertex_positions + c * vertex_stride_float;
|
|
|
|
float ex = vb[0] - va[0], ey = vb[1] - va[1], ez = vb[2] - va[2];
|
|
float fx = vc[0] - va[0], fy = vc[1] - va[1], fz = vc[2] - va[2];
|
|
float gx = vc[0] - vb[0], gy = vc[1] - vb[1], gz = vc[2] - vb[2];
|
|
|
|
// use maximum edge length to establish sampling rate
|
|
// TODO: this is wasteful for thin triangles
|
|
float el = sqrtf(ex * ex + ey * ey + ez * ez);
|
|
float fl = sqrtf(fx * fx + fy * fy + fz * fz);
|
|
float gl = sqrtf(gx * gx + gy * gy + gz * gz);
|
|
|
|
float max_edge = el > fl ? el : fl;
|
|
max_edge = max_edge > gl ? max_edge : gl;
|
|
|
|
// we target 2 samples per voxel edge which should be enough to hit all voxels at any rotation
|
|
int samples = int(max_edge * scale * 2.f);
|
|
samples = samples > 1 ? samples : 1;
|
|
samples = samples < resolution * 2 ? samples : resolution * 2;
|
|
|
|
// normal is used to compute area for weighting as well as for the quadric
|
|
float nx = ey * fz - ez * fy, ny = ez * fx - ex * fz, nz = ex * fy - ey * fx;
|
|
float area = sqrtf(nx * nx + ny * ny + nz * nz);
|
|
|
|
float ns = area == 0.f ? 0.f : 1.f / area;
|
|
nx *= ns;
|
|
ny *= ns;
|
|
nz *= ns;
|
|
|
|
float sx = va[0] - offset[0], sy = va[1] - offset[1], sz = va[2] - offset[2];
|
|
float sr = 1.f / float(samples);
|
|
float weight = area / float((samples + 1) * (samples + 2));
|
|
|
|
// skip degenerate triangles; they don't contribute to voxelization and may result in NaN when computing voxel centroid
|
|
if (weight == 0.f)
|
|
continue;
|
|
|
|
for (int u = 0; u <= samples; ++u)
|
|
for (int v = 0; v <= samples - u; ++v)
|
|
{
|
|
float su = float(u) * sr, sv = float(v) * sr;
|
|
float px = sx + su * ex + sv * fx;
|
|
float py = sy + su * ey + sv * fy;
|
|
float pz = sz + su * ez + sv * fz;
|
|
|
|
// compute coordinates with an extra bit to use in octant mask
|
|
int hx = int(px * (scale * 2));
|
|
int hy = int(py * (scale * 2));
|
|
int hz = int(pz * (scale * 2));
|
|
int x = hx >> 1, y = hy >> 1, z = hz >> 1;
|
|
|
|
// safety: rounding errors and non-finite inputs may produce out of bounds coordinates, so we clamp them
|
|
int cutoff = resolution - 3;
|
|
x = unsigned(x) < unsigned(cutoff) ? x : cutoff;
|
|
y = unsigned(y) < unsigned(cutoff) ? y : cutoff;
|
|
z = unsigned(z) < unsigned(cutoff) ? z : cutoff;
|
|
|
|
size_t row = (y + 1) + size_t(resolution) * (z + 1);
|
|
size_t idx = (x + 1) + size_t(resolution) * row;
|
|
|
|
if (voxels)
|
|
{
|
|
assert(grid[idx] != 0 && grid[idx] != 0xff);
|
|
Voxel& vox = voxels[voxel_rows[row] + (grid[idx] - 1)];
|
|
|
|
vox.coord = (unsigned(x) << 20) | (unsigned(y) << 10) | unsigned(z);
|
|
vox.octants |= 1 << ((hx & 1) | ((hy & 1) << 1) | ((hz & 1) << 2));
|
|
|
|
voxelAccumulate(vox, px, py, pz, weight);
|
|
|
|
if (options & meshopt_RemeshSolve)
|
|
{
|
|
float distance = nx * px + ny * py + nz * pz;
|
|
voxelAccumulateQ(vox, nx, ny, nz, -distance, weight);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
grid[idx] = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
static size_t rowpack(unsigned char* grid, unsigned int* voxel_rows, int resolution)
|
|
{
|
|
size_t result = 0;
|
|
size_t slice = size_t(resolution) * size_t(resolution);
|
|
|
|
for (size_t i = 0; i < slice; ++i)
|
|
{
|
|
unsigned char* data = grid + i * size_t(resolution);
|
|
|
|
int count = 0;
|
|
|
|
for (int x = 0; x < resolution; ++x)
|
|
{
|
|
assert(data[x] <= 1); // voxelize only produces 0/1 values
|
|
|
|
count += data[x];
|
|
data[x] = (data[x] != 0) ? (unsigned char)count : 0;
|
|
}
|
|
|
|
assert(count < 255); // we store offsets in a single byte, with 0 reserved for empty voxels and 0xff reserved for interior voxels
|
|
|
|
voxel_rows[i] = unsigned(result);
|
|
result += count;
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
static void solidifyQueue(unsigned int row, unsigned int* worklist, unsigned char* queued, size_t& pending)
|
|
{
|
|
if (queued[row])
|
|
return;
|
|
|
|
queued[row] = 1;
|
|
worklist[pending++] = row;
|
|
}
|
|
|
|
static void solidify(unsigned char* grid, unsigned int* worklist, unsigned char* queued, int resolution)
|
|
{
|
|
size_t pending = 0;
|
|
memset(queued, 0, size_t(resolution) * size_t(resolution));
|
|
|
|
// mark the interior empty voxels as 'inside'; we will propagate 'empty' state from the boundary
|
|
for (int z = 1; z < resolution - 1; ++z)
|
|
for (int y = 1; y < resolution - 1; ++y)
|
|
{
|
|
unsigned char* data = grid + size_t(resolution) * (y + size_t(resolution) * z);
|
|
|
|
for (int x = 1; x < resolution - 1; ++x)
|
|
data[x] = (data[x] == 0) ? 0xff : data[x];
|
|
}
|
|
|
|
// queue all rows for processing; note that boundary rows are always empty but we need them to propagate into neighboring rows
|
|
for (int z = 0; z < resolution; ++z)
|
|
for (int y = 0; y < resolution; ++y)
|
|
solidifyQueue(y + size_t(resolution) * z, worklist, queued, pending);
|
|
|
|
while (pending)
|
|
{
|
|
unsigned int row = worklist[--pending];
|
|
assert(queued[row]);
|
|
queued[row] = 0;
|
|
|
|
unsigned char* data = grid + size_t(resolution) * row;
|
|
|
|
// propagate outside state to the interior within row
|
|
for (int x = 1; x < resolution - 1; ++x)
|
|
data[x] = (data[x] == 0xff && data[x - 1] == 0) ? 0 : data[x];
|
|
|
|
for (int x = resolution - 2; x >= 1; --x)
|
|
data[x] = (data[x] == 0xff && data[x + 1] == 0) ? 0 : data[x];
|
|
|
|
// propagate outside state to the interior into neighboring rows
|
|
int y = row % resolution, z = row / resolution;
|
|
|
|
for (int k = 0; k < 4; ++k)
|
|
{
|
|
int yn = y + (k == 0 ? -1 : (k == 1 ? 1 : 0));
|
|
int zn = z + (k == 2 ? -1 : (k == 3 ? 1 : 0));
|
|
|
|
if (yn < 1 || yn >= resolution - 1 || zn < 1 || zn >= resolution - 1)
|
|
continue;
|
|
|
|
unsigned char* datan = grid + size_t(resolution) * (yn + size_t(resolution) * zn);
|
|
unsigned char changed = 0;
|
|
|
|
for (int x = 1; x < resolution - 1; ++x)
|
|
{
|
|
unsigned char rep = ((data[x] == 0) & (datan[x] == 0xff)) ? 0 : datan[x];
|
|
changed |= rep ^ datan[x];
|
|
datan[x] = rep;
|
|
}
|
|
|
|
if (changed)
|
|
solidifyQueue(yn + size_t(resolution) * zn, worklist, queued, pending);
|
|
}
|
|
}
|
|
}
|
|
|
|
static void solve(Voxel* voxels, size_t voxel_count, float scale, unsigned int options)
|
|
{
|
|
// regularization factor is a tradeoff between solve stability/uniformity and adherence to the original surface
|
|
float factor = 3e-2f;
|
|
|
|
// cutoff is a squared distance in model space; sqrt(3)/scale is the voxel diagonal
|
|
float cutoff = 3.f / (scale * scale);
|
|
float rscale = 1.f / scale;
|
|
|
|
size_t solved = 0;
|
|
|
|
for (size_t i = 0; i < voxel_count; ++i)
|
|
{
|
|
Voxel& vox = voxels[i];
|
|
float px = vox.px / vox.w, py = vox.py / vox.w, pz = vox.pz / vox.w;
|
|
|
|
// voxel corner, used for clamping and debug visualization
|
|
float cx = int((vox.coord >> 20) & 1023) * rscale;
|
|
float cy = int((vox.coord >> 10) & 1023) * rscale;
|
|
float cz = int((vox.coord >> 0) & 1023) * rscale;
|
|
|
|
if (options & meshopt_RemeshSolve)
|
|
{
|
|
// compute minimizing point; this includes regularization to stabilize the solve
|
|
float sx, sy, sz;
|
|
if (voxelSolve(sx, sy, sz, vox, factor))
|
|
{
|
|
float d = (sx - px) * (sx - px) + (sy - py) * (sy - py) + (sz - pz) * (sz - pz);
|
|
|
|
// reject solutions that move the vertex too far from voxel centroid; this is a safety net in case regularization is insufficient
|
|
if (d < cutoff)
|
|
{
|
|
// use a clamped solution to reduce intersections
|
|
px = sx > cx ? sx : cx;
|
|
px = px < cx + rscale ? px : cx + rscale;
|
|
py = sy > cy ? sy : cy;
|
|
py = py < cy + rscale ? py : cy + rscale;
|
|
pz = sz > cz ? sz : cz;
|
|
pz = pz < cz + rscale ? pz : cz + rscale;
|
|
|
|
solved++;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (options & meshopt_RemeshInternalDebug)
|
|
{
|
|
px = cx + 0.5f * rscale;
|
|
py = cy + 0.5f * rscale;
|
|
pz = cz + 0.5f * rscale;
|
|
}
|
|
|
|
vox.px = px;
|
|
vox.py = py;
|
|
vox.pz = pz;
|
|
}
|
|
|
|
(void)solved;
|
|
|
|
#if TRACE
|
|
if (options & meshopt_RemeshSolve)
|
|
printf("remesher: solved %zu/%zu voxels\n", solved, voxel_count);
|
|
#endif
|
|
}
|
|
|
|
static void emitVertex(float* result, int x, int y, int z, int corner, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution, const float offset[3])
|
|
{
|
|
int ox = corner & 1, oy = (corner >> 1) & 1, oz = (corner >> 2) & 1;
|
|
|
|
size_t row = (y + oy) + size_t(resolution) * (z + oz);
|
|
size_t idx = (x + ox) + size_t(resolution) * row;
|
|
|
|
assert(grid[idx] != 0 && grid[idx] != 0xff);
|
|
const Voxel& vox = voxels[voxel_rows[row] + (grid[idx] - 1)];
|
|
|
|
result[0] = vox.px + offset[0];
|
|
result[1] = vox.py + offset[1];
|
|
result[2] = vox.pz + offset[2];
|
|
}
|
|
|
|
static bool octantDecider(int x, int y, int z, int cube, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution)
|
|
{
|
|
for (int c = 0; c < 8; ++c)
|
|
if (cube & (1 << c))
|
|
{
|
|
int ox = c & 1, oy = (c >> 1) & 1, oz = (c >> 2) & 1;
|
|
|
|
size_t row = (y + oy) + size_t(resolution) * (z + oz);
|
|
size_t idx = (x + ox) + size_t(resolution) * row;
|
|
|
|
assert(grid[idx] != 0);
|
|
|
|
// note: unlike vertex emission, which never consults interior occupied corners, here we need to skip corners without an occupied voxel
|
|
if (grid[idx] == 0xff)
|
|
continue;
|
|
|
|
const Voxel& vox = voxels[voxel_rows[row] + (grid[idx] - 1)];
|
|
|
|
// test octant contained within the cell (opposite of the corner index)
|
|
if (vox.octants & (1 << (7 - c)))
|
|
return false;
|
|
}
|
|
|
|
// select alternate configuration if *all* octants are empty
|
|
return true;
|
|
}
|
|
|
|
static bool quadricDecider(int x, int y, int z, int cube, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution)
|
|
{
|
|
// quads are encoded implicitly as 0xabc 0xcbd
|
|
unsigned int quad = (kTriangleTable[cube][0][0] << 4) | (kTriangleTable[cube][0][1] & 0xf);
|
|
|
|
const Voxel* corner[4];
|
|
|
|
for (int i = 0; i < 4; ++i)
|
|
{
|
|
int c = (quad >> (12 - i * 4)) & 0xf;
|
|
size_t row = (y + ((c >> 1) & 1)) + size_t(resolution) * (z + ((c >> 2) & 1));
|
|
size_t idx = (x + (c & 1)) + size_t(resolution) * row;
|
|
|
|
assert(grid[idx] != 0 && grid[idx] != 0xff);
|
|
corner[i] = &voxels[voxel_rows[row] + (grid[idx] - 1)];
|
|
}
|
|
|
|
// evaluate error for midpoints of primary (bc) and alternate (ad) diagonal
|
|
float mx = (corner[1]->px + corner[2]->px) * 0.5f, my = (corner[1]->py + corner[2]->py) * 0.5f, mz = (corner[1]->pz + corner[2]->pz) * 0.5f;
|
|
float nx = (corner[0]->px + corner[3]->px) * 0.5f, ny = (corner[0]->py + corner[3]->py) * 0.5f, nz = (corner[0]->pz + corner[3]->pz) * 0.5f;
|
|
float error0 = 0, error1 = 0;
|
|
|
|
for (int i = 0; i < 4; ++i)
|
|
{
|
|
error0 += voxelError(*corner[i], mx, my, mz);
|
|
error1 += voxelError(*corner[i], nx, ny, nz);
|
|
}
|
|
|
|
// select alternate configuration if it is clearly better
|
|
return error1 < error0 * 0.9f;
|
|
}
|
|
|
|
static size_t polygonize(float* destination, size_t max_triangle_count, const unsigned char* grid, const Voxel* voxels, const unsigned int* voxel_rows, int resolution, const float offset[3], unsigned int options)
|
|
{
|
|
size_t result = 0;
|
|
size_t slice = size_t(resolution) * size_t(resolution);
|
|
|
|
assert(gRemeshTablesInitialized);
|
|
|
|
for (int z = 0; z < resolution - 1; ++z)
|
|
for (int y = 0; y < resolution - 1; ++y)
|
|
{
|
|
const unsigned char* data = grid + size_t(resolution) * (y + size_t(resolution) * z);
|
|
|
|
// we track each slice as a 8-bit code (matching cube indexing) as we iterate through the row to avoid extra lookups
|
|
int last = (data[0] != 0) | ((data[resolution] != 0) << 2) | ((data[slice] != 0) << 4) | ((data[slice + resolution] != 0) << 6);
|
|
|
|
for (int x = 0; x < resolution - 1; ++x)
|
|
{
|
|
int next = (data[x + 1] != 0) | ((data[x + 1 + resolution] != 0) << 2) | ((data[x + 1 + slice] != 0) << 4) | ((data[x + 1 + slice + resolution] != 0) << 6);
|
|
|
|
// next slice bits are in position for the next iteration so we need to shift them to retrieve current cube code
|
|
int cube = last | (next << 1);
|
|
last = next;
|
|
|
|
if (cube == 0 || cube == 0xff)
|
|
continue;
|
|
|
|
if (!destination)
|
|
{
|
|
// fast path: we can statically determine the upper bound on triangles the loop below will output based on cube configuration
|
|
result += kTriangleCount[cube];
|
|
continue;
|
|
}
|
|
|
|
// deciders are only evaluated for cells with an alternate configuration to reduce overhead
|
|
int alt = 0;
|
|
|
|
if (kTriangleAlt[cube] == 1)
|
|
alt = octantDecider(x, y, z, cube, grid, voxels, voxel_rows, resolution);
|
|
else if (kTriangleAlt[cube] == 2 && (options & meshopt_RemeshSolve))
|
|
alt = quadricDecider(x, y, z, cube, grid, voxels, voxel_rows, resolution);
|
|
|
|
const unsigned short* tris = kTriangleTable[cube][alt];
|
|
|
|
for (int i = 0; tris[i]; ++i)
|
|
{
|
|
// note: we only emit the triangle if we have space for it, but we count it regardless for consistent capacity estimation
|
|
if (result < max_triangle_count)
|
|
{
|
|
unsigned short tri = tris[i];
|
|
|
|
emitVertex(&destination[result * 9 + 0], x, y, z, (tri >> 8) & 0xf, grid, voxels, voxel_rows, resolution, offset);
|
|
emitVertex(&destination[result * 9 + 3], x, y, z, (tri >> 4) & 0xf, grid, voxels, voxel_rows, resolution, offset);
|
|
emitVertex(&destination[result * 9 + 6], x, y, z, (tri >> 0) & 0xf, grid, voxels, voxel_rows, resolution, offset);
|
|
}
|
|
|
|
result++;
|
|
}
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
} // namespace meshopt
|
|
|
|
size_t meshopt_remesh(float* destination, size_t max_triangle_count, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride, int resolution, unsigned int options)
|
|
{
|
|
using namespace meshopt;
|
|
|
|
assert(destination || max_triangle_count == 0);
|
|
assert(index_count % 3 == 0);
|
|
assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
|
|
assert(vertex_positions_stride % sizeof(float) == 0);
|
|
assert(resolution >= 4 && resolution <= 256);
|
|
|
|
meshopt_Allocator allocator;
|
|
|
|
// measure voxel grid to compute position => voxel mapping
|
|
float offset[3] = {};
|
|
float scale = measureGrid(vertex_positions, vertex_count, vertex_positions_stride, resolution, offset);
|
|
|
|
// rasterize triangles into the voxel grid: in the first pass, this simply tags occupied voxels
|
|
unsigned char* grid = allocator.allocate<unsigned char>(size_t(resolution) * size_t(resolution) * size_t(resolution));
|
|
memset(grid, 0, size_t(resolution) * size_t(resolution) * size_t(resolution));
|
|
|
|
voxelize(grid, NULL, NULL, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride, resolution, scale, offset, options);
|
|
|
|
// allocate additional voxel data for each occupied voxel; this can be filled in the second pass to compute positions
|
|
// note that we only do this if we need to compute output triangles; counting runs skip it for performance
|
|
Voxel* voxels = NULL;
|
|
unsigned int* voxel_rows = NULL;
|
|
size_t voxel_count = 0;
|
|
|
|
if (destination)
|
|
{
|
|
voxel_rows = allocator.allocate<unsigned int>(size_t(resolution) * size_t(resolution));
|
|
voxel_count = rowpack(grid, voxel_rows, resolution);
|
|
|
|
voxels = allocator.allocate<Voxel>(voxel_count);
|
|
memset(voxels, 0, voxel_count * sizeof(Voxel));
|
|
|
|
#if TRACE
|
|
printf("remesher: %zu voxels occupied\n", voxel_count);
|
|
#endif
|
|
}
|
|
|
|
// fill in empty voxels that are not reachable from the grid boundary; the inside empty voxels are marked with 0xff
|
|
// note that these voxels do *not* have the associated Voxel data as they were never voxelized
|
|
if ((options & meshopt_RemeshShell) == 0)
|
|
{
|
|
unsigned int* worklist = allocator.allocate<unsigned int>(size_t(resolution) * size_t(resolution));
|
|
unsigned char* queued = allocator.allocate<unsigned char>(size_t(resolution) * size_t(resolution));
|
|
|
|
solidify(grid, worklist, queued, resolution);
|
|
|
|
#if TRACE
|
|
size_t inside_count = 0, occupied_count = 0;
|
|
for (size_t i = 0; i < size_t(resolution) * size_t(resolution) * size_t(resolution); ++i)
|
|
{
|
|
inside_count += (grid[i] == 0xff);
|
|
occupied_count += (grid[i] != 0 && grid[i] != 0xff);
|
|
}
|
|
|
|
printf("remesher: %zu voxels occupied, %zu voxels inside\n", occupied_count, inside_count);
|
|
#endif
|
|
}
|
|
|
|
// accumulate voxel positions: in the second pass, this computes enough data in each voxel to calculate positions
|
|
if (voxels)
|
|
voxelize(grid, voxels, voxel_rows, indices, index_count, vertex_positions, vertex_count, vertex_positions_stride, resolution, scale, offset, options);
|
|
|
|
// compute final voxel positions; each voxel has a single resulting position that will be emitted during polygonization
|
|
if (voxels)
|
|
solve(voxels, voxel_count, scale, options);
|
|
|
|
// output triangles from the voxel grid; if destination is NULL, this still counts the number of triangles that would be generated
|
|
size_t result = polygonize(destination, max_triangle_count, grid, voxels, voxel_rows, resolution, offset, options);
|
|
|
|
#if TRACE
|
|
printf("remesher: %zu triangles (%zu capacity)\n", result, max_triangle_count);
|
|
#endif
|
|
|
|
return result;
|
|
}
|