single threaded cpu based sdf -> octree evaluator

evaluate-cpu.cpp

// a single threaded SDF->octree evaluator on the cpu.

#pragma CFLAGS=-I ../../include -march=native

#include <stdio.h>
#include <glm/glm.hpp>
#include <types.h>
using namespace glm;

#include <vector>
#include <random>
#include <chrono>
#include <fstream>
#include <thread>
using namespace std;

#define MAX_OPERATIONS (1<<17)
#define BIG_BUFFER_SIZE 1000000000

static chrono::high_resolution_clock::time_point program_start_time = chrono::high_resolution_clock::now();
struct Prof {
  double start = 0.0;
  const char *name;
  Prof(const char *name): name(name) {
    start = now();
  }


  f64 now() {
    auto t = chrono::high_resolution_clock::now();
    auto d = chrono::duration_cast<chrono::duration<double>>(
      t - program_start_time
    );
    return d.count();
  }

  ~Prof() {
    printf("%s: %.4fms\n", name, (now() - start) * 1000.0);
  }
};

struct Bounds {
  Bounds(const vec3 &center, const vec3 &radius)
    : center(center), radius(radius)
  {}

  vec3 center;
  vec3 radius;
};

struct Operation {
  enum Shape {
    Sphere,
    Box,
  };

  enum Action {
    Add,
    Cut
  };

  Bounds bounds;

  Shape shape;
  Action action;


  Operation(const vec3 &center, const vec3 &radius, Shape shape = Shape::Sphere, Action action = Action::Add)
    : bounds(center, radius), shape(shape), action(action)
  {}

  void setup(const vec3 &center, const vec3 &radius, Shape shape = Shape::Sphere, Action action = Action::Add) {
    this->bounds.center = center;
    this->bounds.radius = radius;
    this->shape = shape;
    this->action = action;
  }
};

struct OpList {
  u32 count = 0;
  Operation *ops = nullptr;
  OpList() {
    ops = (Operation *)malloc(sizeof(Operation) * MAX_OPERATIONS);
  }
  ~OpList() {
    free(ops);
  }
  Operation *add() {
    u32 i = count;
    count++;
    return &ops[i];
  }
};

f32 sdBox( const vec3 &p, const vec3 &b ) {
  const vec3 q = abs(p) - b;
  return length(glm::max(q, vec3(0.0))) + glm::min(glm::max(q.x, glm::max(q.y, q.z)), 0.0f);
}

f32 sdSphere(const vec3 &p, const float s ) {
  return length(p) - s;
}

// TODO: rotated objects
f32 eval(f32 &d, const vec3 &cell_center, const Operation &op) {
  f32 ld = FLT_MAX;
  const vec3 p = cell_center - op.bounds.center;
  switch (op.shape) {
    case Operation::Shape::Sphere:
      ld = sdSphere(p, op.bounds.radius.x);
      break;
    case Operation::Shape::Box:
      ld = sdBox(p, op.bounds.radius);
      break;
    default:
      return FLT_MAX;
  }

  switch (op.action) {
    case Operation::Action::Add:
      d = glm::min(ld, d);
      break;
    case Operation::Action::Cut:
      d = glm::max(-ld, d);
      break;
    default:
      return FLT_MAX;
  }

  return ld;
}

struct OpIndexList {
  struct Entry {
    u32 start = 0;
    u32 length = 0;
    OpIndexList *list = nullptr;

    void add(u32 idx) {
      length++;
      list->indices[list->count++] = idx;
    }
  };

  u32 count = 0;
  u32 *indices = nullptr;

  OpIndexList() {
    indices = (u32 *)malloc(sizeof(u32) * 100000000);
  }

  ~OpIndexList() {
    free(indices);
  }

  Entry add() {
    Entry e = {};
    e.start = count;
    e.list = this;
    return e;
  }
};

struct ActiveCell {
  OpIndexList::Entry ops;
  Bounds bounds;
  u32 depth;
  u32 child_count = 0;
  ActiveCell *children[8] = {};
  ActiveCell *parent = nullptr;
  ActiveCell(const OpIndexList::Entry &ops, const Bounds &bounds, u32 depth)
    : ops(ops), bounds(bounds), depth(depth)
  {}

  void setup(const OpIndexList::Entry &ops, const Bounds &bounds, u32 depth) {
    this->ops = ops;
    this->bounds = bounds;
    this->depth = depth;
  }
};

struct ActiveCellList {
  u32 count = 0;
  ActiveCell *cells = nullptr;
  ActiveCellList() {
    cells = (ActiveCell *)malloc(sizeof(ActiveCell) * 100000000);
  }

  ~ActiveCellList() {
    free(cells);
  }

  ActiveCell *add() {
    u32 i = count;
    count++;
    return &cells[i];
  }
};

