3D printable models with openFrameworks

For a rethread.art project we wanted to make a sculpture out of a dataset consisting of javascript function calls being executed in the browser when searching the web. Of course my first thought was "that shouldn't be too hard". By documenting the process here, I hope that next time maybe it won't be.

Have a look at the code on GitHub.

Throughout the project, Blender has been immensly helpful for inspecting the meshes created. Especially the 3d-Print Toolbox has been invaluable for finding rare bugs and exactly which vertices are affected.

The Wrong Way

To be honest there's nothing wrong with this technique if all you want is a plane with mountains. Except that you might get huuuuuge files that will eat your memory for breakfast and crash Fusion 360 and Cura (Blender amazingly managed to open even the very hugest files, although it took a few minutes).

My first go at parametrically generating a 3D printable model was to make a plane where the y axis is offset by by a certain force at each vertex. The openFrameworks docs have helpfully provided an example of how to construct the plane in the ofMesh docs. The rest was simple enough: store all the points where you want to offset the plane together with the radius of their influence and the offset at their peaks. Then use a curve to apply to the gravity level within the radius (smoothcurve1 or smoothcurve2 below, inspired by the smoothstep function in GLSL). For each vertex on the plane, get the offset (or the "gravity") from all AttractionPoints and add them together. Easy peasy!

class AttractionPoint {
public:
  glm::vec2 pos;
  float gravity; // the strength of the gravity pull
  float radius; // the maximum radius of influence over which the gravity strength is interpolated
  AttractionPoint(glm::vec2 p, float g, float r): pos(p), gravity(g), radius(r) {}

  float getGravity(glm::vec2 p) {
    float distance = glm::distance(pos, p);
    if(distance > radius) return 0.0f;
    float amount = 1. - (distance/radius);
    // add a nice curve to the gravity influence of the point
    amount = pow(amount, 3.);
    amount = smoothcurve2(amount);
    return gravity * amount;
  }

  // these functions assume a value 0-1
  float smoothcurve1(float x) {
    return x * x * (3 - 2 * x);
  }
  float smoothcurve2(float x) {
    return x * x * x * (x * (x * 6 - 15) + 10);
  }
};

This created a very detailed and smooth plane that felt amazing when 3D printed. Unfortunately a bit too detailed, the model couldn't even be opened in Fusion 360 for editing. The real issue, though, was when we decided we needed holes that went through the mesh to the bottom. Each function in the dataset would be a hole, each script a hill.

Smooth detailed mesh

I tried to adapt the plane idea by pulling down vertices to the base plate, leaving vertices out if they fell into a hole and used all sorts of (what I thought was) clever tricks. Blender became my best friend, as it actually let me open and inspect the mesh, albeit after a few minutes of wait. The 3D printing toolkit let me see just how messed up the model was (hundreds of thousands of non-manifold edges etc.).

Messed up holes with holes in them

After trying to brute force it by raising the resolution for far too long I decided to do it the right way.

The Right Way

For a model to be 3d printable

  • every edge should be connected to exactly 2 faces
    • if an edge is connected to only one face the model has a hole
    • if an edge is connected to more than two faces it's ambiguous
  • no faces should intersect since intersecting faces make it hard to tell what is inside and what is outside

Additionally, the normals, i.e. what side of the triangle faces in or out, should be consistent in the whole model. When making triangles manually, the order of the vertices (the "winding" of the triangle) is what decides its normal.

Delaunay Triangulation

A better way of creating a custom mesh is to use a proper triangulation algorithm that you can feed points and get triangles. Delaunay triangulation is one such algorithm and luckily OpenCV already has an implementation of it, hidden in the Subdiv2D class. With the ofxCv addon, OpenCV is quite easy to use within openFrameworks.

OpenCV Subdiv2D is quite particular about the coordinates you feed it.

  • The max widtand height are non-inclusive so make sure your maximum point is less than the width/height.
  • Feeding it any number outside of its bounding box causes an abrupt crash.
  • It also doesn't like negative coordinates it seems.
  • OpenCV Points are integers which is another pitfall when trying to find points after triangulation since openFrameworks Points (which are really glm::vec3) use floats. The easy solution is to scale up the model until you get the resolution you need. Since we're not doing a vertex per pixel anymore this doesn't increase the complexity or disk size of the model.

