我有一个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;
g->stack_size = new_size;
}
/* Unnecessary check.. */
if (!(g->cell[cell_index] & (CELL_STACKED | CELL_WALKED)))
g->stack[g->stack_used++] = cell_index;
}
static size_t grid_pop(grid *const g)
{
while (g->stack_used > 0 &&
g->cell[g->stack[g->stack_used - 1]] & CELL_WALKED)
g->stack_used--;
if (g->stack_used > 0)
return g->stack[--g->stack_used];
return g->size; /* "none" */
}
验证isosurface通过当前单元格的函数，将这些函数报告给回调函数，然后遍历isosurface，就像是
int isosurface(grid *const g, const double c,
int (*report)(grid *const g,
const size_t x, const size_t y, const size_t z,
const double c,
const double x0y0z0,
const double x1y0z0,
const double x0y1z0,
const double x1y1z0,
const double x0y0z1,
const double x1y0z1,
const double x0y1z1,
const double x1y1z1))
{
const size_t xend = g->xsize - 2; /* Since we examine x+1, too */
const size_t yend = g->ysize - 2; /* Since we examine y+1, too */
const size_t zend = g->zsize - 2; /* Since we examine z+1, too */
const size_t xstride = g->xstride;
const size_t ystride = g->ystride;
const size_t zstride = g->zstride;
unsigned char *const cell = g->cell;
double x0y0z0, x1y0z0, x0y1z0, x1y1z0,
x0y0z1, x1y0z1, x0y1z1, x1y1z1; /* Cell corner samples */
size_t x, y, z, i;
int r;
/* Clear walk stack. */
g->stack_used = 0;
/* Clear walked and stacked flags from the grid cell map. */
i = g->size;
while (i-->0)
g->cell[i] &= ~(CELL_WALKED | CELL_STACKED);
i = grid_find(g, c);
if (i >= g->size)
return errno = ENOENT; /* No isosurface c */
x = (i / g->xstride) % g->xsize;
y = (i / g->ystride) % g->ysize;
z = (i / g->zstride) % g->zsize;
/* We need to limit x,y,z to the valid *cell* coordinates. */
if (x > xend) x = xend;
if (y > yend) y = yend;
if (z > zend) z = zend;
i = x*g->xstride + y*g->ystride + z*g->zstride;
if (x > xend || y > yend || z > zend)
return errno = ENOENT; /* grid_find() returned an edge cell */
grid_push(g, i);
while ((i = grid_pop) < g->size) {
x = (i / g->xstride) % g->xsize;
y = (i / g->ystride) % g->ysize;
z = (i / g->zstride) % g->zsize;
cell[i] |= CELL_WALKED;
x0y0z0 = grid_sample(g, x, y, z);
if (x0y0z0 > c)
continue;
x1y0z0 = grid_sample(g, 1+x, y, z);
x0y1z0 = grid_sample(g, x, 1+y, z);
x1y1z0 = grid_sample(g, 1+x, 1+y, z);
x0y0z1 = grid_sample(g, x, y, 1+z);
x1y0z1 = grid_sample(g, 1+x, y, 1+z);
x0y1z1 = grid_sample(g, x, 1+y, 1+z);
x1y1z1 = grid_sample(g, 1+x, 1+y, 1+z);
/* Isosurface does not pass through this cell?!
* (Note: I think this check is unnecessary.) */
if (x1y0z0 < c && x0y1z0 < c && x1y1z0 < c &&
x0y0z1 < c && x1y0z1 < c && x0y1z1 < c &&
x1y1z1 < c)
continue;
/* Report the cell. */
if (report) {
r = report(g, x, y, z, c, x0y0z0, x1y0z0,
x0y1z0, x1y1z0, x0y0z1, x1y0z1,
x0y1z1, x1y1z1);
if (r) {
errno = 0;
return r;
}
}
/* Could the surface extend to -x? */
if (x > 0 &&
!(cell[i - xstride] & (CELL_WALKED | CELL_STACKED)) &&
( x0y1z0 >= c || x0y0z1 >= c ))
grid_push(g, i - xstride);
/* Could the surface extend to -y? */
if (y > 0 &&
!(cell[i - ystride] & (CELL_WALKED | CELL_STACKED)) &&
( x0y0z1 >= c || x1y0z0 >= c ))
grid_push(g, i - ystride);
/* Could the surface extend to -z? */
if (z > 0 &&
!(cell[i - zstride] & (CELL_WALKED | CELL_STACKED)) &&
( x1y0z0 >= c || x0y1z0 >= c ))
grid_push(g, i - zstride);
/* Could the surface extend to +x? */
if (x < xend &&
!(cell[i + xstride] & (CELL_WALKED | CELL_STACKED)) &&
( x0y1z0 >= c || x0y0z1 >= c ))
grid_push(g, i + xstride);
/* Could the surface extend to +y? */
if (y < xend &&
!(cell[i + ystride] & (CELL_WALKED | CELL_STACKED)) &&
( x1y0z0 >= c || x0y0z1 >= c ))
grid_push(g, i + ystride);
/* Could the surface extend to +z? */
if (z < xend &&
!(cell[i + zstride] & (CELL_WALKED | CELL_STACKED)) &&
( x1y0z0 >= c || x0y1z0 >= c ))
grid_push(g, i + zstride);
}
/* All done. */
errno = 0;
return 0;
}
在这种特殊情况下，我确实相信使用多边形网格可以最佳地可视化/描述等值面，其中单元内的样本是线性插值的。然后，每个report（）调用生成一个多边形（或一个或多个平面三角形）。
请注意，单元格有12个边，并且等值面必须至少穿过其中的三个。假设我们在拐角c0和c1处有两个样本，由边缘跨越，两个角分别具有坐标p0 =（x0，y0，z0）和p1 =（x1，y1，z1）：
if (c0 == c && c1 == c)
/* Entire edge is on the isosurface */
else
if (c0 == c)
/* Isosurface intersects edge at p0 */
else
if (c1 == c)
/* Isosurface intersects edge at p1 */
else
if (c0 < c && c1 > c)
/* Isosurface intersects edge at p0 + (p1-p0)*(c-c0)/(c1-c0) */
else
if (c0 > c && c1 < c)
/* Isosurface intersects edge at p1 + (p0-p1)*(c-c1)/(c0-c1) */
else
/* Isosurface does not intersect the edge */
上述检查对任何类型的连续函数f（x，y，z）都有效;对于非单调函数，问题只是找到相关的单元格。根据本文前面概述的规则，isosurface（）函数需要进行一些更改（检查wrt.x0y0z0..x1y1z1），但它也可以用于任何连续函数f（x，y，z）很少的努力。
如果已知单元角处的样本，尤其是使用线性插值，则构造多边形/三角形非常简单，如您所见。
请注意，通常没有理由担心检查单元格边缘的顺序，因为您几乎肯定会使用矢量微积分和cross product来定位点和多边形。或者，如果你愿意，你可以在点上做Delaunay triangulation（对于任何函数，3到12，尽管超过6个点表示有两个单独的表面，我相信）来构造平面多边形。
有问题吗？评论？
我们在三维中有一个标量场f（x，y，z）。该场的采样/评估成本很高，我们只在整数坐标0≤x，y，z∈ℕ这样做。为了可视化标量场，我们希望使用最小数量的样本/评估来定位一个或多个等值面（具有特定f（x，y，z）值的表面）。
我将在这里尝试描述的方法是fractint中使用的算法的变体，以最小化绘制某些分形所需的迭代次数。一些分形具有相同“值”的大区域，因此不是对区域内的每个点进行采样，而是某些绘制模式跟踪这些区域的边缘。
换句话说，不是定位等值面c的各个点，f（x，y，z）= c，你可以只找到一个点，然后走同构面。步行部分可视化有点复杂，但它实际上只是简单计算机图形中使用的flood fill算法的3D变体。 （实际上，假设场沿着每个维度单调不递减，它实际上主要是2D步行，通常只有几个网格点，而不是与等值面c相关的网格点。这应该是非常有效的。）
我很确定有很好的同行评审论文描述了这种技术（可能在不止一个问题领域），但由于我懒得做一个比谷歌搜索几分钟更好的搜索，我就把它留下来给别人找好参考。道歉。
为简单起见，现在让我们假设该场是连续的并且沿着每个维度单调增加。在面向轴的方框中，大小为N×N×N，该字段在原点（0,0,0）的一个角处具有最小值，在原点的远角处具有最小值，在（N，N，N）处，从（0,0,0）到（N，N，N）的对角线上找到的最小值和最大值之间的所有可能值。换句话说，每个可能的等值面都存在并且是连续的2D表面，不包括点（0,0,0）和（N，N，N），并且每个这样的表面与对角线相交。
如果该字段实际上是非连续的，我们将无法分辨，因为我们的采样方法。在实践中，我们的抽样意味着我们隐含地假设标量场是连续的;无论是否真的如此，我们都将持续对待！
如果函数实际上沿着每个维度单调递增，那么可以将f（x，y，z）= c映射到X（y，z）= x，Y（x，z）= y，Z（x， y）= z，尽管三者中的任何一个足以定义等值面c。这是因为等值面只能跨越最多一点的任何跨越框的线。
如果函数是单调非递减的，则等值面可以与跨越框的任何线相交仍然只有一次，但是交点可以沿着线宽（比一个点）。在实践中，您可以通过仅考虑isovolumes的下表面或上表面（具有静态字段的体积）来处理此问题;也就是说，只有从低于c到c或更大的过渡，或从-c-或 - 低于-c-大于c的过渡。在所有情况下，您并不是真的在寻找等值面值c，而是试图找到一对场样本穿过c的位置。
因为我们在常规网格点处对场进行采样，并且等值面很少（如果有）与这些网格点完全相交，我们将原始框划分为N×N×N个单位大小的立方体，并尝试找到所需的等值面相交的立方体。
这是一个这样的立方体的简单说明，在（x，y，z）到（x + 1，y + 1，z + 1）：
当等值面与立方体相交时，它与标记为X，Y或Z的边和/或标记为D的对角线中的至少一个边相交。特别地，我们将具有f（x，y，z）≤c，并且一个或多个：
我们可以找到一个这样的立方体，而不是对原始体积进行完整的搜索，然后沿着立方体移动以发现等值面相交的立方体。
由于所有等值面必须与（0,0,0）到（N，N，N）的对角线相交，我们可以使用2 + ceil（log2（N））样本找到这样的立方体，使用二进制搜索对角线上的立方体。目标立方体（i，i，i）是f（i，i，i）≤c和f（i + 1，i + 1，i + 1）> c的立方体。 （对于具有等体积的单调非减少场，这表明等体积表面更接近原点作为等值面。）
当我们知道isosurface c与一个立方体相交时，我们基本上可以使用三种方法将这些知识转换为一个点（我们认为isosurface相交）：
由此获得的等值面c的交点可以用于拟合一些表面函数，但我在现实生活中没有看到它。通常，Delaunay triangulation用于将点集转换为多边形网格，然后很容易可视化。
另一种选择是记住每个点与哪个立方体（（x，y，z））和边缘/对角线（X，Y或Z边缘，或对角线D）。然后，您可以直接形成多边形网格。 Voxel技术也可用于快速绘制部分透明的等值面;每个视图光线检查每个立方体一次，如果存在等值面，则可以使用等值面交叉点来插值表面法线向量，使用光线投射/光线跟踪方法生成非常平滑且准确的等值面（不创建任何多边形网格）。
在我看来，这个答案需要编辑 - 至少，一些睡眠和进一步思考，以及澄清。欢迎提出问题，建议，甚至编辑！
如果不仅仅是OP的兴趣，我可以试着看看我是否可以拼凑一个简单的示例C程序。我玩弄了可视化的模拟电子结构，这些领域甚至都不是单调的（虽然采样很便宜）。
您应该查看本文讨论二维案例，并让您深入了解不同的方法：http://leetcode.com/2010/10/searching-2d-sorted-matrix.html 在我看来，逐步线性搜索（在那里的第二部分）对你来说是一个很好的第一步，因为它很容易应用于三维情况，它真的不需要很多经验来理解。 因为这非常简单且仍然非常高效，我会继续使用它，看看它是否符合您在3-d中使用的数据类型的需求。 但是，如果您的唯一目标是性能，则应将二进制分区应用于3-d。这变得有点复杂，因为他所说的“二进制分区”基本上变成了“二进制平面分区”。 因此，您没有将矩阵划分为2个可能较小的矩阵的行。 相反，你有一个平面将你的立方体划分为2个可能的较小的立方体。 要使该平面（或矩阵）中的搜索有效，您首先必须实现他的一种方法:)。 然后用较小的立方体重复所有内容。 请记住，以非常有效的方式实现这一点（即记住内存访问）并非易事。
我将给出这个答案，以尽量减少计算的成本数量。 Matt Ko链接到一个很好的解决方案，但它假设一个廉价的成本函数和基于矩阵的数据，你似乎没有。我给出的方法需要更接近O(log N + k)调用成本函数，其中k是具有所需成本的点数。请注意，这种具有一些性能优化的算法可以在三维矩阵上成为O(N)，几乎没有机会调用性能成本函数，尽管它有点复杂。
psudeocode基于quickselect中使用的技术，如下所示：
While there are still points under considerations:
Find the ideal pivot point and calculate it's cost
Remove the pivot from the point set
If the cost is the desired cost then:
Add the pivot to the solution set
Else:
Separate the points into 3 groups:
G1. Those that are in in the pivot's octant `VII`
G2. Those have the same x, y, or z of the pivot
G3. Those that are not in the pivot's octant `VII`
# Note this can be done in O(N)
If all points are in group 2:
Use 1D binary searches in each dimension to find points with the desired cost
Else:
Compute the cost of the pivot
Keep all points in group 2
If the pivot cost is greater than desired:
Keep only the points in group 1
Else:
Keep only the points in group 3
根据该线的octant VII内外点选择枢轴。如果需要，稍后处理形成八分圆的3条线中的任何一条上的点（G2）。
理想的枢轴点是组1（G1）和组3（G3）中的点数尽可能接近相等。从数学角度看它将沿着最大化两者中较大者的两条线，或者maximize(max(|G1|,|G3|) / min(|G1|,|G3|) )。即使是一个寻找理想枢轴点的相当天真的算法也可以在O(N^2)（可能存在O(N log N)算法）中找到它，但是需要O(N^3)来计算理想枢轴的成本。
在找到理想的枢轴并计算成本之后，每次迭代应该平均看到丢弃的剩余点的大约一半，这再次导致仅对成本函数进行O(log N + k)调用。
最后的说明：
回想起来，我不确定第2组的特殊考虑是否真的需要，因为它可能在第3组，但我不是100％肯定。但是，将它分离似乎并没有改变Big O，所以我没有看到需要改变它，尽管这样做会略微简化算法。
这本身不是答案，只是略微概括的示例C代码。 （代码太长，不能逐字包含。）
基本实现在grid.h (pastebin link)。
在其中，我试图区分网格坐标（0≤x，y，z≤size-1）和单元坐标（0≤x，y，z≤size-2）。特别要注意span类型。每个单元格跨越一系列值：插值，或单元格八角处的离散样本集。因为这个例子使用线性插值来确定每个单元格中等值面与边缘或对角线相交的位置，所以我假设连续跨距。
在实现此示例代码之前，我没有意识到跨越值的单元格对于边缘情况有多么重要。这就是为什么OP和我在评论中讨论了我的另一个答案中的边缘情况，以及为什么单独在我的另一个答案中概述的逻辑不能正确处理边缘情况。
由于OP的特殊情况不是那么常见/有趣，因此这个例子更加通用（因此对于OP的情况来说是非常不优化的）。实际上，这个例子只要求函数没有局部最小值或最大值（允许鞍点和常数区域）;在网格区域内只有一个最小值和一个最大值。最小值和最大值不需要是点状的;他们可以是连续的地区。
因此，在网格生成时，我们不知道哪些单元格包含最小值和最大值。 （在OP的情况下，标量场是单调非递减的并且限于正八分圆，因此最小值为0,0,0，最大值为1，size-1，size-1。）
为了找到最小值和最大值，我实现了两个函数，从网格中的最佳角开始（具有最小或最大的样本值）。 grid_maximum_cell()走非衰减细胞，grid_minimum_cell()走非增加细胞。由于标量场被采样，我们隐含地假设它是连续的。只要没有步行可能停止的局部最大值或最小值，步行将在相对较少的样本中到达正确的单元格。 （但是这种搜索可以进一步优化。考虑这两个函数只是你自己实现的起点。当然，OP根本不需要这些。）
（实际上，对采样标量场的要求是每个等值面是连续的，并且所有等值面都与从使用上述两个函数找到的最小和最大单元格绘制的线相交。）
函数grid_isosurface()可用于定位所需的等值面（字段值）通过的单元格。最后一个参数是一个函数指针。对于等值面通过的每个单元，调用该函数一次。 （注意角落样本的索引顺序，[x][y][z]。）
grid_isosurface()使用二进制搜索（在包含最小样本的单元格的行上，到包含最大样本的单元格）找到所需等值面通过的初始单元格。然后使用我的答案中概述的泛洪填充算法跟踪表面。
例如，grid.c (pastebin link)使用上面的包含文件来评估标量字段
f（x，y，z）= x3 + y3 + z3 + x + y-0.125·（x·y + x·z + y·z + x·y·z）。
在我的Linux机器上，我使用编译并运行了示例
gcc -Wall -std=c99 -Wno-unused -O2 grid.c -o isosurface
./isosurface 50 -1.0 1.0 0.0 > out-0.0
./isosurface 50 -1.0 1.0 0.5 > out-0.5
./isosurface 50 -1.0 1.0 1.0 > out-1.0
并使用Gnuplot绘制出三个等值面：
splot "out-0.0" u 1:2:3 notitle w dots, "out-0.5" u 1:2:3 notitle w dots, "out-1.0" u notitle w dots
这导致了这个相当不错的点云（可在Gnuplot中旋转）：
当最初生成网格时，采集14个样本来定位最大和最小单元。追踪等值面分别需要额外的18024,18199和16953个样本;请注意，如果在同一网格上连续执行第二个和更多等值面，则需要的样本要少得多。
上面的总网格包含51×51×51 = 132651个样本，因此跟踪一个等值面需要约13％的网格点进行采样。对于101×101×101网格，需要的样本下降到约7％; 201×201×201电网，降至3.5％;对于501x501x501网格，为1.4％（125.75M样本中的1.7M）。
这些代码都没有针对OP的情况进行优化，也没有针对OP进行优化。样本缓存用于最小化一般所需的样本数量，但grid_isosurface()等值面行走函数，并且可以修改初始grid_minimum_cell()和grid_maximum_cell()函数以需要略少的样本。对于较大的网格，我不认为优化会产生很大的差异，但对于非常小的网格和非常慢的函数来评估，它可能是值得的。
如果目的是为每个等值面生成多边形网格，我建议在回调函数中生成每个多边形，而不是从整个生成的点云生成。使用上面示例程序中的边/对角交点，可以得到跨越该单元的多边形的所有顶点（不需要高速缓存等）。您所需要的只是正确排序边缘交叉点。
有问题吗？评论？ Bug修复？建议？