Skip to content

Instantly share code, notes, and snippets.

@tmpvar
Created September 27, 2021 07:47
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save tmpvar/fbc4e1ea5e8d2d8a4077c4c1077cae04 to your computer and use it in GitHub Desktop.
Save tmpvar/fbc4e1ea5e8d2d8a4077c4c1077cae04 to your computer and use it in GitHub Desktop.
single threaded cpu based sdf -> octree evaluator
// 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment