1. Although datasets are getting larger, so are computers' internal memory. My latest notebook has 8GB of main memory. My 5-year-old lab machine has 64GB of main memory. Therefore internal algorithms, which do not use external storage, are still important.
2. Although processing the largest databases would require external storage, the cutoff beyond which that is true is higher and higher, currently hundreds of millions of pixels.
3. Space and time efficiency is more important than ever. Reduced space means that larger datasets can be processed internally. That means that reduced asymptotic time growth and multiplicative time constants become more important. Consider the size of {$ (\lg N) $} when {$N=10^{9} $}, the number of voxels that my connected components program can easily process.
4. Speed is also important because a fast implementation becomes a feasible subroutine of a larger algorithm. E.g., hydrography might be used in terrain compression.

1. Everything, terrain and hydrography, is discretized on a regular grid, sometimes called the world.
   1. Water flows off the edge of the world.
   2. The grid has a regular spacing either in meters or in angle in lat/long.
   3. For convenience, we call one grid cell or quantum a pixel (since cell can have alternate meanings, such as a complete DEM-1 array of {$1201^2$} pixels).
   4. A particular pixel is denoted as {$p_{ij}$}.
2. Terrain: elevation,
   1. above a reference surface (sea level, the geoid),
   2. expressed as an array of numbers {$z_{ij}$}.
3. Hydrography: in this context, the water flow through each pixel.
   1. I don't think that 4-connected vs 8-connected makes much difference. Since people prefer 8-connected, that's what I now use.
   2. I don't process non-tree hydrography graphs, like islands, or lakes like Lake George NY with two exits.. They would add to the complexity but not, seemingly, require major new ideas.
4. Each pixel gets one inch of rain.
5. All the water entering a pixel flows out to its lowest neighbor. Ties are broken consistently.

The big problem is that there are a surprisingly large number of basins, or sections of terrain that trap their water. The reason is that all it takes to cause a separate basin is a 1-pixel wide 1 unit hi ridge at the bottom:

1. Binary input is surprisingly faster than ASCII input.
2. Reading {$4097^2$} ASCII unsigned short ints took about 5 CPU seconds.
3. Reading them in binary takes only 0.08 CPU seconds.
4. Perhaps there is something wrong with how ASCII input is implemented; GNU has implemented C++ I/O badly in the past.
5. So, for now, read binary.

1. A 1-pixel pit is a pixel that is lower than (or the same elevation as) all of its 8 neighbors.
2. Removing them is done before adjacencies are computed.
3. It is not strictly necessary, but is very fast, and removes about 60% of the basins, which is a big win. Therefore the connected component computation is sufficiently faster that the total time does not increase. Even better, the basin joining step should be much faster.
4. The process is this: If a pixel {$\le$} than the lowest of its 8 neighbors, then raise its elevation to the mean of its two lowest neighbors.
5. That elevation is probably a fraction, so internal elevations are stored as doubles altho the input elevations are unsigned short ints.
6. Using {$\le$} instead of {$<$} removes many more basins. The reason is that it raises pixels on the margins of flat regions and forestalls many ties that would have occurred in the adjacency matrix computation.
7. This means that when the 2 lowest neighbors are at the same elevation, then this pixel is raised to their elevation. If it were raised higher, it would create a ridge between those pixels, preventing this pixel from joining them into the same basin.
8. This process creates some new 1-pixel pits, from regions that were 2-pixel basins.
9. Therefore, this whole process might be repeated. Each time, it would remove a smaller number of pits.

