#define PROG_NAME "test_integrate_iterative" #define PROG_DESC "test of {pst_integrate_iterative}" #define PROG_VERS "1.0" #define test_integrate_iterative_C_COPYRIGHT "Copyright © 2005 by the State University of Campinas (UNICAMP)" /* Last edited on 2025-03-01 18:32:05 by stolfi */ #define PROG_HELP \ " " PROG_NAME " \\\n" \ " { -slopes {G_FNI_NAME} | -normals {N_FNI_NAME} } \\\n" \ " -reference {REF_Z_FNI_NAME} \\\n" \ " -initial {INIT_OPT} {INIT_NOISE} \\\n" \ " [ -hints {H_FNI_NAME} {WH_MULT} ] \\\n" \ " [ -maxIter {MAX_ITER} ] \\\n" \ " [ -convTol [CONV_TOL} ] \\\n" \ " [ -sortSys [TOPOSORT} ] \\\n" \ " [ -verbose ] \\\n" \ " [ -reportStep {REPORT_STEP} ] \\\n" \ " -outPrefix {PREFIX} \\\n" \ " " argparser_help_info_HELP "" #define PROG_INFO \ "NAME\n" \ " " PROG_NAME " - " PROG_DESC "\n" \ "\n" \ "SYNOPSIS\n" \ PROG_HELP "\n" \ "\n" \ "DESCRIPTION\n" \ PROG_INFO_DESC "\n" \ "\n" \ "OUTPUT FILES\n" \ PROG_INFO_OUT_FILES "\n" \ "\n" \ "OPTIONS\n" \ PROG_INFO_OPTS "\n" \ "\n" \ argparser_help_info_HELP_INFO "\n" \ "\n" \ "SEE ALSO\n" \ " fni_to_pnm(1), pnm_to_fni(1).\n" \ "\n" \ "AUTHOR\n" \ " Created 2005-08-15 by Jorge Stolfi, UNICAMP.\n" \ "MODIFICATION HISTORY\n" \ " By J. Stolfi unless said otherwise.\n" \ " 2008-11-08 by J. Stolfi, IC-UNICAMP: added the weight map option.\n" \ " 2008-11-09 by J. Stolfi, IC-UNICAMP: added the reference height map.\n" \ "\n" \ " 2010-04-28 Computes a {Z} confidence map {OW} from the {G} confidence map {IW}.\n" \ " 2010-04-28 \"-reference\" also writes one-line files with error summaries.\n" \ " 2010-04-28 Explicit \"-slopes\" and \"-weights\" options.\n" \ " 2010-04-28 Removed \"-debugPrefix\" added \"-outPrefix\".\n" \ " 2010-04-28 The debug prefix is {PREFIX} plus \"-dbg\".\n" \ " 2010-04-28 Output {Z} goes to file, not to standard output.\n" \ "\n" \ " 2010-05-04 Added \"-sortSys\".\n" \ "\n" \ " 2025-01-07 Option \"-keepNull\" to keep/excludes zero-weight equations.\n" \ " 2025-01-07 Option \"-initial\" to give optional initial guess.\n" \ "\n" \ " 2025-02-13 Option \"-initial\" replaced by \"-hints\".\n" \ " 2025-02-13 Uses {H} as initial guess. Added the hints weight factor {hintsWeight}.\n" \ " 2025-02-13 Removed separate {Z} weight map {U}.\n" \ "\n" \ " 2025-02-18 Removed \"-keepNull\".\n" \ " 2025-02-18 Made \"-reference\" mandatory and using as initial guess.\n" \ " 2025-02-18 Added \"-initial\" option.\n" \ " 2025-02-18 Added \"-normals\" as alternative to \"-slopes\".\n" \ "WARRANTY\n" \ argparser_help_info_NO_WARRANTY "\n" \ "\n" \ "RIGHTS\n" \ " " test_integrate_iterative_C_COPYRIGHT ".\n" \ "\n" \ argparser_help_info_STANDARD_RIGHTS #define DEFAULT_MAX_ITER 100000 /* Default max iterations per level. */ #define DEFAULT_CONV_TOL 0.0000005 /* Default convergence threshold per level. */ #define stringify(x) strngf(x) #define strngf(x) #x #define PROG_INFO_DESC \ " The program reads a /slope map/ {G}, consisting of" \ " the {X} and {Y} derivatives of some terrain surface {Z(X,Y)}, averaged" \ " in each cell of a 2D integer grid; and" \ " computes the corresponding /height map/ {Z} from that data. Optionally, the" \ " program reads a one- or two-channel /hints map/ {H} with an independent" \ " estimate of the desired height map {Z}. Optionally, the" \ " program may also compare the computed height map {Z} with" \ " a reference height map {R} provided by the user.\n" \ "\n" \ " Instead of a slope map, the program may read a /normal map/ {N} that" \ " specifies the unit normal vector in each cell of the grid. The" \ " program converts {N} internally to a slope map {G}, and proceeds" \ " as in that case.\n" \ "\n" \ " The height map is computed by solving a linear equation system with" \ " the Gauss-Seidel or Gauss-Jordan iterative method.\n" \ "\n" \ " " pst_integrate_build_system_ARGS_INFO("WH_MULT") "\n" \ "\n" \ " The iterative algorithm is" \ " stopped after a specified number of iterations, of after the height changes" \ " between iterations is less than a specified amount. The initial" \ " value for the iteration may be either all zeros or the reference" \ " height map {R}, depending on the \"-initial\" option. In either" \ " case, a specified amount of noise is added to the initial guess.\n" \ "\n" \ " At some selected iterations, the program" \ " will an error map {E}, defined as the difference between" \ " the height map {Z} computed by the program and a reference height" \ " map {R} specified by the \"-reference\" command line" \ " argument. At selected iterations, the error map {E} is written" \ " to the file \"{ERR_Z_FNI_NAME}\". The error map {E} is shifted to" \ " have zero mean. Any {Z} values that are" \ " next to missing or unreliable slope data are excluded from" \ " the comparison. The reference map {R} may also be used" \ " as the initial guess for the iterative system, depending" \ " on the \"-initial\" option.\n" \ "\n" \ " If the \"-reportStep\" option is given, the program" \ " also writes out the current height map at each" \ " iteration. If \"-reportStep {REPORT_STEP}\" is" \ " given and {REPORT_STEP} is nonzero, these diagnostics are written" \ " every that many iterations. The options" \ " \"-debugG\" and \"-debugSys\" cause additional" \ " debugging data to be written---namely, the reduced slope" \ " maps and the linear systems solved at each scale," \ " respectively." #define PROG_INFO_OUT_FILES \ " {PREFIX}-Z.fni\n" \ " The computed height map {Z}.\n" \ "\n" \ " {PREFIX}-E.fni\n" \ " The error map {E} for the cmputed map {Z}.\n" \ "\n" \ " {PREFIX}-E.txt\n" \ " A one-line summary of the error, with the format\n" \ "\n" \ " {NX} {NY} {iter} {change} {sR} {sZ} {sE} {sRelE}.\n" \ "\n" \ " where {NX,NY} is the image size, {iter} is the number of" \ " iterations performed, {change} is the change in" \ " the last iteration, {sR} and {sZ} are the deviations of" \ " the {R} and {Z} maps, {sE} is the RMS value of {E}, and" \ " {sRelE} is the relative error.\n" \ "\n" \ " The following files may be produced at the beginning, during, and" \ " after the end of the iteration, if the \"-reportStep\" option" \ " is given with non-zero alrgument. The {ITER} part is" \ " the number of iterations done, with 9 digits. " \ "\n" \ " {PREFIX}-dbg-{ITER}-Z.fni\n" \ " The computed height map {Z} after {ITER} iterations.\n" \ "\n" \ " {PREFIX}-dbg-{ITER}-E.fni\n" \ " The height error map {E} for {Z} and {R} after {ITER} iterations.\n" \ "\n" \ " {PREFIX}-dbg-{ITER}-E.txt\n" \ " The height error summary after {ITER} iterations." #define PROG_INFO_OPTS \ " -slopes {G_FNI_NAME}\n" \ " -normals {N_FNI_NAME}\n" \ " Exactly one of these two arguments must be given, to" \ " specify the file name containing the input slope map {G}, or" \ " the input normal map {N}. In either case the file must" \ " in the FNI format (see {float_image_read}).\n" \ "\n" \ " The slope map must be a two- or three-channel image, containing" \ " the derivatives {dZ/dX} and {dZ/dY} averaged in each cell of" \ " the image domain grid, and an optional reliability" \ " weight for that data (see pst_slope_map.h).\n" \ "\n" \ " The normal map must be a three- or four-channel, containing" \ " the coordinates {Nx,Ny,Nz} of the unint vector which is the mean" \ " surface normal in each cell, and an optional reliability" \ " weight for the same (see {pst_normal_map.h})." \ "\n" \ " -hints {H_FNI_NAME} {WH_MULT}\n" \ " This optional argument specifies the file name" \ " containing the independent estimated height map {H}, which must" \ " be a one- or two-channel image in FNI format" \ " (see pst_height_map.h). The argument {WH_MULT} must be a number" \ " between 0 and 1 that will be multiplied by the reliability weight" \ " of every height in {H}.\n" \ " -reference {REF_Z_FNI_NAME}\n" \ " This mandatory argument specifies the file that contains the" \ " expecied height map. Ideally this map must have either the same size" \ " as the expected height map, that is, one col and one row more than the" \ " slope map. However, for convenience, the program also accepts a heigh" \ " map with the same size as the slope (or normal) map. In this case, the" \ " reference map will be expanded by one col and one row" \ " with {float_image_expand_by_one}.\n" \ "\n" \ " -initial {INIT_OPT} {INIT_NOISE}\n" \ " This mandatory argument specifies the initial guess of the height" \ " map for the iterative solver. The {INIT_OPT} may be either \"zero\", to" \ " mean all initial heights are zero, or \"reference\", to mean that" \ " the initial guess is the same as the reference height map. In either" \ " case, if {INIT_NOISE} is nonzero, the iitial map will be perturbed before" \ " the iteration starts. The perturbation consists in adding to every" \ " finite height value in the map {Z} a pseudorandom amount in the" \ " interval {[-mag _ +mag]}, where {mag} is {INIT_NOISE} times" \ " {hypot(rad, vmax-vmin)}, {rad} is half of the smallest" \ " dimension (cols or rows) of the map {Z}, and {[vmin _ vmax]} the" \ " original range of the heights in {Z} (ignoring non-finite values).\n" \ "\n" \ " -outPrefix {PREFIX}\n" \ " Specifies the common prefix for all output" \ " file names. Mandatory." \ "\n" \ " -maxIter {MAX_ITER}\n" \ " -convTol {CONV_TOL}\n" \ " These optional parameters specify the stopping criterion for the" \ " iteration. The iteration will stop when the maximuum change in any" \ " height field is less than {CONV_TOL}, or after {MAX_ITER} iterations, whichever" \ " happens first. The defaults are {MAX_ITER = " stringify(DEFAULT_MAX_ITER) "}," \ " {CONV_TOL = " stringify(DEFAULT_CONV_TOL) "}.\n" \ "\n" \ " -sortSys {TOPOSORT}\n" \ " This optional boolean argument specifies the order in which" \ " equations are solved by the program. The {TOPOSORT} may be 'T' or 1 to" \ " solve in order of increasing equation weight, or 'F' of 0 to solve in" \ " scan line orrder. If omitted, the program assumes \"-sortSys F\".\n" \ "\n" \ " -verbose\n" \ " Causes the program to write various diagnostics to stderr, in" \ " particular one line per iteration of the Gauss-Seidel solver.\n" \ "\n" \ " -reportStep {REPORT_STEP}\n" \ " If this optional argument is given and {REPORT_STEP} is" \ " not zero, the heights {Z}, the error maps {E}, and the" \ " error summaries are written out after every {REPORT_STEP} iterations, as" \ " well as at the beginning (0 iterations) and after the end of the" \ " iterative solving. If {REPORT_STEP} is larger than {MAX_ITER}, only" \ " the initial and final states are written. If" \ " {REPORT_STEP} is zero (the default)," \ " no \"...-{ITER}-...\" files are written." #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /* COMMAND-LINE OPTIONS */ typedef struct options_t { char *slopes; /* File name of input slope map file, or NULL. */ char *normals; /* File name of input normal map file, or NULL. */ char *reference; /* File name of input reference height map (NULL if none). */ char *initial_opt; /* Initial option: "zero" or "reference". */ double initial_noise; /* Relative magnitude of perturbations to the initial map. */ char *hints_file; /* File name of external height estimates (NULL if none). */ double hints_weight; /* Weight factor for {H} weights. */ uint32_t maxIter; /* Max iterations. */ double convTol; /* Convergence threshold. */ bool_t sortSys; /* TRUE solves the equations in order of increasing weight. */ char *outPrefix; /* Output file name prefix. */ /* Debugging and analysis */ bool_t verbose; /* TRUE to print various diagnostics. */ uint32_t reportStep; /* Frequency for debugging output during iteration, or 0 if none. */ } options_t; float_image_t *tini_read_fni_file(char *fileName); /* Reads a FNI image from file "{fileName}" (which should include the extension ".fni"). If {fileName} is "-", reads from standard input. */ void tini_write_fni_file(float_image_t *I, char *fileName, int32_t indent); /* Writes the float image {I} in a human-readable format, to a file called "{fileName}" (which should include the extension ".fni"). If {fileName} is "-", writes the image to standard output. Diagnostic images are indented by {indent} spaces. */ options_t *tini_parse_options(int32_t argc, char **argv); /* Parses the command line arguments and packs them as an {options_t}. */ void tini_compute_and_write_height_map ( options_t *o, float_image_t *G, float_image_t *H, float_image_t *Z, float_image_t *R ); /* Computes the height map {Z} from the gradient map {G} and optional hints map {H}. The map {G} must have two or three channels. Channels 0 and 1 are intrepreted as the slopes {dZ/dX} and {dZ/dY}). Channel 2, if it exists, is interpreted as the reliability weight of the corresponding slope data; otherwise the weight is assumed to be 1. The image {H} must have one or two channels, and one row and one col more than {G}. On input, it gives an independent estimate of the heights, and also the initial guess for the iterative system solver. The image {Z} must have two channels, and one row and one col more than {G}. Its contents is ignored on input. On output, sample {Z[0,X,Y]} will be the computed height, and {Z[1,X,Y]} the corresponding reliability weight. If {o.reportStep} is not zero, writes the height map for the first {o.reportStep} iterations, and then every {o.reportStep} iterations thereafter. If {R} is not null, it must be a single-channel image with the same size as {G} or {Z}. Whenever {Z} is reported, the procedure compares {Z} with {R} and writes the error map {E=Z-R} and the error sumary file. */ int32_t main(int32_t argc,char** argv); /* IMPLEMENTATIONS */ int32_t main(int32_t argc, char** argv) { options_t *o = tini_parse_options(argc, argv); float_image_t *G; /* Input gradient map. */ if (o->slopes != NULL) { fprintf(stderr, "reading the slope map {G} ...\n"); G = tini_read_fni_file(o->slopes); } else if (o->normals != NULL) { fprintf(stderr, "reading the normal map {N} ...\n"); float_image_t *N = tini_read_fni_file(o->normals); fprintf(stderr, "converting the normal map to a slope map {G} ...\n"); double maxSlope = 1000.0; /* Should be enough... */ G = pst_normal_map_to_slope_map(N, maxSlope); float_image_free(N); } else { assert(FALSE); } float_image_t *H; /* Input external height estimates. */ if (o->hints_file != NULL) { fprintf(stderr, "reading the hints height map {H} ...\n"); H = tini_read_fni_file(o->hints_file); } int32_t NX_G, NY_G; float_image_get_size(G, NULL, &NX_G, &NY_G); int32_t NC_Z = 2; int32_t NX_Z = NX_G + 1; int32_t NY_Z = NY_G + 1; float_image_t *Z = float_image_new(NC_Z, NX_Z, NY_Z); float_image_t *R; /* Reference height map. */ fprintf(stderr, "reading the reference height map {R} ...\n"); R = tini_read_fni_file(o->reference); fprintf(stderr, "computing the height map {Z} ...\n"); tini_compute_and_write_height_map(o, G, H, Z, R); float_image_free(G); G = NULL; if (H != NULL) { float_image_free(H); H = NULL; } if (Z != NULL) { float_image_free(Z); Z = NULL; } if (R != NULL) { float_image_free(R); R = NULL; } fprintf(stderr, "Done!\n"); return 0; } void tini_compute_and_write_height_map ( options_t *o, float_image_t *G, float_image_t *H, float_image_t *Z, float_image_t *R ) { int32_t NC_G, NX_G, NY_G; float_image_get_size(G, &NC_G, &NX_G, &NY_G); demand((NC_G == 2) || (NC_G == 3), "gradient map {G} must have 2 or 3 channels"); int32_t NX_Z = NX_G + 1; int32_t NY_Z = NY_G + 1; float_image_check_size(Z, 2, NX_Z, NY_Z, "bad height map {Z}, wring size or depth"); if (H != NULL) { int32_t NC_H, NX_H, NY_H; float_image_get_size(H, &NC_H, &NX_H, &NY_H); demand((NC_H == 1) || (NC_H == 2), "hints height map {H} must have 1 or 2 channels"); demand((NX_H == NX_Z) && (NY_H == NY_Z), "bad hints height map {H}, wrong size size"); } float_image_t *IR = NULL; /* Version of {R} with correct size. */ if (R != NULL) { int32_t NC_R, NX_R, NY_R; float_image_get_size(R, &NC_R, &NX_R, &NY_R); demand((NC_R == 1) || (NC_R == 2), "reference height map must have 1 or 2 channels"); if ((NX_R == NX_Z) && (NY_R == NY_Z)) { IR = R; } else if ((NX_R == NX_G) && (NY_R == NY_G)) { /* Ref height map has height per pixel instead of ver corner. */ /* Resize to same size as {Z}: */ fprintf(stderr, "Expanding the reference height map {R} to the correct size ...\n"); IR = float_image_expand_by_one(R, 1); } else { demand(FALSE, "bad reference height map {R}, wrong size"); } } /* Define the initial guess: */ if (strcmp(o->initial_opt, "zero") == 0) { fprintf(stderr, "zeroing the initial solution ...\n"); float_image_fill_channel(Z, 0, 0.0); } else if (strcmp(o->initial_opt, "reference") == 0) { fprintf(stderr, "using the ref height map as the initial solution ...\n"); float_image_assign_channel_rectangle(Z, 0, 0,NX_Z-1, 0,NY_Z-1, IR, 0, 0,0); } else { demand(FALSE, "invalid \"-initial\" option"); } float_image_fill_channel(Z, 1, 1.0); if (o->initial_noise > 0.0) { fprintf(stderr, "Perturbing the initial guess ...\n"); pst_map_perturb_values(Z, 1, o->initial_noise); } { fprintf(stderr, "writing ou the initial guess ...\n"); char *fname_ini = jsprintf("%s-ini-Z.fni", o->outPrefix); tini_write_fni_file(Z, fname_ini, 0); free(fname_ini); } bool_t zeroMean = TRUE; char *debugPrefix = txtcat(o->outPrefix, "-dbg"); /* Prefix for debug file names. */ auto void reportSys(int32_t level, pst_imgsys_t *S); /* This procedure is called by {pst_integrate_iterative} once to report the linear system created from the slope and weight maps. */ auto void reportHeights(int32_t level, int32_t iter, double change, bool_t final, float_image_t *Z); /* This procedure is called by {pst_integrate_iterative} to write the current height map {Z} and the height error map {E}, and the error summary. */ /* Call iterative integrator: */ int32_t level = -1; pst_integrate_iterative ( G, H, o->hints_weight, Z, o->maxIter, o->convTol, o->sortSys, o->verbose, level, &reportSys, o->reportStep, &reportHeights ); /* Free working storage: */ free(debugPrefix); if ((IR != NULL) && (IR != R)) { float_image_free(IR); } return; void reportHeights(int32_t level, int32_t iter, double change, bool_t final, float_image_t *Z) { assert(level == -1); bool_t writeImages = TRUE; bool_t writeError = (IR != NULL); /* Without {H} the solution has indeterminate constant shift, so: */ pst_height_map_level_analyze_and_write ( debugPrefix, level, (final ? -1 : iter), change, Z, IR, writeImages, zeroMean, writeError ); } void reportSys(int32_t level, pst_imgsys_t *S) { char *fileName_sys = float_image_mscale_file_name(debugPrefix, level, -1, "sys", "txt"); FILE* wr_sys = fopen(fileName_sys, "wt"); pst_imgsys_write(wr_sys, S); fclose(wr_sys); free(fileName_sys); float_image_t *SW = pst_imgsys_make_weight_image(S); char *fileName_SW = float_image_mscale_file_name(debugPrefix, level, -1, "SW", "fni"); float_image_write_named(fileName_SW, SW); free(fileName_SW); float_image_free(SW); } } float_image_t *tini_read_fni_file(char *fileName) { demand(fileName != NULL, "file name not given"); fprintf(stderr, "Reading %s ...", fileName); FILE *rd = open_read(fileName, FALSE); float_image_t *I = float_image_read(rd); if (rd != stdin) { fclose(rd); } fprintf(stderr, "\n"); return I; } void tini_write_fni_file(float_image_t *I, char *fileName, int32_t indent) { demand(fileName != NULL, "file name not given"); fprintf(stderr, "%*sWriting %s ...", indent, "", fileName); FILE* wr = open_write(fileName, FALSE); float_image_write(wr, I); if (wr == stdout) { fflush(wr); } else { fclose(wr); } fprintf(stderr, "\n"); } options_t *tini_parse_options(int32_t argc, char **argv) { argparser_t *pp = argparser_new(stderr, argc, argv); argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP); argparser_set_info(pp, PROG_INFO); argparser_process_help_info_options(pp); options_t *o = (options_t *)notnull(malloc(sizeof(options_t)), "no mem"); if (argparser_keyword_present(pp, "-slopes")) { o->slopes = argparser_get_next_non_keyword(pp); } else { o->slopes = NULL; } if (argparser_keyword_present(pp, "-normals")) { o->normals = argparser_get_next_non_keyword(pp); } else { o->normals = NULL; } if ((o->slopes == NULL) == (o->normals == NULL)) { argparser_error(pp, "must specify exaclty one of \"-slopes\" or \"-normals\""); } if (argparser_keyword_present(pp, "-hints")) { o->hints_file = argparser_get_next_non_keyword(pp); o->hints_weight = argparser_get_next_double(pp, 0.0, 1.0); } else { o->hints_file = NULL; o->hints_weight = 0.0; } if (argparser_keyword_present(pp, "-maxIter")) { o->maxIter = (uint32_t)argparser_get_next_int(pp, 0, INT32_MAX); } else { o->maxIter = DEFAULT_MAX_ITER; } if (argparser_keyword_present(pp, "-convTol")) { o->convTol = argparser_get_next_double(pp, 0, DBL_MAX); } else { o->convTol = DEFAULT_CONV_TOL; } if (argparser_keyword_present(pp, "-sortSys")) { o->sortSys = argparser_get_next_bool(pp); } else { o->sortSys = FALSE; } argparser_get_keyword(pp, "-reference"); o->reference = argparser_get_next(pp); argparser_get_keyword(pp, "-initial"); o->initial_opt = argparser_get_next_non_keyword(pp); o->initial_noise = argparser_get_next_double(pp, 0, DBL_MAX); o->verbose = argparser_keyword_present(pp, "-verbose"); if (argparser_keyword_present(pp, "-reportStep")) { o->reportStep = (uint32_t)argparser_get_next_int(pp, 0, INT32_MAX); } else { o->reportStep = 0; } argparser_get_keyword(pp, "-outPrefix"); o->outPrefix = argparser_get_next(pp); argparser_skip_parsed(pp); argparser_finish(pp); return o; }