Also, Delaunay Triangulation looks pretty cool just by itself, doesn't it?

Triangles, lots of triangles

This also means that every AttractionPoint needs generate the points (vertices) to be used for the triangulation, with the added benefit that the resolution of the slope can be fine-tuned:

class AttractionPoint {
public:
  // [...]

  vector<glm::vec2> getPoints() {
    vector<glm::vec2> points;
    points.push_back(pos);
    int numCircles = 10;
    for(int i = 1; i <= numCircles; i++) {
      int numPointsPerCircle = max(4.0, (radius/numCircles) * 0.1 * i);
      for(int k = 0; k < numPointsPerCircle; k++) {
        float angle = (TWO_PI/numPointsPerCircle) * k;
        float r = (radius/numCircles) * i;
        glm::vec2 tp = pos + glm::vec2(cos(angle)*r, sin(angle)*r);
        // round to nearest integer because OpenCV Points are using integers
        tp = glm::round(tp);
        points.push_back(tp);
      }
    }
    return points;
  }
  // [...]
};

The main part of the triangulation is in this snippet:

Subdiv2D subdivTop = Subdiv2D(rect);
for(auto& p : topPoints) {
  // only add points within bounds, points outside the rect can lead to a crash in OpenCV
  if(p.x >= minW && p.x < minW + width && p.y >= minH && p.y < minH + height) {
    subdivTop.insert(toCv(p));
  }
}
// triangulate using Delaunay triangulation and put the results in triangleList which is a vector<Vec6f>
subdivTop.getTriangleList(triangleList);

for( size_t i = 0; i < triangleList.size(); i++ ) {
  // process the triangles one at a time
  Vec6f t = triangleList[i];
  pt[0] = Point(t[0], t[1]);
  pt[1] = Point(t[2], t[3]);
  pt[2] = Point(t[4], t[5]);

  // add triangles to the mesh
  if (rect.contains(pt[0]) && rect.contains(pt[1]) && rect.contains(pt[2])) {
    // used to store the indices of the vertexes used in the mesh (reused old vertexes or newly created)
    int pointIndices[3]; 

    int numPointsInHole = 0;
    size_t lastHoleIndex = -1; // the index of the last point that was a hole point
    // true if points belong to holes, but different holes. only filter out triangles where all points belong to the same hole
    bool differentHoles = false; 
    for(int i = 0; i < 3; i++) {
      auto& p = pt[i];
      int holeMeshIndex = currentIndex;
      // check if the point is already in the topIndexedPoints list of previous points
      IndexedPoint ip = IndexedPoint(toOf(p), currentIndex);
      auto pointIt = std::find(topIndexedPoints.begin(), topIndexedPoints.end(), ip);
      if(pointIt != topIndexedPoints.end()) {
        // point already exists, use it
        pointIndices[i] = pointIt->index;
        holeMeshIndex = pointIt->index;
      } else {
        // point is new, create it
        float offset = getGravityAtPoint(toOf(p));
        mesh.addVertex(ofPoint(p.x, offset, p.y));

        // check if the point is a corner, in that case we save the index for the walls
        for(int i = 0; i < numWallPoints; i++) {
          if(toOf(p) == wallPoints[i]) topCornerIndices.push_back(IndexedPoint(toOf(p), currentIndex));
        }
        pointIndices[i] = currentIndex;
        topIndexedPoints.push_back(ip);
        currentIndex++;
      }
      // check if the point is in a hole's point list. if it is, add it as an index in that hole
      size_t holeIndex = isPointInHoleThenAddIndex(toOf(p), Plane::TOP, holeMeshIndex);
      if(holeIndex != -1) {
        // the point belonged to a hole
        numPointsInHole++;
        // if the holeIndex (which is a unique identifier for a hole, not a point in a hole) is not
        // the same as the last hole found, the hole points belong to different holes
        if(lastHoleIndex == -1) lastHoleIndex = holeIndex;
        else if(holeIndex != lastHoleIndex) {
          differentHoles = true;
        }
      }
    }
    // don't add the triangle if all points are hole points from the same hole
    if(numPointsInHole < 3 || differentHoles) {
      mesh.addIndex(pointIndices[2]);
      mesh.addIndex(pointIndices[1]);
      mesh.addIndex(pointIndices[0]);
    }
  }
}

