Integer labeling for zero set estimation

Everybody knows about marching squares or cubes. These methods are typically used for contour line/surface computation. This problem is actually a particular case of more general zero set problem. In marching squares function f acts from R² → R and zero set problem consists of one equation. Algorithm is fairly trivial in its nature: if there’s a function sign change on the grid edge then we have a contour inside. Can we generalize this approach?

What if f is a vector-function and acts, for example, R³ → R²? Surprisingly, there’s nice family of algorithms to handle this, called labeling methods.

Zero set problem

First, let’s define the problem:

For vector-valued function f: Rⁿ → R^m find x values f(x₁, …, x_n) = (0, …, 0).

Labeling methods

Now, let’s return to the labeling. So it’s labeling algorithms and that means that we have grid/space subdivision and for every vertice of it we’ll compute some label values. Then according to them we can find what grid elements/cells present the zero set of our function. Important thing to say is that function f here is approximated by piecewise linear version of itself f_PL, which is easy to define: values on vertices will be the same and linearly interpolated inside the cell. For brevity let’s call Rⁿ → R^m zero set problems by the tuple (n, m). There’s at least two of labeling strategies known: integer and vector labeling. We focus on the integer labeling, because it’s less obvious and more elegant.

First, suppose we have simplicial space subdivision. Then, we should define the labeling rule for every vertice. For example, according to Introduction to Numerical Continuation Methods by Eugene L. Allgower and Kurt Georg appropriate labeling of simplex vertices can be defined as simple as number of initial coordinates of f which are positive. Simplex is called completely labeled if set of its vertices contains every possible labeling. Every simplicial subdivision under certain conditions on the labeling contain such simplex and this fact is connected with Sperner’s Lemma (somehow).

Then the following statement is true:

For function f : Rⁿ → R^m and simplicial subdivision of Rⁿ define labeling L_s(v) as number of positive initial coordinates of f, then in subdivision completely labeled simplices and only them contain non-empty intersection with f_PL(x) zero set.

But why do we need exactly complete labeling in multidimensional case? For now, I didn’t found corresponding theorem with proof in strict form.

However we can gain some intuition by hand-waving explanation and examples. On the following picture, there’s one example of completely labeled simplex for (3,2)-problem. It has four vertices f1, f2, f3, f4 and labelings corresponding to underlying vector-function component signs ({+, -} means f₁ > 0, f₂ < 0). Then it’s easy to see that there exists a point (0, -) between f1 and f2 and (0, +) between f3 and f4 by linearity of f_PL approximation. So we have one zero point u of f_PL on the line between them. Following the same logic another zero point v can be found. Solution of (3, 2)-problem for linear function f_PL should be a line between u and v.

On this example it’s getting obvious why it’s not enough to use sign change condition (like in marching squares). For example, {{+, -}, {-, +}}-edge with different signs for f coordinates will indeed contain zeros of f_PL in every component, but they will be achieved not simultaneously. So we can see, why completeness of labeling is more forcing and necessary condition.

Interesting notes:

• For (3,2)-problem zero set approximation is generally a polygonal chain. On every subdivision simplex which is completely labeled there’s a door in facet where this chain enters the simplex and door out facet, where it leaves it. Can such facets be found with only labelings given? Actually, for every (n+1, n)-problem zero set approximation is a polygonal chain in Rⁿ⁺¹. Can the general rule be formulated for door simplices facets?

Code

Three years ago in my project I implemented integer labeling. More precisely, simplified version (3,2)-algorithm for rough estimation of zero sets.

Let’s see what is going on in this code and where magic happens. Here’s full version (without extra prints) ~120 lines of code:

#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <vector>

#define COMPLETELY_LABELED(x) ((x) == 0x7)

using namespace std;

typedef struct 
{
   int Nxy;
   int Nz;
   
   double x_min;
   double x_max;
   
   double y_min;
   double y_max;
   
   double z_min;
   double z_max;
   
}grid_t;

typedef struct 
{
    double data[2];
}d_pair;

int label_fast(d_pair f_val){ 
    int map[] = {1, 1, 2, 4};
    int r = ((f_val.data[1] >= 0.0) << 1) | (f_val.data[0] >= 0.0);        
    return map[r];
}

double *f_wrapped_vec(int N, double z, void (*ptr)(int, double, double *), double *slice_data_ptr){
    (*ptr)(N, z, slice_data_ptr);
    return slice_data_ptr;   
}

void label_fast_vec(double *slice_data, char *cache, int size){
    int NN = size*size;
    
    for(int l = 0; l < NN; ++l){
        d_pair f_val;       
        f_val.data[0] = slice_data[l], f_val.data[1] = slice_data[NN+l]; 
        cache[l] = label_fast(f_val);
    }
    
}

