Differential growth

Today i'd like to explore how we can implement a differential growth algorithm in Max 9 using Jitter geometry and JavaScript. Before jumping into the code, let's talk for a second about growth algorithms.

In the context of computer graphics, growth algorithms refer to techniques that simulate the natural growth processes of plants, trees, and other organic structures. These algorithms generate realistic models of biological forms by mimicking the patterns and rules that govern their growth in nature.

Differential growth refers to a process in which different parts of a structure grow at different rates, resulting in complex shapes and patterns. This phenomenon is commonly observed in nature, such as in the growth of leaves, shells, coral reefs, and biological tissues. In computer graphics, differential growth algorithms simulate these natural growth processes to create intricate, organic forms.

So, let's see how differential growth algorithms work:

Algorithm overview

To understand the mechanics of differential growth algorithms, let's start considering a 2D example: Let's drop some points on a plane, somewhat arranged in a circle, and let's connect adjacent points through edges:

Growth is about expansion, so if we want to expand this figure we need to apply some forces to the points to make them move. The first force that we must apply is a repulsive force: within a certain radius, each point repels the others. Let's make an example.

Let's draw a circle of radius r from the center of each point to find the neighboring points (in yellow) that apply a repulsive force on the red one. The yellow arrows are the individual force vectors, and the green one is the sum of the forces. Let's now skip one step ahead and let's apply these forces. Under the repulsive force, the points shift an amount and direction proportional to the force vectors.

You can notice that, due to the repulsive force, some points moved outwards while others moved inwards. So far, so good, but if you think about it, the repulsive force alone can't make the figure grow: if the force is bound within a radius, once the points are far apart from each other of a distance greater than that radius, the figure stops expanding. Here comes into play the next step, that is edge splitting.

If an edge connecting two points is longer than a fixed "split length" (the purple segment), a new point is introduced at the middle of the edge, splitting it in half. This new point is like fuel for the growing mechanism, as now points have a new neighbor to push away. Repeating repulsion and point insertion, the figure starts growing:

In the beginning, we said that points are subjected to some forceS, and, in fact, the repulsive force is not the only one. There's another force that tends to keep the points at the same distance from each other. We could define this force with the somewhat inaccurate term of spring force. Such force can be repulsive or attractive and involves two points connected by an edge.

If the edge length is shorter than half the split length, then the force is repulsive and pushes the two points away in opposite directions along the edge; if the edge length is greater than half the split length, then the force is attractive and pushes the points towards each other.

The repulsive force and the spring force combined make the figure grow while maintaining a "rounded" shape. Without the spring force, the mesh is likely to degenerate under the pressure of the repulsive force alone.

I didn't yet give any formula to chew on, but they'll come later. For now, i'd like to add an extra dimension to our example, and consider the implications of the growth algorithm described above in 3D. In reality, it doesn't change much...if not at all! Every point is now a vertex of a triangle, but the same principles apply to the 3D case as well: each vertex is subjected to a repulsive force and a spring force that develops along the triangle's edges.

Some things become a bit more complicated, but not more complex, for example, the insertion of a new vertex into the mesh. Consider adding one point to a 2-dimensional geometry structure; you only have to add a new point to the list of points, a new edge to the list of edges, and update the points' and edges' pointers to the corresponding edges and points. In 3D, we have more things to care about; inserting a new point means updating the whole geometry structure for the edges, halfedges, and faces connected to the new vertex.

Another thing that becomes trickier is maintaining a well-structured mesh: if, in the 2D example, splitting the edges and introducing new points doesn't affect the structure of the mesh overall, in 3D, we'll have to take care of keeping a good arrangement of the mesh triangles.

That said, let's see how to turn our beloved duck into a Cronenberg-esque nightmare!

Implementation

In this section i'll go through the implementation details. Let's start with the initial setup:

Open the patch differential.growth.maxpat.

This patch can perform a stepped differential growth on a closed mesh.

As you can see, the patch is straightforward, and all the action takes place inside the custom script. So, let's take a look at it -> Double-click on v8 geom.differential.growth.js

The script has two input functions:

function dictionary(dict){

  // Create a reference to the Max dictionary using the dictionary name
  let d = new Dict(dict);

  // Turn this into a JavaScript object
  let fullGeometryDesc = JSON.parse(d.stringify());

  // Get the first geometry
  geom = fullGeometryDesc.geomlist[0];  

  //initialize forces to [0,0,0]
  initForces();

  //output the mesh
  draw();
}

dictionary() takes a Max dictionary containing a Jitter geometry and turns it into a JavaScript object; then, the function initForce() initializes at [0,0,0] an array of size geom.vertices_size that will contain the forces to apply to the vertices.

function initForces(){

  forces = [];
  for(let i = geom.vertices_size - 1; i >= 0; i--) forces.push([0,0,0]);
}

When v8 receives a bang message, the script performs a growth iteration and outputs the vertex position matrix.

function bang(){

  //make it grow!
  grow();

  //output the mesh
  draw();
}

