Meshing a signed distance function using the libfive C/C++ API

[doug-moen]

I am working on adding some new meshing algorithms to Curv. Given all the work that @mkeeter has done in exploring meshing algorithms, I'd like to use libfive as one of my meshing algorithms. It's just not obvious how to do the integration, so I'm looking for suggestions, documentation and examples. (Eg, has someone successfully done this before in open source?)

The only input I have to libfive is a signed distance function that maps [x,y,z] onto a signed distance. That's it. I don't have gradients, I have no ability to do interval arithmetic, etc. I have a vague notion that this should be done by defining a subclass of libfive::Oracle, but that class is full of pure virtual functions that I have no idea how to implement for my use case.

What I actually want is a C++ function that takes a signed distance function, a bounding box, a grid resolution, and a quality parameter, and returns an indexed triangle mesh. So how do I write this function using libfive?

[mkeeter]

Your vague notion is correct!

Here's a sketch of an oracle which wraps a function (float, float, float) -> float:

#include "libfive/oracle/oracle_storage.hpp"

class CurvOracle : public libfive::OracleStorage<LIBFIVE_EVAL_ARRAY_SIZE> {
public:
    CurvOracle(std::function<float(float, float, float)> f)
        : f(f)
    {
        // Nothing to do here
    }

    void evalInterval(libfive::Interval& out) override {
        // Just pick a big ambiguous value.
        out = {-10000.0, 10000.0};
    }

    void evalPoint(float& out, size_t index=0) override {
        const auto pt = points.col(index);
        out = f(pt.x(), pt.y(), pt.z());
    }

    void checkAmbiguous(
            Eigen::Block<Eigen::Array<bool, 1, LIBFIVE_EVAL_ARRAY_SIZE>,
                         1, Eigen::Dynamic> /* out */) override
    {
        // Nothing to do here, because we can only find one derivative
        // per point (points on sharp features may not be handled correctly)
    }

    // Find one derivative with partial differences
    void evalFeatures(boost::container::small_vector<libfive::Feature, 4>& out) override {
        const float EPSILON = 1e-6;
        float center, dx, dy, dz;

        Eigen::Vector3f before = points.col(0);
        evalPoint(center);

        points.col(0) = before + Eigen::Vector3f(EPSILON, 0.0, 0.0);
        evalPoint(dx);

        points.col(0) = before + Eigen::Vector3f(0.0, EPSILON, 0.0);
        evalPoint(dy);

        points.col(0) = before + Eigen::Vector3f(0.0, 0.0, EPSILON);
        evalPoint(dz);

        points.col(0) = before;

        out.push_back(Eigen::Vector3f(
            (dx - center) / EPSILON,
            (dy - center) / EPSILON,
            (dz - center) / EPSILON));
    }

protected:
    std::function<float(float, float, float)> f;
};

In addition, you need an OracleClause which constructs the oracle on demand, which would be something like this:

#include "libfive/oracle/oracle_clause.hpp"

class CurvOracleClause: public libfive::OracleClause
{
public:
    CurvOracleClause(std::function<float(float, float, float)> f)
        : f(f)
    {
        // Nothing to do here
    }

    std::unique_ptr<libfive::Oracle> getOracle() const override {
        return std::unique_ptr<libfive::Oracle>(new CurvOracle(f));
    }

    std::string name() const override { return "CurvOracleClause"; }

protected:
    std::function<float(float, float, float)> f;
};

To put this OracleClause into a Tree then render it, you'd do something like this:

#include "libfive/tree/tree.hpp"
#include "libfive/render/brep/mesh.hpp"

using namespace libfive;

// ...
// Assume you already have std::function<float(float,float,float)> f
auto oc = Tree(std::unique_ptr<OracleClause>(new CurvOracleClause(f));

BRepSettings bs; // edit this to change resolution, quality, etc
Region<3> r({-1, -1, -1}, {1, 1, 1}); // Edit this to change rendering

Mesh::render(oc, r, bs)->saveSTL("out.stl");

(All of this is untested, but should convey the infrastructure to build / use a minimal Oracle)

[doug-moen]

Thanks for the detailed response! I will report back later when this succeeds or fails.

[doug-moen]

I just published a libfive integration to Curv. I'm quite pleased with the mesh generation in my own testing, and another person will start testing it soon.

One issue I haven't addressed. The code

    void evalInterval(libfive::Interval& out) override {
        // Just pick a big ambiguous value.
        out = {-10000.0, 10000.0};
    }

is not scale independent. I'm using this exact code. If I render a huge cube (width 9999999) then the result is not great. I considered using infinity instead of 10,000 but that could cause numerical problems. I suspect that the number 10,000 should be replaced by some function of the shape's bounding box. Any suggestions based on your understanding of the math?

There's one other issue that I resolved. Libfive crashes with a SEGV (null pointer read) unless I change settings.workers = 1. I'm using multi-threaded evaluation with the other mesher, so it's a puzzle. Since there is an easy workaround, it is not a showstopper.

[mkeeter]

If the top-level tree is just this CurvOracle, then returning [-inf, inf] should be fine; however, if you're transforming it further, then I agree that the math could get iffy. You could also return [NaN, NaN], which should cause interval subdivision to always subdivide (down to individual voxels).

I'd have to learn more about the segfault, but that shouldn't be libfive's fault: it should build one Evaluator per thread here. I wonder if there's some shared state in your oracle that's getting messed up if multiple threads access it simultaneously?

[doug-moen]

I tested [-10000,10000], [-inf,inf] and [NaN,NaN] and they all have the same behaviour for various cube tests. A "small" cube is rendered as 12 triangles, while a large cube (cube size=9999999, settings.min_feature=215443) renders as thousands of triangles with rounded edges. I think there must be fixed magic numbers in the libfive mesh rendering in order for the cube to render differently at different scales.

I'm going to use [-inf,inf] since it seems the "least magic" and is not overtly scale dependent.

I will later create an issue for the segfault that doesn't involve any Curv code.

[doug-moen]

This question is now answered. I'll address unexpected behaviour by filing github issues.

[doug-moen]

More feedback for people who try this.

In CurvOracle::evalFeatures, Matt says "Find one derivative with partial differences". Matt's code uses a fixed EPSILON value of 1e-6 for the finite difference.

With certain models, the quality of the result is highly dependent on the value of epsilon. I have a model where 1e-6 produced very bad results, which were improved by increasing the epsilon value. My current code sets epsilon to settings.min_feature / 10, and so far that seems to work well as a default.

Libfive Dual Contouring produces perfect-looking output for CAD-like models with exact distance fields. (Cubes, cylinders, unions, rotations, etc.) The models that have meshing artifacts are the same ones that have "bent" or "twisted" distance fields, where the gradient is different from the surface normal. Since Dual Contouring uses surface normals to position and orient triangles on the surface, I believe that the contract for evalFeatures is to return a surface normal, not to return a distance field gradient. If my theory is correct, then I should get better results if I use Newton's method to find the surface point closest to the query point, then compute an accurate normal, instead of finite differences, for those models where there is a difference. (This won't work perfectly because Newton's method won't provide additional information if the query point is already on the surface.)