You may want to check out the previous post about water. It ended up producing a 2D map of booleans, where true meant there was water in that location of the map. That was great, but we still had a nasty problem ahead. We needed to level the lakes.
It will be easier if we look at some pictures. Imagine this is the tile we need to compute. Initially it is just the terrain height and the sea map (in this case using a sea-level value):
And here you can see the lake and river location phase from the earlier post in pretty pink:
As you can see in the zoomed portion, these lakes often span over very irregular terrain. We knew the shoreline points for the same lake were at similar height. Figuring out the water level from A to B was the problem. We know A-B should be as flat as possible and the same applies to A and any other point in the shoreline. That's a lot of potential points.
If we had all the time in the world, we could try discovering individual lakes in one phase and then do some sort of iterative leveling by relaxation and or filtering. We had only a fraction of a second for this so we knew it had to be hacked.
We tried many different things. (When I say "we", it is not modesty plural. Michael Coates who recently joined Voxel Farm got to implement these desperate things we tried.) Most had to do with rasterization and interpolation, first from shore to shore, then we got looking into water height derivatives. Nothing really worked.
As a last recourse attempt, we tried an approach that involved a Quadtree. The idea was to encode the water bodies in a Quadtree and apply the same leveling algorithms we knew would be too slow if we were operating in individual water samples.
The shorelines produced a lot of small quads, but as we moved inside the lake we started getting larger quads.
We could see that the larger a patch, the more "detached" it was from the shore. That we liked. We could raise or sink those points much more to achieve the level we needed.
Once the water height of the larger patches is decided, we rasterize them into the water heightmap. This is a simple, fast bi-linear interpolation over the quad. Then we do the same for the smaller quads. There is a reason why you need to process all quads of a given size before processing the next level of smaller patches: You need to make sure the small quads match the interpolation values along the edges they share with a larger quad.
This overall process is quite fast. Here is a rendering of the new found water volumes. It is a bit exaggerated, but hopefully you get the idea:
We still need to work out on some of the heuristics, but it looks we are pretty much on our way.
If there is any lesson to be drawn from this, it would be when in despair, apply a Quadtree and everything will be alright.