sdf interpreter

[pca006132]

Some initial implementation for sdf interpreter. Names are quite random, I have no idea how I should name these.

Basically, there are four things here:

1. A general interval arithmetic class (sdf/interval.h). Note that it is not robust against things like NaN and rounding errors (yet).
2. A simple API for expressing symbolic expressions (sdf/value.h), and this file is intended to go to the public API later, probably with another name. It is currently using methods, but I am happy to make it use regular functions.
3. A code gen module (sdf/context.h) that transforms the value into low level opcode, and does some simple optimization. Currently, it only performs constant folding and convert (a * b + c) into FMA, but I plan to do a bit more than these. It is also very limited because I haven't yet added spilling support, and it can only support 256 variables for now.
4. An evaluator (tape.h) that is pretty fast, and does both pointwise evaluation and interval evaluation. Optimizations related to interval evaluation are not implemented yet, but I have idea about how to do them. I will probably add more opcode later as well, e.g. load store for spilling.

And in the test file (sdf_tape_test.cpp), I implemented a simple version of manifold dual contouring (the algorithm that uses an octree and interval arithmetic) but without actually constructing the mesh. I do this just to test the number of evaluations needed. Note that I haven't yet fully debugged my implementation, so there may be bugs that can throw off the result...

Some initial statistics:

Model	Grid Eval Count	Interval Eval Count	Grid Time	Interval Time
Gyroid	143009	111177	7ms	5ms
Blobs	16 363 009	1 474 337	700ms	96ms

Do note that this is run without parallelization. I haven't yet tried to implement the interval evaluation with parallelization yet.

[codecov]

Codecov Report

Attention: Patch coverage is 65.09972% with 245 lines in your changes missing coverage. Please review.

Project coverage is 88.77%. Comparing base (25ce1b1) to head (9716bc3).
Report is 10 commits behind head on master.

Files with missing lines	Patch %	Lines
src/sdf/tape.h	35.46%	111 Missing ⚠️
src/sdf/context.cpp	81.93%	54 Missing ⚠️
src/sdf/interval.h	44.68%	52 Missing ⚠️
src/sdf/value.cpp	81.91%	17 Missing ⚠️
src/sdf/context.h	76.66%	7 Missing ⚠️
src/sdf/affine_value.h	70.00%	3 Missing ⚠️
src/sdf/value.h	66.66%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1122      +/-   ##
==========================================
- Coverage   91.46%   88.77%   -2.70%     
==========================================
  Files          30       37       +7     
  Lines        5908     6617     +709     
==========================================
+ Hits         5404     5874     +470     
- Misses        504      743     +239     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[pca006132]

The API, code generation (with simple optimization) and execution should be done. There are limits on the number of variables and registers that can be used in the code gen, but I will lift that in a day or two, hopefully.

The actual SDF to mesh algorithm is not yet written, will probably take more time to do that.

[pca006132]

This is a pretty general interval arithmetic implementation that covers common math functions. I tried to make things as precise as possible, but things like mod, for example, is hard to make precise and probably subject to changes later.

[pca006132]

This is not the typical std::fmod in c++, this returns something positive. This is nicer to interval arithmetic, and IMO more natural. For negative m however, I am not sure what I should do with that, maybe give something negative?

[pca006132]

This is the bytecode interpreter implementation, which is optimized a bit for faster common arithmetic operations, but can still extend to other domains. It uses operator overloading for common arithmetic operations, but fallback to specialization for other functions (because I can't add a sin method for double, for example).

[pca006132]

Because SDFs typically doesn't have control flow, there is no control flow here as well. Complicated SDFs will require runtime optimization for fast execution and optimize away the unused path, which will be implemented later.

[pca006132]

Simple recursive algorithm to convert the value AST to IR in context.h. context.h handles the IR, performs optimization and eventually output a bytecode tape.

We don't do any optimization here, so it is quite simple (except using an explicit loop instead of recursion).

[pca006132]

AST type for users to build their SDF. Intended to be public.

[pca006132]

We do grid evaluation to check if the SDF is correct.

[pca006132]

ok spilling done

[pca006132]

Originally I copied it from the repository mentioned, but it turns out it was too buggy and misses some required functions, so I changed it quite a bit.

[hrgdavor]

I am following commits and drooling a bit :) ... pls post some pics when you will have some examples to show :)

[pca006132]

I'm still working on the code generation and tape optimization stuff. I'm more comfortable with that and less familiar with the mesh generation code.

@elalish would be great if you can give a level set function that takes two functions, one that outputs an interval and one regular, and I can plug the stuff here in.

[elalish]

@elalish would be great if you can give a level set function that takes two functions, one that outputs an interval and one regular, and I can plug the stuff here in.

I'm not quite sure what you mean; can you give an example?

[pca006132]

@elalish change NearSurface into something like this that is recursive and takes a function that outputs an interval.

[pca006132]

@elalish example here

[pca006132]

Re. Soundness of the interval arithmetic implementation, I think there is a succinct argument why this is sound w.r.t. floating point arithmetic (not real arithmetic, though).

What we want to argue is:

• Unary: Given $Y = f(X)$ for intervals $X, Y$ and some unary operation $f$, we must have $f(x) \in Y$ for all $x in X$.
• Binary: Given $Z = X \ast Y$ for intervals $X, Y, Z$ and some binary operation $\ast$, we must have $x \ast y \in Z$ for all $x \in X$ and $y \in Y$.

The argument basically boils down to monotonicity and case analysis. For monotonic functions, we know that $x \le x' \implies f(x) \le f(x')$. If $x \in X$, i.e. $x_0 \le x \le x_1$ for $X = [x_0, x_1]$, we know that $f(x_0) \le f(x) \le f(x_1)$, so we can have $Y = [f(x_0), f(x_1)]$. This works as long as the rounding mode is consistent, because rounding schemes are monotonic.

For addition, it is also monotonic w.r.t. each of its argument, i.e. $x \le x' \implies x + c \le x' + c$. We can show that $x_0 + y_0 \le x + y_0 \le x + y \le x + y_1 \le x_1 + y_1$. And this get our bound of $Z = [x_0 + y_0, x_1 + y_1]$ as desired.

For multiplication, it is slightly more complicated because we need to split the intervals and union them back, as well as do case analysis. If $X, Y \ge 0$, multiplication is monotonic w.r.t. each argument, so the reasoning is the same as addition. If any of them is negative, we can negate them to get to the $X, Y \ge 0$ case, and negate the final result if only one operand is negative. For cases in which the interval crosses zero, we can split the interval, compute the result, and union them back. For the union result, it must cross zero, so the endpoints of the final range $Z$ must be the product of the endpoints of $X$ and $Y$.

Negation may cause issues when the rounding mode is round up/down (whether it is a real issue requires some more thought), but for the common case of round-to-nearest and ties round-to-even, the result should be the same regardless of sign, so the argument should work for us.