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.

*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.

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.).

After trying to brute force it by raising the resolution for far too long I decided to do it 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.

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
`Point`

s are integers which is another pitfall when trying to find points after triangulation since openFrameworks`Point`

s (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?

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]);
}
}
}
```

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.

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:

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.

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)

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.

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:

- 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.
- Remove any wall points that fall into a hole.
- Filter out wall segments where all points are hole points or where the mid point falls in a hole.
- 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:

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

Resources:

- https://www.hubs.com/knowledge-base/fixing-most-common-stl-file-errors/
- https://pages.mtu.edu/~shene/COURSES/cs3621/SLIDES/Mesh.pdf
- https://en.wikipedia.org/wiki/Delaunay_triangulation

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