c - 在3D网格上以最小的点成本有效地找到等成本点

Question

我有一个 3d网格，其中网格上的每个点(x，y，z) 与 成本值相关联。一点(x，y，z)的成本是，事先不知道。要知道成本，我们需要进行一个复杂的查询，这确实很昂贵。我们对该对象了解的一件事是，成本在所有3个维度内单调不变。

现在给定成本C，我需要找出表面上具有成本C 的点(x，y，z)。这必须通过仅花费最小的来完成。如何解决我的问题？

当我在线搜索时，我得到了与轮廓识别相关的技术，但是所有这些技术都假定所有点的成本都像Marching cubes方法等那样事先知道。在我的情况下，主要度量标准是要计算的点数应最少。

如果有人可以提出一种至少可以获取近似位置的方法，那将很有帮助。

Answer 1

改写说明:
(原始文本(在可能向某人澄清主意的情况下，在该行下保持不变))

我们在三维上有一些函数f(x，y，z)，我们希望找到曲面f(x，y，z)= c。由于该函数产生单个数字，因此它定义了a scalar field，而我们要查找的表面是isosurface c。

在我们的案例中，评估函数f(x，y，z)的成本非常高，因此我们希望将使用它的次数减至最少。不幸的是，大多数等值面算法都相反。

我的建议是使用与Fractint用于二维分形的相似的等值面走动。在代码方面，它很复杂，但是它应该使所需的功能评估的数量减至最少，这恰恰是在Fractint中实现的目的。

背景/历史:

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.

我们在 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(含)。

如果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

扩展到任何标量函数:

The approach shown here actually works for any f(x,y,z) that has only one maximum within the sampled region, say at (x_MAX,y_MAX,z_MAX); and only one minimum, say at (x_MIN,y_MIN,z_MIN); 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 (x_MAX,y_MAX,z_MAX) and (x_MIN,y_MIN,z_MIN), limiting the coordinates to 0 ≤ x_MAX,y_MAX,z_MAX,x_MIN,y_MIN,z_MIN ≤ 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.

在实践中，我将使用两个数组来描述网格:一个用于缓存的样本，一个用于每个网格点的两个标志。一个标志将描述存在高速缓存的值，而另一个标志则表明该行进例程已在该点上行进了网格单元。我用于/构造等值面(用于在规则网格中采样的单调非递减函数)的结构将是

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):很少的努力。

如您所见，在知道单元角处的样本时(尤其是使用线性插值)，构造多边形/三角形非常简单。

请注意，通常无需担心检查单元格边缘的顺序，因为您几乎可以肯定会使用 vector 演算和 cross product特别是确定点和多边形的方向。或者，如果愿意，可以在这些点上执行 Delaunay triangulation(对于任何函数，为3到12，尽管我相信超过6个点表示有两个单独的表面)以构造平面多边形。

问题？注释？



我们在三个维度上都有一个标量场f(x，y，z)。该字段的采样/评估成本很高，我们仅在整数坐标0≤x，y，z∈do时才这样做。为了可视化标量场，我们希望使用最少数量的样本/评估来定位一个或多个等值面(具有特定f(x，y，z)值的面)。

我将在此处尝试描述的方法是 fractint中使用的算法的一种变体，以最大程度地减少绘制某些分形所需的迭代次数。一些分形的大区域具有相同的“值”，因此某些绘图模式不是对区域内的每个点进行采样，而是跟踪了这些区域的边缘。

换句话说，除了定位等值面c的单个点(f(x，y，z)= c)外，您还可以只定位一个点，然后沿等值面走动。步行部分的可视化有些复杂，但实际上只是简单计算机图形中使用的 flood fill算法的3D变体。 (实际上，鉴于该字段沿每个维度单调递减，实际上实际上将是一个2D行走，通常只有几个网格点，而不是与等值面c相关的网格点。这应该非常有效。)

我敢肯定，有很好的经过同行评审的论文都描述了这种技术(可能在多个问题领域)，但是由于我懒于进行比几分钟的Google搜索更好的搜索，因此我将其保留给别人找好的引用。道歉。

为了简单起见，现在让我们假设该字段是连续的，并且沿每个维度单调递增。在大小为N×N×N的轴心框中，该字段在原点(0,0,0)的一个角处具有最小值，在原点最远的角(N，N，N)处具有最大值，所有可能的值都在从(0,0,0)到(N，N，N)的对角线中找到的最小值和最大值之间。换句话说，除了点(0,0,0)和(N，N，N)，每个可能的等值面都存在并且是连续的2D曲面，并且每个这样的曲面都与对角线相交。

如果该字段实际上是不连续的，则由于我们的采样方法，我们将无法分辨。在实践中，我们的采样意味着我们隐式地假设标量场是连续的。无论是否真的，我们都会将其视为连续的!

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.

因为我们在规则的网格点处对场进行采样，并且等值面很少(如果有的话)与这些网格点精确相交，所以我们将原始框划分为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(log2(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是另一端的样本，则沿该直线进行线性插值的交点出现在(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 triangulation用于将点集转换为多边形网格，然后可以轻松地对其进行可视化。

另一个选择是记住每个点与哪个立方体((x，y，z))和边缘/对角线(X，Y或Z边缘，或D对角线)相关。然后，您可以直接形成多边形网格。 Voxel技术还可以用于快速绘制部分透明的等值面。每个视线都会检查每个立方体一次，如果存在等值面，则可以使用等值面相交点对表面法线 vector 进行插值，从而通过射线转换/光线追踪方法(不创建任何多边形网格)生成非常平滑且外观精确的等值面。

在我看来，这个答案需要修改-至少需要一些 sleep ，进一步思考和澄清。欢迎提出问题，建议甚至进行编辑!

如果不仅仅是OP引起了人们的兴趣，我可以尝试看看是否可以为此整理一个简单的示例C程序。我玩弄了可视化的模拟电子结构，这些字段甚至都不是单调的(尽管采样很便宜)。