Rewritten explanation: (the source text, if it can clarify the idea to someone, remains unchanged under the line)
We have some function f (x, y, z) in three dimensions, and we want to find the surface f (x, y, z) = c. Since the function gives a singular, it defines a scalar field , and the surface we are looking for is isosurface c.
In our case, estimating the function f (x, y, z) is very expensive, so we want to minimize the amount of time it uses. Unfortunately, most isosurface algorithms accept the opposite.
My suggestion is to use a similar caterpillar view as Fractint can use for two-dimensional fractals. From a code point of view, this is complicated, but this should minimize the number of required function evaluations - this is what was implemented in Fractint.
History / History:
In the late 1980s and early 1990s, I included a fractal set of fractint drawings. Computers were much slower, and each point was painfully slow. Fractint made a lot of effort to display fractals as quickly as possible, but still for sure. (Some of you may remember the color cycle that he could make by rotating the colors in the palette used. It was hypnotic; here is a Youtube clip from the 1995 documentary, Infinity Colors, which includes both color cycles and scaling: Calculating a full-screen fractal can take several hours (with high scaling factors close to the actual fractal set), but then you can (save it as an image and) use color-cycling to “bring it to life”.)
Some of these fractals were or had regions where the number of required iterations did not monotonously decrease with respect to the actual fractal fractal set, i.e. no “islands” stuck out, just a steady periodic increase in iterative steps - - in one quick estimation mode, edge tracing was used to find the border where the number of iterations changed: in other words, areas filled with the same color. After the region is closed, it is traced to the center of this region to find the next edge of the iteration; after it was closed, it could just fill the donut or C-shape area between these borders with the correct constant color without evaluating the function for these pixels!
Here we have a very similar situation, with the exception of three dimensions instead of two. Each isosurface is also two-dimensional by definition, so really, all that changes is how we pass the boundary.
The walk itself is similar to flood fill algorithms, except that we pass in three dimensions, and our boundary is an isosurface, we are tracing.
We select the original function in the correct grid , say, the grid N × N × N. (This is not the only possibility, but this is the simplest and most useful case, and what the OP does).
In the general case, isosurfaces will not pass exactly through the grid points, but between the grid points. Therefore, our task is to find the mesh cells that the isosurface passes through .
In a regular N × N × N grid, there are (N-1) × (N-1) × (N-1) cubic cells: 
Each cell has eight angles in ( 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 ∈ ℕ - coordinates of the whole grids, 0 ≤ x , y , z ≤ N -2 - coordinates of the whole grids.
Carefully pay attention to the limits of the integer grid coordinates. If you think about it, you will realize that the N × N × N grid has only ( N -1) × ( N -1) × ( N -1), and since we use the grid coordinates for the angle closest to the origin, the valid coordinate range for this angle is 0 to N -2 inclusive.
If f (x, y, z) monotonically increases in each dimension, then the isosurface c passes through the 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, i.e. its partial derivatives are either zero or positive at all points, then in the above case two-dimensional isosurfaces are determined, and the outer surface for isovolums (volumes where f (x, y, z) is constant). The inner surface for isovolume c is 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 of any scalar function:
The approach given here really works for any f (x, y, z) that has only one maximum within the sampling area, for example, ( x MAX , <i> y <sub> MAXsub>, <i> g <sub> MAXsub >); and only one minimum, for example, ( x MIN , y MIN , z <sub> MINsub>); without local highs or lows in the sample area.
In this 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 + ( X, y, z), f (x, y + 1, z), f (x + 1, y + 1, z) must be lower than or equal to c, and at least one higher or equal to c, and not all equal to c.
In addition, the initial cell can be found first, and the isosource c can 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 , <i> y <sub> MAXsub>, <i> g <sub> MAXsub>, <i> x <sub> MIN, y MIN , z MIN ≤ N -2 (only count valid cells, in other words).
If the function is not monotonic, then the placement of the original cell that passes through the isosurface c is more complicated. In this case, you need a different approach. (If you can find the grid coordinates for all local maxima and minima, then you can perform binary searches from the global minimum to local maxima above c and from local minima below c to the global maximum.)
Since we select the function f (x, y, z) with intervals, we implicitly assume it to be continuous. If this is not the case, and you need to show gaps as well, you can enlarge the grid with information about the gap at each point (seven Boolean flags or bits for each grid point, for a “gap from (x, y, z) to (x +, y +, r +) "). Surface walking then must also respect (not cross) such breaks.
In practice, I would use two arrays to describe the grid: one for cached samples and one for two flags for each grid point. One flag will describe that the cached value exists, and another that the running procedure has already passed the grid cell at this point. The structure that I will use / need for walking and constructing isosurfaces (for a monotonically non-decreasing function selected in a regular grid) will be
typedef struct { size_t xsize; size_t ysize; size_t zsize; size_t size; size_t xstride; size_t ystride; size_t zstride; double xorigin; double yorigin; double zorigin; double xunit; double yunit; double zunit; void *data; double *sample(void *data, double x, double y, double z); size_t stack_size; size_t stack_used; size_t *stack; unsigned char *cell; double *cache; } 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 a function to find a cell to start walking using a binary search on the diagonal of the grid (assuming a 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 along an isosurface, we cannot rely on recursion; call chains will be too long. Instead, we have a stack stack for cell indices that we need to examine:
static void grid_push(grid *const g, const size_t cell_index) { 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; } 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) { } g->stack = new_stack; g->stack_size = new_size; } 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; }
A function that checks that the isosurface passes through the current cell, reports the callback function and walks along the isosurface, will 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; const size_t yend = g->ysize - 2; const size_t zend = g->zsize - 2; 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; size_t x, y, z, i; int r; g->stack_used = 0; 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; x = (i / g->xstride) % g->xsize; y = (i / g->ystride) % g->ysize; z = (i / g->zstride) % g->zsize; 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_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); if (x1y0z0 < c && x0y1z0 < c && x1y1z0 < c && x0y0z1 < c && x1y0z1 < c && x0y1z1 < c && x1y1z1 < c) continue; if (report) { r = report(g, x, y, z, c, x0y0z0, x1y0z0, x0y1z0, x1y1z0, x0y0z1, x1y0z1, x0y1z1, x1y1z1); if (r) { errno = 0; return r; } } if (x > 0 && !(cell[i - xstride] & (CELL_WALKED | CELL_STACKED)) && ( x0y1z0 >= c || x0y0z1 >= c )) grid_push(g, i - xstride); if (y > 0 && !(cell[i - ystride] & (CELL_WALKED | CELL_STACKED)) && ( x0y0z1 >= c || x1y0z0 >= c )) grid_push(g, i - ystride); if (z > 0 && !(cell[i - zstride] & (CELL_WALKED | CELL_STACKED)) && ( x1y0z0 >= c || x0y1z0 >= c )) grid_push(g, i - zstride); if (x < xend && !(cell[i + xstride] & (CELL_WALKED | CELL_STACKED)) && ( x0y1z0 >= c || x0y0z1 >= c )) grid_push(g, i + xstride); if (y < xend && !(cell[i + ystride] & (CELL_WALKED | CELL_STACKED)) && ( x1y0z0 >= c || x0y0z1 >= c )) grid_push(g, i + ystride); if (z < xend && !(cell[i + zstride] & (CELL_WALKED | CELL_STACKED)) && ( x1y0z0 >= c || x0y1z0 >= c )) grid_push(g, i + zstride); } errno = 0; return 0; }
In this particular case, I believe that isosurfaces are best visualized / described using a polygonal mesh, with patterns inside the cell being linearly interpolated. Then each call to report () creates one polygon (or one or more flat triangles).
Note that the cell has 12 edges, and at least three of them must intersect the isosurface. Suppose we have two samples at angles c0 and c1 stretched over the edges, with two angles having the coordinates p0 = (x0, y0, z0) and p1 = (x1, y1, z1), respectively:
if (c0 == c && c1 == c) else if (c0 == c) else if (c1 == c) else if (c0 < c && c1 > c) else if (c0 > c && c1 < c) else
The above check is true for any continuous function f (x, y, z); for nonmonotonic functions, the problem is simply finding the corresponding cells. The isosurface () function needs some changes (checks on x0y0z0..x1y1z1) in accordance with the rules described earlier in this post, but this can also be done to work for any continuous function f (x, y, z) with little effort .
Building the polygon / triangle (s) when the patterns at the corners of the cell 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 the cell are checked, since you almost certainly use vector calculus and cross product , in particular, to orient points and polygons. Or, if you want, you can do Delaunay triangulation at points (from 3 to 12 for any function, although more than 6 points indicates that there are two separate surfaces, I suppose) to create flat polygons.
Questions? Comments?
We have a scalar field f (x, y, z) in three dimensions. The field is expensive for sampling / estimation, and we do this only for integer coordinates 0 ≤ x, y, z ∈ ℕ. To visualize a scalar field, we want to find one or more isosurfaces (surfaces with a certain value f (x, y, z)) using the minimum number of samples / estimates.
The approach that I will try to describe here is a variant of the algorithm used in fractint , in order to minimize the number of iterations, it was necessary to draw certain fractals. Some fractals have large areas with the same "value", so instead of selectively selecting each point within the area, a certain drawing mode traces the edges of these areas.
In other words, instead of finding individual points of the isosurface
c ,
f (x, y, z) = c , you can only find one point and then go along the isosurface. Part of the walk is a little difficult to visualize, but in reality it is just a three-dimensional version of the
fill fill algorithm, which is used in simple computer graphics. (In fact, if the field does not monotonically decrease along each dimension, it will be basically a 2D walk, usually just a few grid points, different from those related to the isosurface
c . Sample, should be really effective.)
I'm sure there are good peer-reviewed articles describing this technique (perhaps in more than one problem area), but since I'm too lazy to do a better search than after a couple of minutes of Google searches, I leave it to others to find good links . Apologies
For simplicity, for now, suppose that the field is continuously and monotonically increasing along each dimension. Within an axis-oriented block of size N × N × N, the field will have at least one corner at the beginning (0,0,0), a maximum at the far corner from the origin, at (N, N, N), with all possible values between the minimum and maximum, found diagonally from (0,0,0) to (N, N, N). In other words, there is every possible isosurface and is a continuous two-dimensional surface, excluding points (0,0,0) and (N, N, N), and each such surface intersects the diagonal.
If the field is not actually continuous, we will not be able to say because of our sampling method. In practice, our sample means that we implicitly assume that the scalar field is continuous; we will see it as continuous, whether or not it is!
If the function actually increases monotonically along each dimension, then we can map f (x, y, z) = c to X (y, z) = x , Y (x, z) = y , Z (x, y) = z , although any of the three is sufficient to determine the isosurface c . This is due to the fact that the isosurface can intersect any line covering the field of at most one point.
If the function is monotonically non-decreasing instead, the isosurface can intersect any line that spans the field, still only once, but the intersection can be wider (than a point) along the line. In practice, you can cope with this by considering only the lower or upper surfaces of isovolutes (volumes with a static field); those. only the transition from-below- c -to- c -or-more or the transition from c -or- lower to more-than- <i> s. In all cases, you are not really looking for the value of isosurface c , but trying to find where a pair of field patterns intersect c .
Since we select the field at regular grid points, and the isosurface rarely (if ever) exactly intersects these grid points, we divide the original block by N × N × N and try to find cubes that intersect the desired isosurface.
Here is a simple illustration of one such cube with (x, y, z) k (x + 1, y + 1, z + 1) :
When an 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 will have f (x, y, z) ≤ c and one or more of:
- f (x + 1, y, z)> c (the isosurface c intersects the edge of the cube marked with X) (Note: in this case we want to go along the sizes y and z)
- f (x, y + 1, z)> c (the isosurface c intersects the edge of the cube marked with Y) (Note: in this case we want to go along the sizes x and z)
- f (x, y, z + 1)> c (the isosurface c intersects the edge of the cube marked Z) (Note: in this case we want to go along the sizes x and y)
- f (x + 1, y + 1, z + 1)> c (the isosurface c intersects the diagonal of the cube denoted by D) (Note: in this case, we may need to study all directly connected grid points to see which direction we need go.)
Instead of doing a full search for the source volume, we can just find one such cube and go through the cubes to find cubes intersecting the isosurface.
Since all isosurfaces must intersect the diagonal from (0,0,0) to (N, N, N) , we can find such a cube using only 2 + ceil (log 2 (N)) , using a binary search on the cubes along the diagonal . The target cube (i, i, i) is one for which f (i, i, i) ≤ c and f (i + 1, i + 1, i + 1) > c . (For monotonically nondecreasing fields with isovolume, this shows an isovolutionary surface closer to origin as an isosurface.)
When we know that an isosurface c intersects a cube, we can mainly use three approaches to transform this knowledge into a point (which we consider an isosurface to intersect):
- The cube has eight angles, each at a grid point. We can select a 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 additional samples, since we already know the patterns 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 this line occurs at (cu) / (vu) , and 0 is in (x, y, z ) , and 1 is at the other end of the edge / diagonal (for (x + 1, y, z) , (x, y + 1, z) , (x, y, z + 1) or <i> (x + 1 , y + 1, r + 1)).
- You can use binary edge / diagonal search to find the intersection point. This requires n additional samples per front / diagonal to get the intersection point at n- bit accuracy along the edge / diagonal. (Since the original mesh cannot be too rough compared to the details in the field or the details will not be visible anyway, you usually use something like n = 2 , n = 3 , n = 4 , or n = 5 no more. )
The intersection points for the isosurface c obtained in this way can be used to fit some surface function, but I have not seen this in real life. Typically, Delaunay triangulation is used to transform a point into a polygonal mesh, which is then easily visualized.
Another option is to remember which cube ( (x, y, z) ) and edge / diagonal ( X , Y , or Z edge, or D for the diagonal) each point is connected to. Then you can create a mesh with a polygon. Voxel methods can also be used to quickly isolate partially transparent isosurfaces; each line of sight examines each cube once, and if an isosurface is present, the intersection points of the isosurface can be used to interpolate the surface normal vector, creating very smooth and accurate perspective isosurfaces with raycasting / raytracing methods (without creating a polygon mesh).
It seems to me that this answer needs to be edited - at least, some sleep and further thought, as well as clarifications. Questions, suggestions and even corrections are welcome!
If you have interest not only from OP, I can try and see if I can combine a simple C program example for this. I played with the visualization of simulated electronic structures, and these fields are not even monotonous (although the selection is cheap).