以最小的点成本有效地在3d网格上找到等成本点 [英] find iso-cost points on a 3d grid efficiently with minimum costing of points
问题描述
我有一个 3d网格，其中网格上的每个点(x，y，z) 与一个 成本值.任何点(x，y，z)的成本事先未知.要知道成本，我们需要进行一个复杂的查询，这确实很昂贵.我们对该对象了解的一件事是，成本在所有3个维度上均单调不减.
现在给定成本C，我需要在表面上找到成本C的点(x，y，z).这必须通过仅计算最低费用来完成.如何解决我的问题?
当我在线搜索时，我得到了与轮廓识别相关的技术，但是所有这些技术都假设所有点的成本都像Marching cubes方法等那样是预先知道的.在我的情况下，主要度量标准是要计算的点数最少. >
如果有人可以提出一种至少可以获取近似位置的方法，将会很有帮助.
重写的说明: (原始文本(如果可能会向某人说明这一点，请在该行下保持不变)
我们在三维中有一些函数 f(x，y，z)，我们希望找到曲面 f(x，y，z) = c .由于该函数产生单个数字，因此它定义了标量字段，并且我们正在寻找表面为等值表面 c .
在我们的例子中，评估函数 f(x，y，z)的成本很高，因此我们希望最大程度地减少使用它的次数.不幸的是，大多数等值面算法都假设相反.
我的建议是使用与Fractint用于二维分形的相似的等值面走动.在代码方面，它很复杂，但是它应该减少所需的功能评估数量，而这恰恰是在Fractint中实现的目的.
背景/历史记录:
在1980年代末和1990年代初，我设计了一个分形绘图套件 Fractint .那时的计算机要慢得多，评估每个点的速度也很慢.在Fractint中进行了很多努力，以使其尽可能快地显示分形，但仍能准确显示. (有些人可能还记得通过旋转使用的调色板中的颜色可以实现的颜色循环.这很催眠；
这些分形中的一些是或具有某些区域，其中所需的迭代次数相对于实际的分形集合分形单调非减少-也就是说，没有孤岛"突出，只是迭代步骤中偶尔出现了稳定的增加- -，一种快速评估模式使用 edge tracing (边缘跟踪)来定位迭代次数发生变化的边界:换句话说，用单色填充的区域.关闭一个区域后，它便向该区域的中心追溯，以找到下一个迭代边缘.同样，在关闭后，它可以用正确的恒定颜色填充边界之间的甜甜圈形或C形区域，而无需评估这些像素的功能！
在这里，我们有一个非常相似的情况，只是在三个维度而不是两个维度上.按照定义，每个等值面也是二维的，因此，实际上所有变化都是我们行走边界的方式.
步行本身类似于洪水填充算法，除了我们在三个维度上行走，而且边界是我们要追踪的等值面.
我们在常规网格中对原始函数进行采样，例如 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(含).
如果 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) em>， f(x + 1，y，z)， f(x，y + 1，z)， f(x + 1，y + 1，z)必须小于或等于 c ，并且至少一个大于或等于 c ，并且不都等于 c .
此外，等值面 c 穿过的初始单元格始终可以使用( x MAX ， y MAX ， z MAX )和( x MIN ， y MIN ， z MIN )，将坐标限制为0≤ x MAX ， y MAX ， z MAX ， x MIN ， y MIN ， z MIN ≤ 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)
上面的三个宏显示了如何将网格索引转换回网格坐标.
要走等值面，我们不能依赖递归；通话链太长.取而代之的是，我们有一个应该检查单元格索引的遍历堆栈:
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" */
}
用于验证等值面是否通过当前单元格，将其报告给回调函数并遍历等值面的函数将类似于
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条边，等值面必须与其中至少3条相交.假设我们在角 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 ，以最大程度地减少绘制某些分形所需的迭代次数.有些分形具有值"相同的大面积，因此某些绘图模式不是对区域内的每个点进行采样，而是跟踪了这些区域的 edges .
换句话说，代替定位等值面 c 的单个点， f(x，y，z)= c ，您只能定位一个点，并且然后行走等值面.步行部分的可视化有些复杂，但实际上只是洪水填充的3D变体洪水填充在简单计算机图形学中使用的算法. (实际上，鉴于该磁场沿每个维度单调递减，实际上实际上是一个2D漫游，除了与等值面 c 有关的那些栅格点外，通常只有几个栅格点.这应该真的很有效.)
我敢肯定，有很好的经过同行评审的论文描述了这种技术(可能在多个问题领域)，但是由于我懒于进行比几分钟的Google搜索更好的搜索，我将其留给其他人以找到好的参考资料.抱歉.
为简单起见，现在让我们假设该字段是连续的，并且沿每个维度单调递增.在大小为 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 .这是因为等值面最多只能跨框的任何一条线.
如果函数是单调非递减的，则等值面可以仅与跨越该框的任何线相交一次，但沿该线的交点可以更宽(比一个点).在实践中，您可以通过仅考虑等容腔(具有静态场的体积)的下表面或上表面来解决此问题.即，仅是从低于 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 ，以及以下一项或多项:
- 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 . (对于具有等容孔的单调非递减场，这表示等体积面更接近于原点为等值面.)
当我们知道等值面 c 与一个立方体相交时，基本上可以使用三种方法将该知识转换为一个点(我们认为等值面要相交):
- 该多维数据集具有八个角，每个角位于一个网格点.我们可以选择字段值最接近 c 的角/网格点.
- 我们可以插值-选择一个近似点-等值面 c 与边线/对角线相交.我们可以进行线性插值而无需任何额外的样本，因为我们已经知道交叉边缘/对角线末端的样本. 如果 u = f(x，y，z)< c ，而 v > c 是另一端的样本，沿着该线的线性插值交点出现在(cu)/(vu)，其中0位于(x，y，z)，1位于边缘/对角线的另一端(位于 (x + 1，y，z)，(x，y + 1，z)，(x，y，z + 1)或(x + 1，y + 1，z + 1)).
- 您可以沿边/对角线使用二进制搜索来找到相交点. 每个边缘/对角线需要 n 个额外的样本，才能沿边缘/对角线以 n 位精度获得交点. (由于原始网格与该字段中的细节相比不能太粗糙，否则这些细节将始终不可见，因此通常使用 n = 2 ， n = 3 ， n = 4 或 n = 5 .)
这样获得的等值面 c 的交点可以用于拟合某些表面函数，但我在现实生活中还没有看到.通常， Delaunay三角剖分用于将点集转换为多边形网格，这样很容易可视化.
另一种选择是记住哪个立方体((x，y，z))和边/对角线( X ， Y 或 Z 边缘，或 D 对角线)，每个点都与之相关.然后，您可以直接形成多边形网格. Voxel 技术也可以用于快速绘制部分透明的等值面；每个视线都会检查每个立方体一次，如果存在等值面，则可以使用等值面相交点对表面法线向量进行插值，从而通过射线投射/光线追踪方法(不创建任何多边形网格)生成非常平滑且外观精确的等值面.
在我看来，这个答案需要编辑-至少需要一些睡眠，进一步思考和澄清.欢迎提出问题，建议甚至进行编辑！
如果不仅仅是OP有兴趣，我可以尝试看看是否可以为此整理一个简单的示例C程序.我玩弄了可视化的模拟电子结构，这些字段甚至都不是单调的(尽管采样很便宜).
I have a 3d grid where in each point (x,y,z) on the grid is associated with a cost value. The cost of any point (x,y,z) is not known in advance. To know the cost, we need to make a complex query which is really expensive. One thing we know about the object is that cost is monotonically non-decreasing in all 3 dimensions.
Now given a cost C, I need to find the points (x,y,z) on the surface which have cost C. This has to be done by costing only bare minimum. How to solve my problem?
When I searched online, I am getting contour identification related techniques but all these techniques assume all point's cost is known in advance like Marching cubes method etc. In my case major metric is the number of points costed should be minimum.
It would be helpful if some one can suggest a way to get approximate locations at least if not exact.
Rewritten explanation: (original text, in case it might clarify the idea to someone, is kept unchanged below the line)
We have some function f(x,y,z) in three dimensions, and we wish to find the surface f(x,y,z) = c. Since the function yields a single number, it defines a scalar field, and the surface we are looking for is the isosurface c.
In our case, evaluating the function f(x,y,z) is very costly, so we wish to minimize the number of times we use it. Unfortunately, most isosurface algorithms assume the opposite.
My suggestion is to use a similar isosurface walk as Fractint could use for two-dimensional fractals. Code-wise, it is complicated, but it should minimize the amount of function evaluations needed -- that was exactly the purpose it was implemented in Fractint.
Background / History:
In the late 1980s and early 1990s, I encoutered a fractal drawing suite Fractint. Computers were much slower then, and evaluating each point was painfully slow. A lot of effort was made in Fractint to make it display the fractals as fast as possible, but still accurately. (Some of you might remember the color-cycling it could do, by rotating the colors in the palette used. It was hypnotic; here is a Youtube clip from the 1995 documentary "Colors of Infinity", which both color-cycles and zooms in. Calculating a full-screen fractal could take hours (at high zoom factors, close to the actual fractal set), but then you could (save it as an image and) use the color-cycling to "animate" it.)
Some of those fractals were, or had regions, where the number of iterations needed was monotonically non-decreasing toward the actual fractal set fractal -- that is, no "islands" sticking out, just steady occasional increase in iteration steps --, one fast evaluation mode used edge tracing to locate the boundary where the number of iterations changed: in other words, the regions filled with a single color. After closing a region, it then traced towards the center of that region to find the next iteration edge; after that was closed too, it could just fill the donut- or C-shaped region between those boundaries with the correct constant color, without evaluating the function for those pixels!
Here, we have a very similar situation, except in three dimensions instead of two. Each isosurface is also two-dimensional by definition, so really, all that changes, is how we walk the boundary.
The walk itself is similar to flood fill algorithms, except that we walk in three dimensions, and our boundary is the isosurface we're tracing.
We sample the original function in a regular grid, say an N×N×N grid. (This is not the only possibility, but it is the easiest and most useful case, and what the OP is doing.)
In general, the isosurfaces will not pass through the grid points exactly, but between the grid points. Therefore, our task is to find the grid cells the isosurface passes through.
In an N×N×N regular grid, there are (N-1)×(N-1)×(N-1) cubic cells:
Each cell has eight corners at (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), and (x+1,y+1,z+1), where x,y,Z ∈ ℕ are the integer grid coordinates, 0 ≤ x,y,z ≤ N-2 are the integer grid coordinates.
Carefully note the integer grid coordinate limits. If you think about it, you'll realize that an N×N×N grid has only (N-1)×(N-1)×(N-1) cells, and since we use the grid coordinates for the corner closest to origin, the valid coordinate range for that corner is 0 to N-2, inclusive.
If f(x,y,z) increases monotonically in each dimension, then isosurface c passes through cell (x,y,z) if
f(x,y,z) ≤ c
and at least one of
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
If f(x,y,z) is monotonically non-decreasing -- that is, its partial derivatives are either zero or positive at all points --, then the above locates two-dimensional isosurfaces, and the outer surface for isovolumes (volumes where f(x,y,z) is constant). The inner surface for isovolumes c are then those cells (x,y,z) for which
f(x,y,z) < c
and at least one of
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
Extension to any scalar function:
The approach shown here actually works for any f(x,y,z) that has only one maximum within the sampled region, say at (xMAX,yMAX,zMAX); and only one minimum, say at (xMIN,yMIN,zMIN); with no local maxima or minima within the sampled region.
In that case, the rule is that at least one of 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) must be below or equal to c, and at least one above or equal to c, and not all equal to c.
Also, an initial cell an isosurface c passes through can then always be found using a binary search between (xMAX,yMAX,zMAX) and (xMIN,yMIN,zMIN), limiting the coordinates to 0 ≤ xMAX,yMAX,zMAX,xMIN,yMIN,zMIN ≤ N-2 (to only consider valid cells, in other words).
If the function is not monotonic, locating an initial cell the isosurface c passes through is more complicated. In that case, you need a different approach. (If you can find the grid coordinates for all local maxima and minima, then you can do binary searches from global minimum to local maxima above c, and from local minima below c to global maximum.)
Because we sample the function f(x,y,z) at intervals, we implicitly assume it to be continous. If that is not true -- and you need to show also the discontinuities -- you can augment the grid with discontinuity information at each point (seven boolean flags or bits per grid point, for "discontinuity from (x,y,z) to (x+,y+,z+)"). The surface walking then must also respect (not cross) such discontinuities.
In practice, I would use two arrays to describe the grid: one for cached samples, and one for two flags per grid point. One flag would describe that the cached value exists, and another that the walking routine has already walked the grid cell at that point. The structure I'd use/need for walking and constructing isosurfaces (for a monotonically non-decreasing function sampled in a regular grid) would be
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];
}
and the function to find the cell to start the walk on, using a binary search along the grid diagonal (assuming non-decreasing monotonic function, so all isosurfaces must cross the diagonal):
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)
The three macros above show how to convert the grid index back to grid coordinates.
To walk the isosurface, we cannot rely on recursion; the call chains would be too long. Instead, we have a walk stack for cell indexes we should examine:
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" */
}
The function that verifies that the isosurface passes through the current cell, reports those to a callback function, and walks the isosurface, would be something like
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;
}
In this particular case, I do believe the isosurfaces are best visualized/described using a polygon mesh, with samples within a cell linearly interpolated. Then, each report() call produces one polygon (or one or more flat triangles).
Note that the cell has 12 edges, and the isosurface must cross at least three of these. Let's assume we have two samples at corners c0 and c1, spanned by an edges, with the two corners having coordinates p0=(x0,y0,z0) and p1=(x1,y1,z1) respectively:
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 */
The above check is valid for any kind of continuous function f(x,y,z); for non-monotonic functions the problem is just finding the relevant cells. The isosurface() function needs some changes (the checks wrt. x0y0z0..x1y1z1), according to the rules outlined earlier in this post, but it too can be made to work for any continuous function f(x,y,z) with little effort.
Constructing the polygon/triangle(s) when the samples at the cell corners are known, especially using linear interpolation, is very simple as you can see.
Note that there is usually no reason to worry about the order in which the edges of a cell are checked, as you will almost certainly use vector calculus and cross product in particular to orient the points and polygons. Or, if you like, you can do Delaunay triangulation on the points (3 to 12 for any function, although more than 6 points indicates there are two separate surfaces, I believe) to construct flat polygons.
Questions? Comments?
We have a scalar field f(x,y,z) in three dimensions. The field is costly to sample/evaluate, and we do so only at integer coordinates 0 ≤ x,y,z ∈ ℕ. To visualize the scalar field, we wish to locate one or more isosurfaces (surfaces with a specific f(x,y,z) value), using the minimum number of samples/evaluations.
The approach I'll try to describe here is a variant of the algorithm used in fractint, to minimize the number of iterations needed to draw certain fractals. Some fractals have large areas with the same "value", so instead of sampling every point within the area, certain drawing mode traced the edges of those areas.
In other words, instead of locating individual points of the isosurface c, f(x,y,z) = c, you can locate just one point, and then walk the isosurface. The walk part is a bit complicated to visualize, but it really is just a 3D variant of the flood fill algorithm used in simple computer graphics. (Actually, given the field is monotonically non-decreasing along each dimension, it'll actually be a mostly 2D walk, with typically just a few grid points other than those relevant to the isosurface c sampled. This should be really efficient.)
I'm pretty sure there are good peer-reviewed papers describing this very technique (probably in more than one problem domain), but since I'm too lazy to do a better search than a couple of minutes of Google searches, I leave it to others to find good references. Apologies.
For simplicity, for now, let's assume that the field is continuous and monotonically increasing along each dimension. Within an axis-oriented box of size N×N×N, the field will have a minimum at one corner at origin (0,0,0), a maximum at the far corner from origin, at (N,N,N), with all possible values between the minimum and maximum found along the diagonal from (0,0,0) to (N,N,N). In other words, that every possible isosurface exists and is a continuous 2D surface, excluding points (0,0,0) and (N,N,N), and every such surface intersects the diagonal.
If the field is actually non-continuous, we won't be able to tell, because of our sampling method. In practice, our sampling means we implicitly assume the scalar field is continuous; we will treat is as continuous, whether or not it really is!
If the function is actually monotonically increasing along each dimension, then it is possible to map f(x,y,z)=c to X(y,z)=x, Y(x,z)=y, Z(x,y)=z, although any one of the three is sufficient to define the isosurface c. This is because the isosurface can only cross any line spanning the box in at most one point.
If the function is monotonically non-decreasing instead, the isosurface can intersect any line spanning the box still only once, but the intersection can be wider (than a point) along the line. In practice, you can handle this by considering only the lower or upper surfaces of the isovolumes (volumes with a static field); i.e. only the transition from-lower-than-c-to-c-or-greater, or the transition from-c-or-lower-to-greater-than-c. In all cases, you're not really looking for the isosurface value c, but trying to locate where a pair of the field samples crosses c.
Because we sample the field at regular grid points, and the isosurface rarely (if ever) intersects those grid points exactly, we divide the original box into N×N×N unit-sized cubes, and try to find the cubes the desired isosurface intersects.
Here is a simple illustration of one such cube, at (x,y,z) to (x+1,y+1,z+1):
When the isosurface intersects a cube, it intersects at least one of the edges marked X, Y, or Z, and/or the diagonal marked D. In particular, we'll have f(x,y,z) ≤ c, and one or more of:
- f(x+1,y,z) > c (isosurface c crosses the cube edge marked with X) (Note: In this case, we wish to walk along the y and z dimensions)
- f(x,y+1,z) > c (isosurface c crosses the cube edge marked with Y) (Note: In this case, we wish to walk along the x and z dimensions)
- f(x,y,z+1) > c (isosurface c crosses the cube edge marked with Z) (Note: In this case, we wish to walk along the x and y dimensions)
- f(x+1,y+1,z+1) > c (isosurface c crosses the cube diagonal, marked with D) (Note: In this case, we may need to examine all directly connected grid points, to see which direction we need to walk to.)
Instead of doing a complete search of the original volume, we can just find one such cube, and walk along the cubes to discover the cubes the isosurface intersects.
Since all isosurfaces have to intersect the diagonal from (0,0,0) to (N,N,N), we can find such a cube using just 2+ceil(log2(N)) samples, using a binary search over the cubes on the diagonal. The target cube (i,i,i) is the one for which f(i,i,i) ≤ c and f(i+1,i+1,i+1) > c. (For monotonically non-decreasing fields with isovolumes, this shows the isovolume surface closer to origin as the isosurface.)
When we know that the isosurface c intersects a cube, we can use basically three approaches to convert that knowledge to a point (that we consider the isosurface to intersect):
- The cube has eight corners, each at a grid point. We can pick the corner/grid point with the field value closest to c.
- We can interpolate -- choose an approximate point -- where the isosurface c intersects the edge/diagonal. We can do linear interpolation without any extra samples, since we already know the samples at the ends of the crossed edge/diagonal. If u = f(x,y,z) < c, and v > c is the sample at the other end, the linearly interpolated intersection point along that line occurs at (c-u)/(v-u), with 0 being at (x,y,z), and 1 being at the other end of the edge/diagonal (at (x+1,y,z), (x,y+1,z), (x,y,z+1), or (x+1,y+1,z+1)).
- You can use a binary search along the edge/diagonal, to find the intersection point. This needs n extra samples per edge/diagonal, to get the intersection point at n-bit accuracy along the edge/diagonal. (As the original grid cannot be too coarse compared to the details in the field, or the details will not be visible anyway, you normally use something like n=2, n=3, n=4, or n=5 at most.)
The intersection points for the isosurface c thus obtained can be used for fitting some surface function, but I have not seen that in real life. Typically, Delaunay triangulation is used to convert the point set to a polygon mesh, which is then easy to visualize.
Another option is to remember which cube ((x,y,z)) and edge/diagonal (X, Y, or Z edge, or D for diagonal) each point is related to. Then, you can form a polygon mesh directly. Voxel techniques can also be used to quickly draw partially transparent isosurfaces; each view ray examines each cube once, and if the isosurface is present, the isosurface intersection points can be used to interpolate a surface normal vector, producing very smooth and accurate-looking isosurfaces with raycasting/raytracing methods (without creating any polygon mesh).
It seems to me I this answer is in need of editing -- at minimum, some sleep and further thought, and clarifications. Questions, suggestions, and even edits are welcome!
If there is interest from more than just the OP, I could try and see if I can cobble together a simple example C program for this. I've toyed with visualizing simulated electronic structures, and those fields are not even monotonic (although sampling is cheap).
这篇关于以最小的点成本有效地在3d网格上找到等成本点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
