code: def eval_grid_octree(
      coords,
      eval_func,
      init_resolution=64,
      threshold=0.05,
      num_samples=512 * 512 * 512,
  ):
      resolution = coords.shape[1:4]
  
      sdf = np.zeros(resolution)
  
      notprocessed = np.zeros(resolution, dtype=np.bool)
      notprocessed[:-1, :-1, :-1] = True
      grid_mask = np.zeros(resolution, dtype=np.bool)
  
      reso = resolution[0] // init_resolution
  
      while reso > 0:
          # subdivide the grid
          grid_mask[
              0 : resolution[0] : reso,
              0 : resolution[1] : reso,
              0 : resolution[2] : reso,
          ] = True
          # test samples in this iteration
          test_mask = np.logical_and(grid_mask, notprocessed)
          # print('step size:', reso, 'test sample size:', test_mask.sum())
          points = coords[:, test_mask]
  
          sdf[test_mask] = batch_eval(
              points, eval_func, num_samples=num_samples
          )
          notprocessed[test_mask] = False
  
          # do interpolation
          if reso <= 1:
              break
          x_grid = np.arange(0, resolution[0], reso)
          y_grid = np.arange(0, resolution[1], reso)
          z_grid = np.arange(0, resolution[2], reso)
  
          v = sdf[
              tuple(np.meshgrid(x_grid, y_grid, z_grid, indexing="ij"))
          ]
  
          v0 = v[:-1, :-1, :-1]
          v1 = v[:-1, :-1, 1:]
          v2 = v[:-1, 1:, :-1]
          v3 = v[:-1, 1:, 1:]
          v4 = v[1:, :-1, :-1]
          v5 = v[1:, :-1, 1:]
          v6 = v[1:, 1:, :-1]
          v7 = v[1:, 1:, 1:]
  
          x_grid = x_grid[:-1] + reso // 2
          y_grid = y_grid[:-1] + reso // 2
          z_grid = z_grid[:-1] + reso // 2
  
          nonprocessed_grid = notprocessed[
              tuple(np.meshgrid(x_grid, y_grid, z_grid, indexing="ij"))
          ]
  
          v = np.stack([v0, v1, v2, v3, v4, v5, v6, v7], 0)
          v_min = v.min(0)
          v_max = v.max(0)
          v = 0.5 * (v_min + v_max)
  
          skip_grid = np.logical_and(
              ((v_max - v_min) < threshold), nonprocessed_grid
          )
  
          n_x = resolution[0] // reso
          n_y = resolution[1] // reso
          n_z = resolution[2] // reso
  
          xs, ys, zs = np.where(skip_grid)
          for x, y, z in zip(xs * reso, ys * reso, zs * reso):
              sdf[
                  x : (x + reso + 1),
                  y : (y + reso + 1),
                  z : (z + reso + 1),
              ] = v[x // reso, y // reso, z // reso]
              notprocessed[
                  x : (x + reso + 1),
                  y : (y + reso + 1),
                  z : (z + reso + 1),
              ] = False
          reso //= 2
  
      return sdf.reshape(resolution)
  
imports: import numpy
  import time
label_desc: Use image node_id 5 for calling a global function with example usage: eval_grid_octree(coords, eval_func, init_resolution, threshold, num_samples) and returns: sdf