title: How to smooth the blocks of a 3D voxel world?
body: <p>In my (Minecraft-like) 3D voxel world, I want to smooth the shapes for more natural visuals. Let's look at this example in 2D first.</p>
  
  <p><a href="https://i.stack.imgur.com/aG0PU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/aG0PU.png" alt="No smoothing, circular smoothing, and bezier smoothing"></a></p>
  
  <p>Left is how the world looks without any smoothing. The terrain data is binary and each voxel is rendered as a unit size cube.</p>
  
  <p>In the center you can see a naive circular smoothing. It only takes the four directly adjacent blocks into account. It is still not very natural looking. Moreover, I'd like to have flat 45-degree slopes emerge.</p>
  
  <p>On the right you can see a smoothing algorithm I came up with. It takes the eight direct and diagonal neighbors into account in order to come up with the shape of a block. I have <a href="https://github.com/danijar/voxel-smoothing-2d" rel="nofollow noreferrer">the C++ code</a> online. Here is the code that comes up with the control points that the bezier curve is drawn along.</p>
  
  <pre><code>#include &lt;iostream&gt;
  
  using namespace std;
  using namespace glm;
  
  
  list&lt;list&lt;dvec2&gt;&gt; Points::find(ivec2 block)
  {
      // Control points
      list&lt;list&lt;ivec2&gt;&gt; lines;
      list&lt;ivec2&gt; *line = nullptr;
  
      // Fetch blocks, neighbours start top left and count
      // around the center block clock wise
      int center = m_blocks-&gt;get(block);
      int neighs[8];
      for (int i = 0; i &lt; 8; i++) {
          auto coord = blockFromIndex(i);
          neighs[i] = m_blocks-&gt;get(block + coord);
      }
  
      // Iterate over neighbour blocks
      for (int i = 0; i &lt; 8; i++) {
          int current = neighs[i];
          int next = neighs[(i + 1) % 8];
          bool is_side   = (((i + 1) % 2) == 0);
          bool is_corner = (((i + 1) % 2) == 1);
  
          if (line) {
              // Border between air and ground needs a line
              if (current != center) {
                  // Sides are cool, but corners get skipped when they don't
                  // stop a line
                  if (is_side || next == center)
                      line-&gt;push_back(blockFromIndex(i));
              } else if (center || is_side || next == center) {
                  // Stop line since we found an end of the border. Always
                  // stop for ground blocks here, since they connect over
                  // corners so there must be open docking sites
                  line = nullptr;
              }
          } else {
              // Start a new line for the border between air and ground that
              // just appeared. However, corners get skipped if they don't
              // end a line.
              if (current != center) {
                  lines.emplace_back();
                  line = &amp;lines.back();
                  line-&gt;push_back(blockFromIndex(i));
              }
          }
      }
  
      // Merge last line with first if touching. Only close around a differing corner for air
      // blocks.
      if (neighs[7] != center &amp;&amp; (neighs[0] != center || (!center &amp;&amp; neighs[1] != center))) {
          // Skip first corner if enclosed
          if (neighs[0] != center &amp;&amp; neighs[1] != center)
              lines.front().pop_front();
          if (lines.size() == 1) {
              // Close circle
              auto first_point = lines.front().front();
              lines.front().push_back(first_point);
          } else {
              // Insert last line into first one
              lines.front().insert(lines.front().begin(), line-&gt;begin(), line-&gt;end());
              lines.pop_back();
          }
      }
  
      // Discard lines with too few points
      auto i = lines.begin();
      while (i != lines.end()) {
          if (i-&gt;size() &lt; 2)
              lines.erase(i++);
          else
              ++i;
      }
  
      // Convert to concrete points for output
      list&lt;list&lt;dvec2&gt;&gt; points;
      for (auto &amp;line : lines) {
          points.emplace_back();
          for (auto &amp;neighbour : line)
              points.back().push_back(pointTowards(neighbour));
      }
      return points;
  }
  
  glm::ivec2 Points::blockFromIndex(int i)
  {
      // Returns first positive representant, we need this so that the
      // conditions below "wrap around"
      auto modulo = [](int i, int n) { return (i % n + n) % n; };
  
      ivec2 block(0, 0);
      // For two indices, zero is right so skip
      if (modulo(i - 1, 4))
          // The others are either 1 or -1
          block.x = modulo(i - 1, 8) / 4 ? -1 : 1;
      // Other axis is same sequence but shifted
      if (modulo(i - 3, 4))
          block.y = modulo(i - 3, 8) / 4 ? -1 : 1;
      return block;
  }
  
  dvec2 Points::pointTowards(ivec2 neighbour)
  {
      dvec2 point;
      point.x = static_cast&lt;double&gt;(neighbour.x);
      point.y = static_cast&lt;double&gt;(neighbour.y);
  
      // Convert from neighbour space into
      // drawing space of the block
      point *= 0.5;
      point += dvec2(.5);
  
      return point;
  }
  </code></pre>
  
  <p>However, this is still in 2D. How to translate this algorithm into three dimensions?</p>