Now let's dive into grow(), the core function of the script:

function grow(){

  findMinMax();
  initGrid(); 
  fillGrid();   

  calcRepulsiveForces();  
  calcSpringForces();
  
  applyForces();

  insertNewVertices();

  remesh();
  smooth();
}

Judging by the functions' names, you can identify parts of the algorithm described in the introductory section with a few extra things. The remaining part of the current section has been subdivided into paragraphs to give you an exhaustive and organized explanation of how these functions work.

1) The parameters

The growth algorithm is controlled by a few attributes:

• growth_speed is a factor used to reduce or increase the applied forces easily
• split_len determines the edge length over which a new vertex is inserted into the mesh; it also affects the spring force and the search radius for the repulsive force
• radius_ratio determines the extent of the force radius as a function of split_len ($radius = radiusratio*splitlen$)
• remesh_enable activates/deactivates re-meshing (see paragraph 7)
• smooth_iterations sets the number of smoothing iterations (see paragraph 7)
• smooth_amt controls the amount of smoothing (see paragraph 7)
• inertia controls the amount of inertia (see paragraph 5)

2) The acceleration structure

Let's think about the forces that we have to apply to the mesh, and in particular, let's focus on the repulsive force. We said that each vertex exerts a repulsive force over the other vertices if the distance separating them is less than a fixed radius; in other words, we should evaluate each vertex against each other vertex to know, at least, if they're close enough to interact. Such an operation is veeeeery slow because for $N$ vertices, we have to perform $(N^2)/2$ distance checks. Is there anything we can do to reduce the problem's complexity?

We know that two vertices repel each other only if their relative distance is less than a radius; knowing this, we could create a grid of $cellsize=radius$ and insert in each cell of the grid the indexes of the vertices living inside it. If we do so, we then have to compute the repulsive force only between the vertices living in neighboring grid cells.

Open the JavaScript file geom.differential.growth.include.gridfunctions.js

function findMinMax(){

  //initialize meshMin and meshMax
  meshMin = [100000, 100000, 100000];
  meshMax = [-100000, -100000, -100000];

  let p;
  //iterate over the vertices to find the minimum and maximum
  for(let i = geom.vertices_size - 1; i >= 0; i--){

    p = geom.vertices[i].point;

    meshMin[0] = Math.min(meshMin[0], p[0]);
    meshMin[1] = Math.min(meshMin[1], p[1]);
    meshMin[2] = Math.min(meshMin[2], p[2]);
    meshMax[0] = Math.max(meshMax[0], p[0]);
    meshMax[1] = Math.max(meshMax[1], p[1]);
    meshMax[2] = Math.max(meshMax[2], p[2]);
  }

  //compute the extend of the mesh
  shape_size = [  meshMax[0] - meshMin[0], 
                  meshMax[1] - meshMin[1], 
                  meshMax[2] - meshMin[2]];
}

The function findMinMax() computes the extent of the mesh by iterating through vertices and computing the minimum and maximum vertex positions. The mesh's extent is used to determine the acceleration grid's size.

function initGrid(){

  //compute the number of cells in each dimension of the grid
  grid_size = new Array(3);
  grid_size[0] = Math.floor(shape_size[0] / radius);
  grid_size[1] = Math.floor(shape_size[1] / radius);
  grid_size[2] = Math.floor(shape_size[2] / radius);

  //fill the grid with an empty array
  grid = new Array(grid_size[0]);
  for(let x = grid_size[0]-1; x >= 0; x--){
    grid[x] = new Array(grid_size[1]);
    for(let y = grid_size[1]-1; y >= 0; y--){
      grid[x][y] = new Array(grid_size[2]);
      for(let z = grid_size[2]-1; z >= 0; z--){
        grid[x][y][z] = {vertices: []};
      }
    }
  }
}

The function initGrid() creates and initializes the grid. In theory, one could allocate a grid of arbitrary size, but to save memory and computing power, we make it fit the mesh tightly by computing the number of cells in each dimension as $gridsize[n] = floor(shapesize[n] / radius)$. Each cell of the grid is initialized with an empty array.

function pos2grid(p){

  //convert the vertex position into a grid position
  p = sub(p, meshMin);
  p = div(p, shape_size);
  p = floor(mul(p, [grid_size[0]-1, grid_size[1]-1, grid_size[2]-1] ));
  return p;
}

function fillGrid(){

  let g;

  //iterate over the vertices
  for(let i = geom.vertices_size - 1; i >= 0; i--){

    //calc grid position
    g = pos2grid(geom.vertices[i].point);

    //push vertex i in this grid cell
    grid[g[0]][g[1]][g[2]].vertices.push(i);
  }
}

The function fillGrid() iterates over the mesh's vertices and pushes their indexes into the corresponding grid cell.

3) Repulsive force

Let's get to the juicy part. As said before, vertices exert a repulsive force on the other vertices closer than radius. We already prepared an acceleration structure to speed up the process, so let's see how the repulsive force computation is implemented:

