Castle Story looks like this due to technical constraints: Were there to be a heightmap per each voxel in the entire volume, rather than only a heightmap per each surface voxel, storage cost would be vastly greater, on the order of O(n^3) which can be prohibitive, as opposed to a more favourable O(n^2), where n is the side length of a cubic voxel space representing your world. Bear in mind that heightmap information for subterranean voxels is implicit in the grid structure, so heightmap info need only be explicitly stored for those voxels which lie on the surface. So the Castle Story guys have snapped vertices to grid interstices to save on storage cost and mesh-construction complexity... read on.
First off, let's look at your options:
(1) is going to make things difficult for you, since in a single voxel column you only want mini-heightmap information for the topmost voxel -- refer back to the first paragraph for the reasons why. Looking at (1), the column on the right has heightmap information for both the topmost and the second-from-top (and this could apply to the third-from-top and so on, if the slope were extreme enough.). This is no good.
(2) Is probably a better option then; that is, ensuring that corner vertices snap to voxel grid interstices. But how then do we address the extreme slope issue? Well, we need to pick some gradient where we simply snap to a vertical column and thus ensure we have no gradients than can crosscut through n topmost voxels. A 45 degree gradient is the natural cutoff for reasons I explain below. So instead of (3), we'd have (4):
(4) Snapping the voxel's corner vertex to the nearest grid interstice is the solution. The visual effect -- diagonal aliasing -- can be seen in your screenshot of Castle Story. The cutoff slope for the voxels is a 1:1 gradient, or 45 degrees (as viewed orthogonally). Working backward from the solution, let's look at the reasons:
- The only way we could have an unbroken, extreme slope, is if a voxel were vertically elongated...
- .... But no voxel's mesh can exceed it's bounding area, irrespective of it's actual smoothed shape; a voxel needs to sit in a defined space in a 3D grid which accounts, if not for it's exact form, then at least for it's axis-aligned bounding box.
Another reason to approach it this way is that, as you've already discovered, not using snapping (discretization) leads to a range of geometrically complex surface smoothing scenarios, which are best avoided altogether... games of this sort generally do not require the degree of accuracy that a proper CSG algorithm would provide, which is the whole reason we are using voxels in the first place: Voxels make it far easier to work incrementally with volumes, than floating-point (continuous) polygon intersection / CSG algorithms do.