From f8ed18ed73902af0c27395c9f2da3e85fa185b6d Mon Sep 17 00:00:00 2001
From: Ilhan Özgen Xian <iozgen@lbl.gov>
Date: Mon, 19 Jul 2021 15:41:13 -0700
Subject: add elevation interpolation

---
 header.h         | 43 +++++++++++++++++++++++++++++++++++++++++--
 input/test.input |  2 +-
 main.c           | 25 ++++++++++++++++++++-----
 3 files changed, 62 insertions(+), 8 deletions(-)

diff --git a/header.h b/header.h
index 1391b5b..cfb4162 100644
--- a/header.h
+++ b/header.h
@@ -227,7 +227,7 @@ int get_border(double *data, double nodata, int nrow, int ncol)
   return 0;
 }
 
-void write_avs(struct triangulateio *io, int counter)
+void write_avs(struct triangulateio *io, int counter, double *elevations, int ncol, int nrow, double xllcorner, double yllcorner, double cellsize, double nodata)
 {
 
   char fname[1024];
@@ -239,7 +239,46 @@ void write_avs(struct triangulateio *io, int counter)
   fprintf(fp, "%d %d %d %d %d\n", io->numberofpoints, io->numberoftriangles, 0, 1, 0);
 
   for (int i = 0; i < io->numberofpoints; ++ i) {
-    fprintf(fp, "%d %f %f %f\n", i+1, io->pointlist[2*i], io->pointlist[2*i+1], 0.0);
+
+    /* map elevation to the grid */
+    double x = io->pointlist[2*i];
+    double y = io->pointlist[2*i+1];
+
+    int ix = (int) floor((x - xllcorner) / cellsize);
+    int iy = nrow - (int) floor((y - yllcorner) / cellsize);
+
+    double z = elevations[ix + ncol * iy];
+    if (z == nodata) {
+      int ix0 = (ix + 1) < ncol ? ix + 1 : ix;
+      int ix1 = (ix - 1) > 0 ? ix - 1 : ix;
+      int iy0 = (iy + 1) < nrow ? iy + 1 : iy;
+      int iy1 = (iy - 1) > 0 ? iy - 1 : iy;
+
+      /* indices of von neumann neighborhood */
+      int neumann[4] = {ix0 + ncol * iy, ix1 + ncol * iy, ix + ncol * iy0, ix + ncol * iy1};
+
+      for (int j = 0; j < 4; j ++) {
+	double znew = elevations[neumann[j]];
+	if (znew != nodata) {
+	  z = znew;
+	  break;
+	}
+      }
+
+      int moore[4] = {ix0 + ncol * iy0, ix0 + ncol * iy1, ix1 + ncol * iy0, ix1 + ncol * iy1};
+
+      for (int j = 0; j < 4; j ++) {
+	double znew = elevations[moore[j]];
+	if (znew != nodata) {
+	  z = znew;
+	  break;
+	}
+      }
+      
+    }
+    
+    fprintf(fp, "%d %f %f %f\n", i+1, io->pointlist[2*i], io->pointlist[2*i+1], z);
+
   }
 
   for (int i = 0; i < io->numberoftriangles; ++ i) {
diff --git a/input/test.input b/input/test.input
index de91967..2141065 100644
--- a/input/test.input
+++ b/input/test.input
@@ -1,3 +1,3 @@
 datafile  input/dem.asc
 N         8
-threshold 0.001
\ No newline at end of file
+threshold 0.001
diff --git a/main.c b/main.c
index 1a1a9d3..82ac8f7 100644
--- a/main.c
+++ b/main.c
@@ -49,10 +49,11 @@ int main(int argc, char** argv)
     } else if (argv[1][1] == 'm') {
       printf("[OK] Generate mesh\n");
       multiresolution_meshing(argc, argv);
-    } else {
-      fprintf(stderr, "[ERR] Flag not understood: %s\n", argv[1]);
     }
     
+  } else {
+    fprintf(stderr, "[ERR] Flag not understood: %s\n", argv[1]);
+    fprintf(stderr, "[ERR] Usage: ./mammoth -[a|m] <input/file>\n");
   }
   
   return 0;
@@ -81,6 +82,7 @@ int multiresolution_meshing(int argc, char** argv)
   int n_bounds;
   int *levels;
   double *boundary_coords;
+  double *elevations;
   
   if (argc < 2) {
     fprintf(stderr, "[ERR] No input file provided.\n");
@@ -121,6 +123,16 @@ int multiresolution_meshing(int argc, char** argv)
 
   printf("[OK] Metadata of %s: nrow: %d, ncol: %d, xll: %lf, yll: %lf, dx: %lf, nodata: %lf\n", datafilename, nrow, ncol, xllcorner, yllcorner, cellsize, nodata);
 
+  /*
+   * ----------------------------------------------------------------------
+   * read in dem data
+   * ----------------------------------------------------------------------
+   */
+
+  elevations = calloc(nrow * ncol, sizeof(double));
+  get_data(datafilename, nrow, ncol, ncol, elevations);
+  printf("[OK] DEM data read\n");
+
   /*
    * ----------------------------------------------------------------------
    * read in boundary
@@ -305,7 +317,7 @@ int multiresolution_meshing(int argc, char** argv)
 
   printf(" done\n");
 
-  write_avs(&out[0], 1);
+  write_avs(&out[0], 1, elevations, ncol, nrow, xllcorner, yllcorner, cellsize, nodata);
 
   for (int current_level = 1; current_level < nlevels; ++ current_level) {
 
@@ -389,7 +401,7 @@ int multiresolution_meshing(int argc, char** argv)
     
     printf(" done\n");
 
-    write_avs(&out[current_level], current_level+1);
+    write_avs(&out[current_level], current_level+1, elevations, ncol, nrow, xllcorner, yllcorner, cellsize, nodata);
 
   }
   
@@ -409,6 +421,7 @@ int multiresolution_meshing(int argc, char** argv)
   free(out[0].neighborlist);
   
   free(levels);
+  free(elevations);
   free(boundary_coords);
   
   return 0;
@@ -949,8 +962,10 @@ int get_data(char *datafilename, int nrow, int ncol, int ncol_padded, double *da
     if (strcmp(key, "nodata_value") == 0) {
 
       fscanf(dp, "%lf", &nodata);
+#if VERBOSE
       printf("[OK] Starting to read data\n");
-
+#endif
+      
       for (int i = 0; i < nrow; i ++) {
 	for (int j = 0; j < ncol; j ++) {
 	  fscanf(dp, "%lf", &data[j + ncol_padded * i]);
-- 
cgit v1.2.3-13-gbd6f