3d von Koch by Chris M. Thomasson

I start with the following vertexs for a triangular regular pyramid, in the unit cube.

glm::vec3 t[] = {
   glm::vec3(-1.0, -1.0, 1.0),
   glm::vec3(1.0, -1.0, 1.0),
   glm::vec3(0.0, -1.0, -1.0),
   glm::vec3(0.0, 1.0, 0.0)
};


I then index into these to create four triangles.

One for each side.

g.push_points3(g.m_vindex, 0, 1, 3);
g.push_points3(g.m_vindex, 1, 2, 3);
g.push_points3(g.m_vindex, 2, 0, 3);
g.push_points3(g.m_vindex, 0, 1, 2);

This creates a tetrahedron. Keep in mind that this scales to any triangle, so a cube would work with two triangles per face.


I then iterate each triangle in the tetrahedron:

g.iterate_triangle(0, imax, length, t[0], t[1], t[3]);
g.iterate_triangle(0, imax, length, t[1], t[2], t[3]);
g.iterate_triangle(0, imax, length, t[2], t[0], t[3]);
g.iterate_triangle(0, imax, length, t[1], t[0], t[2]);


The code for the recursive function (iterate_triangle) is:

_____________________________________
glm::vec3 get_scale_point(glm::vec3 p0, glm::vec3 p1, double scale)
{
  glm::vec3 diff = p1 - p0;

  return glm::vec3(
    p0.x + diff.x * scale,
    p0.y + diff.y * scale,
    p0.z + diff.z * scale
  );
}

glm::vec3
get_tip_from_triangle(
  glm::vec3 p0, // bottom left
  glm::vec3 p1, // bottom right
  glm::vec3 p2,  // top,
  glm::vec3 topb, // bottom of the new top,
  double length
){
  glm::vec3 tcross = glm::cross(p1 - p0, p2 - p0);

  double tlen = glm::length(tcross);
  double s = length / tlen;

  glm::vec3 ret(
    topb.x + tcross.x * s,
    topb.y + tcross.y * s,
    topb.z + tcross.z * s
  );

  double d = glm::distance(ret, topb);

  std::printf("length: %lf - %lf\r\n", length, d);

  return ret;
}


/* iterate a triangle */
void iterate_triangle(
  unsigned int irecur,
  unsigned int irecur_max,
  double length,
  glm::vec3 p0, // bottom left
  glm::vec3 p1, // top bottom right
  glm::vec3 p2  // top
){

  if (irecur >= irecur_max) return;

  unsigned int vbase = (m_verts.size() / 3);


  glm::vec3 bmid = get_scale_point(p0, p1, 0.5);
  glm::vec3 topb = get_scale_point(bmid, p2, 0.5);

  double scale = 1.0 / 2.0;
  glm::vec3 top_left = get_scale_point(p0, p2, scale);
  glm::vec3 top_right = get_scale_point(p1, p2, scale);
  glm::vec3 bottom_mid = bmid;
  glm::vec3 ttop = get_tip_from_triangle(p0, p1, p2, topb, length);


  push_points3(m_verts, top_left);
  push_points4(m_colors, glm::vec4(1.0, 0.0, 0.0, 1.0));

  push_points3(m_verts, top_right);
  push_points4(m_colors, glm::vec4(0.0, 1.0, 0.0, 1.0));

  push_points3(m_verts, bottom_mid);
  push_points4(m_colors, glm::vec4(0.0, 0.0, 1.0, 1.0));

  push_points3(m_verts, ttop);
  push_points4(m_colors, glm::vec4(1.0, 1.0, 0.0, 1.0));

  push_points3(m_vindex, vbase + 0, vbase + 1, vbase + 3);
  push_points3(m_vindex, vbase + 0, vbase + 2, vbase + 3);
  push_points3(m_vindex, vbase + 2, vbase + 1, vbase + 3);



  double blength = length * 0.33333;
  double tlength = length * 0.3333333;

  iterate_triangle(irecur + 1, irecur_max, blength, top_right, top_left, ttop);
  iterate_triangle(irecur + 1, irecur_max, blength, top_left, bottom_mid, ttop);
  iterate_triangle(irecur + 1, irecur_max, blength, bottom_mid, top_right, ttop);

  iterate_triangle(irecur + 1, irecur_max, tlength, top_left, top_right, p2);
  iterate_triangle(irecur + 1, irecur_max, tlength, p0, bmid, top_left);
  iterate_triangle(irecur + 1, irecur_max, tlength, bmid, p1, top_right);
}
_____________________________________