function calcRepulsiveForces(){

  let gx,gy,gz,gtx,gty,gtz,x,y,z;

  //for each cell of the grid
  for(gx = grid_size[0]-1; gx >= 0; gx--){
    for(gy = grid_size[1]-1; gy >= 0; gy--){
      for(gz = grid_size[2]-1; gz >= 0; gz--){


        //explore the neighborood of the grid cell of coordinates [gx][gy][gz]
        for(x = 0; x <= 1; x++){
          gtx = gx + x;
          if(gtx >= grid_size[0]) continue; //skip if outside the grid
          for(y = 0; y <= 1; y++){
            gty = gy + y;
            if(gty >= grid_size[1]) continue; //skip if outside the grid
            for(z = 0; z <= 1; z++){
              gtz = gz + z;
              if(gtz >= grid_size[2]) continue; //skip if outside the grid

              if(gx == gtx && gy == gty && gz == gtz){
                calcRepulsionSameCell(gx,gy,gz);
              } else {
                calcRepulsionDifferentCell(gx,gy,gz,gtx,gty,gtz);
              }
            } 
          }
        }
      }
    }
  }
}

First of all, we set up a nester for loop to iterate over the grid cells.

  //for each cell of the grid
  for(gx = grid_size[0]-1; gx >= 0; gx--){
    for(gy = grid_size[1]-1; gy >= 0; gy--){
      for(gz = grid_size[2]-1; gz >= 0; gz--){

        [...]
      }
    }
  }

For each grid cell, we run a new nested cycle to access one of these grid cell neighbors and the cell itself.

        //explore the neighborood of the grid cell of coordinates [gx][gy][gz]
        for(x = 0; x <= 1; x++){
          gtx = gx + x;
          if(gtx >= grid_size[0]) continue; //skip if outside the grid
          for(y = 0; y <= 1; y++){
            gty = gy + y;
            if(gty >= grid_size[1]) continue; //skip if outside the grid
            for(z = 0; z <= 1; z++){
              gtz = gz + z;
              if(gtz >= grid_size[2]) continue; //skip if outside the grid

            [...]
            }
          } 
        }

This image shows the first few iterations:

• the red grid cell is the current cell in the iteration (of coordinates $gx$, $gy$, $gz$), and it's the cell to process.
• the yellow cells are 7 neighboring cells (of coordinates $gtx$, $gty$, $gtz$).
• in green, the cells that have been processed already.
• the blue cells are the ones still to be processed.

The vertices in the red cell must be tested against the vertices contained in the yellow cells and in itself.

              if(gx == 0 gtx && gy == gty && gz == gtz){
                calcRepulsionSameCell(gx,gy,gz);
              } else {
                calcRepulsionDifferentCell(gx,gy,gz,gtx,gty,gtz);
              }

The first step is to compute the repulsive interactions between the vertices in the red cell. To do so, we run the function calcRepulsionSameCell()

function calcRepulsionSameCell(gx, gy, gz){

  //declare the variables
  let v0,v1,p0,p1,diff,dist,dist2,dir,force,k,j;

  //for each vertex of the grid cell of coordinates [gx][gy][gz]
  for(k = grid[gx][gy][gz].vertices.length-1; k >= 0; k--){

    //get the vertex index
    v0 = grid[gx][gy][gz].vertices[k];

    //get the vertex position
    p0 = geom.vertices[v0].point;

    //for each other vertex of the grid cell of coordinates [gx][gy][gz]
    for(let j = grid[gx][gy][gz].vertices.length-1; j >= k+1; j--){

      //get the vertex index
      v1 = grid[gx][gy][gz].vertices[j];

      //get the vertex position
      p1 = geom.vertices[v1].point;

      calcDistanceAndSetForces(v0,v1,p0,p1);
    }
  }
}

The function iterates over the vertices living in the red cell, and for each possible couple of vertices, it runs the function calcDistanceAndSetForces()

function calcDistanceAndSetForces(v0,v1,p0,p1){

      //compute the squared distance between the two vertices
      diff = sub(p0, p1);
      dist2 = dot(diff, diff);

      //if the squared distance is greater than the squared radius, 
      //skip this vertex
      if(dist2 > radius2) return; 

      //compute the distance
      dist = Math.sqrt(dist2);

      //compute the force direction
      dir = div1(diff, dist);

      //compute the force: force_dir / (p0p1_dist^2)
      force = div1(dir, dist2+0.001);

      //add this force to the forces pushing v0
      forces[v0] = sum(forces[v0], force);

      //add an equal and opposite force to v1
      forces[v1] = sum(forces[v1], mul1(force, -1));

}

The function calcDistanceAndSetForces(), as the name suggests, computes the distance between the vertices $v0$ and $v1$. If such distance is less than radius, it computes the amount of repulsive force to apply on $v0$ and $v1$.

The force is obtained by normalizing the difference between the position of the two vertices and dividing it by the squared distance separating them. This makes the force stronger when the vertices are close and weaker when they're far apart. In the code, a small offset of $0.001$ is added to the squared distance to avoid numerical issues, as we could potentially divide by a value very close to $0$.

