Small bug in excesstopography code

[wschwanghart]

This issue could also be submitted to the TTLEM3D repository, but it is probably better placed here. I ran following code in MATLAB using excesstopo:

G = gradient8(DEM);
MASK = DEM;
MASK.Z(2:end-1,2:end-1) = inf;
DEM2 = excesstopo(MASK,G);
imageschs(DEM2)

The resulting DEM2 contains numerous pixels with inf (any(isinf(DEM2.Z),'all')). I suspect that some pixels at the end of the queue do not get updated. The same is true, if I set all mask pixels to 5000.

MASK.Z(2:end-1,2:end-1) = 5000;

Then, any(DEM2 == 5000) is true.

[wkearn]

This is what you normally do when you solve the eikonal equation: fix the boundaries, set the interior to some high values, and march inward from the boundary. The two differences between excesstopography_* and a traditional eikonal solver are 1. the solution is constrained by the original DEM and 2. the boundaries of the DEM are handled a little differently to ensure they also satisfy the slope constraints. Basically they aren't fixed to their original values in the DEM, they are determined by marching from the lowest elevations in the domain. If you fix the boundaries, but you have a depression in your DEM that lies below the boundary elevation, there might be some inconsistency between the value fixed on the boundary and the value determined by marching out of the depression, so instead I let the boundary elevations "float." In your case, the lowest fixed elevations should be on the DEM boundary. My guess is that this leads to a situation where there aren't enough constraints to propagate to all of the interior points. I'll need to do some further investigation to figure out how to handle this situation.

[wkearn]

@wschwanghart: #54 should fix this. Let me know if you're still seeing similar errors.

You'll need to recompile libtopotoolbox and your excesstopography_* MEX files. The CMakeLists.txt file in TTLEM3D should automatically pull the updated main branch of libtopotoolbox when you run it again.