void evalCell(
  const vec3 &cell_center,
  const vec3 &cell_radius,
  const f32 h,
  const OpList *op_list,
  const OpIndexList::Entry &active_ops,
  ActiveCellList *output,
  u32 depth = 0
) {
  f32 d = FLT_MAX;

  const u32 op_count = active_ops.length;
  auto filtered_ops = active_ops.list->add();
  for (u32 i=0; i<op_count; i++) {
    u32 op_idx = active_ops.list->indices[active_ops.start + i];
    const auto &op = op_list->ops[op_idx];
    f32 op_d = eval(
      d,
      cell_center,
      op
    );

    if (op_d < h) {
      filtered_ops.add(op_idx);
    }
  }

  if (abs(d) < h) {
    // printf("subdivide d(%f) h(%f) depth(%u) ops(%u -> %u)\n", d, h, depth, active_ops.size(), filtered_ops.size());
    vec3 r = cell_radius * 0.5f;
    if (r.x < 1.0f) {
      output->add()->setup(
        filtered_ops,
        Bounds(cell_center, cell_radius),
        depth
      );
      return;
    }

    vec3 c = cell_center;
    f32 lower_h = length(r);
    evalCell(c + vec3(-r.x, -r.y, -r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3( r.x, -r.y, -r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3(-r.x,  r.y, -r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3( r.x,  r.y, -r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3(-r.x, -r.y, r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3( r.x, -r.y, r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3(-r.x,  r.y, r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
    evalCell(c + vec3( r.x,  r.y, r.z), r, lower_h, op_list, filtered_ops, output, depth+1);
  }
}

int main() {
  OpList *op_list = new OpList;
  ActiveCellList *output = new ActiveCellList;
  OpIndexList *op_index_list = new OpIndexList;
  OpIndexList::Entry active_ops = op_index_list->add();

  // build a list of ops
  {
    auto _ = Prof("setup");

    // wiffle cube
    if (1) {
      op_list->add()->setup(
        vec3(1024),
        vec3(512),
        Operation::Shape::Box
      );

      op_list->add()->setup(
        vec3(1024),
        vec3(600),
        Operation::Shape::Sphere,
        Operation::Action::Cut
      );
    }

    // random distribution
    if (1) {
      std::random_device rd;
      std::mt19937 gen(rd());
      std::uniform_real_distribution<> dis(0.0, 1.0);

      for (u32 i=0; i<1000; i++) {
        vec3 radius = 5.0f + vec3(dis(gen) * 50.0);
        vec3 center = vec3(dis(gen), dis(gen), dis(gen)) * 2048.0f;
        op_list->add()->setup(center, radius, Operation::Shape::Box);
        op_list->add()->setup(center, radius * 1.25f, Operation::Shape::Sphere, Operation::Action::Cut);
      }
    }

    {
      u32 c = op_list->count;

      for (u32 i=0; i<c; i++) {
        active_ops.add(i);
      }
    }
  }

  {
    auto _ = Prof("eval");
    evalCell(vec3(1024), vec3(1024), length(vec3(1024)), op_list, active_ops, output);
  }

  printf("summary:\n");
  printf("  total leaves: %u\n", output->count);

  u64 total_op_indices = 0;
  u64 total_leaves = 0;
  f64 mean_ops = 0;
  u32 max_depth = 0;

  const u32 cell_count = output->count;
  for (u32 i=0; i<cell_count; i++) {
    total_op_indices += output->cells[i].ops.length;
  }

  printf("  total ops: %u\n", op_list->count);
  printf("  total op index entries: %lu\n", total_op_indices);
  printf("  max depth: %u\n", max_depth);
  mean_ops = static_cast<f64>(total_op_indices) / static_cast<f64>(output->count);
  printf("  mean: %.2f\n", mean_ops);

  f64 sdev_sum = 0;
  for (u32 i=0; i<cell_count; i++) {
    sdev_sum += glm::pow(
      static_cast<f64>(output->cells[i].ops.length) - mean_ops,
      2.0
    );
  }

  f64 sample_variance = sdev_sum / static_cast<f64>(total_leaves - 1ULL);

  printf("  stddev: %.4f\n", glm::sqrt(sample_variance));

  printf("  total index size: %.2fMB\n", static_cast<f64>(total_op_indices * 4ULL) / 1024.0 / 1024.0);

  #ifndef PROFILE
    ofstream f;
    f.open("out.xyz");
    for (u32 i=0; i<cell_count; i++) {
      const auto &cell = output->cells[i];
      f << cell.bounds.center.x << " ";
      f << cell.bounds.center.y << " ";
      f << cell.bounds.center.z << " ";
      f << "1.0\n";
    }
    f.flush();
    f.close();
    printf("done\n");
  #endif


  return 0;
}