diff options
Diffstat (limited to 'main.c')
| -rw-r--r-- | main.c | 319 |
1 files changed, 227 insertions, 92 deletions
@@ -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"); |