Remember that every action has an equal and opposite reaction; in other words, if $v0$ pushes $v1$, $v1$ pushes $v0$ in the opposite direction. For this reason, the force is applied both on $v0$ and $v1$ equally, but with flipped signs.

One may ask why we still have to test the distance between two vertices if we managed to put them in a grid of $cellsize = radius$; the distance test is still relevant because the grid cells are cubes of side radius, so potentially two vertices in the same cell could be separated by a distance greater than radius, as the diagonal of a cell has $length = radius*\sqrt{3}$.

Note: the distance check is performed on the squared distance against the squared radius; $dist > radius$ is mathematically equivalent to $dist^2 > radius^2$, but in case the two vertices are too far apart to interact, the costly sqrt() function is skipped. That's a minor optimization.

Then, we perform a similar computation between the vertices populating the red cell and the neighboring cells (the ones colored yellow).

function calcRepulsionDifferentCell(gx, gy, gz, gtx, gty, gtz){

  //declare the variables
  let v0,v1,p0,p1,diff,dist,dist2,dir,force,k,j;

  //for each vertex of the grid cell of coordinates [gx][gy][gz]
  for(k = grid[gx][gy][gz].vertices.length-1; k >= 0; k--){

    //get the vertex index
    v0 = grid[gx][gy][gz].vertices[k];

    //get the vertex position
    p0 = geom.vertices[v0].point;

    //for each vertex of the neighboring cell of coordinates [gtx][gty][gtz]
    for(let j = grid[gtx][gty][gtz].vertices.length-1; j >= 0; j--){

      //get the vertex index
      v1 = grid[gtx][gty][gtz].vertices[j];

      //get the vertex position
      p1 = geom.vertices[v1].point;

      calcDistanceAndSetForces(v0,v1,p0,p1);
    }
  }
}

The function calcRepulsionDifferentCell() iterates over the vertices in the red cell and couples them with all the vertices in a yellow cell. It then performs the distance check and the repulsive force computation just as above.

Once we managed to successfully compute all the repulsive forces, we can move on to the next force to compute.

4) Spring force

function calcSpringForces(){

  //compute the length of half split distance
  let halfSplit_len = split_len*0.5;

  //for each edge in the geometry
  for(let i = geom.edges_size - 1; i >= 0; i--){

    //get the first halfedge of this edge
    let he0 = geom.edges[i].halfedges[0];

    //get the endpoints of the edge
    let v0 = geom.halfedges[he0].from;
    let v1 = geom.halfedges[he0].to;
    let p0 = geom.vertices[v0].point;
    let p1 = geom.vertices[v1].point;

    //compute their distance
    let diff = sub(p0, p1);
    let dist2 = dot(diff, diff);
    let edge_len = Math.sqrt(dist2);

    //compute the direction of this edge
    let edge_dir = div1(diff, edge_len);

    //compute the force: edge_dir / (dist^2)
    let force = div1(edge_dir, dist2+0.001);

    //if the length of the edge is greater than half the split length, make the force attractive
    if(edge_len > halfSplit_len) force = mul1(force, -1);
    
    //apply the force to v0
    forces[v0] = sum( forces[v0], force );

    //and an equal and opposite force to v1
    forces[v1] = sum( forces[v1], mul1(force, -1) );
  }
}

The spring force is way easier to manage than the repulsive force, as there's no acceleration structure involved in the computation. The calcSpingForces() function iterates over the edges of the Jitter geometry and computes the spring force between the endpoints of the edge.

The computation of the force is exactly like the repulsive force with the only difference that when $edgelen > splitlen/2$, the sign of the force vector is flipped, turning the force from repulsive to attractive.

    //if the length of the edge is greater than half the split length, make the force attractive
    if(edge_len > halfSplit_len) force = mul1(force, -1);

5) Applying forces

Finally, we computed all the forces to which each vertex is subjected. It's now time to use them to move the vertices. The way forces are applied to the vertices is effortless:

function applyForces(){

  //for each vertex
  for(let i = geom.vertices_size - 1; i >= 0; i--){

    //sum to the current vertex position the force, dampened by a fixed 
    //scaling factor (0.00003) and the param growth_speed
    geom.vertices[i].point = sum(geom.vertices[i].point, mul1(forces[i], 0.00003*growth_speed));

    //multiply this force by the inertia factor
    forces[i] = mul1(forces[i], inertia);
  }
}

The function applyForces() iterates over the vertices and sums to their positions the corresponding force, scaled by a small arbitrary value and the growth_speed attribute.

Since our final goal is to simulate a natural behavior, we should consider that bodies with mass tend to conserve their state of motion. For this reason, the forces we just computed aren't discarded at the next growth iteration, but they're used as the initial forces[] state for the next iteration. This creates an inertia of the growing movement. The last line of the function scales the forces by inertia, preparing the forces[] array for the next growth iteration.

Once the forces are all applied and the vertices have moved, we can pass to the next step of the algorithm: edge splitting.

