Interpolation on marching cubes algorithm
I worked it out on my own. I'm trying to write it down correctly so it's readable for newcomers.

Introduction

First of all, interpolation itself requires that there are float values and not binary values. This gives the ability to create different edge-crossing-vertices. How it's going I'm explaining now.

The interpolation which is needed here is linear interpolation.
For understanding purposes I'll explain it here too.

General Interpolation

After all, what is interpolation? Well good question.
Interpolation is described as name for algorithms which smooth a given input of data e.g. points (with n-dimensions). Simpler there are points searched between the points for a given interpolation value. First it's required to give the points a base so that's clear where points are searched for in between. After that a specific algorithm for the interpolation has to be chosen. There are different kinds of interpolation algorithm out there each for a specific application.
Best known is linear interpolation which results for a value one point on a line which crosses two input points.

The formula for linear interpolation (2D here) is as following (simply maths):
$$
y = y_{0} + (value - x_{0} ) * \frac{y_{1} - y_{0} }{x_{1} - x_{0}}
$$
with $$ value \in   \mathbb{R} \wedge 0 &lt;= value &lt;= 1 $$

After all with this formula it's possible to apply interpolation on the vertices of the marching cubes algorithm which requires following rules.


The voxels must have a floating value instead of a binary value.
The value of the voxels must be only between zero and one
There must be a boundary value be declared which describes whether a voxel is supposed to be in surface. 


Looking back to the blog the boundary value is called the iso surface value. 
This value describes as explained whether a voxel is being asumed to be in the surface so for acquiring the edge crossing. This is important to know because without it you wouldn't find out where even a edge-crossing computation is needed.

The interpolation is done on each coordinate of the two voxels position where and crossing-edge vertex is.
The implementation is quite easy and looks like this (extracted from here):

/*
   Linearly interpolate the position where an isosurface cuts
   an edge between two vertices, each with their own scalar value
*/
XYZ VertexInterp(isolevel,p1,p2,valp1,valp2)
double isolevel;
XYZ p1,p2;
double valp1,valp2;
{
   double mu;
   XYZ p;

   if (ABS(isolevel-valp1) &lt; 0.00001)
      return(p1);
   if (ABS(isolevel-valp2) &lt; 0.00001)
      return(p2);
   if (ABS(valp1-valp2) &lt; 0.00001)
      return(p1);
   mu = (isolevel - valp1) / (valp2 - valp1);
   p.x = p1.x + mu * (p2.x - p1.x);
   p.y = p1.y + mu * (p2.y - p1.y);
   p.z = p1.z + mu * (p2.z - p1.z);

   return(p);
}


Conclusion

The interpolation of the edge-crossing vertices is a nice feature for smoothing the presentation of voxels.However the interpolation relies on that there are already floating values which means that it cannot smooth your existing binary data. For such cases different algorithms and approaches are needed.