1. A pair of pixels is adjacenct if one flows into the other, or if both flow into the same third pixel.
2. As of April 4, 2011, we assume 8-connectivity.
3. In this step, pixels on the terrain boundary are assumed to have a neighbor outside the edge of the world with a hi elevation. This keeps two basins that both flow off the edge of the world from joining.
4. We assume that the water flowing thru a pixel all flows out to that pixel's lowest neighbor. (That would be easy to change.)
5. When comparing elevations, ties are broken in a consistent manner. E.g., following the ideas of Simulation of Simplicity, instead of comparing just {$z_{i_1,j_1} $}, we lexicographically compare the triples {$ (z_{i_1,j_1}, i_1, j_1) $}.
6. Two pixels {$p_{i_1,j_1}$} and {$p_{i_2,j_2}$} are geometrically adjacent iff the 2nd is one pixel away from the 1st in any of the 8 directions.
7. Two geometrically adjacent pixels {$p_{i_1,j_1}$} and {$p_{i_2,j_2}$} are flow-adjacent if the water flowing out from one flows to the other. That is, either one of the two pixels is the other's lowest neighbor.
8. A basin is a connected set in the transitive completion of flow adjacency.
9. In what follows, we index arrays from 0.
10. The basic idea is:
    1. Create a graph where the nodes are the pixels and there is an edge between 2 nodes that are flow-adjacent, then
    2. Find its connected components.
11. connect2 works on a Euclidean bitmap, where 2 1-pixels are in the same component if they are adjacent. There are no explicit edges.
12. Therefore we will create such a bitmap from the terrain, and then find its connected components.
13. From the {$n\times n$} elevation array {$z$}, create a {$(3n-2)\times(3n-2)$} adjacency matrix {$\alpha$} (called adj in the program) as follows ¹.
    1. {$\alpha_{3i,3j}$} corresponds to {$p_{ij}$}.
    2. Set {$\alpha_{i,j}=0$}.
    3. Then set {$\alpha_{3i,3j}=1$}.
    4. Set {$\alpha_{3i+1,3j}=\alpha_{3i+2,3j}=1 $} iff {$p_{ij}$} is flow-adjacent to {$p_{i+1,j}$}.
    5. Set {$\alpha_{3i+1,3j+1}=\alpha_{3i+2,3j+2}=1 $} iff {$p_{ij}$} is flow-adjacent to {$p_{i+1,j+1}$}.
    6. Proceed analogously for the other 7 directions.

1. Now we process {$\alpha$} to find its connected components, or basins.
2. I am gradually changing the terminology in the program over from component, which was natural when connected components were the topic to basin, which is more natural here.
3. First we pad {$\alpha$} on all 4 sides with a row of pixels of a high elevation. This is to prevent basins that both flow off the edge of the world from joining.
4. The connection program is an extension of connect2, my fast Connected Components program. One change is that I am now connecting the 1 pixels, not the 0s that were more natural for connect2's original application (connecting voids).
5. The result is a matrix basinid giving the component (basin) number of each pixel.

1. Each basin has exactly one local minimum. (If there were a 2nd local min, the pixels flowing into it would not be flow connected to the others.)
2. That local min may or may not be at the edge of the basin.
3. In the normal case, a basin's local min is at the edge of the basin.
4. However, many basins have interior local mins, sometimes genuine geographic features, but usually caused by data errors or too coarse quantization.
5. These are unacceptable because they trap rivers and prevent any long rivers from forming. (The effect on the computed hydrography is quite dramatic.)
6. We will fix this by raising the elevations of certain pixels to remove all interior local mins, before calculating the flow.
7. This will result in all basins with interior local mins being destroyed by being joined to other basins.
8. Let the unique lowest point of basin {${\cal B}_i$} be {$d_i$}. {$d_i$} is lower than any other point in {${\cal B}_i$} and also lower than all of {$d_i$}'s 8 neighbors, even if they be in a different basin. (If {$d_i$} had a lower neighbor, then {${\cal B}_i$} wouldn't be a separate basin.)
9. That lowest point is unique because ties are broken consistently.
10. All the water in {${\cal B}_i$} flows to {$d_i$}.
11. The problem is to combine the many small basins into fewer larger ones.
12. A sill is a pixel {$ (x_1,y_1) $} in a basin {$ {\cal B}_1 $}, which is adjacent to a pixel {$ (x_2,y_2) $} in basin {$ {\cal B}_2 $}. Any other pixel in {$b_2$} that is adjacent to a pixel in {$ {\cal B}_1 $} is higher than pixel {$(x_2,y_2)$}. IOW, if {$ {\cal B}_1 $} were to overflow into {$ {\cal B}_2 $}, it would overflow from {$ (x_1,y_1) $} to {$ (x_2,y_2) $}.
13. {$ (x_1,y_1) $} may be lower, the same elevation as, or higher than {$(x_2,y_2)$}. Here are 1-D examples. In each case, red pixels are {${\cal B}_1$} and blue pixels {${\cal B}_2$}.
    lower: 1 2 3 1, same: 1 2 2 1, higher 1 3 2 1
