The good news: your Marching Cubes algorithm looks just fine to me! That 3d surface reconstruction looks gorgeous. If you're committed to a voxel-based approach with isosurface visualization, you're off to a fine start.
The problem is that 3d noise in this form really isn't suitable for use with a Marching Cubes-type algorithm for terrain. If you want to be able to do the sort of terrain modification that the demo you point to does, what I would suggest is to build the heightmap first and then build your volumetric data based on the heightmap.
Building the heightmap in the first place is relatively straightforward; I'd just use one of the classic fractal methods. You can use an approach like square-diamond midpoint displacement (probably the most common), but I would encourage a Fourier-based approach instead because it should allow you a lot more control over how 'rolling' your terrain is (and you can always 'postfilter' by applying a gamma-type effect, e.g. setting h(x,y) -> c*h(x,y)^3 or h(x,y) -> c*sqrt(h(x,y)), to make the terrain more or less jagged - if you look at the early articles on fractal terrain you'll see a lot of this).
Once you've got your heightmap data, there are a bunch of different procedures for converting that information into volume data. For simplicity's sake I'll assume that your heightmap data has been stored at higher resolution than your voxel data has; this makes a lot of sense, since a 3d array with the same XY resolution would (obviously) have to be a lot larger than the 2d heightmap. The easiest way to do it is as simple as
for (0 < vx < VOXEL_X_RES) {
for (0 < vy < VOXEL_Y_RES) {
pick the closest point (x,y) on the heightmap corresponding to (vx,vy)
(e.g. x = HEIGHTMAP_X_RES*vx/VOXEL_X_RES, etc)
find the z value of the heightmap at (x,y) (and convert it to a voxel-space value vz)
fill every value below (vx, vy, vz) with 1 to indicate it's in the terrain
fill every value above (vx, vy, vz) with 0
}
}
Because this fills the voxel array with only 0s and 1s, it can mean that the marching cubes interpolation along edges looks a little too 'regular' and results in surfaces with too many identical slopes on them. Instead, you may want to 'subsample' the heightmap to build your voxel array, using something like:
for (0 < vx < VOXEL_X_RES) {
for (0 < vy < VOXEL_Y_RES) {
find all 'pixels' on the heightmap corresponding to (vx,vy)
for ( 0 < vz < VOXEL_Z_RES) {
compute how many of the pixels in the set above have z > vz
set the value of the voxel array at (vx, vy, vz) to the proportion of pixels with z > vz
}
}
}
The set of pixels on the heightmap for this second version can be anything from a straight subsampled grid (i.e., find hxmin = HEIGHTMAP_X_RES*vx/VOXEL_X_RES
and hxmax = HEIGHTMAP_X_RES*(vx+1)/VOXEL_X_RES
, similarly for hymin
and hymax
, and consider all pixels (hx, hy)
with (hxmin <= hx < hxmax)
and (hymin <= hy < hymax)
) to an 'overlapping samples' approach where you look at every heightmap-pixel within a circle centered around wherever (vx, vy) maps to on the heightmap; you could even subsample your heightmap and interpolate between heightmap pixels for finer detail.
One major caveat with this whole approach is that fractal terrain isn't really very 'terrain-aware' - it knows height, but it doesn't know anything about features : it doesn't really know what a river is or the effect it has on the surrounding geography; it can't distinguish easily between 'old' and 'new' geography; etc. If you want truly realistic terrain then you should consider some of the more simulationist approaches - but these tend to directly conflict with the kind of terrain modifications that the linked video shows, so if realism is your goal then it may be worth rethinking the voxel approach entirely.