Manifold edges

With the Subdiv2D algorithm you feed it with points and you get triangles back. That means that the same point will be given back to you several times when it is part of different triangles. Because we don't want to add the same point several times to the same mesh we need to manually keep track of them and their index in the mesh. That way, whenever we get a vertex from the triangulation we can check if it's already in the mesh and if so what index it has. This is especially crucial for applying filters to the holes and walls as outlined below.

Holes

With the Delaunay triangulation in place, I could just add points for the hole and the algorithm would make triangles around them. Adding the same points for both the top and bottom planes meant I could just connect those with triangles like a cylinder. This came with the added benefit of being able to decide on the resolution for the holes.

Here are the interesting parts of the HolePoint class:

class HolePoint {
public:
  glm::vec2 pos;
  // if points are added too close to a hole the Delaunay triangulation can create triangles
  // covering the hole and creating non-manifold edges. Therefore a clearence radius is needed.
  float clearRadius;
  float holeRadius;
  vector<IndexedPoint> topPlaneHoleIndices;
  vector<IndexedPoint> bottomPlaneHoleIndices;
  vector<glm::vec2> cirPoints; // contains cylinder points and wall points if there are any
  vector<glm::vec2> wallPoints; // only contains the wall points

  HolePoint(glm::vec2 p, float cr, float hr): pos(p), clearRadius(cr), holeRadius(hr) {}

  HoleRelation insideHole(glm::vec2 p) {
    float distance = glm::distance(pos, p);
    if(distance < holeRadius) return HoleRelation::HOLE;
    else if (distance < clearRadius) return HoleRelation::CLEARANCE;
    return HoleRelation::NONE;
  }

  void generateCircumferencePoints() {
    cirPoints.clear();
    int numPoints = 6;
    for(int i = 0; i < numPoints; i++) {
      float angle = (TWO_PI/numPoints) * i;
      glm::vec2 tp = glm::vec2(pos.x + cos(angle) * holeRadius, pos.y + sin(angle) * holeRadius);
      // round to nearest integer because OpenCV Points are using integers
      tp = glm::round(tp);
      cirPoints.push_back(tp);
    }
  }

  void addIndex(glm::vec2 p, Plane plane, size_t index) {
    IndexedPoint ip = IndexedPoint(p, index);
    if(plane == Plane::TOP) {
      // add the point only if it does not already exist
      // this would perhaps be cleaner with an std::set
      if(find(topPlaneHoleIndices.begin(), topPlaneHoleIndices.end(), ip) == topPlaneHoleIndices.end())
        topPlaneHoleIndices.push_back(ip);
    } else if (plane == Plane::BOTTOM) {
      if(find(bottomPlaneHoleIndices.begin(), bottomPlaneHoleIndices.end(), ip) == bottomPlaneHoleIndices.end())
        bottomPlaneHoleIndices.push_back(ip);
    }
  }

