diff --git a/include/topotoolbox.h b/include/topotoolbox.h index 801e2637..901c50b4 100644 --- a/include/topotoolbox.h +++ b/include/topotoolbox.h @@ -2651,4 +2651,499 @@ ptrdiff_t simplify_line(float *out_i, float *out_j, const float *track_i, const float *track_j, ptrdiff_t n_points, float tolerance, int method); +/** + @brief Count edges in a bitmap edge set + + @param[in] bitmap The input bitmap + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + @endparblock + + @param[in] dims The dimensions of the bitmap + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @return Number of edges in the edge set +*/ +TOPOTOOLBOX_API +ptrdiff_t edgeset_count(uint8_t *bitmap, ptrdiff_t dims[2]); + +/** + @brief Determine edge weight offsets in a bitmap edge set + + @details A bitmap edge set stores the edge weights in a packed + array of floats sorted by the source pixel. Since each pixel may + have a different number of outgoing edges, the offset into the + weights array corresponding to the first edge leaving a source + pixel depends on the number of edges leaving the pixels preceding + the source. This is computed with a prefix sum (a scan) of the + bitmap array. + + @param[out] scan The edge weight offsets + @parblock + + An array of ptrdiff_t of the same size as the bitmap. The value at + pixel (i, j) is the offset into the weights array at which the + first outgoing edge of pixel (i, j) can be found if it is + present. Subsequent edges are stored contiguously. + + @endparblock + + @param[in] bitmap The input bitmap + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + @endparblock + + @param[in] dims The dimensions of the bitmap + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @return Number of edges in the edge set + */ +TOPOTOOLBOX_API +ptrdiff_t edgeset_scan(ptrdiff_t *scan, uint8_t *bitmap, ptrdiff_t dims[2]); + +/** + @brief Count edges in a merged bitmap edge set + + @details To merge two bitmap edge sets one takes the bitwise or of + corresponding pixels in the two bitmaps. This function computes the + number of edges that will result from merging two input + bitmaps. This is necessary to preallocate the edge weight array + before merging. + + @param[in] b1 The first input bitmap + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + @endparblock + + @param[in] b2 The second input bitmap + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + @endparblock + + @param[in] dims The dimensions of the bitmap + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @return Number of edges in the edge set +*/ +TOPOTOOLBOX_API +ptrdiff_t edgeset_count_merged(uint8_t *b1, uint8_t *b2, ptrdiff_t dims[2]); + +/** + @brief Merge two bitmap edge sets + + @details To merge two bitmap edge sets one takes the bitwise or of + corresponding pixels in the two bitmaps and sorts the edge weights + appropriately. + + @param[out] w0 The output edge weight array + @parblock + + An array of floats containing the edge weights of the merged edge + set. This array should be preallocated with a size of at least the + number of edges returned by edgeset_count_merged. + + @endparblock + + @param[out] s0 The output offset array + @parblock + + An array of ptrdiff_t containing the offsets in the edge weight + array of the first outgoing edge for each pixel, if it exists. This + can also be computed after the merge with edgeset_scan. + + @endparblock + + @param[inout] b1 The first input bitmap + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + The first input bitmap will be updated with edges from the second + input bitmap. + + @endparblock + + @param[inout] w1 The first input weight array + @parblock + + An array of floats stored in the packed weight format corresponding + to the first input bitmap. + + @endparblock + + @param[in] b2 The second input bitmap + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + @endparblock + + @param[in] w2 The second input weight array + @parblock + + An array of floats stored in the packed weight format corresponding + to the second input bitmap. + + @endparblock + + @param[in] dims The dimensions of the bitmap + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @return Number of edges in the edge set +*/ +TOPOTOOLBOX_API +ptrdiff_t edgeset_merge(float *w0, ptrdiff_t *s0, uint8_t *b1, float *w1, + uint8_t *b2, float *w2, ptrdiff_t dims[2]); + +/** + @brief Topologically sort a set of edges + + @details flow_routing_tsort performs a topological sort of the + edges stored in a bitmap edge set. This is a necessary step in + creation of a TopoToolbox Flow object. This function can be used + with any flow routing algorithm that outputs an edge set. + + @param[out] stream The topologically sorted list of nodes + @parblock + + This is an array of ptrdiff_t that will be filled with the linear + indices of nodes in topological order. It should be preallocated + with a size equal to the number of pixels in the DEM. + + @endparblock + + @param[out] source The topologically sorted list of source pixels + @parblock + + An array of ptrdiff_t that will be filled with the source indices + for each edge in the flow network. Topologically sorted so that all + incoming edges to a pixel occur before any outgoing pixels. This + array should be preallocated with a size equal to the number of + edges in the edge set as returned by any of the edgeset functions. + + @endparblock + + @param[out] target The topologically sorted list of target pixels + @parblock + + An array of ptrdiff_t that will be filled with the target indices + for each edge in the flow network. Topologically sorted so that all + incoming edges to a pixel occur before any outgoing pixels. This + array should be preallocated with a size equal to the number of + edges in the edge set as returned by any of the edgeset functions. + + @endparblock + + @param[out] sorted_weight The topologically sorted list of edge weights + @parblock + + An array of floats that will be filled with the edge weights for + each edge in the flow network. The order of the weights in this + array corresponds to the source and target arrays. This array + should be preallocated with a size equal to the number of edges in + the edge set as returned by any of the edgeset functions. + + @endparblock + + @param[in] stack A stack used for the graph traversal + @parblock + + An array of ptrdiff_t that is used as a stack of node indices + during the graph traversal. This is a temporary array only needed + by this function. It should be preallocated with a size equal to + the number of pixels in the DEM. + + @endparblock + + @param[in] stackdir A stack used for the graph traversal + @parblock + + An array of uint8_t that is used as a stack of edge directions + during the graph traversal. This is a temporary array only needed + by this function. It should be preallocated with a size equal to + the number of pixels in the DEM. + + @endparblock + + @param[in] direction The input bitmap edge set + @parblock + + The input bitmap is an array of uint8_t. Each of the 8 bits + corresponds to an outgoing edge from that pixel starting with the + least significant bit representing the edge to the right and moving + counterclockwise. + + @endparblock + + @param[in] weight The edge weight array + @parblock + + An array of floats containing the edge weights for each edge. The + edge weights are densely packed and ordered by the source pixel + then by the edge index as described in edgeset_scan. + + @endparblock + + @param[in] weightscan The array of offsets into the edge weight array + @parblock + + An array of ptrdiff_t containing the offsets of the weight of the + first outgoing edge for each pixel in the weight array. This should + be generated by edgeset_scan or edgeset_merge. + + @endparblock + + @param[in] visited A marker array for the graph traversal + @parblock + + An array of uint8_t that is used to mark visited pixels during the + graph traversal. This is a temporary array only needed by this + function. It should be preallocated with a size equal to the number + of pixels in the DEM. + + @endparblock + + @param[in] edge_count The number of edges in the edge set + @parblock + + This can be obtained as the return value from any of the edgeset + functions. + + @endparblock + + @param[in] dims The dimensions of the bitmap + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @param[in] order The memory order of the underlying arrays + @parblock + 0 for column-major, 1 for row-major + @endparblock + */ +TOPOTOOLBOX_API +void flow_routing_tsort(ptrdiff_t *stream, ptrdiff_t *source, ptrdiff_t *target, + float *sorted_weight, ptrdiff_t *stack, + uint8_t *stackdir, uint8_t *direction, float *weight, + ptrdiff_t *weightscan, uint8_t *visited, + ptrdiff_t edge_count, ptrdiff_t dims[2], int order); + +/** + @brief Compute flow directions using the D-inf method of Tarboton (1997) + + @details D-inf chooses up to two downstream neighbors based on the + steepest downward slope of 8 triangular facets surrounding a center + pixel and apportions flow between the two pixels according to the + angle of the downward slope vector. + + Tarboton, D. G. (1997), A new method for the determination of flow + directions and upslope areas in grid digital elevation models, + Water Resour. Res., 33(2), 309–319, doi:10.1029/96WR03137. + + @param[out] direction A bitmap edge set + @parblock + + A uint8_t array of the same size as the DEM. Each of the 8 bits of + each pixel corresponds to an outgoing edge. + + @endparblock + + @param[out] prop Flow proportion towards the first outgoing edge + @parblock + + A float array of the same size as the DEM. This is the proportion + of flow that is routed towards the first outgoing edge from each + pixel. + + @endparblock + + @param[in] dem The input DEM + @parblock + + A float array with size dims[0] x dims[1]. + + @endparblock + + @param[in] dims The dimensions of the DEM + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @param[in] order The memory order of the underlying arrays + @parblock + 0 for column-major, 1 for row-major + @endparblock + */ +TOPOTOOLBOX_API +void flow_routing_dinf_directions(uint8_t *direction, float *prop, float *dem, + ptrdiff_t dims[2], int order); + +/** + @brief Compute and store the D-inf edge weights + + @details Once the flow directions and flow proportions have been + computed, the edge weights need to be stored in an array in an + order that will be used by functions like flow_routing_tsort. + + @param[out] weight The weight array + @parblock + + An array of floats that will be filled with the edge weights. This + should be preallocated with a size equal to the number of edges in + the direction bitmap. This size can be computed using + edgeset_count. + + @endparblock + + @param[in] direction The input direction bitmap + @parblock + + An array of uint8_t that should come from the direction argument of + flow_routing_dinf_directions. This is a bitmap edge set encoded as + described elsewhere. + + @endparblock + + @param[in] prop The input flow proportions + @parblock + + An array of floats that should come from the prop argument of + flow_routing_dinf_directions. This is the proportion of flow routed + to the first outgoing edge of each pixel. + + @endparblock + + @param[in] dims The dimensions of the DEM + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + */ +TOPOTOOLBOX_API +void flow_routing_dinf_weights(float *weight, uint8_t *direction, float *prop, + ptrdiff_t dims[2]); + +/** + @brief Resolve flats using a shortest path algorithm + + @details This flat resolution algorithm flows flats along the + shortest path to sills. Sills are resolved pixels with the same + elevation as neighboring flats. The flat pixels neighboring sills + are resolved to flow into their neighboring sill. The remaining + flats are then resolved to flow along the shortest path to a sill. + + Uses Dijkstra's algorithm with a preallocated priority queue. + + @param[out] direction The output direction bitmap + @param[in] resolved nonzero for pixels that have resolved directions + @parblock + + An array of uint8_t flagging unresolved pixels with a zero. Only + pixels with a zero in this array will be resolved by the shortest + path algorithm. + + @endparblock + + @param path_distance A temporary array used by the shortest path algorithm + @parblock + + An array of floats of size dims[0] x dims[1]. + @endparblock + + @param heap A temporary array used by the shortest path algorithm + @parblock + + An array of ptrdiff_t of size dims[0] x dims[1]. + + @endparblock + + @param back A temporary array used by the shortest path algorithm + @parblock + + An array of ptrdiff_t of size dims[0] x dims[1]. + + @endparblock + + @param[in] dem The input dem + @parblock + + This is used to make sure that the pixels only flow to neighbors of + the same elevation. + + @endparblock + + @param[in] dims The dimensions of the DEM + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @param[in] order The memory order of the underlying arrays + @parblock + 0 for column-major, 1 for row-major + @endparblock + */ +TOPOTOOLBOX_API +void resolve_flats_shortest_path(uint8_t *direction, uint8_t *resolved, + float *path_distance, ptrdiff_t *heap, + ptrdiff_t *back, float *dem, ptrdiff_t dims[2], + int order); + +/** + @brief Compute the weights of the shortest path flat resolution + algorithm + + @details The weights are all 1. This function is provided for + consistency with other flat resolution implementations. + + @param[out] weight The array of edge weights + @parblock + + An array of floats preallocated with a size count. + + @endparblock + + @param[in] count The number of edges + @parblock + + This can be computed from the output direction bitmap of + resolve_flats_shortest_path using edgeset_count. + + @endparblock + */ +TOPOTOOLBOX_API +void resolve_flats_shortest_path_weights(float *weight, ptrdiff_t count); + #endif // TOPOTOOLBOX_H diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4b844ae1..55289aa5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -36,6 +36,8 @@ add_library(topotoolbox helpers/stat_func.h helpers/deque.c helpers/deque.h + helpers/edgeset.c + dinf.c ) # Define the include directory diff --git a/src/dinf.c b/src/dinf.c new file mode 100644 index 00000000..2e4e0669 --- /dev/null +++ b/src/dinf.c @@ -0,0 +1,165 @@ +#define TOPOTOOLBOX_BUILD + +#include +#include +#include +#include + +#include "topotoolbox.h" + +#define SQRT2 1.41421356237309504880f +#define PI 3.14159265358979323846f + +static float replicate_boundaries(float *dem, ptrdiff_t i, ptrdiff_t j, + ptrdiff_t ei, ptrdiff_t ej, + ptrdiff_t dims[2]) { + float z = dem[j * dims[0] + i]; + if (i + ei < 0 || i + ei >= dims[0]) { + if (j + ej < 0 || j + ej >= dims[1]) { + // Corners: z = dem[i, j] + } else { + // Top and bottom boundaries + z = dem[(j + ej) * dims[0] + i]; + } + } else if (j + ej < 0 || j + ej >= dims[1]) { + // Left and right boundaries. The ei value is valid because we + // skipped the previous branch. + z = dem[j * dims[0] + i + ei]; + } else { + // Interior + z = dem[(j + ej) * dims[0] + i + ei]; + } + return z; +} + +typedef struct { + int n; + float a; +} facet; + +static facet dinf_flow_prop(float v1, float v2, int order) { + facet out = {-1, 0}; + if (v1 == 0.0 && v2 == 0.0) { + return out; + } + // Compute the angle + // + // Despite appearances, this is correct. One must rotate the + // coordinate axes to line up the (v1, v2) vector with the (ei, ej) + // vectors. This depends on the memory order of the array. Then we + // work in units of PI/4 so that we only have to do this division + // once. We shift the angles from [-PI,PI] to [0,2PI] with the mod(x + // + 8, 8) (8 is 2PI in units of PI/4). + if (order) { + out.a = fmodf(atan2f(-v2, v1) / (PI / 4) + 8, 8); + } else { + out.a = fmodf(atan2f(-v1, v2) / (PI / 4) + 8, 8); + } + + // The result is a number between 0 and 8. The integer part is the + // facet number. The fractional part is 1 - a, the correct flow + // proportion proportion. + out.n = (int)floorf(out.a); + out.a = (out.n + 1) - out.a; + + return out; +} + +TOPOTOOLBOX_API +void flow_routing_dinf_directions(uint8_t *direction, float *prop, float *dem, + ptrdiff_t dims[2], int order) { + // Compute and store the direction bitmap and the flow proportions + // rather than the 2D flow vectors. + ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1}, + {1, 1, 0, -1, -1, -1, 0, 1}}; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + ptrdiff_t idx0 = j * dims[0] + i; + float z0 = dem[idx0]; + + direction[idx0] = 0; + prop[idx0] = 0.0; + + if (isnan(z0)) { + prop[idx0] = NAN; + continue; + } + + float maxslope = 0.0; + for (int n = 0; n < 8; n++) { + ptrdiff_t i1 = e[order & 1][(n + (order & 1)) & 7]; + ptrdiff_t j1 = e[(order ^ 1) & 1][(n + (order & 1)) & 7]; + float z1 = replicate_boundaries(dem, i, j, i1, j1, dims); + + ptrdiff_t i2 = e[order & 1][(n + 1 - (order & 1)) & 7]; + ptrdiff_t j2 = e[(order ^ 1) & 1][(n + 1 - (order & 1)) & 7]; + float z2 = replicate_boundaries(dem, i, j, i2, j2, dims); + + float s1 = j1 * (z2 - z0) - j2 * (z1 - z0); + float s2 = i1 * (z0 - z2) - i2 * (z0 - z1); + + float mag = sqrtf(s1 * s1 + s2 * s2); + float d1 = sqrtf((float)(i1 * i1 + j1 * j1)); + float d2 = sqrtf((float)(i2 * i2 + j2 * j2)); + + if (i1 * s2 - j1 * s1 < 0) { + mag = (s1 * i1 + s2 * j1) / d1; + s1 = mag * i1 / d1; + s2 = mag * j1 / d1; + } else if (j2 * s1 - i2 * s2 < 0) { + mag = (s1 * i2 + s2 * j2) / d2; + s1 = mag * i2 / d2; + s2 = mag * j2 / d2; + } + + if (mag > maxslope) { + facet flowdir = dinf_flow_prop(s1, s2, order); + + prop[idx0] = flowdir.a; + + direction[idx0] = 0; + if (flowdir.a > 0) { + direction[idx0] |= 1 << (flowdir.n & 7); + } + + if (flowdir.a < 1) { + direction[idx0] |= 1 << ((flowdir.n + 1) & 7); + } + maxslope = mag; + } + } + } + } +} + +TOPOTOOLBOX_API +void flow_routing_dinf_weights(float *weight, uint8_t *direction, float *prop, + ptrdiff_t dims[2]) { + ptrdiff_t e = 0; + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + ptrdiff_t idx = j * dims[0] + i; + float a = prop[idx]; + + if (direction[idx] == 129) { + // Swap the order of the edges + if (a < 1) { + weight[e++] = 1 - a; + } + + if (a > 0) { + weight[e++] = a; + } + } else if (direction[idx] > 0) { + if (a > 0) { + weight[e++] = a; + } + + if (a < 1) { + weight[e++] = 1 - a; + } + } + } + } +} diff --git a/src/flow_routing.c b/src/flow_routing.c index 76434b14..9978b32c 100644 --- a/src/flow_routing.c +++ b/src/flow_routing.c @@ -5,8 +5,11 @@ #include #include +#include "helpers/priority_queue.h" #include "topotoolbox.h" +#define SQRT2 1.41421356237309504880f + // Compute the steepest descent flow direction from the DEM with flow // routed over flat regions by the auxiliary topography in dist. uint8_t compute_flowdirection(ptrdiff_t i, ptrdiff_t j, float *dem, float *dist, @@ -307,3 +310,206 @@ ptrdiff_t flow_routing_d8_edgelist(ptrdiff_t *source, ptrdiff_t *target, } return edge_count; } + +/////////////////////////////////// +// Topological sorting of edge sets + +static uint8_t next_neighbor(uint8_t d) { + uint8_t n = 0; + while (d >>= 1) { + n++; + } + return n; +} + +TOPOTOOLBOX_API +void flow_routing_tsort(ptrdiff_t *stream, ptrdiff_t *source, ptrdiff_t *target, + float *sorted_weight, ptrdiff_t *stack, + uint8_t *stackdir, uint8_t *direction, float *weight, + ptrdiff_t *weightscan, uint8_t *visited, + ptrdiff_t edge_count, ptrdiff_t dims[2], int order) { + ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1}, + {1, 1, 0, -1, -1, -1, 0, 1}}; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + visited[j * dims[0] + i] = 0; + } + } + + ptrdiff_t stack_top = 0; + ptrdiff_t next_node = dims[0] * dims[1]; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + ptrdiff_t idx = j * dims[0] + i; + + if (visited[idx]) continue; + + visited[idx] = 1; + + stack[stack_top] = idx; + stackdir[stack_top++] = direction[idx]; + + while (stack_top > 0) { + while (stackdir[stack_top - 1]) { + // We still have neighbors left to check + ptrdiff_t node = stack[stack_top - 1]; + ptrdiff_t iu = node % dims[0]; + ptrdiff_t ju = node / dims[0]; + + uint8_t n = next_neighbor(stackdir[stack_top - 1]); + stackdir[stack_top - 1] ^= (1 << n); + + ptrdiff_t in = iu + e[order & 1][n]; + ptrdiff_t jn = ju + e[(order ^ 1) & 1][n]; + + if (!visited[jn * dims[0] + in]) { + visited[jn * dims[0] + in] = 1; + + // Push the neighbor and its outgoing edges on the stack + stack[stack_top] = jn * dims[0] + in; + stackdir[stack_top++] = direction[jn * dims[0] + in]; + } else if (visited[jn * dims[0] + in] == 1) { + // This is a cycle. The neighbor is already on the stack + assert(0 && "flow_routing_tsort detected a cycle. This is a bug."); + } + } + // Pop the stack. + ptrdiff_t node = stack[--stack_top]; + ptrdiff_t iu = node % dims[0]; + ptrdiff_t ju = node / dims[0]; + + visited[node] = 2; + + if (next_node <= 0) { + assert(0 && "flow_routing_tsort ran out of nodes. This is a bug."); + } + stream[--next_node] = node; + + ptrdiff_t offset = weightscan[node]; + + // Add its edges to source and target + for (int n = 0; n < 8; n++) { + if (direction[node] & (1 << n)) { + ptrdiff_t in = iu + e[order & 1][n]; + ptrdiff_t jn = ju + e[(order ^ 1) & 1][n]; + + if (edge_count <= 0) { + assert(0 && + "flow_routing_tsort ran out of nodes. This is a bug."); + } + source[--edge_count] = node; + target[edge_count] = jn * dims[0] + in; + + sorted_weight[edge_count] = weight[offset++]; + } + } + } + } + } +} + +static uint8_t get_direction(ptrdiff_t i, ptrdiff_t j, int order) { + uint8_t e[2][9] = {{8, 16, 32, 4, 0, 64, 2, 1, 128}, + {8, 4, 2, 16, 0, 1, 32, 64, 128}}; + + return e[order][(j + 1) * 3 + (i + 1)]; +} + +TOPOTOOLBOX_API +void resolve_flats_shortest_path(uint8_t *direction, uint8_t *resolved, + float *path_distance, ptrdiff_t *heap, + ptrdiff_t *back, float *dem, ptrdiff_t dims[2], + int order) { + ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1}, + {1, 1, 0, -1, -1, -1, 0, 1}}; + float chamfer[8] = {1.0f, SQRT2, 1.0f, SQRT2, 1.0f, SQRT2, 1.0f, SQRT2}; + + PriorityQueue pq = pq_create(dims[0] * dims[1], heap, back, path_distance, 0); + + // Initialize distances based on the resolved array. + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + ptrdiff_t idx = j * dims[0] + i; + if (resolved[idx] == 0) { + pq_insert(&pq, idx, INFINITY); + } else { + path_distance[idx] = -INFINITY; + } + } + } + // Identify presills, flats that border a resolved pixel with the same + // elevation + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + ptrdiff_t idx = j * dims[0] + i; + if (resolved[idx]) { + continue; + } + float z = dem[idx]; + for (ptrdiff_t n = 0; n < 8; n++) { + ptrdiff_t in = i + e[order & 1][n]; + ptrdiff_t jn = j + e[(order ^ 1) & 1][n]; + + if (in < 0 || in >= dims[0] || jn < 0 || jn >= dims[1]) { + // Pixel is on the boundary, call it a presill + pq_decrease_key(&pq, idx, 0.0); + break; + } + + ptrdiff_t idxn = jn * dims[0] + in; + + if (resolved[idxn] && z == dem[idxn]) { + // These are the presills, so their distance is 0 + pq_decrease_key(&pq, idx, 0.0); + + // Flow the presill towards the neighbor + direction[idx] = get_direction(in - i, jn - j, order); + break; + } + } + } + } + + while (!pq_isempty(&pq)) { + ptrdiff_t idx = pq_deletemin(&pq); + float d = pq_get_priority(&pq, idx); + + ptrdiff_t j = idx / dims[0]; + ptrdiff_t i = idx % dims[0]; + + float z = dem[idx]; + + for (ptrdiff_t n = 0; n < 8; n++) { + ptrdiff_t in = i + e[order & 1][n]; + ptrdiff_t jn = j + e[(order ^ 1) & 1][n]; + + // Skip any neighboring pixels that are outside the domain, are + // resolved or are not equal in elevation. The weird elevation + // comparison is in case the neighboring pixel is a NaN, in + // which case !(neighbor_z <= z) is true and we skip the pixel. + if (in < 0 || in >= dims[0] || jn < 0 || jn >= dims[1] || + resolved[jn * dims[0] + in] || !(dem[jn * dims[0] + in] == z)) { + continue; + } + float proposal = d + chamfer[n]; + + if (proposal < pq_get_priority(&pq, jn * dims[0] + in)) { + pq_decrease_key(&pq, jn * dims[0] + in, proposal); + + // Save the flow direction to the current pixel as the neighbor's flow + // direction + direction[jn * dims[0] + in] = get_direction(i - in, j - jn, order); + } + } + } +} + +// Weights are 1 for shortest path flat resolution +TOPOTOOLBOX_API +void resolve_flats_shortest_path_weights(float *weight, ptrdiff_t count) { + for (ptrdiff_t e = 0; e < count; e++) { + weight[e] = 1.0; + } +} diff --git a/src/helpers/edgeset.c b/src/helpers/edgeset.c new file mode 100644 index 00000000..d004d4cd --- /dev/null +++ b/src/helpers/edgeset.c @@ -0,0 +1,118 @@ +#define TOPOTOOLBOX_BUILD + +#include "topotoolbox.h" + +///////////////////////////////////////// +// Bit maps for storing weighted edge sets +// +// Bit maps are a combination of a uint8_t array to identify outgoing +// edges of each pixel and a float array of corresponding edge +// weights. The k-th of byte [i, j] in the array is set if the k-th +// outgoing edge is present. The edges are numbered like +// +// 3 | 2 | 1 +// 4 | x | 0 +// 5 | 6 | 7 +// +// This is the ordering used by the Dinf algorithm. Note that in +// column-major, i increases down and j increases right. In row-major +// i increases right and j increases down. The bitmaps don't care +// about row or column major in the usual way of libtopotoolbox +// functions: reverse the order of the dims array if you are using +// row-major arrays. +// +// The edge weights are stored in an appropriately sized float +// array. The offset of each pixel's edge weights within this array is +// given by the running sum of the number of edges as you scan through +// the bitmap. The edges for each pixel are stored in the order of the +// edge numbering (counterclockwise from right). + +// Count the number of bits set in a byte. +static int bitcount(uint8_t v) { + int c; + for (c = 0; v; c++) v &= v - 1; + return c; +} + +// Count all the edges in the bitmap. This is used to preallocate +// arrays of edges. +TOPOTOOLBOX_API +ptrdiff_t edgeset_count(uint8_t *bitmap, ptrdiff_t dims[2]) { + ptrdiff_t n = 0; + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + n += bitcount(bitmap[j * dims[0] + i]); + } + } + return n; +} + +// Prefix sum (scan) of the number of edges. The scan array contains +// offsets into the list of weights for each pixel. Returns the number +// of edges, so you can use scan_edges to prepare the scan array and +// count edges for allocating weights. +TOPOTOOLBOX_API +ptrdiff_t edgeset_scan(ptrdiff_t *scan, uint8_t *bitmap, ptrdiff_t dims[2]) { + ptrdiff_t n = 0; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + scan[j * dims[0] + i] = n; + n += bitcount(bitmap[j * dims[0] + i]); + } + } + return n; +} + +// Actually storing the weights in the edge array is the +// responsibility of the caller because they may need to compute the +// weights or other quantities. + +// Count the total number of edges in the set merge of b1 and b2 This +// is used to preallocate the array of weights for the merged +// set. Duplicated edges are counted only once: this is not +// a multiset. +TOPOTOOLBOX_API +ptrdiff_t edgeset_count_merged(uint8_t *b1, uint8_t *b2, ptrdiff_t dims[2]) { + ptrdiff_t n = 0; + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + n += bitcount(b1[j * dims[0] + i] | b2[j * dims[0] + i]); + } + } + return n; +} + +// Merge the two edge sets (b1, w1) and (b2, w2). The edge weights are +// stored in the array w0 at the offsets in s0. If an edge exists in +// both b1 and b2, the weight from b1 is taken. The offset arrays of +// b1 and b2 are not necessary because the edge indices are computed +// along the way. The b1 array is updated with the new edges. +TOPOTOOLBOX_API +ptrdiff_t edgeset_merge(float *w0, ptrdiff_t *s0, uint8_t *b1, float *w1, + uint8_t *b2, float *w2, ptrdiff_t dims[2]) { + ptrdiff_t n = 0; + ptrdiff_t n1 = 0; + ptrdiff_t n2 = 0; + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + s0[j * dims[0] + i] = n; + + // Now loop over the neighbors + for (ptrdiff_t k = 0; k < 8; k++) { + if (b1[j * dims[0] + i] & (1 << k)) { + w0[n++] = w1[n1++]; + + // If the bit is also set in b2, skip that weight. + if (b2[j * dims[0] + i] & (1 << k)) n2++; + } else if (b2[j * dims[0] + i] & (1 << k)) { + w0[n++] = w2[n2++]; + + // Set the bit in b1 + b1[j * dims[0] + i] |= (1 << k); + } + } + } + } + return n; +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c32392c9..e7507c8e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -97,7 +97,8 @@ endif() target_link_libraries(polyline PRIVATE topotoolbox) add_test(NAME polyline COMMAND polyline) set_tests_properties(polyline PROPERTIES ENVIRONMENT_MODIFICATION - "PATH=path_list_prepend:$<$:$>") + "PATH=path_list_prepend:$<$:$>") + # TEST : snapshots # @@ -150,9 +151,23 @@ if(GDAL_FOUND AND TT_SNAPSHOT_DIR) target_link_libraries(snapshot PRIVATE topotoolbox GDAL::GDAL) add_test(NAME snapshot COMMAND snapshot ${TT_SNAPSHOT_DIR}) set_tests_properties(snapshot PROPERTIES ENVIRONMENT_MODIFICATION - "PATH=path_list_prepend:$<$:$>") + "PATH=path_list_prepend:$<$:$>") + + # TEST : dinf + # + add_executable(dinf dinf.c utils.c grid.c flow.c) + if(TT_SANITIZE AND NOT MSVC) + target_compile_options(dinf PRIVATE "$<$:-fsanitize=address>") + target_link_options(dinf PRIVATE "$<$:-fsanitize=address>") + endif() + target_link_libraries(dinf PRIVATE topotoolbox GDAL::GDAL) + add_test(NAME dinf COMMAND dinf ${TT_SNAPSHOT_DIR}) + set_tests_properties(dinf PROPERTIES ENVIRONMENT_MODIFICATION + "PATH=path_list_prepend:$<$:$>") endif() + + # TEST: clang-format # # Runs clang format to verify formatting of code diff --git a/test/dinf.c b/test/dinf.c new file mode 100644 index 00000000..80ffb666 --- /dev/null +++ b/test/dinf.c @@ -0,0 +1,763 @@ +#undef NDEBUG + +#include +#include +#include +#include +#include +#include + +#include "flow.h" +#include "grid.h" +#include "topotoolbox.h" +#include "utils.h" + +#define test(desc) \ + for (int _res = 0, _break = 1; _break; _break = 0, \ + _res ? printf("not ok %d - %s\n", ++_testcount, desc) \ + : printf("ok %d - %s\n", ++_testcount, desc), \ + _success |= _res) +#define test_fail _res = 1 +#define test_init() \ + for (int _success = 0, _testcount = (printf("TAP version 14\n"), 0), \ + _break = 1; \ + _break; _break = 0, printf("1..%d\n", _testcount), exit(_success)) + +ptrdiff_t ei[8] = {0, -1, -1, -1, 0, 1, 1, 1}; +ptrdiff_t ej[8] = {1, 1, 0, -1, -1, -1, 0, 1}; + +void print_flow(Flow f) { + printf("%ld edges\n", f.count); + for (ptrdiff_t e = 0; e < f.count; e++) { + printf("%ld %ld %f\n", f.source[e], f.target[e], f.fraction[e]); + } +} + +int test_node_once(Flow f) { + // Has every node been assigned exactly once? + Grid count = GridCreate(GridF32, NULL, f.cellsize, f.dims); + for (ptrdiff_t j = 0; j < f.dims[1]; j++) { + for (ptrdiff_t i = 0; i < f.dims[0]; i++) { + ((float *)count.data)[f.stream[j * f.dims[0] + i]] += 1.0; + } + } + + Grid ones = GridFill(1.0, count.cellsize, count.dims); + + int res = GridEq(count, ones); + + GridFree(&count); + GridFree(&ones); + + return res; +} + +int test_tsort_node(Flow fd) { + // Is the node list topologically sorted? + // For every edge u => v, u comes before v in the node list + + int res = 0; + + // idxs maps every cell to its position in the node list + Grid idxs = GridCreate(GridIdx, NULL, fd.cellsize, fd.dims); + for (ptrdiff_t j = 0; j < fd.dims[1]; j++) { + for (ptrdiff_t i = 0; i < fd.dims[0]; i++) { + ((ptrdiff_t *)idxs.data)[fd.stream[j * fd.dims[0] + i]] = + j * fd.dims[0] + i; + } + } + + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; + + res = ((ptrdiff_t *)idxs.data)[u] < ((ptrdiff_t *)idxs.data)[v] ? 1 : 0; + } + + GridFree(&idxs); + return res; +} + +int test_tsort(Flow fd) { + // Are the edge lists topologically sorted? + // All incoming edges to a vertex v occur before any outgoing edges + + int tsort_success = 1; + + Grid outdegree = GridCreate(GridU8, NULL, fd.cellsize, fd.dims); + + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; + + ((uint8_t *)outdegree.data)[u] += 1; + if (((uint8_t *)outdegree.data)[v] > 0) { + tsort_success = 0; + } + } + if (!tsort_success) { + printf("# Edge list was not topologically sorted\n"); + } + GridFree(&outdegree); + return tsort_success; +} + +int test_nanflow(Flow fd, Grid dem) { + // Missing data in the DEM should have no flow into or out of them. + int res = 1; + + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; + if (isnan(((float *)dem.data)[u])) { + res = 0; + printf("# NaN pixel of DEM is a source: %ld\n", u); + } + + if (isnan(((float *)dem.data)[v])) { + res = 0; + printf("# NaN pixel of DEM is a target: %ld\n", v); + } + } + + return res; +} + +int test_no_boundary_edges(Flow fd) { + // No edges should originate from outside the domain or leave the + // domain + int res = 1; + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; + + ptrdiff_t ui = u % fd.dims[0]; + ptrdiff_t uj = u / fd.dims[0]; + + ptrdiff_t vi = v % fd.dims[0]; + ptrdiff_t vj = v / fd.dims[0]; + + if (llabs(ui - vi) > 1 || llabs(uj - vj) > 1) { + printf("# Indices (%ld, %ld) and (%ld,%ld) are not neighbors.\n", ui, uj, + vi, vj); + } + + if (ui < 0 || ui >= fd.dims[0] || uj < 0 || uj >= fd.dims[1]) { + printf("# Edge %ld => %ld originates from outside the domain.\n", u, v); + res = 0; + } + + if (vi < 0 || vi >= fd.dims[0] || vj < 0 || vj >= fd.dims[1]) { + printf("# Edge %ld => %ld leaves the domain.\n", u, v); + res = 0; + } + } + return res; +} + +int test_outlets(Flow fd) { + // There should be at least one outlet, a pixel with incoming edges + // but no outgoing edges. + Grid outdegree = GridCreate(GridU8, NULL, fd.cellsize, fd.dims); + Grid indegree = GridCreate(GridU8, NULL, fd.cellsize, fd.dims); + + uint8_t *outdegree_data = outdegree.data; + uint8_t *indegree_data = indegree.data; + + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; + + outdegree_data[u] += 1; + indegree_data[v] += 1; + } + + int res = 0; + for (ptrdiff_t j = 0; j < fd.dims[1]; j++) { + for (ptrdiff_t i = 0; i < fd.dims[0]; i++) { + if (indegree_data[j * fd.dims[0] + i] > 0 && + outdegree_data[j * fd.dims[0] + i] == 0) { + // This is an outlet + res = 1; + } + } + } + + GridFree(&outdegree); + GridFree(&indegree); + + return res; +} + +int test_sum_edge_weights(Flow fd) { + // Outgoing edge weights should sum to one + int res = 1; + + Grid edge_sums = GridCreate(GridF32, NULL, fd.cellsize, fd.dims); + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + + ((float *)edge_sums.data)[u] += fd.fraction[e]; + } + + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + + if (((float *)edge_sums.data)[u] != 1.0) { + printf("# Outgoing edges from pixel %ld do not sum to 1.\n", u); + res = 0; + } + } + GridFree(&edge_sums); + + return res; +} + +Flow dinf_flow(Grid dem, int order) { + Flow fd = {0}; + + Grid direction = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + Grid prop = GridCreate(GridF32, NULL, dem.cellsize, dem.dims); + + flow_routing_dinf_directions(direction.data, prop.data, dem.data, dem.dims, + order); + + ptrdiff_t edge_count = edgeset_count(direction.data, dem.dims); + + float *weight = calloc(edge_count, sizeof *weight); + + flow_routing_dinf_weights(weight, direction.data, prop.data, direction.dims); + + // Resolve flats + Grid rdirection = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + Grid resolved = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + + // Initialize resolved according to the flow directions + for (ptrdiff_t j = 0; j < dem.dims[1]; j++) { + for (ptrdiff_t i = 0; i < dem.dims[0]; i++) { + if (isnan(((float *)dem.data)[j * dem.dims[0] + i])) { + // If DEM is a NaN, treat the flow direction as resolved + ((uint8_t *)resolved.data)[j * dem.dims[0] + i] = 0xff; + } else { + ((uint8_t *)resolved.data)[j * dem.dims[0] + i] = + ((uint8_t *)direction.data)[j * dem.dims[0] + i] ? 0xff : 0; + } + } + } + + Grid path_distance = GridCreate(GridF32, NULL, dem.cellsize, dem.dims); + Grid heap = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + Grid back = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + + resolve_flats_shortest_path(rdirection.data, resolved.data, + path_distance.data, heap.data, back.data, + dem.data, dem.dims, order); + GridFree(&resolved); + GridFree(&path_distance); + GridFree(&heap); + GridFree(&back); + + // Compute weights for shortest path resolution + ptrdiff_t rcount = edgeset_count(rdirection.data, rdirection.dims); + float *rweight = calloc(rcount, sizeof *rweight); + + resolve_flats_shortest_path_weights(rweight, rcount); + + // Merge edge sets + edge_count = + edgeset_count_merged(direction.data, rdirection.data, direction.dims); + + float *merged_weights = calloc(edge_count, sizeof *merged_weights); + Grid mergedscan = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + + ptrdiff_t edge_count2 = + edgeset_merge(merged_weights, mergedscan.data, direction.data, weight, + rdirection.data, rweight, direction.dims); + + if (edge_count != edge_count2) { + printf("# Merged edge counts not equal: %ld != %ld", edge_count, + edge_count2); + } + + // Compute the topological sort of the merged edge set. + Grid stream = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + ptrdiff_t *source = calloc(edge_count, sizeof *source); + ptrdiff_t *target = calloc(edge_count, sizeof *target); + float *sorted_weight = calloc(edge_count, sizeof *sorted_weight); + + Grid stack = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + Grid stackdir = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + Grid visited = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + + flow_routing_tsort(stream.data, source, target, sorted_weight, stack.data, + stackdir.data, direction.data, merged_weights, + mergedscan.data, visited.data, edge_count, dem.dims, + order); + + fd.count = edge_count; + fd.stream = stream.data; + fd.source = source; + fd.target = target; + fd.fraction = sorted_weight; + fd.cellsize = dem.cellsize; + fd.dims[0] = dem.dims[0]; + fd.dims[1] = dem.dims[1]; + + GridFree(&stack); + GridFree(&stackdir); + GridFree(&visited); + GridFree(&mergedscan); + GridFree(&direction); + GridFree(&prop); + free(weight); + free(merged_weights); + free(rweight); + GridFree(&rdirection); + + return fd; +} + +int main(int argc, char *argv[]) { + GDALAllRegister(); + + const char *snapshot_directory = ""; + + if (argc >= 2) { + snapshot_directory = argv[1]; + } + + test_init() { + test("Slope") { + ptrdiff_t dims[2] = {3, 3}; + float data[9] = {2, 2, 2, 1, 1, 1, 0, 0, 0}; + Grid dem = GridCreate(GridF32, data, 1.0, dims); + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + } + + test("Larger slope") { + ptrdiff_t dims[2] = {5, 5}; + float data[25] = {4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, + 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0}; + Grid dem = GridCreate(GridF32, data, 1.0, dims); + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + } + + test("Diagonal") { + ptrdiff_t dims[2] = {3, 3}; + float data[9] = {2, 3, 4, 1, 2, 3, 0, 1, 2}; + Grid dem = GridCreate(GridF32, data, 1.0, dims); + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + } + + test("Sink") { + ptrdiff_t dims[2] = {3, 3}; + float data[9] = {2, 2, 2, 2, 1, 2, 2, 2, 2}; + Grid dem = GridCreate(GridF32, data, 1.0, dims); + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + } + + test("Outward sloping cone") { + ptrdiff_t dims[2] = {16, 16}; + float cellsize = 150.0 / (dims[0] - 1); + Grid dem = GridCreate(GridF32, NULL, cellsize, dims); + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + float x = cellsize * i - 75.0; + float y = cellsize * j - 75.0; + ((float *)dem.data)[j * dims[0] + i] = 200 - sqrtf(x * x + y * y); + } + } + + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&dem); + } + + test("Outward sloping cone with flats") { + ptrdiff_t dims[2] = {16, 16}; + float cellsize = 150.0 / (dims[0] - 1); + Grid dem = GridCreate(GridF32, NULL, cellsize, dims); + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + float x = cellsize * i - 75.0; + float y = cellsize * j - 75.0; + ((float *)dem.data)[j * dims[0] + i] = + fminf(150.0, 200 - sqrtf(x * x + y * y)); + } + } + + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&dem); + } + + if (*snapshot_directory) { + test("bigtujunga") { + const char *filepath = "/bigtujunga_30m/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + + test("perfectworld") { + const char *filepath = "/perfectworld/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + + test("tibet") { + const char *filepath = "/tibet/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + + test("kedarnath") { + const char *filepath = "/kedarnath/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + + test("kunashiri") { + const char *filepath = "/kunashiri/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + + test("taiwan") { + const char *filepath = "/taiwan/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + + test("taalvolcano") { + const char *filepath = "/taalvolcano/dem.tif"; + char *path = calloc(strlen(snapshot_directory) + strlen(filepath) + 1, + sizeof *path); + strcpy(path, snapshot_directory); + strncat(path, filepath, strlen(filepath) + 1); + printf("# %s\n", path); + + Grid dem = GridFromFile(path); + free(path); + + Grid demf = GridFillsinks(dem); + Flow fd = dinf_flow(demf, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&demf); + GridFree(&dem); + } + } + + test("nans") { + ptrdiff_t dims[2] = {5, 5}; + ptrdiff_t cellsize = 1.0; + Grid dem = GridCreate(GridF32, NULL, cellsize, dims); + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + float x = cellsize * i; + float y = cellsize * j; + ((float *)dem.data)[j * dims[0] + i] = x + y; + } + } + ((float *)dem.data)[2 * dims[0] + 2] = NAN; + + Flow fd = dinf_flow(dem, 0); + + if (!test_node_once(fd)) test_fail; + if (!test_tsort_node(fd)) test_fail; + if (!test_tsort(fd)) test_fail; + if (!test_nanflow(fd, dem)) test_fail; + if (!test_no_boundary_edges(fd)) test_fail; + if (!test_outlets(fd)) test_fail; + if (!test_sum_edge_weights(fd)) test_fail; + + FlowFree(&fd); + GridFree(&dem); + } + + test("Memory order slope") { + ptrdiff_t dims[2] = {5, 7}; + + Grid fdem = GridCreate(GridF32, NULL, 1.0, dims); + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + float z = j * dims[0] + i; + ((float *)fdem.data)[j * dims[0] + i] = z; + } + } + + Grid cdem = GridTranspose(fdem); + + Grid fdemf = GridFillsinks(fdem); + GridFree(&fdem); + + Grid cdemf = GridFillsinks(cdem); + GridFree(&cdem); + + Flow cfd = dinf_flow(cdemf, 1); + Flow ffd = dinf_flow(fdemf, 0); + + Grid cacc = FlowAccumulation(cfd); + Grid facc = FlowAccumulation(ffd); + Grid caccT = GridTranspose(cacc); + + if (!GridAllClose(caccT, facc, 1e-5)) test_fail; + + GridFree(&caccT); + GridFree(&cacc); + GridFree(&facc); + GridFree(&fdemf); + GridFree(&cdemf); + + FlowFree(&ffd); + FlowFree(&cfd); + } + + test("Memory order saddle (ties)") { + ptrdiff_t dims[2] = {3, 3}; + ptrdiff_t fdims[2] = {dims[0], dims[1]}; + + float fdem_buffer[9] = {2, 2, 0, 2, 1, 2, 0, 2, 2}; + + Grid fdem = GridCreate(GridF32, fdem_buffer, 1.0, fdims); + Grid cdem = GridTranspose(fdem); + + Flow cfd = dinf_flow(cdem, 1); + Flow ffd = dinf_flow(fdem, 0); + + Grid cacc = FlowAccumulation(cfd); + Grid facc = FlowAccumulation(ffd); + Grid caccT = GridTranspose(cacc); + + if (!GridAllClose(caccT, facc, 1e-5)) test_fail; + + GridFree(&caccT); + GridFree(&cacc); + GridFree(&facc); + GridFree(&cdem); + + FlowFree(&ffd); + FlowFree(&cfd); + } + + test("Memory order random") { + ptrdiff_t dims[2] = {19, 17}; + + Grid fdem = GridRandom(0x746f706f7f6f6f6c, 1.0, dims); + + Grid cdem = GridTranspose(fdem); + + Grid fdemf = GridFillsinks(fdem); + GridFree(&fdem); + Grid cdemf = GridFillsinks(cdem); + GridFree(&cdem); + + Flow cfd = dinf_flow(cdemf, 1); + Flow ffd = dinf_flow(fdemf, 0); + + Grid cacc = FlowAccumulation(cfd); + Grid facc = FlowAccumulation(ffd); + + Grid caccT = GridTranspose(cacc); + + if (!GridAllClose(caccT, facc, 1e-5)) test_fail; + + GridFree(&caccT); + GridFree(&cacc); + GridFree(&facc); + + FlowFree(&ffd); + FlowFree(&cfd); + + GridFree(&fdemf); + GridFree(&cdemf); + } + } +} diff --git a/test/filters.cpp b/test/filters.cpp index b4389963..157458c1 100644 --- a/test/filters.cpp +++ b/test/filters.cpp @@ -7,6 +7,8 @@ // TODO: no tests with real 3d SEs performed! +#undef NDEBUG + #include #include #include diff --git a/test/flow.c b/test/flow.c index 7ea50c3f..8417a28e 100644 --- a/test/flow.c +++ b/test/flow.c @@ -79,6 +79,7 @@ Grid FlowAccumulation(Flow f) { Grid Outlets(Flow f) { Grid outlets = GridCreate(GridU8, NULL, f.cellsize, f.dims); + if (outlets.data && f.count > 0) { uint8_t *outletdata = outlets.data; @@ -101,6 +102,7 @@ Grid Outlets(Flow f) { GridFree(&outdegree); } } + return outlets; } @@ -139,5 +141,6 @@ Grid FlowDrainageBasins(Flow f) { free(weights); } } + return basins; } diff --git a/test/grid.c b/test/grid.c index 6ffb780a..1772a411 100644 --- a/test/grid.c +++ b/test/grid.c @@ -13,7 +13,7 @@ void GridFree(Grid *g) { free(g->data); g->data = NULL; g->dims[0] = 0; - g->dims[0] = 0; + g->dims[1] = 0; g->cellsize = 0.0; }