From b883a9085fb483c5631ac8d2cd768ae9877d2204 Mon Sep 17 00:00:00 2001 From: Ilhan Özgen Xian Date: Mon, 26 Apr 2021 13:45:33 -0700 Subject: add triangle mesh generation -- not working --- Makefile | 4 +- header.h | 113 +++++++++++++++++- main.c | 396 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---- 3 files changed, 487 insertions(+), 26 deletions(-) diff --git a/Makefile b/Makefile index 343d134..7174c85 100644 --- a/Makefile +++ b/Makefile @@ -7,8 +7,10 @@ SRC=main.c CC=/usr/local/bin/gcc-10 +TRIANGLE_DIR=/Users/IOzgen/Documents/work/coding/tools/triangle + CFLAGS=-g -Wall -Wextra -pedantic -std=c99 -LFLAGS= +LFLAGS=-I$(TRIANGLE_DIR) CPU= RUNFLAGS= INCLUDES= diff --git a/header.h b/header.h index 2a67546..3988abb 100644 --- a/header.h +++ b/header.h @@ -3,10 +3,22 @@ #include #include +#ifdef SINGLE +#define REAL float +#else +#define REAL double +#endif + +#include "triangle.h" + #ifndef DEBUG_IO #define DEBUG_IO 0 #endif +#ifndef VERBOSE +#define VERBOSE 0 +#endif + /* * ---------------------------------------------------------------------- * wavelet functions @@ -23,6 +35,10 @@ int decompress(double *m, int nrow, int ncol); int hard_threshold(double *m, int *flag, int n, int nlevels, int nrow, int ncol, double threshold); +int multiresolution_analysis(int argc, char** argv); + +int multiresolution_meshing(int argc, char** argv); + /* * ---------------------------------------------------------------------- */ @@ -100,6 +116,101 @@ 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 @@ -186,6 +297,6 @@ int get_border(double *data, double nodata, int nrow, int ncol) } fclose(fp); - + return 0; } diff --git a/main.c b/main.c index fad6af2..f59bdad 100644 --- a/main.c +++ b/main.c @@ -10,28 +10,11 @@ /* Main logic */ -int main(int argc, char** argv) +int main(int argc, char**argv) { - FILE *fp; - - char datafilename[255]; - int nlevels; - double threshold; - - int ncol; - int nrow; - - double xllcorner; - double yllcorner; - double cellsize; - double nodata; - - double *data; - printf("\n\n"); - printf(" __ __ _ __ __ __ __ ___ _____ _ _ \n"); printf(" | \\/ | / \\ | \\/ | \\/ |/ _ \\_ _| | | |\n"); printf(" | |\\/| | / _ \\ | |\\/| | |\\/| | | | || | | |_| |\n"); @@ -56,17 +39,334 @@ int main(int argc, char** argv) printf(" simulations, Journal of Hydroinformatics, 22, 1059-1077, 2020.\n"); printf(" doi: 10.2166/hydro.2020.198\n\n"); + + + if (argv[1][0] == '-') { + + if (argv[1][1] == 'a') { + printf("[OK] Carry out multiresolution analysis\n"); + multiresolution_analysis(argc, 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]); + } + + } + + return 0; + +} + +int multiresolution_meshing(int argc, char** argv) +{ + + FILE *fp; + FILE *lp; + FILE *bp; + + char datafilename[255]; + int nlevels; + double threshold; + + int ncol; + int nrow; + + double xllcorner; + double yllcorner; + double cellsize; + double nodata; + + int n_bounds; + int *levels; + double *boundary_coords; + + if (argc < 2) { + fprintf(stderr, "[ERR] No input file provided.\n"); + return -1; + } + + printf("[OK] Reading in: %s\n", argv[2]); + + fp = fopen(argv[2], "r"); + + if (fp == NULL) { + fprintf(stderr, "[ERR] Couldn't read input file: %s\n", argv[2]); + return -1; + } + + char key[1024]; + + while (fscanf(fp, "%s", key) != EOF) { + + if (strcmp(key, "datafile") == 0) { + fscanf(fp, "%s", datafilename); + } else if (strcmp(key, "N") == 0) { + fscanf(fp, "%d", &nlevels); + } else if (strcmp(key, "threshold") == 0) { + fscanf(fp, "%lf", &threshold); + } else { + printf("[!!] Key not recognised: %s\n", key); + printf("[!!] Continue until known key is encountered\n"); + } + + } + + fclose(fp); + + if (get_metadata(datafilename, &nrow, &ncol, &xllcorner, &yllcorner, &cellsize, &nodata)) { + return -1; + } + + 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 boundary + * ---------------------------------------------------------------------- + */ + + printf("[OK] Reading boundary points {\n"); + + bp = fopen("output/boundary.asc", "r"); + + if (bp == NULL) { + fprintf(stderr, "[ERR] No boundary file found. Try running multiresolution analysis first.\n"); + return -1; + } + + int index_x; + int index_y; + double z; + + 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); + + 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); + +#if DEBUG_IO + FILE *boundary_debug_file; + boundary_debug_file = fopen("output/_tmp_boundary.asc", "w"); + for (int i = 0; i < n_bounds; i ++) { + fprintf(boundary_debug_file, "%lf %lf %lf \n", boundary_coords[3*i], boundary_coords[3*i+1], boundary_coords[3*i+2]); + } + fclose(boundary_debug_file); +#endif + + printf("[OK] } done\n"); + + /* + * ---------------------------------------------------------------------- + * read in level data + * ---------------------------------------------------------------------- + */ + + int powx = (int) ceil(log2(ncol)); + int powy = (int) ceil(log2(nrow)); + + int ncol_padded = ipow(2, powx); + int nrow_padded = ipow(2, powy); + + 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; + } + } + + int *_tmp_levels = calloc(ncol_padded * nrow_padded, sizeof(int)); + + lp = fopen("output/l.asc", "r"); + + if (lp == NULL) { + fprintf(stderr, "[ERR] No level file found. Try running multiresolution analysis first.\n"); + return -1; + } + + while(fscanf(lp, "%s", key) != EOF) { + + for (int i = 0; i < nrow_padded; i ++) { + for (int j = 0; j < ncol_padded; j ++) { + fscanf(lp, "%d", &_tmp_levels[j + ncol_padded * i]); + } + } + + } + + 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]; + printf(" %d ", levels[j + ncol * i]); + } + printf("\n"); + } +#endif + + /* + * ---------------------------------------------------------------------- + * build the mesh + * ---------------------------------------------------------------------- + */ + + double a0 = cellsize * cellsize; + + struct triangulateio in; + struct triangulateio mid; + struct triangulateio out; + struct triangulateio vorout; + + in.numberofpoints = n_bounds; + in.numberofpointattributes = 1; + + 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.numberofsegments = 0; + in.numberofholes = 0; + in.numberofregions = 1; + + for (int i = 0; i < n_bounds; i ++) { + + 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; + + } + + in.regionlist = (double *) calloc(in.numberofregions * 4, sizeof(double)); + in.regionlist[3] = a0; + + /* + * Make necessary initializations so that Triangle can return a + * triangulation in `mid' and a voronoi diagram in `vorout'. + */ + + 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; + + vorout.pointlist = (double *) NULL; + vorout.pointattributelist = (double *) NULL; + vorout.edgelist = (int *) NULL; + vorout.normlist = (double *) NULL; + + triangulate("pczevn", &in, &mid, &vorout); + + 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); + + mid.trianglearealist = (double *) calloc(mid.numberoftriangles, sizeof(double)); + + for (int i = 0; i < mid.numberoftriangles; i ++) { + mid.trianglearealist[i] = a0; + } + + out.pointlist = (double *) NULL; + out.pointattributelist = (double *) NULL; + out.trianglelist = (int *) NULL; + out.triangleattributelist = (double *) NULL; + + 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(levels); + free(boundary_coords); + + return 0; + +} + +int multiresolution_analysis(int argc, char** argv) +{ + + FILE *fp; + + char datafilename[255]; + int nlevels; + double threshold; + + int ncol; + int nrow; + + double xllcorner; + double yllcorner; + double cellsize; + double nodata; + + double *data; + if (argc < 2) { fprintf(stderr, "[ERR] No input file provided.\n"); return -1; } - printf("[OK] Reading in: %s\n", argv[1]); + printf("[OK] Reading in: %s\n", argv[2]); - fp = fopen(argv[1], "r"); + fp = fopen(argv[2], "r"); if (fp == NULL) { - fprintf(stderr, "[ERR] Couldn't read input file: %s\n", argv[1]); + fprintf(stderr, "[ERR] Couldn't read input file: %s\n", argv[2]); return -1; } @@ -367,6 +667,7 @@ int main(int argc, char** argv) } } + #if DEBUG_IO FILE *_tmp; _tmp = fopen("_tmp.asc", "w"); @@ -376,10 +677,53 @@ int main(int argc, char** argv) } fprintf(_tmp, "\n"); } + fclose(_tmp); + #endif + printf("[OK] Get boundary using square tracing... "); get_border(_data, nodata, nrow, ncol); + printf("done\n"); + + //FILE *dp; + // + //dp = fopen("output/boundary.asc", "r"); + // + //int index_x; + //int index_y; + //double _tmp; + // + //while (!feof(dp)) { + // + // fscanf(dp, "%d %d %lf", &index_x, &index_y, &_tmp); + // struct node p = {.x = index_x * cellsize + xllcorner, .y = index_y * cellsize + yllcorner, .next = NULL}; + // push(&bcpoints, &p); + // + //} + // + //fclose(dp); + // + //printf("[OK] Boundary polygon constructed with %ld points\n", bcpoints.len); + // + //struct node *tmp = pop(&bcpoints); + // + //printf("[OK] %lf %lf\n", tmp->x, tmp->y); + // + /* size_t len_xybc = bcpoints.len; */ + /* xybc = (double *) calloc(2 * len_xybc, sizeof(double)); */ + + /* printf("[OK] Processing points: "); */ + /* for (size_t i = 0; i < len_xybc; i ++) { */ + /* printf(" %ld ", i); */ + /* struct node *p = pop(&bcpoints); */ + /* xybc[2*i] = p->x * cellsize + xllcorner; */ + /* xybc[2*i+1] = p->y * cellsize + yllcorner; */ + /* } */ + /* printf("\n"); */ + + /* */ + /* * ---------------------------------------------------------------------- * Housekeeping @@ -392,6 +736,7 @@ int main(int argc, char** argv) fclose(lptr); + /* free(xybc); */ free(levels); free(flag); free(_data); @@ -450,7 +795,8 @@ int get_metadata(char *datafilename, int *nrow, int *ncol, double *xllcorner, do return 0; } -int get_data(char *datafilename, int nrow, int ncol, int ncol_padded, double *data) { +int get_data(char *datafilename, int nrow, int ncol, int ncol_padded, double *data) +{ FILE *dp; dp = fopen(datafilename, "r"); @@ -545,7 +891,8 @@ int compress(double* data, double *m, int nrow, int ncol) return 0; } -int decompress(double *data, int nrow, int ncol) { +int decompress(double *data, int nrow, int ncol) +{ printf("[OK] Decompressing a %d-by-%d matrix\n", nrow, ncol); @@ -579,7 +926,8 @@ int decompress(double *data, int nrow, int ncol) { return 0; } -int hard_threshold(double *m, int *flag, int n, int nlevels, int nrow, int ncol, double threshold) { +int hard_threshold(double *m, int *flag, int n, int nlevels, int nrow, int ncol, double threshold) +{ /* * Hard thresholding as explained in (Caviedes-Voullieme and * Kesserwani, 2015: Benchmarking a multi-resolution discontinuous -- cgit v1.2.3-13-gbd6f