14. Note the effect of the tie-breaking rule: 2 1 3 3 3 2.
15. The spill elevation of a sill of {$ {\cal B}_1 $} is the greater of the that sill's elevation and the elevation of the pixel in {$ {\cal B}_2$} that it would spill over into.
16. {$ {\cal B}_1 $} usuall has many sills, each one to a different basin.
17. We need to remember the higher sills because the lowest sill will be deleted should {$ {\cal B}_1 $} and {$ {\cal B}_2 $} be joined. That would cause the sill between them to be deleted.
18. Conceptually, gradually fill a basin until it flows over a sill into a lower neighboring basin, and join the 2 basins.
19. One complication is this case: 4 4 4 1 2 3 2
    We want to fill (red) {$ {\cal B}_1$} until it overflows into (blue) {$ {\cal B}_2$}. However, raising {$ {\cal B}_1$} until the level of the adjacent pixel in {$ {\cal B}_2$}, creates a new low region.

1. This routine processes one basin, {$ {\cal B}_1$}, each time it's called.
2. It handles the case that some pixels in {$ {\cal B}_1$} might be lower than {$ {\cal B}_1$}'s spill elevation (because the spill elevation might be higher because the adjacent pixel in the adjacent basin might be higher).
3. The process of raising certain pixels in {$ {\cal B}_1$} goes as follows.
   1. Create an empty deque {$\cal D$} of pixels in {$ {\cal B}_1$} waiting to be processed. (Using a deque causes the pixels to be processed in breadth-first order.)
   2. Create an empty set {$\cal S$} of pixels in {$ {\cal B}_1$} that have been processed.
   3. Assume a constant small delta elevation {$dz$}. It should be small enough that, even when added many times to a legal pixel elevation, we do not get as high as the next legal elevation. It should be large enough that it is different from zero when added.
   4. Initialize a current elevation {$z$} to the spill elevation plus {$dz$}.
   5. Insert the sill into {$\cal D$} and {$\cal S$}.
   6. Repeat the following until {$\cal D$} is empty.
      1. Pull the 1st pixel off of {$\cal D$}. Call it {$z_{ij}$}.
      2. If its elevation is greater than {$z$}, then this pixel was already elevated, so go back to pull another pixel. It is possible that its elevation is between the spill elevation and {$z$}, if it was earlier elevated by virtue of being in the depression of another sill. If so, it needs to be elevated again.
      3. Insert some or all of its 8 neighbors into {$\cal D$}. To be inserted, a neighbor must -
         1. not be in {$\cal S$},
         2. be in {$ {\cal B}_1$}, and
         3. have elevation less than {$z$}.
      4. Increase {$z$} by {$dz$}.
      5. Set {$z_{ij}$} to {$z$}.
      6. Note: This can cause a pixel that would have flowed to another pixel in this basin, now to flow in a different direction into another basin. I don't think that this is a problem.
      7. This also means that the adjacent pixel to a spill pixel might now be higher. Therefore the spill elevation must be raised.

The process of joining {$ {\cal B}_1$} into {$ {\cal B}_2$} goes as follows.

