Use cmake, tail recursions, fix missing header

.gitignore

@@ -0,0 +1,3 @@
+perf.data*
+build/
+.vscode/

CMakeLists.txt

@@ -0,0 +1,63 @@
+cmake_minimum_required (VERSION 3.9)
+cmake_policy(SET CMP0077 NEW)
+
+project(meshed
+  VERSION 1.0.0
+  LANGUAGES C
+)
+
+option(USE_ASAN "Enable or disable address sanitizer" OFF)
+option(USE_USAN "Enable or disable undefined sanitizer" OFF)
+
+find_package(GDAL REQUIRED)
+find_package(OpenMP REQUIRED)
+
+set(CMAKE_C_STANDARD 99)
+
+set(WARNING_FLAGS
+  -Waddress
+  -Wall
+  -Wdeprecated
+  -Wextra
+  -Wimplicit-fallthrough
+  -Wmisleading-indentation
+  -Wmissing-prototypes
+  -Wshadow
+  -Wshift-sign-overflow
+  -Wswitch
+  -Wuninitialized
+  -Wunreachable-code
+  -Wunused
+  -Wunused-but-set-variable
+  -Wunused-function
+  # -Werror
+)
+
+add_executable(meshed.exe
+  delineate.c
+  delineate_lessmem.c
+  delineate_moremem.c
+  gettimeofday.c
+  hierarchy.c
+  main.c
+  outlet_list.c
+  outlets.c
+  raster.c
+  timeval_diff.c
+)
+
+target_link_libraries(meshed.exe PRIVATE GDAL::GDAL OpenMP::OpenMP_C m)
+target_compile_options(meshed.exe PRIVATE ${WARNING_FLAGS})
+set_property(TARGET meshed.exe PROPERTY INTERPROCEDURAL_OPTIMIZATION TRUE)
+if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
+  target_compile_options(meshed.exe PRIVATE -foptimize-sibling-calls)
+endif()
+if(USE_ASAN)
+  target_compile_options(meshed.exe PRIVATE -fsanitize=address)
+  target_link_options(meshed.exe PRIVATE -fsanitize=address)
+endif()
+
+if(USE_USAN)
+  target_compile_options(meshed.exe PRIVATE -fsanitize=undefined)
+  target_link_options(meshed.exe PRIVATE -fsanitize=undefined)
+endif()

README.md

@@ -13,12 +13,17 @@ Flow direction encoding in GeoTIFF:<br>
 
 * C compiler with [OpenMP](https://www.openmp.org/) support
 * [GDAL](https://gdal.org/)
+* [Cmake](https://cmake.org/)
 
 For Windows, use [MSYS2](https://www.msys2.org/) and [OSGeo4W](https://trac.osgeo.org/osgeo4w/) to install [GCC](https://gcc.gnu.org/) and [GDAL](https://gdal.org/), respectively.
 
 ## How to compile MESHED
 
+Use the standard cmake sequence:
 ```bash
+mkdir build
+cd build
+cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo ..
 make
 ```
 

delineate_funcs.h

@@ -163,7 +163,7 @@ static void trace_up(struct raster_map *dir_map, int row, int col, int id,
 
     /* use gcc -O2 or -O3 flags for tail-call optimization
      * (-foptimize-sibling-calls) */
-    trace_up(dir_map, row_next, col_next, id, up_stack);
+    MUST_TAIL return trace_up(dir_map, row_next, col_next, id, up_stack);
 }
 
 static void init_up_stack(struct cell_stack *up_stack)

global.h

@@ -83,4 +83,10 @@ struct hierarchy *analyze_hierarchy(struct raster_map *,
                                     struct outlet_list *);
 int write_hierarchy(const char *, struct hierarchy *);
 
+#if defined(__clang_major__) && __clang_major__ >= 13
+#define MUST_TAIL __attribute__((musttail))
+#else
+#define MUST_TAIL
+#endif
+
 #endif

raster.c

@@ -228,7 +228,7 @@ struct raster_map *read_raster(const char *path, int type, int get_stats)
 #pragma omp parallel
     {
 #pragma omp single
-        datasets = malloc(sizeof *datasets * omp_get_num_threads());
+        datasets = malloc(sizeof(GDALDatasetH*) * omp_get_num_threads());
         datasets[omp_get_thread_num()] = GDALOpen(path, GA_ReadOnly);
     }
 
@@ -246,7 +246,7 @@ struct raster_map *read_raster(const char *path, int type, int get_stats)
 #pragma omp parallel
     {
 #pragma omp single
-        bands = malloc(sizeof *bands * omp_get_num_threads());
+        bands = malloc(sizeof(GDALRasterBandH*) * omp_get_num_threads());
         bands[omp_get_thread_num()] =
             GDALGetRasterBand(datasets[omp_get_thread_num()], 1);
     }
@@ -321,6 +321,9 @@ struct raster_map *read_raster(const char *path, int type, int get_stats)
     if (error)
         return NULL;
 
+    free(bands);
+    free(datasets);
+
     return rast_map;
 }
 
@@ -331,8 +334,6 @@ int write_raster(const char *path, struct raster_map *rast_map, int type)
     GDALDatasetH dataset;
     GDALRasterBandH band;
     GDALDataType data_type, gdt_type;
-    size_t row_size;
-    int row;
 
     if (!driver)
         return 1;
@@ -342,9 +343,7 @@ int write_raster(const char *path, struct raster_map *rast_map, int type)
         return 2;
 
     if (rast_map->compress)
-        options = CSLSetNameValue(options, "COMPRESS", "LZW");
-
-    row_size = rast_map->ncols;
+        options = CSLSetNameValue(options, "COMPRESS", "ZSTD");
 
     /* actual data size */
     switch (rast_map->type) {
@@ -364,7 +363,6 @@ int write_raster(const char *path, struct raster_map *rast_map, int type)
         data_type = GDT_Byte;
         break;
     }
-    row_size *= GDALGetDataTypeSizeBytes(data_type);
 
     /* requested data type */
     gdt_type = data_type;
@@ -399,13 +397,13 @@ int write_raster(const char *path, struct raster_map *rast_map, int type)
     band = GDALGetRasterBand(dataset, 1);
     GDALSetRasterNoDataValue(band, rast_map->null_value);
 
-    for (row = 0; row < rast_map->nrows; row++) {
-        if (GDALRasterIO
-            (band, GF_Write, 0, row, rast_map->ncols, 1,
-             (char *)rast_map->cells.v + row * row_size, rast_map->ncols,
-             1, data_type, 0, 0) != CE_None)
-            return 4;
-    }
+    const int ret = GDALRasterIO(
+        band, GF_Write, 0, 0, rast_map->ncols, rast_map->nrows,
+        (char *)rast_map->cells.v, rast_map->ncols,
+        rast_map->nrows, data_type, 0, 0);
+
+    if (ret != CE_None)
+        return 4;
 
     GDALClose(dataset);
 

timeval_diff.c

@@ -1,4 +1,6 @@
 #include <stdlib.h>
+#include "global.h"
+
 #ifdef _MSC_VER
 #include <winsock2.h>
 #else