summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIlhan Özgen Xian <iozgen@lbl.gov>2021-06-09 18:56:14 -0700
committerIlhan Özgen Xian <iozgen@lbl.gov>2021-06-09 18:56:14 -0700
commit5c4926bca4e6479ae49521683abe0b0f4d3cfcda (patch)
tree2e25aba227a359a4c64c1daf1ca76ccb755a7846
parentb883a9085fb483c5631ac8d2cd768ae9877d2204 (diff)
add meshing functionality
-rw-r--r--.gitignore4
-rw-r--r--Makefile27
-rw-r--r--header.h229
-rw-r--r--main.c319
4 files changed, 386 insertions, 193 deletions
diff --git a/.gitignore b/.gitignore
index ae1fdfa..5ce281c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,6 @@
output/
*~
mammoth
-*.o \ No newline at end of file
+*.o
+.DS_STORE
+TAGS \ No newline at end of file
diff --git a/Makefile b/Makefile
index 7174c85..6840954 100644
--- a/Makefile
+++ b/Makefile
@@ -1,11 +1,20 @@
#
-# Makefile
+# Makefile for mammoth
#
+# Author:
+# Dr. Ilhan Ozgen-Xian
+# Lawrence Berkeley National Laboratory 2021
+#
+# Type "make" to compile mammoth
+#
+
PRG=mammoth
SRC=main.c
+HDR=header.h
-CC=/usr/local/bin/gcc-10
+CC=/usr/local/bin/gcc-11
+CTAGS=uctags
TRIANGLE_DIR=/Users/IOzgen/Documents/work/coding/tools/triangle
@@ -14,7 +23,7 @@ LFLAGS=-I$(TRIANGLE_DIR)
CPU=
RUNFLAGS=
INCLUDES=
-LIB=
+LIB=$(TRIANGLE_DIR)/triangle.o
#
#
@@ -23,6 +32,7 @@ LIB=
OBJ=$(SRC:.c=.o)
all: $(PRG)
+ make tags
run: $(PRG)
$(RUNFLAGS) ./$(PRG)
@@ -34,9 +44,16 @@ $(PRG): $(OBJ)
$(CC) $(LFLAGS) $(CFLAGS) -c $< -o $@
clean:
- rm -vf $(PRG) $(OBJ) output/*
+ rm -vf $(PRG) $(OBJ)
rm -vf $(P)
-.PHONY: all clean
+cleanall:
+ rm -vf output/*
+ make clean
+
+tags: $(SRC)
+ $(CTAGS) -Re $(SRC) $(HDR)
+
+.PHONY: all clean cleanall tags
# eof
diff --git a/header.h b/header.h
index 3988abb..1391b5b 100644
--- a/header.h
+++ b/header.h
@@ -43,6 +43,27 @@ int multiresolution_meshing(int argc, char** argv);
* ----------------------------------------------------------------------
*/
+int sign(double x)
+{
+ int flag = (x > 0) ? 1 : -1;
+ if (x == 0) {
+ flag = 0;
+ }
+ return flag;
+}
+
+double max(double a, double b)
+{
+ double results = (a > b) ? a : b;
+ return results;
+}
+
+double min(double a, double b)
+{
+ double results = (a < b) ? a : b;
+ return results;
+}
+
int ipow(int base, int exponent)
{
int result = 1;
@@ -116,101 +137,6 @@ int transpose(double* matrix, int nrow, int ncol) {
return 0;
}
-void report(io, markers, reporttriangles, reportneighbors, reportsegments,
- reportedges, reportnorms)
-struct triangulateio *io;
-int markers;
-int reporttriangles;
-int reportneighbors;
-int reportsegments;
-int reportedges;
-int reportnorms;
-{
- int i, j;
-
- for (i = 0; i < io->numberofpoints; i++) {
- printf("Point %4d:", i);
- for (j = 0; j < 2; j++) {
- printf(" %.6g", io->pointlist[i * 2 + j]);
- }
- if (io->numberofpointattributes > 0) {
- printf(" attributes");
- }
- for (j = 0; j < io->numberofpointattributes; j++) {
- printf(" %.6g",
- io->pointattributelist[i * io->numberofpointattributes + j]);
- }
- if (markers) {
- printf(" marker %d\n", io->pointmarkerlist[i]);
- } else {
- printf("\n");
- }
- }
- printf("\n");
-
- if (reporttriangles || reportneighbors) {
- for (i = 0; i < io->numberoftriangles; i++) {
- if (reporttriangles) {
- printf("Triangle %4d points:", i);
- for (j = 0; j < io->numberofcorners; j++) {
- printf(" %4d", io->trianglelist[i * io->numberofcorners + j]);
- }
- if (io->numberoftriangleattributes > 0) {
- printf(" attributes");
- }
- for (j = 0; j < io->numberoftriangleattributes; j++) {
- printf(" %.6g", io->triangleattributelist[i *
- io->numberoftriangleattributes + j]);
- }
- printf("\n");
- }
- if (reportneighbors) {
- printf("Triangle %4d neighbors:", i);
- for (j = 0; j < 3; j++) {
- printf(" %4d", io->neighborlist[i * 3 + j]);
- }
- printf("\n");
- }
- }
- printf("\n");
- }
-
- if (reportsegments) {
- for (i = 0; i < io->numberofsegments; i++) {
- printf("Segment %4d points:", i);
- for (j = 0; j < 2; j++) {
- printf(" %4d", io->segmentlist[i * 2 + j]);
- }
- if (markers) {
- printf(" marker %d\n", io->segmentmarkerlist[i]);
- } else {
- printf("\n");
- }
- }
- printf("\n");
- }
-
- if (reportedges) {
- for (i = 0; i < io->numberofedges; i++) {
- printf("Edge %4d points:", i);
- for (j = 0; j < 2; j++) {
- printf(" %4d", io->edgelist[i * 2 + j]);
- }
- if (reportnorms && (io->edgelist[i * 2 + 1] == -1)) {
- for (j = 0; j < 2; j++) {
- printf(" %.6g", io->normlist[i * 2 + j]);
- }
- }
- if (markers) {
- printf(" marker %d\n", io->edgemarkerlist[i]);
- } else {
- printf("\n");
- }
- }
- printf("\n");
- }
-}
-
/*
* ----------------------------------------------------------------------
* meshing functions
@@ -300,3 +226,116 @@ int get_border(double *data, double nodata, int nrow, int ncol)
return 0;
}
+
+void write_avs(struct triangulateio *io, int counter)
+{
+
+ char fname[1024];
+ sprintf(fname, "output/mesh%d.inp", counter);
+
+ FILE *fp;
+ fp = fopen(fname, "w");
+
+ 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);
+ }
+
+ for (int i = 0; i < io->numberoftriangles; ++ i) {
+ fprintf(fp, "%d %d tri %d %d %d\n", i+1, 1, io->trianglelist[3*i], io->trianglelist[3*i+1], io->trianglelist[3*i+2]);
+ }
+
+ fclose(fp);
+
+}
+
+void report(struct triangulateio *io, int markers, int reporttriangles,
+ int reportneighbors, int reportsegments, int reportedges,
+ int reportnorms)
+{
+
+ int i, j;
+
+ for (i = 0; i < io->numberofpoints; i++) {
+ printf("Point %4d:", i);
+ for (j = 0; j < 2; j++) {
+ printf(" %.6g", io->pointlist[i * 2 + j]);
+ }
+ if (io->numberofpointattributes > 0) {
+ printf(" attributes");
+ }
+ for (j = 0; j < io->numberofpointattributes; j++) {
+ printf(" %.6g",
+ io->pointattributelist[i * io->numberofpointattributes + j]);
+ }
+ if (markers) {
+ printf(" marker %d\n", io->pointmarkerlist[i]);
+ } else {
+ printf("\n");
+ }
+ }
+ printf("\n");
+
+ if (reporttriangles || reportneighbors) {
+ for (i = 0; i < io->numberoftriangles; i++) {
+ if (reporttriangles) {
+ printf("Triangle %4d points:", i);
+ for (j = 0; j < io->numberofcorners; j++) {
+ printf(" %4d", io->trianglelist[i * io->numberofcorners + j]);
+ }
+ if (io->numberoftriangleattributes > 0) {
+ printf(" attributes");
+ }
+ for (j = 0; j < io->numberoftriangleattributes; j++) {
+ printf(" %.6g", io->triangleattributelist[i *
+ io->numberoftriangleattributes + j]);
+ }
+ printf("\n");
+ }
+ if (reportneighbors) {
+ printf("Triangle %4d neighbors:", i);
+ for (j = 0; j < 3; j++) {
+ printf(" %4d", io->neighborlist[i * 3 + j]);
+ }
+ printf("\n");
+ }
+ }
+ printf("\n");
+ }
+
+ if (reportsegments) {
+ for (i = 0; i < io->numberofsegments; i++) {
+ printf("Segment %4d points:", i);
+ for (j = 0; j < 2; j++) {
+ printf(" %4d", io->segmentlist[i * 2 + j]);
+ }
+ if (markers) {
+ printf(" marker %d\n", io->segmentmarkerlist[i]);
+ } else {
+ printf("\n");
+ }
+ }
+ printf("\n");
+ }
+
+ if (reportedges) {
+ for (i = 0; i < io->numberofedges; i++) {
+ printf("Edge %4d points:", i);
+ for (j = 0; j < 2; j++) {
+ printf(" %4d", io->edgelist[i * 2 + j]);
+ }
+ if (reportnorms && (io->edgelist[i * 2 + 1] == -1)) {
+ for (j = 0; j < 2; j++) {
+ printf(" %.6g", io->normlist[i * 2 + j]);
+ }
+ }
+ if (markers) {
+ printf(" marker %d\n", io->edgemarkerlist[i]);
+ } else {
+ printf("\n");
+ }
+ }
+ printf("\n");
+ }
+}
diff --git a/main.c b/main.c
index f59bdad..1a1a9d3 100644
--- a/main.c
+++ b/main.c
@@ -1,4 +1,4 @@
-/* -*- mode: c; style: linux -*- */
+/* -*- mode: c -*- */
/* Purpose: Given input data in the form of a 2D array, perform a
* discrete Haar wavelet transform.
@@ -10,7 +10,7 @@
/* Main logic */
-int main(int argc, char**argv)
+int main(int argc, char** argv)
{
printf("\n\n");
@@ -21,23 +21,23 @@ int main(int argc, char**argv)
printf(" | | | |/ ___ \\| | | | | | | |_| || | | _ |\n");
printf(" |_| |_/_/ \\_\\_| |_|_| |_|\\___/ |_| |_| |_|\n\n");
- printf(" +----------------------------+\n");
- printf(" | |\n");
- printf(" | . ........... |\n");
- printf(" | ................... |\n");
- printf(" | .... .. ........... .. |\n");
- printf(" | ... . .............. . |\n");
- printf(" | .. .............. |\n");
- printf(" | .. .. .. |\n");
- printf(" | .. .. |\n");
- printf(" | |\n");
- printf(" +----------------------------+\n\n");
-
- printf(" I. Ozgen-Xian, G. Kesserwani, D. Caviedes-Voullieme,\n");
- printf(" S. Molins, Z. Xu, D. Dwivedi, J. D. Moulton, C. I. Steefel,\n");
- printf(" Wavelet-based local mesh refinement for rainfall-runoff\n");
- printf(" simulations, Journal of Hydroinformatics, 22, 1059-1077, 2020.\n");
- printf(" doi: 10.2166/hydro.2020.198\n\n");
+ /* printf(" +----------------------------+\n"); */
+ /* printf(" | |\n"); */
+ /* printf(" | . ........... |\n"); */
+ /* printf(" | ................... |\n"); */
+ /* printf(" | .... .. ........... .. |\n"); */
+ /* printf(" | ... . .............. . |\n"); */
+ /* printf(" | .. .............. |\n"); */
+ /* printf(" | .. .. .. |\n"); */
+ /* printf(" | .. .. |\n"); */
+ /* printf(" | |\n"); */
+ /* printf(" +----------------------------+\n\n"); */
+
+ printf(" I. Ozgen-Xian, G. Kesserwani, D. Caviedes-Voullieme,\n");
+ printf(" S. Molins, Z. Xu, D. Dwivedi, J. D. Moulton, C. I. Steefel,\n");
+ printf(" Wavelet-based local mesh refinement for rainfall-runoff\n");
+ printf(" simulations, Journal of Hydroinformatics, 22, 1059-1077, 2020.\n");
+ printf(" doi: 10.2166/hydro.2020.198\n\n");
@@ -161,9 +161,9 @@ int multiresolution_meshing(int argc, char** argv)
boundary_coords[3*counter] = index_x * cellsize + xllcorner;
boundary_coords[3*counter+1] = (nrow - index_y) * cellsize + yllcorner;
boundary_coords[3*counter+2] = z;
- #if VERBOSE
+#if VERBOSE
printf("[ok] Read point %d/%d\n", counter + 1, n_bounds);
- #endif
+#endif
}
fclose(bp);
@@ -193,6 +193,7 @@ int multiresolution_meshing(int argc, char** argv)
levels = calloc(ncol * nrow, sizeof(int));
+
for (int i = 0; i < nrow; i ++) {
for (int j = 0; j < ncol; j ++) {
levels[j + ncol * i] = (int) nodata;
@@ -220,15 +221,19 @@ int multiresolution_meshing(int argc, char** argv)
fclose(lp);
-#if DEBUG_IO
for (int i = 0; i < nrow; i ++) {
for (int j = 0; j < ncol; j ++) {
- levels[j + ncol * i] = _tmp_levels[j + ncol * i];
+ levels[j + ncol * i] = _tmp_levels[j + ncol_padded * i];
+#if DEBUG_IO
printf(" %d ", levels[j + ncol * i]);
+#endif
}
+#if DEBUG_IO
printf("\n");
- }
#endif
+ }
+
+ printf("[OK] Level data read\n");
/*
* ----------------------------------------------------------------------
@@ -236,99 +241,172 @@ int multiresolution_meshing(int argc, char** argv)
* ----------------------------------------------------------------------
*/
- double a0 = cellsize * cellsize;
-
+ double a0 = cellsize * cellsize * ipow(2, nlevels-1);
+
+ printf("[OK] Building mesh with largest cell area %.2f m2 {\n\n", a0);
+
struct triangulateio in;
- struct triangulateio mid;
- struct triangulateio out;
- struct triangulateio vorout;
in.numberofpoints = n_bounds;
in.numberofpointattributes = 1;
-
+ in.numberofsegments = n_bounds;
+ in.numberofholes = 0;
+ in.numberofregions = 0;
+
in.pointlist = (double *) calloc(in.numberofpoints * 2, sizeof(double));
in.pointattributelist = (double *) calloc(in.numberofpoints * in.numberofpointattributes, sizeof(double));
in.pointmarkerlist = (int *) calloc(in.numberofpoints, sizeof(int));
+ in.segmentlist = (int *) calloc(in.numberofsegments * 2, sizeof(int));
+ in.segmentmarkerlist = (int *) calloc(in.numberofsegments, sizeof(int));
+
+ /* Each of the columns of this array contains the constraint's x and
+ y coordinates are at indices [1] and [2], followed by the
+ regional attribute at index [3], followed by the maximum area at
+ index [4]. So we have size(regionlist,1)==4. Note that each
+ regional attribute is used only if you select the 'A' switch, and
+ each area constraint is used only if you select the 'a' switch
+ (with no number following), but omitting one of these switches
+ does not change the memory layout. Input only, although the
+ pointer is copied to the output structure for convenience. */
+ in.regionlist = (double *) calloc(in.numberofregions * 4, sizeof(double));
- in.numberofsegments = 0;
- in.numberofholes = 0;
- in.numberofregions = 1;
-
- for (int i = 0; i < n_bounds; i ++) {
+ for (int i = 0; i < in.numberofpoints; ++ i) {
+ in.pointmarkerlist[i] = 1;
+ //in.pointattributelist[i] = boundary_coords[3 * i + 2]; // set z coordinates
+ for (int j = 0; j < 2; ++ j) {
+ in.pointlist[2 * i + j] = boundary_coords[3 * i + j];
+ }
+ }
- in.pointlist[2*i] = boundary_coords[3*i];
- in.pointlist[2*i+1] = boundary_coords[3*i+1];
- in.pointattributelist[i] = 1.0;
- in.pointmarkerlist[i] = 0;
-
+ for (int i = 0; i < in.numberofsegments-1; ++ i) {
+ in.segmentlist[2*i] = in.numberofsegments-i;
+ in.segmentlist[2*i+1] = in.numberofsegments-1-i;
}
- in.regionlist = (double *) calloc(in.numberofregions * 4, sizeof(double));
- in.regionlist[3] = a0;
+ in.segmentlist[2*in.numberofpoints-2] = 1;
+ in.segmentlist[2*in.numberofpoints-1] = in.numberofsegments;
- /*
- * Make necessary initializations so that Triangle can return a
- * triangulation in `mid' and a voronoi diagram in `vorout'.
- */
+ struct triangulateio out[nlevels];
+
+ out[0].pointlist = (double *) NULL;
+ out[0].pointattributelist = (double *) NULL;
+ out[0].pointmarkerlist = (int *) NULL;
+ out[0].trianglelist = (int *) NULL;
+ out[0].triangleattributelist = (double *) NULL;
+ out[0].neighborlist = (int *) NULL;
+ out[0].segmentlist = (int *) NULL;
+ out[0].segmentmarkerlist = (int *) NULL;
- mid.pointlist = (double *) NULL;
- mid.pointattributelist = (double *) NULL;
- mid.pointmarkerlist = (int *) NULL;
- mid.trianglelist = (int *) NULL;
- mid.triangleattributelist = (double *) NULL;
- mid.segmentlist = (int *) NULL;
- mid.segmentmarkerlist = (int *) NULL;
- mid.edgelist = (int *) NULL;
- mid.edgemarkerlist = (int *) NULL;
+ printf("[ok] Initial triangulation...");
+
+ char triangleargs[1024];
+ sprintf(triangleargs, "QDpa%f", a0);
+ triangulate(triangleargs, &in, &out[0], (struct triangulateio *) NULL);
- vorout.pointlist = (double *) NULL;
- vorout.pointattributelist = (double *) NULL;
- vorout.edgelist = (int *) NULL;
- vorout.normlist = (double *) NULL;
+ printf(" done\n");
- triangulate("pczevn", &in, &mid, &vorout);
+ write_avs(&out[0], 1);
- printf("Initial triangulation:\n\n");
- report(&mid, 1, 1, 1, 1, 1, 0);
- printf("Initial Voronoi diagram:\n\n");
- report(&vorout, 0, 0, 0, 0, 1, 1);
+ for (int current_level = 1; current_level < nlevels; ++ current_level) {
- mid.trianglearealist = (double *) calloc(mid.numberoftriangles, sizeof(double));
+ out[current_level-1].trianglearealist = (double *) calloc(out[current_level-1].numberoftriangles, sizeof(double));
- for (int i = 0; i < mid.numberoftriangles; i ++) {
- mid.trianglearealist[i] = a0;
- }
+ // tag the stuff
+ for (int i = 0; i < out[current_level-1].numberoftriangles; ++ i) {
- out.pointlist = (double *) NULL;
- out.pointattributelist = (double *) NULL;
- out.trianglelist = (int *) NULL;
- out.triangleattributelist = (double *) NULL;
+ int p1 = out[current_level-1].trianglelist[3*i] - 1;
+ int p2 = out[current_level-1].trianglelist[3*i+1] - 1;
+ int p3 = out[current_level-1].trianglelist[3*i+2] - 1;
+
+ double x1 = out[current_level-1].pointlist[2*p1] - xllcorner;
+ double y1 = out[current_level-1].pointlist[2*p1+1] - yllcorner;
+
+ double x2 = out[current_level-1].pointlist[2*p2] - xllcorner;
+ double y2 = out[current_level-1].pointlist[2*p2+1] - yllcorner;
+
+ double x3 = out[current_level-1].pointlist[2*p3] - xllcorner;
+ double y3 = out[current_level-1].pointlist[2*p3+1] - yllcorner;
+
+ double xmin = min(x1, min(x2, x3));
+ double ymin = min(y1, min(y2, y3));
+ double xmax = max(x1, max(x2, x3));
+ double ymax = max(y1, max(y2, y3));
+
+ int ixmin = (int) floor(xmin / cellsize);
+ int iymin = (int) floor(ymin / cellsize);
+ int ixmax = (int) floor(xmax / cellsize);
+ int iymax = (int) floor(ymax / cellsize);
+
+ /* double xc = (x1 + x2 + x3) / 3.0; */
+ /* double yc = (y1 + y2 + y3) / 3.0; */
+ /* double a = 0.5 * (- y2 * x3 + y1 * (-x2 + x3) + x1 * (y2 - y3) + x2 * y3); */
+ /* double c = 1.0 / (2.0 * a); */
+
+ /* double f1 = y1 * x2; */
+ /* double f2 = x1 * y3; */
+ /* double f3 = y3 - y1; */
+ /* double f4 = x1 - x3; */
+
+ /* double g1 = x1 * y2; */
+ /* double g2 = y1 * x2; */
+ /* double g3 = y1 - y2; */
+ /* double g4 = x2 - x1; */
+
+ int flag = -1;
+
+ /* printf("[ok] Triangle #%d/%d\n", i+1, mid.numberoftriangles); */
+ /* printf("\t %d %d %d / %d\n", p1, p2, p3, mid.numberofpoints); */
+ /* printf("\t %d %d : %d %d\n", ixmin, ixmax, iymin, iymax); */
+ /* printf("\t %f %f : %f %f\n", xmin, xmax, ymin, ymax); */
+
+ for (int iy = iymin; iy < iymax; ++ iy) {
+ for (int ix = ixmin; ix < ixmax; ++ ix) {
+ /* printf("%d < %d?\n", flag, levels[ix + ncol * iy]); */
+ flag = (flag < levels[ix + ncol * (nrow - iy)]) ? levels[ix + ncol * (nrow -iy)] : flag;
+ }
+ }
+
+ /* printf("\t flag: %d...", flag); */
+
+ if (flag >= current_level) {
+ /* printf(" tagged\n"); */
+ out[current_level-1].trianglearealist[i] = a0 / ipow(2, current_level);
+ }
+
+ }
+
+ out[current_level].pointlist = (double *) NULL;
+ out[current_level].pointattributelist = (double *) NULL;
+ out[current_level].pointmarkerlist = (int *) NULL;
+ out[current_level].trianglelist = (int *) NULL;
+ out[current_level].triangleattributelist = (double *) NULL;
+ out[current_level].segmentlist = (int *) NULL;
+ out[current_level].segmentmarkerlist = (int *) NULL;
+
+ printf("[ok] Triangulation #%d ...", current_level + 1);
+
+ triangulate("QDarp", &out[current_level-1], &out[current_level], (struct triangulateio *) NULL);
+
+ printf(" done\n");
+
+ write_avs(&out[current_level], current_level+1);
+
+ }
+
+ printf("[OK] } done\n");
- triangulate("prazB", &mid, &out, (struct triangulateio *) NULL);
free(in.pointlist);
free(in.pointattributelist);
free(in.pointmarkerlist);
free(in.regionlist);
- free(mid.pointlist);
- free(mid.pointattributelist);
- free(mid.pointmarkerlist);
- free(mid.trianglelist);
- free(mid.triangleattributelist);
- free(mid.trianglearealist);
- free(mid.neighborlist);
- free(mid.segmentlist);
- free(mid.segmentmarkerlist);
- free(mid.edgelist);
- free(mid.edgemarkerlist);
- free(vorout.pointlist);
- free(vorout.pointattributelist);
- free(vorout.edgelist);
- free(vorout.normlist);
- free(out.pointlist);
- free(out.pointattributelist);
- free(out.trianglelist);
- free(out.triangleattributelist);
+
+ free(out[0].pointlist);
+ free(out[0].pointattributelist);
+ free(out[0].pointmarkerlist);
+ free(out[0].trianglelist);
+ free(out[0].triangleattributelist);
+ free(out[0].neighborlist);
free(levels);
free(boundary_coords);
@@ -686,6 +764,63 @@ int multiresolution_analysis(int argc, char** argv)
printf("done\n");
+/* FILE *bp; */
+/* FILE *op; */
+
+/* bp = fopen("output/boundary.asc", "r"); */
+
+/* int index_x; */
+/* int index_y; */
+/* double z; */
+
+/* int n_bounds = 0; */
+/* char _tmp_char; */
+
+/* while (!feof(bp)) { */
+/* _tmp_char = fgetc(bp); */
+/* if (_tmp_char == '\n') { */
+/* n_bounds ++; */
+/* } */
+/* } */
+
+/* printf("[ok] Found %d points\n", n_bounds); */
+
+/* rewind(bp); */
+
+/* double *boundary_coords = calloc(3 * n_bounds, sizeof(double)); */
+
+/* for (int counter = 0; counter < n_bounds; counter ++) { */
+/* fscanf(bp, "%d %d %lf", &index_x, &index_y, &z); */
+/* boundary_coords[3*counter] = index_x * cellsize + xllcorner; */
+/* boundary_coords[3*counter+1] = (nrow - index_y) * cellsize + yllcorner; */
+/* boundary_coords[3*counter+2] = z; */
+/* #if VERBOSE */
+/* printf("[ok] Read point %d/%d\n", counter + 1, n_bounds); */
+/* #endif */
+/* } */
+
+/* fclose(bp); */
+
+/* op = fopen("output/boundary.poly", "w"); */
+
+/* fprintf(op, "%d %d %d %d\n", n_bounds, 2, 1, 0); */
+
+/* for (int i = 0; i < n_bounds; ++ i) { */
+/* fprintf(op, "%d %f %f %f\n", i+1, boundary_coords[3*i], boundary_coords[3*i+1], boundary_coords[3*i+2]); */
+/* } */
+
+/* fprintf(op, "%d %d\n", n_bounds, 0); */
+
+/* for (int i = 0; i < n_bounds-1; ++ i) { */
+/* fprintf(op, "%d %d %d\n", i+1, n_bounds-i, n_bounds-1-i); */
+/* } */
+
+/* fprintf(op, "%d %d %d\n0 0\n", n_bounds, 1, n_bounds); */
+
+/* fclose(op); */
+
+/* free(boundary_coords); */
+
//FILE *dp;
//
//dp = fopen("output/boundary.asc", "r");