6) Edge splitting

When the length of an edge is greater than a user-defined split length, such edge is split in half, and a new vertex is inserted into the mesh at the midpoint of the edge.

function insertNewVertices(){

  //For each edge of the mesh
  for(let i = geom.edges_size - 1; i >= 0; i--){

    //het the first halfedge of this edge
    let he0 = geom.edges[i].halfedges[0];

    //get the edges endpoints
    let v0 = geom.halfedges[he0].from;
    let v1 = geom.halfedges[he0].to;
    let p0 = geom.vertices[v0].point;
    let p1 = geom.vertices[v1].point;

    //calc the edge length
    let diff = sub(p0, p1);
    let edge_len2 = dot(diff, diff);

    //if the edge's length is greater than split length, insert a new vertex on this edge
    if(edge_len2 > split_len2) pushVertex(i);

  }
}

The function insertNewVertices() iterates over the edges and computes the edge's length. When it's greater than split length, the function pushVertex() is called:

function pushVertex(e0){

  //current state
  let h0 = geom.edges[e0].halfedges[0]
  let h1 = geom.halfedges[h0].next;
  let h2 = geom.halfedges[h1].next;
  let h3 = geom.edges[e0].halfedges[1];
  let h4 = geom.halfedges[h3].next;
  let h5 = geom.halfedges[h4].next;

  let f0 = geom.halfedges[h0].face;
  let f1 = geom.halfedges[h3].face;

  let v0 = geom.halfedges[h0].from;
  let v1 = geom.halfedges[h3].from;
  let v2 = geom.halfedges[h2].from;
  let v3 = geom.halfedges[h5].from;

  //compute the position of the new vertex
  let p0 = geom.vertices[v0].point;
  let p1 = geom.vertices[v1].point;
  let pm = mul1(sum(p0, p1), 0.5);

  //push the new vertex
  let vm = geom.vertices_size++;
  geom.vertices.push({  point: pm,
              halfedges: []
            });

  //add a new element to the force array that's the sum 
  //of the forces of the surrounding vertices
  forces.push([0,0,0]);
  forces[vm] = sum(forces[vm], forces[v0]);
  forces[vm] = sum(forces[vm], forces[v1]);
  forces[vm] = mul1(forces[vm], 0.5);

  
  //push new faces
  let f2 = geom.faces_size++;
  geom.faces.push({halfedge: h5});

  let f3 = geom.faces_size++;
  geom.faces.push({halfedge: h1});

  //update existing faces
  geom.faces[f0].halfedge = h2;
  geom.faces[f1].halfedge = h4;

  //add new halfedges
  let h6 = geom.halfedges_size++;
  geom.halfedges.push({   from: vm,
              to: v2,
              next: h2,
              prev: h0,
              face: f0,
              opposite: h6+1 /*h7*/
            });

  let h7 = geom.halfedges_size++;
  geom.halfedges.push({   from: v2,
              to: vm,
              next: h7+1 /*h8*/,
              prev: h1,
              face: f3,
              opposite: h6
            });

  let h8 = geom.halfedges_size++;
  geom.halfedges.push({   from: vm,
              to: v1,
              next: h1,
              prev: h7,
              face: f3,
              opposite: h3
            });

  let h9 = geom.halfedges_size++;
  geom.halfedges.push({   from: vm,
              to: v3,
              next: h5,
              prev: h4,
              face: f2,
              opposite: h9+1 /*h10*/
            });

  let h10 = geom.halfedges_size++;
  geom.halfedges.push({   from: v3,
              to: vm,
              next: h10+1 /*h11*/,
              prev: h4,
              face: f1,
              opposite: h9
            });

  let h11 = geom.halfedges_size++;
  geom.halfedges.push({   from: vm,
              to: v0,
              next: h4,
              prev: h10,
              face: f1,
              opposite: h0
            });

  //update existing halfedges
  geom.halfedges[h0].from = v0;
  geom.halfedges[h0].to = vm;
  geom.halfedges[h0].next = h6;
  geom.halfedges[h0].prev = h2;
  geom.halfedges[h0].face = f0;
  geom.halfedges[h0].opposite = h11;

  geom.halfedges[h1].from = v1;
  geom.halfedges[h1].to = v2;
  geom.halfedges[h1].next = h7;
  geom.halfedges[h1].prev = h8;
  geom.halfedges[h1].face = f3;

  geom.halfedges[h2].from = v2;
  geom.halfedges[h2].to = v0;
  geom.halfedges[h2].next = h0;
  geom.halfedges[h2].prev = h6;
  geom.halfedges[h2].face = f0;

  geom.halfedges[h3].from = v1;
  geom.halfedges[h3].to = vm;
  geom.halfedges[h3].next = h9;
  geom.halfedges[h3].prev = h5;
  geom.halfedges[h3].face = f2;
  geom.halfedges[h3].opposite = h8;

  geom.halfedges[h4].from = v0;
  geom.halfedges[h4].to = v3;
  geom.halfedges[h4].next = h10;
  geom.halfedges[h4].prev = h11;
  geom.halfedges[h4].face = f1;

  geom.halfedges[h5].from = v3;
  geom.halfedges[h5].to = v1;
  geom.halfedges[h5].next = h3;
  geom.halfedges[h5].prev = h9;
  geom.halfedges[h5].face = f2;

  //add the new edges
  let e1 = geom.edges_size++;
  geom.edges.push({ halfedges: [h8, h3] });

  let e2 = geom.edges_size++;
  geom.edges.push({ halfedges: [h6, h7] });

  let e3 = geom.edges_size++;
  geom.edges.push({ halfedges: [h9, h10] });

  //update edges
  geom.edges[e0].halfedges = [h0, h11];

  //update valences
  geom.vertices[v2].halfedges.push(h7);
  geom.vertices[v3].halfedges.push(h10);
  geom.vertices[vm].halfedges.push(h6);
  geom.vertices[vm].halfedges.push(h11);
  geom.vertices[vm].halfedges.push(h9);
  geom.vertices[vm].halfedges.push(h8);
  
}