1. Call fill_depression on {$ {\cal B}_1$}.
2. Merge the sill list of {$ {\cal B}_1$} into {$ {\cal B}_2$}, keeping the list sorted by increasing spill elevation.
3. Make a note in the data structure for {$ {\cal B}_1$}, that it has been joined into {$ {\cal B}_2$}. However, if {$ {\cal B}_2$} was earlier joined into {$ {\cal B}_3$}, then point {$ {\cal B}_1$} directly to {$ {\cal B}_3$}.
4. If {$ {\cal B}_2$} now has two or more sills with the same spill basin, then delete the higher one(s).
5. If there is a sill whose spill basin is {$ {\cal B}_2$}, which includes {$ {\cal B}_1$}, then delete it.

1. There may arise chains of joined basins, which is inefficient and makes it hard to determine if two basins are ultimately joined to each other. Therefore I frequently follow these chains and collapse them so that each basin points straight to the root of the tree of joined basins.
2. The root basin id is the id of the original basin with the lowest spill elevation.
3. The process of joining all joinable basins goes as follows.
   1. Create a heap (priority queue) {$\cal H$} of all the basins in increasing order of spill height.
   2. Repeat until {$\cal H$} is empty.
      1. Pop the first basin {$ {\cal B}_1$} off of {$\cal H$}.
      2. If {$ {\cal B}_1$}'s spill elevation has changed since it was originally pushed onto {$\cal H$}, then push it onto {$\cal H$} again with its new elevation and continue to the next iteration.
      3. Its spill elevation could have increased because two basins mutually spilled into each other. When they were joined, the combined basin's new spill elevation is now higher. This figure shows the problem.
         
      4. If {$ {\cal B}_1$} does not spill off the edge of the world, (indicated by a basinid of -1), then join it to the first basin on its sill list, {$ {\cal B}_2$}.
      5. If it does spill off the edge of the world, then call fill_depression because it might have a depression next to its sill.
4. After joining, the remaining basins will be those that spill off the edge of the world.

Follow Eddie's suggestion to use the Terraflow idea, since, it's fast enough for the initial version here.

1. Sort the pixels in descending order of elevation. Break ties as usual.
2. Initialize each element of the flow array to 1, for the incoming rain.
3. Traversing the sorted pixel array, add each pixel's flow to the flow of its lowest neighbor, provided that the neighbor is lower than the pixel, or is of equal elevation, breaking ties the same as during the sorting.
   If there is no lower neighbor, then we've found a surviving interior sink, which should have been removed by now.
4. If the flow stops with this pixel, then sum its flow into totflow as a check.

1. Function ind uses a variable that is a pair of indices to index into a 2D array. It would be more obvious to overload operator[] but I can't see an easy way to do that.
2. Telev is double because a lot of significant digits are needed to handle the small increments added in fill_depression.