extern "C" int integer_labeling_3d_vec(grid_t GridDesc, void (*ptr)(int, double, double *), 
                                   double *result, int *size){
    vector <double> x, y, z;
    
    int N = GridDesc.Nxy, Nz = GridDesc.Nz;   
    int NN = N*N, NN_sub = (N-1)*(N-1);

    unsigned long long cnt = 0, iter_cnt = 0, grid_size = Nz - 1;
    
    char *cache_prev = new char[NN], *cache_next = new char[NN], *tmp; 
    double *slice_data_ptr = new double[2*NN];
    
    grid_size *= NN_sub;
    
    x.resize(N), y.resize(N), z.resize(Nz);
        
    for(int i = 0; i < N; ++i){
        x[i] = GridDesc.x_min + (GridDesc.x_max - GridDesc.x_min)*i/(N-1);
        y[i] = GridDesc.y_min + (GridDesc.y_max - GridDesc.y_min)*i/(N-1);
    }
        
    for(int i = 0; i < Nz; ++i){
        z[i] = GridDesc.z_min + (GridDesc.z_max - GridDesc.z_min)*i/(Nz-1);
    }
        
    label_fast_vec(f_wrapped_vec(N, z[0], ptr, slice_data_ptr), cache_prev, N);
       
    for(int k = 0; k < Nz-1; ++k){  
        label_fast_vec(f_wrapped_vec(N, z[k+1], ptr, slice_data_ptr), cache_next, N);
               
        for(int l = 0; l < NN; )
            if((l / N) != N-1 && (l % N) != N-1){
            
                char label_flags = (cache_prev[l] | cache_prev[l + N] | cache_next[l] | cache_next[l + N]);
                ++l;
                label_flags |= (cache_prev[l] | cache_prev[l + N] | cache_next[l] | cache_next[l + N]);
            
                iter_cnt++;
 
                if(iter_cnt % 100000 == 0)
                    cout << iter_cnt << " cubes checked from " << grid_size << endl;                
            
                if(COMPLETELY_LABELED(label_flags)){
                    result[3*cnt    ] = x[(l - 1) / N];
                    result[3*cnt + 1] = y[(l - 1) % N];
                    result[3*cnt + 2] = z[k];
                    
                    ++cnt;
                }
            }else
                 ++l;
                 
        tmp        = cache_prev;
        cache_prev = cache_next;
        cache_next = tmp;
    }
    
    *size = cnt * 3;
    
    cout << cnt << "  completely labeled cubes found, " << grid_size << " checked!" << endl;
    
    delete[] cache_prev;
    delete[] cache_next;
    delete[] slice_data_ptr;
        
    return 0;    
}

Notice that instead of simplicial space subdivision we have cubic grid. Here’s a little trick: let’s call cube completely labeled if set of its vertices contains every possible labeling.

It’s easy to see that completely labeled cube contains completely labeled simplex. If we find such cube we can return necessary simplex, compute line segment where or … just return corner coordinates which were enough for my purposes back then. Because for small enough grids completely labeled cubes also can be used as rough, “voxel” approximation to zero set.

Details

Labeling data of cube layer stored in cache_prev, cache_next variables. On every iteration we load two consecutive function data slices. Cube labeling set is encoded by bits in label_flags variable.

Basically, here we sum up some integer values for every cube of our grid in some speed-optimized manner.

   char label_flags = (cache_prev[l] | cache_prev[l + N] | cache_next[l] | cache_next[l + N]);
   ++l;
   label_flags |= (cache_prev[l] | cache_prev[l + N] | cache_next[l] | cache_next[l + N]);  
   if(COMPLETELY_LABELED(label_flags)){
      /* cube is found */
   }

Array int map[] = {1, 1, 2, 4} is a map between function sign combinations and bit flag values {11 → 100, 10 → 010, 01 → 001, 00 → 001}. Macros COMPLETELY_LABELED(…) is defined simply as three 1’s binary number

#define COMPLETELY_LABELED(x) ((x) == 0x7)

and that’s it.

Note, that computations for each data slice (and actually for single grid cell) are independent. Such code can be easily parallelized with standard frameworks like OMP or even ported on GPU.

Instead of conclusion

Now we have extended picture of algorithms for zero set computation. Moreover, one can view marching algorithms as simplified versions of general simplex labeling.

Full version of the code you can find here.

Repository contains corresponding python package (warning: this version is not well-tested) and where you can find examples and try it by your own.