Now, the function pushVertex() may seem pretty long and scary, but in reality, its job is relatively simple: if we want to add a new vertex to the mesh, we must reconstruct a coherent Jitter geometry around the new vertex. This means updating all the parts of the Jitter geometry structure that are connected or somehow related to the new vertex. Take a look at the following image:

The image shows what is happening: starting from the two triangles on the left, we introduce a new vertex $vm$ at the middle of edge $e0$. To do so, we have to add new elements to the Jitter geometry (vertex: $vm$, faces: $f2$, $f3$, edges: $e1$, $e2$, $e3$, halfedges: $h6$, $h7$, $h8$, $h9$, $h10$, $h11$) and update the rest of the structure (vertex: $v0$, $v1$, $v2$, $v3$, faces: $f0$, $f1$, edge $e0$, halfedges: $h0$, $h1$, $h2$, $h3$, $h4$, $h5$).

Moreover, we must add a new element to the forces[] array for the newborn vertex $vm$. To do it, we compute the average force of the neighboring vertices $v0$ and $v1$, which are the endpoints of the original edge $e0$:

  //add a new element to the force array that's the sum 
  //of the forces of the surrounding vertices
  forces.push([0,0,0]);
  forces[vm] = sum(forces[vm], forces[v0]);
  forces[vm] = sum(forces[vm], forces[v1]);
  forces[vm] = mul1(forces[vm], 0.5);

7) Remeshing

The mesh can finally grow, and we can add new vertices where needed. The problem is that such new vertices can potentially pop up anywhere in the mesh and ruin the mesh quality. What, mesh quality? By "quality" i don't mean any aesthetic consideration, nor i'm referring to the amount of details that a mesh can have; The mesh quality is something that can be measured and responds to precise criteria. A mesh is considered "good" when it exhibits the following features:

• equal-length edges
• equilateral triangles
• valence close to 6

Since our growth algorithm uses a sort of spring force to keep the vertices at the same distance from each other, we know for sure that the first criterion is met; the same applies to the second criterion as equal-length edges form equilateral triangles. Regarding the last, we don't have any assurance.

But first of all, what is the valence? It is simply the number of edges connected to a vertex. A "sane" mesh has a valence of 6 for the internal vertices and 4 for vertices at the edge. The image below shows a "bad" and a "good" mesh, reporting the valence of each vertex:

What can we do then to improve the mesh quality? We can perform an edge flip where needed.

Edge flipping consists of taking two adjacent triangles and flipping the internal edge by swapping the vertices it connects. This can improve or worsen the mesh mean valence, so we don't perform it on every edge but only where it's worth it.

function remesh(){

  //if remesh is disabled, return
  if(remesh_enable == 0) return; 

  //declare the variables
  let flipValence, noFlipValence, tv, h0,h3,h1,h4,v0,v1,v2,v3;

  //for each edge of the mesh
  for(let e0 = geom.edges_size - 1; e0 >= 0; e0--){

    //get the halfedges of this edge
    h0 = geom.edges[e0].halfedges[0];
    h3 = geom.edges[e0].halfedges[1];
    h1 = geom.halfedges[h0].next;
    h4 = geom.halfedges[h3].next;

    //get the endpoints of this edge and the opposite vertices
    v0 = geom.halfedges[h0].from;
    v1 = geom.halfedges[h4].to;
    v2 = geom.halfedges[h3].from;
    v3 = geom.halfedges[h1].to;

    //calc if it's worth to flip this edge
    tv = geom.vertices[v0].halfedges.length - 6; 
    noFlipValence = tv*tv;
    tv = geom.vertices[v1].halfedges.length - 6; 
    noFlipValence += tv*tv;
    tv = geom.vertices[v2].halfedges.length - 6;
    noFlipValence += tv*tv;
    tv = geom.vertices[v3].halfedges.length - 6;
    noFlipValence += tv*tv;

    tv = geom.vertices[v0].halfedges.length - 7;
    flipValence = tv*tv;
    tv = geom.vertices[v1].halfedges.length - 5;
    flipValence += tv*tv;
    tv = geom.vertices[v2].halfedges.length - 7;
    flipValence += tv*tv;
    tv = geom.vertices[v3].halfedges.length - 5;
    flipValence += tv*tv;

    //if the valence improves (gets closer to 6)
    if(flipValence < noFlipValence){

      //get the rest of the geometry
      let h2 = geom.halfedges[h1].next;
      let h5 = geom.halfedges[h4].next;
      let f0 = geom.halfedges[h0].face;
      let f1 = geom.halfedges[h3].face;

      //and flip this edge
      flipEdge(e0, h0, h1, h2, h3, h4, h5, v0, v1, v2, v3, f0, f1);
    } 
  }
}

