在三维网格上有效地找到等成本点，并且点数成本最低
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了在三维网格上有效地找到等成本点，并且点数成本最低相关的知识，希望对你有一定的参考价值。
我有一个3d网格，其中网格上的每个点（x，y，z）与成本值相关联。任何点（x，y，z）的成本都不是预先知道的。要知道成本，我们需要进行一个非常昂贵的复杂查询。我们对该目标了解的一件事是，在所有3个维度中，成本是单调不减少的。
现在给出成本C，我需要在表面上找到具有成本C的点（x，y，z）。这必须通过仅花费最小成本来完成。如何解决我的问题？
当我在线搜索时，我正在获得与轮廓识别相关的技术，但所有这些技术都假设所有点的成本都是预先知道的，比如Marching cubes方法等。在我的情况下，主要指标是计算的点数应该是最小的。
如果某人能够建议一种获得近似位置的方法，至少如果不准确的话会很有帮助。
重写的解释:(原文，如果它可能向某人澄清这个想法，在线下保持不变）
我们在三维中有一些函数f（x，y，z），我们希望找到表面f（x，y，z）= c。由于函数产生一个数字，它定义了a scalar field，我们正在寻找的表面是isosurface c。
在我们的例子中，评估函数f（x，y，z）非常昂贵，因此我们希望尽量减少使用它的次数。不幸的是，大多数等值面算法都假设相反。
我的建议是使用类似的等值面行走，因为Fractint可以用于二维分形。代码方面，它很复杂，但它应该最小化所需的功能评估量 - 这正是它在Fractint中实现的目的。
背景/历史：
在20世纪80年代末和90年代初期，我遇到了一个分形绘图套件Fractint。计算机速度要慢得多，评估每个点的速度非常慢。在Fractint中做了很多努力，使其尽可能快地显示分形，但仍然准确。 （你们中的一些人可能会记住它可以做的颜色循环，通过旋转所用调色板中的颜色。它是催眠的; here是1995年纪录片“无限的颜色”的Youtube剪辑，它的颜色循环和放大计算全屏分形可能需要数小时（在高变焦因子下，接近实际的分形集），但是你可以（保存为图像）并使用颜色循环来“动画”它。）
这些分形中的一些是或者有区域，其中所需的迭代次数朝着实际的分形集分形单调非减少 - 也就是说，没有“岛”伸出，只是在迭代步骤中偶尔稳定增加 - ，快速评估模式使用边缘跟踪来定位迭代次数改变的边界：换句话说，区域填充单一颜色。关闭一个区域后，它会跟踪该区域的中心以找到下一个迭代边缘;在它关闭之后，它可以用正确的恒定颜色填充这些边界之间的环形或C形区域，而不评估这些像素的功能！
在这里，我们有一个非常相似的情况，除了三个维度而不是两个维度。根据定义，每个等值面也是二维的，所以真正改变的是我们如何走边界。
步行本身类似于flood fill算法，除了我们走三维，我们的边界是我们追踪的等值面。
我们在regular grid中采样原始函数，比如N×N×N网格。 （这不是唯一的可能性，但它是最容易和最有用的情况，也是OP正在做的事情。）
通常，等值面不会精确地穿过网格点，而是在网格点之间。因此，我们的任务是找到等值面通过的网格单元。
在N×N×N规则网格中，存在（N-1）×（N-1）×（N-1）个立方单元：
每个单元格在（x，y，z），（x + 1，y，z），（x，y + 1，z），（x + 1，y + 1，z），（x，y）处有八个角，z + 1），（x + 1，y，z + 1），（x，y + 1，z + 1）和（x + 1，y + 1，z + 1），其中x，y， Z∈ℕ是整数网格坐标，0≤x，y，z≤N-2是整数网格坐标。
仔细注意整数网格坐标限制。如果你考虑一下，你会发现N×N×N网格只有（N-1）×（N-1）×（N-1）个单元格，因为我们使用最接近角落的网格坐标到原点，该角的有效坐标范围是0到N-2，包括0和N-2。
如果f（x，y，z）在每个维度上单调增加，则等值面c通过单元格（x，y，z），如果
f(x,y,z) ≤ c
至少有一个
f(x+1, y, z) > c
f(x, y+1, z) > c
f(x+1, y+1, z) > c
f(x, y, z+1) > c
f(x+1, y, z+1) > c
f(x, y+1, z+1) > c
f(x+1, y+1, z+1) > c
如果f（x，y，z）是单调非递减的 - 也就是说，它的偏导数在所有点都是零或正 - 然后上面定位二维等值面，而外层表面是等体积（体积）其中f（x，y，z）是常数）。等体积c的内表面是那些细胞（x，y，z）
f(x,y,z) < c
至少有一个
f(x+1, y, z) ≥ c
f(x, y+1, z) ≥ c
f(x+1, y+1, z) ≥ c
f(x, y, z+1) ≥ c
f(x+1, y, z+1) ≥ c
f(x, y+1, z+1) ≥ c
f(x+1, y+1, z+1) ≥ c
扩展到任何标量函数：
这里显示的方法实际上适用于在采样区域内只有一个最大值的任何f（x，y，z），例如（xMAX，yMAX，zMAX）;并且只有一个最小值，例如（xMIN，yMIN，zMIN）;在采样区域内没有局部最大值或最小值。
在这种情况下，规则是f（x，y，z），f（x + 1，y，z），f（x，y + 1，z），f（x + 1，y）中的至少一个+ 1，z），f（x，y，z），f（x + 1，y，z），f（x，y + 1，z），f（x + 1，y + 1，z）必须低于或等于c，并且至少一个高于或等于c，并且不等于c。
此外，可以始终使用（xMAX，yMAX，zMAX）和（xMIN，yMIN，zMIN）之间的二进制搜索找到初始单元和等值面c通过，将坐标限制为0≤xMAX，yMAX，zMAX，xMIN， yMIN，zMIN≤N-2（换句话说，仅考虑有效单元）。
如果函数不是单调的，则定位等值面c通过的初始单元更复杂。在这种情况下，您需要一种不同的方法。 （如果您可以找到所有局部最大值和最小值的网格坐标，那么您可以从c以上的全局最小值到局部最大值进行二进制搜索，并从c以下的局部最小值到全局最大值进行二进制搜索。）
因为我们以间隔对函数f（x，y，z）进行采样，所以我们隐含地假设它是连续的。如果不是这样 - 并且您还需要显示不连续性 - 您可以使用每个点处的不连续信息来增加网格（七个布尔标记或每个网格点的位数，用于“从（x，y，z）到”的不连续性（X +，Y +，Z +）“）。然后，表面行走也必须尊重（不交叉）这种不连续性。
在实践中，我将使用两个数组来描述网格：一个用于缓存样本，一个用于每个网格点两个标记。一个标志将描述缓存的值存在，另一个标志表示步行例程已经在该点处走过网格单元。我使用/需要用于行走和构造等值面的结构（对于在规则网格中采样的单调非递减函数）将是
typedef struct {
size_t xsize;
size_t ysize;
size_t zsize;
size_t size; /* xsize * ysize * zsize */
size_t xstride; /* [z][y][x] array = 1 */
size_t ystride; /* [z][y][x] array = xsize */
size_t zstride; /* [z][y][x] array = xsize * ysize */
double xorigin; /* Function x for grid coordinate x = 0 */
double yorigin; /* Function y for grid coordinate y = 0 */
double zorigin; /* Function z for grid coordinate z = 0 */
double xunit; /* Function x for grid coordinate x = 1 */
double yunit; /* Function y for grid coordinate y = 1 */
double zunit; /* Function z for grid coordinate z = 1 */
/* Function to obtain a new sample */
void *data;
double *sample(void *data, double x, double y, double z);
/* Walking stack */
size_t stack_size;
size_t stack_used;
size_t *stack;
unsigned char *cell; /* CELL_ flags */
double *cache; /* Cached samples */
} grid;
#define CELL_UNKNOWN (0U)
#define CELL_SAMPLED (1U)
#define CELL_STACKED (2U)
#define CELL_WALKED (4U)
double grid_sample(const grid *const g, const size_t gx, const size_t gy, const size_t gz)
{
const size_t i = gx * g->xstride + gy * g->ystride + gz * g->zstride;
if (!(g->cell[i] & CELL_SAMPLED)) {
g->cell[i] |= CELL_SAMPLED;
g->cache[i] = g->sample(g->data, g->xorigin + (double)gx * g->xunit,
g->yorigin + (double)gy * g->yunit,
g->zorigin + (double)gz * g->zunit);
}
return g->cache[i];
}
以及找到要开始行走的单元格的功能，使用沿网格对角线的二分搜索（假设非减少单调函数，因此所有等值面必须穿过对角线）：
size_t grid_find(const grid *const g, const double c)
{
const size_t none = g->size;
size_t xmin = 0;
size_t ymin = 0;
size_t zmin = 0;
size_t xmax = g->xsize - 2;
size_t ymax = g->ysize - 2;
size_t zmax = g->zsize - 2;
double s;
s = grid_sample(g, xmin, ymin, zmin);
if (s > c) {
return none;
}
if (s == c)
return xmin*g->xstride + ymin*g->ystride + zmin*g->zstride;
s = grid_sample(g, xmax, ymax, zmax);
if (s < c)
return none;
if (s == c)
return xmax*g->xstride + ymax*g->ystride + zmax*g->zstride;
while (1) {
const size_t x = xmin + (xmax - xmin) / 2;
const size_t y = ymin + (ymax - ymin) / 2;
const size_t z = zmin + (zmax - zmin) / 2;
if (x == xmin && y == ymin && z == zmin)
return x*g->xstride + y*g->ystride + z*g->zstride;
s = grid_sample(g, x, y, z);
if (s < c) {
xmin = x;
ymin = y;
zmin = z;
} else
if (s > c) {
xmax = x;
ymax = y;
zmax = z;
} else
return x*g->xstride + y*g->ystride + z*g->zstride;
}
}
#define GRID_X(grid, index) (((index) / (grid)->xstride)) % (grid)->xsize)
#define GRID_Y(grid, index) (((index) / (grid)->ystride)) % (grid)->ysize)
#define GRID_Z(grid, index) (((index) / (grid)->zstride)) % (grid)->zsize)
上面的三个宏显示了如何将网格索引转换回网格坐标。
走同构面，我们不能依靠递归;调用链太长了。相反，我们应该检查单元索引的walk堆栈：
static void grid_push(grid *const g, const size_t cell_index)
{
/* If the stack is full, remove cells already walked. */
if (g->stack_used >= g->stack_size) {
const size_t n = g->stack_used;
size_t *const s = g->stack;
unsigned char *const c = g->cell;
size_t i = 0;
size_t o = 0;
while (i < n)
if (c[s[i]] & CELL_WALKED)
i++;
else
s[o++] = s[i++];
g->stack_used = o;
}
/* Grow stack if still necessary. */
if (g->stack_used >= g->stack_size) {
size_t *new_stack;
size_t new_size;
if (g->stack_used < 1024)
new_size = 1024;
else
if (g->stack_used < 1048576)
new_size = g->stack_used * 2;
else
new_size = (g->stack_used | 1048575) + 1048448;
new_stack = realloc(g->stack, new_size * sizeof g->stack[0]);
if (new_stack == NULL) {
/* FATAL ERROR, out of memory */
}
g->stack = new_stack;
以上是关于在三维网格上有效地找到等成本点，并且点数成本最低的主要内容，如果未能解决你的问题，请参考以下文章