在三维网格上有效地找到等成本点，并且点数成本最低
我有一个3d网格 ，其中网格上的每个点（x，y，z） 与 成本值相关联 。 任何点（x，y，z）的成本都不是事先知道的 。 要知道成本，我们需要进行一个非常昂贵的复杂查询。 我们对这个目标知道的一件事是， 所有三个维度的成本都是单调不减少的 。
现在给出成本C，我需要在表面上找到成本为C 的点（x，y，z） 。 这必须通过仅花费最低成本来完成。 如何解决我的问题？
当我在网上搜索时，我得到了与轮廓识别相关的技术，但是所有这些技术都假设所有点的成本都是预先知道的，比如Marching cubes方法等。在我的例子中，主要指标是成本计算的点数应该是最小的。
如果有人能够建议一种获得近似位置的方法，至少如果不准确的话会很有帮助。
重写的解释:(原文，如果它可能向某人澄清这个想法，保持在线下不变）
我们在三维中有一些函数f（x，y，z） ，我们希望找到表面f（x，y，z） = c 。 由于函数产生一个数字，它定义了一个标量场 ，我们要寻找的表面是等值面 c 。
在我们的例子中，评估函数f（x，y，z）非常昂贵，因此我们希望尽量减少使用它的次数。 不幸的是，大多数等值面算法都假设相反。
我的建议是使用类似的等值面走路，因为Fractint可以用于二维分形。 代码方面，它很复杂，但它应该最小化所需的function评估量 – 这正是它在Fractint中实现的目的。
背景/历史：
在20世纪80年代末和90年代初期，我遇到了分形绘图套件Fractint 。 计算机速度慢得多，评估每个点的速度非常慢。 在Fractint中做了很多努力，使其尽可能快地显示分形，但仍然准确。 （你们中的一些人可能会记得它可以做的颜色循环，通过旋转所用调色板中的颜色。它是催眠的; 这是一个来自1995年纪录片“无限的颜色”的Youtube剪辑，它的颜色循环和放大计算全屏分形可能需要数小时（在高变焦因子下，接近实际的分形集），但是你可以（保存为图像）并使用颜色循环来“动画”它。）
这些分形中的一些是或者有区域，其中所需的迭代次数朝着实际的分形集分形单调非减少 – 即，没有“岛”伸出，只是迭代步骤中的稳定偶尔增加 – ，一个快速评估模式使用边缘跟踪来定位迭代次数改变的边界：换句话说，区域填充单一颜色。 关闭一个区域后，它会跟踪该区域的中心以找到下一个迭代边缘; 在关闭之后，它可以用正确的恒定颜色填充这些边界之间的环形或C形区域，而不评估这些像素的function！
在这里，我们有一个非常相似的情况，除了三个维度而不是两个维度。 根据定义，每个等值面也是二维的，所以真正发生的一切，就是我们如何走边界。
漫步本身类似于泛洪填充算法，除了我们走三维，我们的边界是我们正在追踪的等值面。
我们在常规网格中对原始函数进行采样，比如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） ，例如（ x MAX ， y MAX ， z MAX ）; 并且只有一个最小值，比如（ x MIN ， y MIN ， z MIN ）; 在采样区域内没有局部最大值或最小值。
在这种情况下，规则是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 。
此外，通过（ x MAX ， y MAX ， z MAX ）和（ x MIN ， y MIN ， z MIN ）之间的二进制搜索，可以始终找到通过等值面c通过的初始单元，将坐标限制为0≤xMAX ， y MAX ， z MAX ， x MIN ， y MIN ， 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]; }
以及找到要开始行走的单元格的function，使用沿网格对角线的二分搜索（假设非减少单调函数，因此所有等值面必须穿过对角线）：
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" */ }
validationisosurface通过当前单元格的函数，将这些函数报告给回调函数，然后遍历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）很少的努力。
当单元角落处的样本已知时，特别是使用线性插值，构造多边形/三角形非常简单，如您所见。
请注意，通常没有理由担心检查单元格边缘的顺序，因为您几乎肯定会使用向量计算和交叉乘积来定位点和多边形。 或者，如果你愿意，你可以对点进行Delaunay三角测量 （对于任何函数，3到12，尽管超过6个点表示有两个单独的表面，我相信）构造平面多边形。
有问题吗？ 评论？
我们在三维中有一个标量场f（x，y，z） 。 该场的采样/评估成本很高，我们只在整数坐标0≤x，y， z∈ℕ这样做。 为了可视化标量场，我们希望使用最小数量的样本/评估来定位一个或多个等值面（具有特定f（x，y，z）值的表面）。
我将在这里尝试描述的方法是fractint中使用的算法的变体，以最小化绘制某些分形所需的迭代次数。 一些分形具有相同“值”的大区域，因此不是对区域内的每个点进行采样，而是某些绘制模式跟踪这些区域的边缘 。
换句话说，不是定位等值面c的各个点， f（x，y，z）= c ，你只能定位一个点，然后走同构面 。 步行部分可视化有点复杂，但它实际上只是简单计算机图形中使用的泛洪填充算法的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 -to- c或更大的过渡，或从-c-或-low-更大于c的过渡。 在所有情况下，您并不是真的在寻找等值面值c ，而是试图找到一对场样本穿过 c的位置 。
因为我们在常规网格点处对场进行采样，并且等值面很少（如果有）与这些网格点完全相交，我们将原始框划分为N×N×N个单位大小的立方体，并尝试找到所需的等值面相交的立方体。
这是一个这样的立方体的简单图示，在（x，y，z）到（x + 1，y + 1，z + 1） ： 
当等值面与立方体相交时，它与标记为X ， Y或Z的边和/或标记为D的对角线中的至少一个相交。 特别是，我们将有f（x，y，z） ≤c ，以及以下一个或多个：
- f（x + 1，y，z）> c （等值面c穿过用X标记的立方体边缘）（注意：在这种情况下，我们希望沿y和z尺寸行走）
- f（x，y + 1，z）> c （等值面c穿过用Y标记的立方体边缘）（注意：在这种情况下，我们希望沿x和z尺寸行走）
- f（x，y，z + 1）> c （等值面c穿过用Z标记的立方体边缘）（注意：在这种情况下，我们希望沿x和y尺寸行走）
- f（x + 1，y + 1，z + 1）> c （等值面c穿过立方体对角线，用D标记）（注意：在这种情况下，我们可能需要检查所有直接连接的网格点，看看哪个方向我们需要走路。）
我们可以找到一个这样的立方体，而不是对原始体积进行完整的搜索，然后沿着立方体走动以发现等值面相交的立方体。
由于所有等值面必须从（0,0,0）到（N，N，N）的对角线相交，我们可以使用2 + ceil（log 2 （N））样本找到这样的立方体，使用二进制搜索对角线上的立方体。 目标立方体（i，i，i）是f（i，i，i） ≤c和f（i + 1，i + 1，i + 1） > c的立方体。 （对于具有等体积的单调非减少场，这表明等体积表面更靠近原点作为等值面。）
当我们知道isosurface c与一个立方体相交时，我们基本上可以使用三种方法将这些知识转换为一个点（我们认为isosurface相交）：
- 立方体有八个角，每个角都在一个网格点。 我们可以选择最接近c的字段值的角点/网格点。
- 我们可以插值 - 选择一个近似点 - 等值面c与边/对角线相交。 我们可以在没有任何额外样本的情况下进行线性插值，因为我们已经知道交叉边缘/对角线末端的样本。 如果u = f（x，y，z）
- 您可以沿边/对角线使用二分搜索来查找交点。 这需要每边/对角线n个额外样本，以获得沿边/对角线的n位精度的交点。 （由于原始网格与字段中的细节相比不能太粗，或者无论如何都不会显示细节，通常使用n = 2 ， n = 3 ， n = 4或n = 5的内容。 ）
由此获得的等值面c的交点可以用于拟合一些表面函数，但是我在现实生活中没有看到它。 通常， Delaunay三角剖分用于将点集转换为多边形网格，然后易于可视化。
另一个选择是记住每个点与哪个立方体（ （x，y，z） ）和边缘/对角线（ X ， Y或Z边缘，或对角线D ）。 然后，您可以直接形成多边形网格。 体素技术也可用于快速绘制部分透明的等值面; 每个视图光线检查每个立方体一次，如果存在等值面，则可以使用等值面交叉点来插值表面法线向量，使用光线投射/光线跟踪方法生成非常平滑和准确的等值面（不创建任何多边形网格）。
在我看来，这个答案需要编辑 - 至少，一些睡眠和进一步思考，以及澄清。 欢迎提出问题，建议，甚至编辑！
如果不仅仅是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是具有期望成本的点的数量。 请注意，这种具有一些性能优化的算法可以在3D矩阵上做成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^2)找到它，但是需要O(N^3)来计算理想的成本发现后转动。
在找到理想的枢轴并计算成本之后，每次迭代应该平均看到丢弃的剩余点的大约一半，这再次导致仅对成本函数的O(log N + k)调用。
最后的说明：
回想起来，我不确定第2组的特殊考虑是否真的需要，因为它可能在第3组，但我不是100％肯定。 然而，将它分离似乎并没有改变Big O，所以我没有看到需要改变它，尽管这样做会略微简化算法。
这本身不是答案，只是略微概括的示例C代码。 （代码太长，不能逐字包含。）
基本实现在grid.h（pastebin链接）中 。
在其中，我试图区分网格坐标 （ 0≤x ， y ，z≤size-1）和单元坐标 （ 0≤x ， y ，z≤size-2）。 特别要注意span类型。 每个单元格跨越一系列值：插值，或单元格八角处的离散样本集。 因为这个例子使用线性插值来确定每个单元格中等值面与边缘或对角线相交的位置，所以我假设连续跨距。
在实现此示例代码之前，我没有意识到跨越值的单元格对于边缘情况有多么重要。 这就是为什么OP和我在评论中讨论了我的另一个答案中的边缘情况，以及为什么单独在我的另一个答案中概述的逻辑不能正确处理边缘情况。
由于OP的特殊情况不是那么常见/有趣，因此这个例子更通用（因此对于OP来说是非常不优化的）。 实际上，这个例子只要求函数没有局部最小值或最大值（允许鞍点和常数区域）; 在网格区域内只有一个最小值和一个最大值。 最小值和最大值不需要是点状的; 他们可以是连续的地区。
因此，在网格生成时，我们不知道哪些单元格包含最小值和最大值。 （在OP的情况下，标量场是单调非递减的并且仅限于正八分圆，因此最小值为0,0,0，最大值为-1， 大小为 -1， 大小为 -1。）
为了找到最小值和最大值，我实现了两个函数，从网格中的最佳角开始（具有最小或最大的样本值）。 grid_maximum_cell()非递减的单元格， grid_minimum_cell()非递增的单元格。 由于标量场被采样，我们隐含地假设它是连续的。 只要没有步行可能停止的局部最大值或最小值，步行将在相对较少的样本中到达正确的单元格。 （但是这种搜索可以进一步优化。考虑这两个函数只是你自己实现的起点。当然，OP根本不需要这些。）
（实际上，对采样标量场的要求是每个等值面是连续的，并且所有等值面都与从使用上述两个函数找到的最小和最大单元格绘制的线相交。）
函数grid_isosurface()可用于定位所需的等值面（字段值）通过的单元格。 最后一个参数是一个函数指针。 对于等值面通过的每个单元，调用该函数一次。 （注意角样本的索引顺序， [x][y][z] 。）
grid_isosurface()使用二进制搜索（从包含最小样本的单元格到包含最大样本的单元格的行grid_isosurface()定位所需等值面通过的初始单元格。 It then traces the surface, using the flood-fill-like algorithm outlined in my answer.
For an example, grid.c (pastebin link) uses the above include file, to evaluate the scalar field
f ( x , y , z ) = x 3 + y 3 + z 3 + x + y – 0.125·( x · y + x · z + y · z + x · y · z ).
On my Linux machine, I compiled and ran the example using
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
and used Gnuplot to plot out the three isosurfaces:
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
which leads to this pretty nice point cloud (rotatable in Gnuplot): 
When the grid is initially generated, 14 samples are taken to locate the maximum and minimum cells. Tracing the isosurfaces required additional 18024, 18199, and 16953 samples, respectively; note that much fewer samples are needed for the second and further isosurfaces, if you do them consecutively on the same grid.
The total grid above contains 51×51×51 = 132651 samples, so tracing one isosurface required about 13% of the grid points to be sampled. For a 101×101×101 grid, the samples needed drops down to about 7%; for a 201×201×201 grid, down to 3.5%; for a 501x501x501 grid, to 1.4% (1.7M out of 125.75M samples).
None of this code is optimized for OP’s case, nor optimized in general. A sample cache is used to minimize the number of samples needed in general, but the grid_isosurface() isosurface walking function, and the initial grid_minimum_cell() and grid_maximum_cell() functions can be modified to require slightly fewer samples. For larger grids, I don’t expect the optimizations to make much of a difference, but for very small grids and very slow functions to evaluate, it might be worthwhile.
If the intent is to generate a polygon mesh for each isosurface, I recommend generating each polygon in the callback function, not from the overall generated point cloud. Using the edge/diagonal intersections like in the above example program, you get all the vertices for the polygon spanning that cell (no caches or such are needed). All you need is to order the edge intersection points correctly.
有问题吗？ 评论？ Bug fixes? 建议？