在三维网格上有效地找到等成本点,并且点数成本最低

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;
    

以上是关于在三维网格上有效地找到等成本点,并且点数成本最低的主要内容,如果未能解决你的问题,请参考以下文章

在 C++ 中通过网格/矩阵找到成本优化路径

如何找到成本最高的路径

经管-8

Java 的最低成本斯坦纳树库

寻找最小阶梯成本的动态规划问题的错误答案

搭建一个网站最低需要多少钱?如何以最低成本做一个网站?