1. hydrog.cc is the major program.
2. Usage: hydrog NROWS < ELEVATION-FILE.elevb
3. Wrap function function don() { time hydrog $2 < data/$1.elevb convert flow.pgm $1-flow.png convert flow.pgm -gamma 5 -normalize $1-flow2.png rm -f flow.pgm xv $1-flow.png&
4. Output image file:
   1. FOO-flow.pgm The water flow thru each pixel, scaled with sqrt to reduce the range and to keep the max flow under 65535, which image manipulation programs prefer. Converting this to png saves space. Converting it to jpg would distort it.
5. Optional or obsolete output image files. The optional outputs are enabled by setting conditional compilation macros in the program. Enabling them significantly slows the program, since I/O is slow. Some of them might not be correct if that code hasn't been updated.
   1. FOO-pits.pbm The 1-pixel pits (that were then removed).
   2. FOO-comps.pgm The basins, before they were joined.
   3. FOO-edges.pbm The edges between the basins.
   4. FOO-edgesover.png Image of the edges overlaid on the original terrain.

1. float2int.cc Reads floats from cin, write ints to cout. Used to process ASCII data files written by Matlab: float2int < t10.dat >! t10.eleva cat head10.pgm t10.eleva >! t10.pgm
2. int2bin.cc Reads ascii unsigned short ints from cin and write in 2-byte binary to cout.
3. To convert 16-bit png to 16-bit ascii pgm:
   convert -compress None ps_height_1k.png ps_height_1k.pgm
4. To delete 1st 3 lines of a file:
   tail --lines=+4 in > out
5. To convert eleva file to pgm file: Prepend a header like head400.pgm
6. To cause ints to be written with commas thus: 12,345,678, set a locale in the program.

1. t10.eleva 10x10 test with a saddle, from t10.dat from Matlab. A subset of hill2.
2. {hill,mtn}{1,2,3}.eleva The GeoStar 400x400 datasets.
3. mtn1a.eleva 100x100 top left corner of mtn1
4. ps{1,4,16}k.eleva Puget Sound, from http://www.cc.gatech.edu/projects/large_models/ps.html . {$1025^2,4097^2,16387^2 $}. Originally called ps_height_1k.eleva etc.
5. FOO.elevb The above elevation files in binary format.

1. fns zsh functions to run the programs.
2. head10.pgm etc. Head for a 10x10 ASCII PGM file with a particular max value for a pixel.

1. .comp Output from connect2
2. .bita ASCII bitmap, containing only '0' and '1' characters, of a 2D or 3D rectangular region intended to have the connected components of regions of 0s (or is it the 1s?) computed. If the region is {$ (n_x,n_y,n_z) $} voxels, then this file's size is {$ n_xn_yn_z $} bytes.
3. .bitb Binary bitmap. File size is {$ \lceil n_xn_yn_z/8 \rceil $} bytes.
4. .svg Drawings created in inkscape. (However, Firefox doesn't display svg files in an IMG tag.)
5. .png Derivative drawings for display.
6. .eleva ASCII file of square matrix of elevations, containing {$ N^2 $} blank-separated positive ints.
7. .elevb binary file of square matrix of elevations, containing {$ N^2 $} 2-byte binary unsigned ints.

1. hydrography/ - the code and smaller data files.
2. Larger data files (Puget Sound):
   1. ps1k.eleva.bz2
   2. ps4k.eleva.bz2
   3. ps16k.eleva.bz2

1. When xv reads an ascii pgm it automatically scales stated max pixel value up or down to be 255. Therefore, sometimes set stated max to 255.
2. Ditto ImageMagick display apparently.
3. Dunno how to handle max over 255.
4. The display programs display pgm file pixels left to right, top to bottom.
5. However if you probe a pixel in xv, it lists its coords in the order (col,row), e.g., (y,x).

1. The current version is hydrography-20110424.tar.gz.
2. Hydrog apparently works, and very quickly, taking slightly worse than linear time per pixel.
3. It is economical with storage. The {$4097^2$} case uses 690MB of memory (including libraries, determined by viewing the process with top). That is only 43B per pixel. It might be reduceable.
4. Since the flow is so uneven (most pixels are 1; the biggest flow may have much of all the water in the world), you need to gamma scale the image to see a lot.

If you want to make hydrog even faster or smaller, start here:

1. Replace the STL sort routine, which is the slowest part of the whole program, with a radix sort.
2. Consider replacing the STL set with a boost set or a specially coded version.
3. Represent elevations with 4-byte fixed point numbers instead of 8-byte ints.

1. hydrog is fast enough that you can call it thousands of times. This introduces new research possibilities.
2. Compare to other methods, in speed, space, and in quality of output.
3. Use it to postprocess terrain produced by other methods.
4. Try on random terrain.
5. Study effects on hydrography of lower res terrain, or terrain with noise added.
6. Simple things:
   1. Overlay the flow.pgm image on the elevation image.
   2. Write new elevated terrain as a separate file.

1. http://terrain.cs.duke.edu/watershed/index.php
2. grass: r.terraflow
3. http://portal.acm.org/citation.cfm?id=940218