  void addCylinderToMesh(ofMesh& mesh) {
    if(topPlaneHoleIndices.size() != bottomPlaneHoleIndices.size()) {
      ofLogError("HolePoint::addCylinderToMesh") << "top and bottom hole index vectors have different sizes";
      ofLogError("HolePoint::addCylinderToMesh") << "top: " << topPlaneHoleIndices.size();
      for(auto& ip : topPlaneHoleIndices) {
        ofLogError("HolePoint::addCylinderToMesh") << ip.index << " " << ip.pos.x << ", " << ip.pos.y;
      }
      ofLogError("HolePoint::addCylinderToMesh") << "bottom: " << bottomPlaneHoleIndices.size();
      for(auto& ip : bottomPlaneHoleIndices) {
        ofLogError("HolePoint::addCylinderToMesh") << ip.index << " " << ip.pos.x << ", " << ip.pos.y;
      }
    }
    // sort the points in circular order
    // 1. calculate the angles of the points in polar coordinates relative to the center of the hole
    for(auto& i : topPlaneHoleIndices) {
      i.calculateAngle(pos);
    }
    for(auto& i : bottomPlaneHoleIndices) {
      i.calculateAngle(pos);
    }
    // 2. sort the points
    std::sort(topPlaneHoleIndices.begin(), topPlaneHoleIndices.end());
    std::sort(bottomPlaneHoleIndices.begin(), bottomPlaneHoleIndices.end());

    // create the cylinder
    for (int i = 0; i < topPlaneHoleIndices.size(); i++){
      // for the last wall we want to wrap around to the first point
      int nextIndex = (i+1)%(topPlaneHoleIndices.size());
      // filter out pieces belonging to a wall
      bool isWallPiece = false;
      // if both positions exist in the wallPoints vector the segment is a wall segment that should not be
      if(find(wallPoints.begin(), wallPoints.end(), topPlaneHoleIndices[i].pos) != wallPoints.end()
        && find(wallPoints.begin(), wallPoints.end(), topPlaneHoleIndices[nextIndex].pos) != wallPoints.end()) 
      {
        isWallPiece = true;
      }
      if(!isWallPiece) {
        mesh.addIndex(bottomPlaneHoleIndices[i].index);           // 10
        mesh.addIndex(topPlaneHoleIndices[nextIndex].index);           // 1
        mesh.addIndex(topPlaneHoleIndices[i].index);               // 0

        mesh.addIndex(bottomPlaneHoleIndices[i].index);           // 10
        mesh.addIndex(bottomPlaneHoleIndices[nextIndex].index);       // 11
        mesh.addIndex(topPlaneHoleIndices[nextIndex].index);           // 1
      }
    }
  }
};

There were a few more hiccups on my quest for a completely manifold mesh with holes. Because of how the Delaunay triangulation works, a point very close to a hole can skip the closest point when forming triangles leading to non-manifold edges over the hole:

Hole covered by a face

The simple solution to this was to filter out non-hole points that were too close to a hole. This, with a few more tiny tweaks, was enough to create a proper working mesh if the whole thing could be made in one piece.

Whole mesh working

But we wanted it much bigger than the 3d printer could manage in one go so we're back to splitting it into a grid. The sides of the mesh had problems with non-manifold edges: the triangulation algorithm doesn't want to make triangles that are too skinny (where one angle is much larger than the others) so a point close to the edge would not connect with the wall points if that would result in a skinny triangle. Solution: filter out points too close to the wall and use more wall points per walls to create less skinny triangles that connect to the wall. Success! (with slightly lower resolution right at the edge)

Grid pieces with clear edges

At this point the mesh had no non-manifold edges, but the holes would still not go all the way through when the model was sliced in Cura. Turns out the normals of the walls were flipped compared to all the other normals so just winding the wall triangles the other way flipped the normals right and let the holes pass all the way through.

Normals corrected

Cutting this mesh into a grid of pieces caused two new problems: the indices from the last piece have to be cleared between pieces (easy) and walls intersecting with a hole have to be broken up (harder). The solution to the latter issue was in four steps:

  1. For each hole, check if the hole intersects with a wall and if it does, add the two points where the hole intersects the wall to the hole and to the wall.
  2. Remove any wall points that fall into a hole.
  3. Filter out wall segments where all points are hole points or where the mid point falls in a hole.
  4. When making the hole cylinders, don't create a segment between two points that are part of the wall.

Doing all of this resulted in this cleanly sliced mesh:

Holes through walls

Hope this is useful for someone, even if that someone is myself in a year!

Resources:

rethread.art is a KTH project supported by CASTOR, WASP and IoD.

Previous Post