The function remesh() iterates over the edges of the mesh and calculates if it's worth to flip this edge or not. The benefit of flipping this edge is computed comparing the squared mean valence of the no-flip scenario against the flip route.

    //calc if it's worth to flip this edge
    tv = geom.vertices[v0].halfedges.length - 6; 
    noFlipValence = tv*tv;
    tv = geom.vertices[v1].halfedges.length - 6; 
    noFlipValence += tv*tv;
    tv = geom.vertices[v2].halfedges.length - 6;
    noFlipValence += tv*tv;
    tv = geom.vertices[v3].halfedges.length - 6;
    noFlipValence += tv*tv;

    tv = geom.vertices[v0].halfedges.length - 7;
    flipValence = tv*tv;
    tv = geom.vertices[v1].halfedges.length - 5;
    flipValence += tv*tv;
    tv = geom.vertices[v2].halfedges.length - 7;
    flipValence += tv*tv;
    tv = geom.vertices[v3].halfedges.length - 5;
    flipValence += tv*tv;

if the mean valence improves, then flip this edge:

    //if the valence improves (gets closer to 6)
    if(flipValence < noFlipValence){

      //get the rest of the geometry
      let h2 = geom.halfedges[h1].next;
      let h5 = geom.halfedges[h4].next;
      let f0 = geom.halfedges[h0].face;
      let f1 = geom.halfedges[h3].face;

      //and flip this edge
      flipEdge(e0, h0, h1, h2, h3, h4, h5, v0, v1, v2, v3, f0, f1);
    } 

Like we did for pushing a new vertex into the mesh, we can flip an edge by changing the elements of the Jitter geometry structure connected to that edge:

function flipEdge(e0, h0, h1, h2, h3, h4, h5, v0, v1, v2, v3, f0, f1){

  //flip the edge by updating the halfedge structure
  geom.halfedges[h0].from = v1; 
  geom.halfedges[h0].to = v3; 
  geom.halfedges[h0].next = h2; 
  geom.halfedges[h0].prev = h4; 
  geom.halfedges[h0].face = f0;

  geom.halfedges[h1].from = v2; 
  geom.halfedges[h1].to = v3; 
  geom.halfedges[h1].next = h3; 
  geom.halfedges[h1].prev = h5; 
  geom.halfedges[h1].face = f1;

  geom.halfedges[h2].from = v3; 
  geom.halfedges[h2].to = v0; 
  geom.halfedges[h2].next = h4; 
  geom.halfedges[h2].prev = h0; 
  geom.halfedges[h2].face = f0;

  geom.halfedges[h3].from = v3; 
  geom.halfedges[h3].to = v1; 
  geom.halfedges[h3].next = h5; 
  geom.halfedges[h3].prev = h1; 
  geom.halfedges[h3].face = f1;

  geom.halfedges[h4].from = v0; 
  geom.halfedges[h4].to = v1; 
  geom.halfedges[h4].next = h0; 
  geom.halfedges[h4].prev = h2; 
  geom.halfedges[h4].face = f0;

  geom.halfedges[h5].from = v1; 
  geom.halfedges[h5].to = v2; 
  geom.halfedges[h5].next = h1; 
  geom.halfedges[h5].prev = h3; 
  geom.halfedges[h5].face = f1;

  geom.faces[f0].halfedge = h0;
  geom.faces[f1].halfedge = h3; 

  //increase the valence of v1 and v3
  geom.vertices[v1].halfedges.push(h0);
  geom.vertices[v3].halfedges.push(h3);

  //decrease the valence of v0 and v2
  geom.vertices[v0].halfedges.splice(geom.vertices[v0].halfedges.indexOf(h0), 1);
  geom.vertices[v2].halfedges.splice(geom.vertices[v2].halfedges.indexOf(h3), 1);

}

The function flipEdge() re-assigns the elements of the Jitter geometry of the two triangles divided by the edge to flip:

this is the difference with and without remeshing:

• without remeshing

• with remeshing

The growth pattern looks completely different. I decided to leave the possibility to enable and disable remeshing because, although not ideal for the quality of the mesh, the shapes that the mesh could assume without remeshing may be interesting and desirable.

