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

CRM*_*CRM 6 c algorithm math contour data-structures

我有一个3d网格,其中网格上的每个点(x,y,z) 成本相关联.任何点(x,y,z)的成本都不是预先知道的.要知道成本,我们需要进行一个非常昂贵的复杂查询.我们对该目标了解的一件事是,在所有3个维度中,成本是单调不减少的.

现在给出成本C,我需要在表面上找到成本为C 的点(x,y,z).这必须通过仅花费最低成本来完成.如何解决我的问题?

当我在线搜索时,我正在获得与轮廓识别相关的技术,但所有这些技术都假设所有点的成本都是预先知道的,比如Marching cubes方法等.在我的情况下,主要指标是成本计算的点数应该是最小的.

如果某人能够建议一种获得近似位置的方法,至少如果不准确的话会很有帮助.

Nom*_*mal 6

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.

我的建议是使用类似的等值面行走,因为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.

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.)

通常,等值面不会精确地穿过网格点,而是在网格点之间.因此,我们的任务是找到等值面通过的网格单元.

N × N × N规则网格中,存在(N -1)×(N -1)×(N -1)个立方单元: 网格中的一个立方体单元格

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
Run Code Online (Sandbox Code Playgroud)

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
Run Code Online (Sandbox Code Playgroud)

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
Run Code Online (Sandbox Code Playgroud)

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
Run Code Online (Sandbox Code Playgroud)

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];
}
Run Code Online (Sandbox Code Playgroud)

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) 
Run Code Online (Sandbox Code Playgroud)

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" */
}
Run Code Online (Sandbox Code Playgroud)

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;
}
Run Code Online (Sandbox Code Playgroud)

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 */
Run Code Online (Sandbox Code Playgroud)

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):

  1. The cube has eight corners, each at a grid point. We can pick the corner/grid point with the field value closest to c.
  2. 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)).
  3. 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).


在我看来,这个答案需要编辑 - 至少,一些睡眠和进一步思考,以及澄清.欢迎提出问题,建议,甚至编辑!

如果不仅仅是OP的兴趣,我可以尝试看看我是否可以拼凑一个简单的示例C程序.我玩弄了可视化的模拟电子结构,这些领域甚至都不是单调的(虽然采样很便宜).


小智 2

您应该查看这篇文章,其中讨论了二维情况,并让您深入了解不同的方法: http://leetcode.com/2010/10/searching-2d-sorted-matrix.html

在我看来,逐步线性搜索(在第二部分中)对您来说将是一个很好的第一步,因为它很容易应用于 3-d 情况,并且确实不需要很多经验来理解。
因为这非常简单而且仍然非常高效,所以我会继续使用它,看看它是否符合您在 3-d 中处理的数据类型的需求。

但是,如果您的唯一目标是性能,那么您应该将二进制分区应用于 3-d。这变得有点复杂,因为他所说的“二进制分区”本质上变成了“二进制平面分区”。
所以你没有一条线将你的矩阵划分成两个可能的更小的矩阵。
相反,您有一个平面将立方体划分为 2 个可能的更小的立方体。
为了使该平面(或矩阵)中的搜索有效,您首先必须实现他的方法之一:)。
然后你用较小的立方体重复一切。
请记住,以非常有效的方式实现这一点(即记住内存访问)并非易事。