There's still something we can do to improve the structural quality of the mesh; under the pressure of the forces, the vertices are pushed here and there, creating an uneven mesh surface. While this may be desirable, i thought to include a controllable smoothing function to create a sort of "surface tension":

function smooth(){

  //if smoothing is disabled, quit
  if(smooth_iterations <= 0) return;
  let newPos;

  //compute mixing factors
  let filterA = smooth_amt;
  let filterB = 1 - filterA;

  //compute "smooth_iterations" of smoothing
  for(let s = 0; s < smooth_iterations; s++){
    
    //prepare an array for the smoothed positions
    newPos = new Array(geom.vertices_size);
    
    //for each vertex
    for(let i = geom.vertices_size - 1; i >= 0; i--){

      //initialize the mean position to [0,0,0]
      let meanPos = [0,0,0];

      //for each halfedge starting from this edge
      for(let n = geom.vertices[i].halfedges.length - 1; n >= 0; n--){

        //get the position of this neighboring vertex
        let h = geom.vertices[i].halfedges[n];
        let v = geom.halfedges[h].to;
        let p = geom.vertices[v].point;

        //add its position to the running mean of the positions
        meanPos = sum(meanPos, p);
      }

      //compute the mean position dividing by the number of neighboring vertices
      newPos[i] = div1(meanPos, geom.vertices[i].halfedges.length);
    }

    //for each vertex
    for(let i = geom.vertices_size - 1; i >= 0; i--){

      //update the position of this vertex blending the smoothed and non-smoothed position
      geom.vertices[i].point = mul1(geom.vertices[i].point, filterB);
      geom.vertices[i].point = sum(geom.vertices[i].point, mul1(newPos[i], filterA));
    }   
  }
}

The function smooth() iterates over the vertices of the mesh, and, for each vertex, it accesses the neighboring vertices, which are the ones connected to this vertex by an edge. The function then computes the mean position of the neighboring vertices and updates the vertex positions by blending the current vertex position with the mean position of the neighboring vertices. The smoothing is controlled by the smooth_iterations attribute, which decides how many iterations of the smoothing algorithm to perform, and by the attribute smooth_amt, which determines the blending factor between the current vertex position and the mean position of the surrounding vertices.

The smoothing operation seems to improve the cohesion of the mesh surface.

Next, we just need to triangulate the vertices and output the mesh.

8) Output the mesh

The last stage of our algorithm consists of triangulating the vertices and outputting the geometry as a Jitter matrix.

function draw(){

  //set the dimension of the output matrix
  mOut.dim = geom.faces_size*3;

  //declare the variables
  let he, v, p;
  let points = new Float32Array(geom.faces_size*9);
  let count = 0;

  //for each face, get the vertices of this face, and push them into the array
  for(let i = geom.faces_size - 1; i >= 0; i--){

    he = geom.faces[i].halfedge;
    v = geom.halfedges[he].from;
    p = geom.vertices[v].point;
    points[count++] = p[0];
    points[count++] = p[1];
    points[count++] = p[2];

    he = geom.halfedges[he].next;
    v = geom.halfedges[he].from;
    p = geom.vertices[v].point;
    points[count++] = p[0];
    points[count++] = p[1];
    points[count++] = p[2];

    v = geom.halfedges[he].to;
    p = geom.vertices[v].point;
    points[count++] = p[0];
    points[count++] = p[1];
    points[count++] = p[2];
  }

  //cooy the array into the matrix
  mOut.copyarraytomatrix(points);

  //output the matrix
  outlet(0, "jit_matrix", mOut.name);
}

The function draw() iterates over the faces of the mesh and triangulates their vertices by pushing their positions into an array in the right order. The array gets then copied into a Jitter matrix at the end of the for loop, and the matrix is then sent out from the script.

Final thoughts

This concludes the tutorial about implementing a differential growth algorithm using Jitter geometry and JavaScript. Remember that you could take the result of this growth process and turn it into a Jitter geometry again for further manipulation, like subdivision, normals generation, and to create texture coordinates. This is not the only growth algorithm that can be implemented using these tools, and that's not the only way to implement this specific one.

I give you some quick inputs to hopefully tickle your curiosity:

• what about generating new random vertices instead of waiting for an edge to overcome the split_length?
• what if radius, split_len, and halfsplit_len were independent from each other?
• what if the growth_speed wasn't constant but more prominent for the vertices pointing up?
• what about introducing a flow field made with jit.bfg to affect the growth speed and direction?
• What about adding a surrounding mesh as a container and computing its repulsive force? Will the grown mesh take the shape of the container?
• what about introducing an age parameter? At each growth iteration, vertices get older and less reactive, but the newborn vertices start at age 0

I hope you can use this patch and this script to explore the fascinating world of growth algorithms. Check out also other growth strategies, such as Lindenmayer Systems, Fractal-Based Growth, Space Colonization Algorithms, Diffusion-Limited Aggregation, and Reaction-Diffusion Systems.

Thanks much for following me to this point, and have a tremendous patching day!

MM