#define PROG_NAME "make_test_slope_maps" #define PROG_DESC "generate gradient maps for tests of slope_to_height" #define PROG_VERS "1.0" /* Copyright � 2005 by the State University of Campinas (UNICAMP). */ /* See the copyright, authorship, and warranty notice at end of file. */ /* Last edited on 2011-02-28 00:38:16 by stolfi */ #define PROG_HELP \ " " PROG_NAME " \\\n" \ " -function {ZFUNC} \\\n" \ " -size {NX} {NY} \\\n" \ " [ -smoothZ {LZ} {NZ} ] \\\n" \ " [ -smoothG {LG} {NG} ] \\\n" \ " [ -numGrad ] \\\n" \ " -noiseG {SIGMA_G} \\\n" \ " -noiseW {SIGMA_W} \\\n" \ " -maxGDiff {MAXGDIFF} \\\n" \ " -outPrefix {OUT_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" \ " Generates a height map {Z(X,Y)}, its gradient (slope) map {G(X,Y)}, and" \ " a corresponding weight map, in the format required" \ " by {slope_to_height}. The height map is written to the" \ " single-channel float image file \"{OUT_PREFIX}-Z.fni\"." \ " The gradient is written as" \ " channels 0 (X derivative) and channel 1 (Y derivative) of the two-channel float image" \ " file \"{OUT_PREFIX}-G.fni\". The weight" \ " map is written to \"{OUT_PREFIX}-W.fni\"." \ " The normal map is written as the three-channel float image" \ " file \"{OUT_PREFIX}-N.fni\".\n" \ "\n" \ " The height map is defined by the {ZFUNC} parameter, an integer" \ " that selects one of the program's built-in functions" \ " ({function_00}, {function_01}, ...). The parameters" \ " {NX} and {NY} are the width and height of the slope" \ " map, in pixels.\n" \ "\n" \ " The height function is nominally evaluated" \ " at the corners of a rectangular grid with" \ " {NX} by {NY} square cells. Note that the height map" \ " has {NX+1} columns and {NY+1} rows. If \"-smoothZ\" is not" \ " specified, the height is sampled at that point only, otherwise it" \ " is an average of {NZ} by {NZ} samples in a" \ " window {LZ} by {LZ} pixels centered at that point.\n" \ "\n" \ " The slope maps are nominally evaluated at the center of each pixel" \ " of the grid, so it has {NX} columns and {NY} rows. If \"-numGrad\" is" \ " specified, the average gradient in each cell is computed by finite" \ " differences from the height values at the four corners of the" \ " cell. If \"-numGrad\" is not specified, the gradient is computed" \ " by sampling the analytic gradient of the function. In any case," \ " if \"-smoothG\" is not specified, the analytic gradient is sampled" \ " at the center of the pixel, otherwise it is an average of {NG} by" \ " {NG} samples in a window {LG} by {LG} pixels centered at that" \ " point around that point. Note that the analytic gradient may be" \ " widely inconsistent with the height field, especially if the function" \ " is discontinuous or highly random inside the pixel.\n" \ "\n" \ " The allowed values for {LZ} and {LG} are 0, 1, and 2. If {NZ} is 1, a" \ " single sample is taken at the nominal position; this is the default" \ " behavior. One of the samples will fall at the nominal position if" \ " and only if {NZ} is odd. If {LZ} is zero, then {NZ} must be 1. \n" \ "\n" \ " In any case, the reliability weight of each pixel is set ot 1.0 if the" \ " \"-numGrad\" option is given. If \"-numGrad\" is not given, the weight" \ " is 1 if the average analytic gradient coincides with the numeric gradient, and" \ " decreases to zero as the Euclidean norm of the difference between the two" \ " gradients increases to or beyond {MAXGDIFF}.\n" \ "\n" \ " If {SIGMA_G} is positive, the program also randomly perturbs each" \ " component of the gradient value, to simulate noisy data. The" \ " perturbation consists in adding a random number {E} with" \ " Gaussian distribution, zero mean, and standard" \ " deviation {SIGMA_G}. The perturbation affects only" \ " the gradient map, not the height map.\n" \ "\n" \ " Similarly, if {SIGMA_W} is positive, each weight gets mutiplied" \ " by {exp(E)} where {E} is a random number with Gaussian distribution,"\ " zero mean, and standard deviation {SIGMA_W}.\n" \ "\n" \ " The normal map is computed from the (perturbed) gradient map. The Z"\ " component is always positive.\n" \ "\n" \ "SEE ALSO\n" \ " fni_to_pnm(1), pnm_to_fni(1), slope_to_height(1)\n" \ "\n" \ "AUTHOR\n" \ " Created 2005-08-15 by Jorge Stolfi, UNICAMP.\n" \ "MODIFICATION HISTORY\n" \ " 2010-05-04 by J. Stolfi, IC-UNICAMP:\n" \ " * Weight is always 1 with \"-numGrad\"." #define _GNU_SOURCE #include #include #include #include #include #include #include #include #include #include #include #include #include #include typedef struct sampling_opts_t { int N; /* Number of sampling points per axis. */ int L; /* Width of window in pixels. */ } sampling_opts_t; /* A sampling specification. The window width must be an integer to ensure the partition-of-unit property. */ typedef struct options_t { int ZFunc; /* Integer code of height map. */ int NX, NY; /* Image size. */ sampling_opts_t smoothZ; /* Sampling options for the heights. */ sampling_opts_t smoothG; /* Sampling options for the analytic gradient. */ bool_t numGrad; /* Output numeric gradient instead of analytic one. */ double maxGDiff; /* Max discrepancy between analyic and numeric gradient. */ double sigmaG; /* Deviation of gaussian gradient noise. */ double sigmaW; /* Deviation of log-gaussian weight noise. */ char* outPrefix; /* Output name prefix. */ } options_t; typedef struct sampling_t { int N; /* Number of sampling points per axis. */ double *d; /* Sampling positions relative to center point (indexed {0..NS-1}). */ double *w; /* Sampling weights (indexed {0..NS-1}). */ } sampling_t; /* A sampling parameters and tables. */ /* INTERNAL PROTOTYPES */ typedef void zfunc_t (r2_t p, double *z, r2_t *dz); /* A function that returns a height field {Z} and its gradient at a given point {p} of the plane. A function of this type should should return in {*z} the height {Z(p)}, which should be a number in the range {[-1 _ +1]}. If {dz} is not NULL, it should also return in {*dz} the gradient of {Z} with respect to {p}. */ void make_test_images ( zfunc_t *func, sampling_opts_t *smoothZ, sampling_opts_t *smoothG, bool_t numGrad, double maxGDiff, double sigmaG, double sigmaW, float_image_t *IZ, float_image_t *IG, float_image_t *IW, float_image_t *IN ); /* Create a height map {Z(X,Y)}, the associated gradient map {G(X,Y) = dZ/dX(X,Y),dZ/dY(X,Y)}, a reliability weight map {W(X,Y)}, and a normal map {N(X,Y)}, from the given function {func}. All images refer to the same grid {G} with unit square cells, spanning a rectangle {[0 _ NX]�[0 _ NY]}. The height {Z(X,Y)} is nominally computed at the grid vertex {(X,Y)} and stored in pixel {[X,Y]} of the grid of image {IZ}, which must have {NX+1} columns and {NY+1} rows. The sampling is controlled by {smpZ}. The analytic gradient is nominally computed at pixel centers by sampling the analytic derivatives as computed by {func}. The sampling is controlled by {smpG}. The numeric gradient is computed by averaging the height differences at the four corners of the pixel. The {numGrad} parameter selects which gradient will be returned as the gradient map {G(X,Y)}. In any case the gradient map is stored into the two-channel float image {IG}, which must have {NX} columns and {NY} rows If {numGrad} is true, the weights {IW[X,Y]} are 1 everywhere. If {numGrad} is false, each weight {IW[X,Y]} is computed by comparing the analytic gradient {IG[X,Y]} for the same pixel with the gradient computed by numerical differentiation of the {IZ} values at the corners of that pixel. The image {IW} must have {NX} columns and {NY} rows. The normal map, computed from the slope map, is stored in the three-channel float image {IN}, which must have {NX} columns and {NY} rows. The parameters {sigmaG} and {sigmaW} control the amount of noise to be introduced in the gradient and weight maps, respectively. When calling the {func} procedure, the image domain coordinates {P=(X,Y)} are scaled down and shifted so that the {func} argument {p=(x,y)} spans the signed unit square {[-1_+1]�[-1_+1]} of the plane while {P} spans a square that fits snugly in the grid {G}, centered. The function values returned by {func} are scaled up so that the computed gradient remains consistent with the image coordinate system. */ sampling_t compute_sampling_tables(int L, int N, double pixels_per_unit); /* Returns a sampling table {smp} for the given sampling parameters. Allocates {smp->d} and {smp->w} and fills them with proper values as specified by {L} and {N}. If {N} is 1 then {L} is ignored. Otherwise, if {L} is 1 the weights are uniform. Otherwise, Hann cosine weights are used. The displacement are in function domain units. The {pixels_per_unit} parameter is the number of pixels in one function domain unit. */ double compute_height_value ( zfunc_t *func, r2_t p, sampling_t *smp ); /* Computes the height field defined by {func} at the grid corner point {p}. Evaluates {f} at a grid of {NW}�{NW} subsampling points near {p} with displacements {smp->d[0..NW-1]} and averages them with the weights {smp->w[0..NW-1]}. The count {NW} must be odd. */ r2_t compute_analytic_pixel_gradient ( zfunc_t *func, r2_t p, sampling_t *smp ); /* Computes the mean gradient of the height field defined by {func} at the point {p}. Evaluates the gradient of {f} at a grid of {NW}�{NW} subsampling points near {p} with displacements {r*smp->d[0..NW-1]} and averages them with the weights {smp->w[0..NW-1]}. The count {NW} must be odd. */ r2_t compute_numeric_pixel_gradient ( double f00, double f10, double f01, double f11 ); /* Estimates the gradient of the height in a pixel by numerical differences from the four height values at the corners of the pixel: {f00} at bottom left, {f10} at bottom right, {f01} at top left, and {f11} at top right. */ double compute_pixel_weight ( zfunc_t *func, r2_t *dza, r2_t *dzn, double maxGDiff ); /* Computes the reliability weight of a pixel given the average gradient {*dza} in the pixel and the numerical gradient {*dzn}. The weight depends on the Euclidean norm {e} of the difference between the two gradients, and is zero if {e >= maxGDiff}. */ void perturb_gradient(r2_t *dz, double sigma); /* Adds to each component of {*dz} a random value with Gaussian distribution, mean 0 and deviation {sigma}. */ void perturb_weight(double *w, double sigma); /* Multiplies {*w} by {max(0,1-S)} where {S} is a random value with log-Gaussian distribution with mean of log 0 and deviation of log {sigma}. */ void write_test_image ( char *pref, char *tag, float_image_t *I ); /* Writes the image {I} to file named "{pref}-{tag}.fni", in FNI format (see float_image.h}). */ void normalize(double v[]); /* Normalizes the three-vector {v[0..2]} to unit Euclidean length. */ void write_fni_image(char *fileName, float_image_t *I); /* Writes the image {I} to file "{fileName}" as a FNI image, with {float_image_write}. */ void function_00(r2_t p, double *z, r2_t *dz); void function_01(r2_t p, double *z, r2_t *dz); void function_02(r2_t p, double *z, r2_t *dz); void function_03(r2_t p, double *z, r2_t *dz); void function_04(r2_t p, double *z, r2_t *dz); void function_05(r2_t p, double *z, r2_t *dz); void function_06(r2_t p, double *z, r2_t *dz); void function_07(r2_t p, double *z, r2_t *dz); void function_08(r2_t p, double *z, r2_t *dz); void function_09(r2_t p, double *z, r2_t *dz); void function_10(r2_t p, double *z, r2_t *dz); void function_11(r2_t p, double *z, r2_t *dz); void function_12(r2_t p, double *z, r2_t *dz); void function_13(r2_t p, double *z, r2_t *dz); void function_14(r2_t p, double *z, r2_t *dz); void function_15(r2_t p, double *z, r2_t *dz); void function_16(r2_t p, double *z, r2_t *dz); void function_17(r2_t p, double *z, r2_t *dz); void function_18(r2_t p, double *z, r2_t *dz); void function_19(r2_t p, double *z, r2_t *dz); void function_20(r2_t p, double *z, r2_t *dz); void function_21(r2_t p, double *z, r2_t *dz); void function_22(r2_t p, double *z, r2_t *dz); void function_23(r2_t p, double *z, r2_t *dz); /* The built-in height map functions. They assume that the signed unit square {[-1_+1]�[-1_+1]} fits snugly inside the image's domain. !!! Each function should define a weight too. !!! !!! Either that, or let {z=NAN} to get weight 0. !!! !!! Each function should be told the resolution of the map? !!! */ #define MIN_ZFUNC 0 #define MAX_ZFUNC 23 /* Min and max values of the {ZFUNC} command line parameter. */ void function_wave(r2_t *p, r2_t *f, double phase, double *z, r2_t *dz); /* Computes the generic wave function {z = 0.5*sin(2*PI*dot(p,f) + phase)}. Also computes its gradient {dz}. */ void function_babel(r2_t p, double RI, double RO, double N, double EF, double *z, r2_t *dz); /* Computes the generic Tower of Babel function {z} with inner radius {RI}, outer radius {RO}, a ramp with {N} full turns. Also computes its gradient {dz}. If {EF = 0}, the ramp is bounded by a vertical cliff on both sides. If {EF > 0} a filling is added adjacent to the cliff, to smooth it out. The filling width is {EF} times the bare ramp width, so the max value of {EF} is 1. */ options_t *parse_options(int argc, char **argv); /* Parses the command line arguments and packs them as an {options_t}. */ sampling_opts_t parse_sampling_options(argparser_t *pp, char *key); /* Parses an optional smoothing option consisting of the keyword {key} followed by the window width {L} and the number of samples per axis {N}. */ int main(int argc, char** argv); /* IMPLEMENTATIONS */ int main(int argc, char** argv) { options_t *o = parse_options(argc, argv); float_image_t *IZ = float_image_new(1, o->NX+1, o->NY+1); float_image_t *IG = float_image_new(2, o->NX, o->NY); float_image_t *IW = float_image_new(1, o->NX, o->NY); float_image_t *IN = float_image_new(3, o->NX, o->NY); zfunc_t *func = NULL; switch(o->ZFunc) { case 0: func = &function_00; break; case 1: func = &function_01; break; case 2: func = &function_02; break; case 3: func = &function_03; break; case 4: func = &function_04; break; case 5: func = &function_05; break; case 6: func = &function_06; break; case 7: func = &function_07; break; case 8: func = &function_08; break; case 9: func = &function_09; break; case 10: func = &function_10; break; case 11: func = &function_11; break; case 12: func = &function_12; break; case 13: func = &function_13; break; case 14: func = &function_14; break; case 15: func = &function_15; break; case 16: func = &function_16; break; case 17: func = &function_17; break; case 18: func = &function_18; break; case 19: func = &function_19; break; case 20: func = &function_20; break; case 21: func = &function_21; break; case 22: func = &function_22; break; case 23: func = &function_23; break; default: fprintf(stderr, "bad ZFunc = %d\n", o->ZFunc); exit(1); } /* Compute height and gradient maps: */ make_test_images(func, &(o->smoothZ), &(o->smoothG), o->numGrad, o->maxGDiff, o->sigmaG, o->sigmaW, IZ, IG, IW, IN); float_image_t* IEW = pst_weight_map_heights_from_slopes(IZ->sz[1], IZ->sz[2],IW); /* Write images: */ write_test_image(o->outPrefix, "Z", IZ); write_test_image(o->outPrefix, "G", IG); write_test_image(o->outPrefix, "W", IW); write_test_image(o->outPrefix, "N", IN); write_test_image(o->outPrefix, "EW", IEW); return 0; } void make_test_images ( zfunc_t *func, sampling_opts_t *smoothZ, sampling_opts_t *smoothG, bool_t numGrad, double maxGDiff, double sigmaG, double sigmaW, float_image_t *IZ, float_image_t *IG, float_image_t *IW, float_image_t *IN ) { assert(IZ->sz[0] == 1); assert(IG->sz[0] == 2); assert(IW->sz[0] == 1); int NX = IG->sz[1]; assert(IW->sz[1] == NX); assert(IZ->sz[1] == NX+1); int NY = IG->sz[2]; assert(IW->sz[2] == NY); assert(IZ->sz[2] == NY+1); /* Center of image domain (in pixels, from bottom left corner): */ double cX = 0.5*NX; double cY = 0.5*NY; int t = (NX < NY ? NX : NY); /* Smallest dimension of image. */ double pixels_per_unit = t/2.0; /* Number of grid pixels per {func} domain unit. */ /* Compute the sample displacements and weights: */ sampling_t smpZ = compute_sampling_tables(smoothZ->L, smoothZ->N, pixels_per_unit); sampling_t smpG = compute_sampling_tables(smoothG->L, smoothG->N, pixels_per_unit); /* Compute height at each grid corner: */ int X, Y, c; for (Y = 0; Y <= NY; Y++) { for (X = 0; X <= NX; X++) { /* Compute function domain coordinates {(x,y)} of grid corner: */ double x = (X - cX)/pixels_per_unit; double y = (Y - cY)/pixels_per_unit; r2_t p = (r2_t){{ x, y }}; /* Compute height at grid corner: */ double za = compute_height_value(func, p, &smpZ); /* Scale height from {func} units to pixels: */ float_image_set_sample(IZ, 0, X, Y, za * pixels_per_unit); } } /* Ensure that the mean Z value is zero: */ double mean_z; float_image_compute_sample_avg_dev(IZ, 0, &mean_z, NULL); for (Y = 0; Y <= NY; Y++) { for (X = 0; X <= NX; X++) { double z = float_image_get_sample(IZ, 0, X, Y); float_image_set_sample(IZ, 0, X, Y, z - mean_z); } } /* Compute slope, weight, and normal map at each grid pixel: */ for (Y = 0; Y < NY; Y++) { for (X = 0; X < NX; X++) { /* Compute function domain coordinates {(x,y)} of pixel center: */ double x = (X + 0.5 - cX)/pixels_per_unit; double y = (Y + 0.5 - cY)/pixels_per_unit; r2_t p = (r2_t){{ x, y }}; /* Compute analytic gradient at pixel center: */ r2_t dza = compute_analytic_pixel_gradient(func, p, &smpG); /* Compute the numeric gradient by simple differences of heights: */ double f00 = float_image_get_sample(IZ, 0, X, Y); double f10 = float_image_get_sample(IZ, 0, X+1, Y); double f01 = float_image_get_sample(IZ, 0, X, Y+1); double f11 = float_image_get_sample(IZ, 0, X+1, Y+1); r2_t dzn = compute_numeric_pixel_gradient(f00, f10, f01, f11); /* Compute weight by comparing the analytic and numeric gradients: */ double wt = (numGrad ? 1.0 : compute_pixel_weight(func, &dza, &dzn, maxGDiff)); /* Choose final gradient: */ r2_t dz = (numGrad ? dzn : dza); /* Add perturbations to the slopes and weights: */ if (sigmaG > 0.0) { perturb_gradient(&dz, sigmaG); } if (sigmaW > 0.0) { perturb_weight(&wt, sigmaW); } /* Compute average normal from the average gradient: */ r3_t nrm = pst_normal_map_normal_from_slope(&dz); /* Save the results in the images: */ float_image_set_sample(IG, 0, X, Y, dz.c[0]); float_image_set_sample(IG, 1, X, Y, dz.c[1]); float_image_set_sample(IW, 0, X, Y, wt); for (c = 0; c < 3; c++) { float_image_set_sample(IN, c, X, Y, nrm.c[c]); } } } } sampling_t compute_sampling_tables(int L, int N, double pixels_per_unit) { assert(L <= 2); /* For now. */ /* Compute sample positions and weights for height sampling: */ sampling_t smp; smp.N = N; smp.d = notnull(malloc(N*sizeof(double)), "no mem"); smp.w = notnull(malloc(N*sizeof(double)), "no mem"); int k; for (k = 0; k < N; k++) { double rk = (k + 0.5)/N - 0.5; /* Relative sample position in {[-1/2 _ +1/2]}. */ smp.d[k] = rk*L/pixels_per_unit; /* Relative sample position in function domain units. */ smp.w[k] = (L <= 1 ? 1.0 : 0.5*(1 + cos(2*M_PI*rk))); fprintf(stderr, " sample %2d position = %12.8f weight = %12.8f\n", k, smp.d[k], smp.w[k]); } return smp; } double compute_height_value ( zfunc_t *func, r2_t p, sampling_t *smp ) { int NS = smp->N; double sum_w = 0; double sum_wz = 0; int xs, ys; for (ys = 0; ys < NS; ys++) for (xs = 0; xs < NS; xs++) { /* Generate a sample point {(xk,yk)} inside cell {[X,Y]}: */ double xk = p.c[0] + smp->d[xs]; double yk = p.c[1] + smp->d[ys]; r2_t pk = (r2_t){{ xk, yk }}; /* Get height value, ignore gradient: */ double zk; func(pk, &zk, NULL); assert(isfinite(zk)); /* Neither infinity nor NAN. */ /* Get sample weight: */ double wk = smp->w[xs]*smp->w[ys]; /* Accumulate height value: */ sum_wz += wk*zk; sum_w += wk; } /* Compute average: */ assert(sum_w > 0); double za = sum_wz/sum_w; return za; } r2_t compute_analytic_pixel_gradient ( zfunc_t *func, r2_t p, sampling_t *smp ) { int NS = smp->N; double sum_w = 0; r2_t sum_wdz = (r2_t){{ 0, 0 }}; int xs, ys; for (ys = 0; ys < NS; ys++) for (xs = 0; xs < NS; xs++) { /* Generate a sample point {(xk,yk)} inside cell {[X,Y]}: */ double xk = p.c[0] + smp->d[xs]; double yk = p.c[1] + smp->d[ys]; r2_t pk = (r2_t){{ xk, yk }}; /* Get height value, derivatives: */ double zk; r2_t dzk; func(pk, &zk, &dzk); assert(isfinite(dzk.c[0])); /* Neither infinity nor NAN. */ assert(isfinite(dzk.c[1])); /* Neither infinity nor NAN. */ /* Get sample weight: */ double wk = smp->w[xs]*smp->w[ys]; /* Accumulate gradient: */ int c; for (c=0; c < 2; c++) { sum_wdz.c[c] += wk*dzk.c[c]; } sum_w += wk; } /* Compute mean gradient from sum of gradients: */ assert(sum_w > 0); r2_t dza; r2_scale(1/sum_w, &sum_wdz, &dza); return dza; } r2_t compute_numeric_pixel_gradient ( double f00, double f10, double f01, double f11 ) { r2_t dzn; dzn.c[0] = 0.5*(f10 - f00 + f11 - f01); dzn.c[1] = 0.5*(f01 - f00 + f11 - f10); return dzn; } double compute_pixel_weight ( zfunc_t *func, r2_t *dza, r2_t *dzn, double maxGDiff ) { bool_t debug = FALSE; double diff = r2_dist(dza, dzn); double rerr = diff/maxGDiff; double wa = (rerr > 1.0 ? 0.0 : (1.0 - rerr)*(1.0 - rerr)); if (debug && (wa < 0.95)) { fprintf(stderr, "compute_pixel_weight:"); fprintf(stderr, "dza = %e %e dzn = %e %e\n", dza->c[0], dza->c[1], dzn->c[0], dzn->c[1]); fprintf(stderr, "diff = %e rerr = %e wa = %e\n", diff, rerr, wa); } return wa; } void write_test_image(char *pref, char *tag, float_image_t *I) char *fileName = jsprintf("%s-%s.fni", pref, tag); write_fni_image(fileName, I); free(fileName); } void perturb_gradient(r2_t *dz, double sigma) { int c; for (c = 0; c < 2; c++) { dz->c[c] += sigma * dgaussrand(); } } void perturb_weight(double *w, double sigma) { double fac = exp(sigma * dgaussrand()); (*w) *= fac; } void normalize(double v[]) { int c; double m2 = 0; for (c = 0; c < 3; c++) { double vax = v[c]; m2 += vax*vax; } double m = sqrt(m2); for (c = 0; c < 3; c++) { if (m > 0) { /* Normalize vector to unit length: */ v[c] /= m; } else { /* Messy spot, give +Z as the normal: */ v[c] = (c == 2); } } } void write_fni_image(char *fileName, float_image_t *I) { FILE *wr = open_write(fileName, TRUE); float_image_write(wr, I); fclose(wr); } void function_00(r2_t p, double *z, r2_t *dz) { /* Constant function: */ double zval = 0.50; (*z) = zval; if (dz != NULL) { dz->c[0] = dz->c[1] = 0; } } void function_01(r2_t p, double *z, r2_t *dz) { /* Linear ramp along X: */ double x = p.c[0], y = p.c[1]; double Cx = 0.5, Cy = 0.0; (*z) = Cx*x + Cy*y; if (dz != NULL) { dz->c[0] = Cx; dz->c[1] = Cy; } } void function_02(r2_t p, double *z, r2_t *dz) { /* Linear ramp along Y: */ double x = p.c[0], y = p.c[1]; double Cx = 0.0, Cy = 0.5; (*z) = Cx*x + Cy*y; if (dz != NULL) { dz->c[0] = Cx; dz->c[1] = Cy; } } void function_03(r2_t p, double *z, r2_t *dz) { /* Linear ramp along diagonal: */ double x = p.c[0], y = p.c[1]; double Cx = 1/4.0, Cy = 1/3.0; (*z) = Cx*x + Cy*y; if (dz != NULL) { dz->c[0] = Cx; dz->c[1] = Cy; } } void function_04(r2_t p, double *z, r2_t *dz) { /* Slanted parabolic hump: */ double x = p.c[0], y = p.c[1]; double A = 0.10; double B = 0.05; double Cx = 1.0, Cy = 1.5; double S = Cx*x + Cy*y; double T = -Cy*x + Cx*y; (*z) = 1.0 - A*S*S - B*T*T; if (dz != NULL) { dz->c[0] = -2.0*(A*S*Cx - B*T*Cy); dz->c[1] = -2.0*(A*S*Cy + B*T*Cx); } } void function_05(r2_t p, double *z, r2_t *dz) { /* Buried sphere: */ double x = p.c[0], y = p.c[1]; double zmin = 0.05; /* Ground level. */ /* Initialize {*z,*dz} with the ground plane: */ (*z) = zmin; if (dz != NULL) { dz->c[0] = dz->c[1] = 0.00; } /* Where the sphere is above the ground, set {*z,*dz} to it: */ double R = 0.8; /* Sphere radius. */ double D2 = x*x + y*y; if (D2 < R*R) { double S = R*R - D2; double dSdx = -2*x; double dSdy = -2*y; double F = sqrt(S); if (F > zmin) { double dFdS = 0.5/F; (*z) = F; if (dz != NULL) { dz->c[0] = dFdS*dSdx; dz->c[1] = dFdS*dSdy; } } } } void function_06(r2_t p, double *z, r2_t *dz) { /* Pentagonal pyramid: */ double x = p.c[0], y = p.c[1]; double zmin = 0.1; /* Ground level. */ /* Initialize {*z,*dz} with the ground plane: */ (*z) = zmin; if (dz != NULL){ dz->c[0] = dz->c[1] = 0.00; } /* Where the pyramid is above the ground, set {*z,*dz} to it: */ double R = 0.75; /* Pyramid radius at z = 0. */ double H = 0.80; /* Pyramid height. */ double D2 = x*x + y*y; if (D2 == 0) { /* At the tip, assume slope = 0: */ (*z) = H; } else { double ang = atan2(y,x); double dang = 2*M_PI/5; double rang = dang*((int)floor((ang + 2*M_PI)/dang + 0.83) - 0.33); double rx = cos(rang), ry = sin(rang); double S = rx*x + ry*y; if (S < R) { double dSdx = rx; double dSdy = ry; double F = H*(1 - S/R); if (F > zmin) { double dFdS = -H/R; (*z) = F; if (dz != NULL) { dz->c[0] = dFdS*dSdx; dz->c[1] = dFdS*dSdy; } } } } } void function_07(r2_t p, double *z, r2_t *dz) { /* Conical mound: */ double x = p.c[0], y = p.c[1]; double zmin = 0.1; /* Ground level. */ /* Initialize {*z,*dz} with the ground plane: */ (*z) = zmin; if (dz != NULL){ dz->c[0] = dz->c[1] = 0.00; } double R = 0.75; /* Cone radius at z = 0. */ double H = 0.80; /* Cone height. */ double D2 = x*x + y*y; if (D2 == 0) { /* At the tip, assume slope = 0: */ (*z) = H; } else { /* Where the cone is above the ground, set {*z,*dz} to it: */ if (D2 < R*R) { double S = sqrt(D2); double dSdx = x/S; double dSdy = y/S; double F = H*(1 - S/R); if (F > zmin) { double dFdS = -H/R; (*z) = F; if (dz != NULL) { dz->c[0] = dFdS*dSdx; dz->c[1] = dFdS*dSdy; } } } } } void function_08(r2_t p, double *z, r2_t *dz) { /* Ripples: */ double x = p.c[0], y = p.c[1]; double A = 0.450; /* Ripple amplitude (one half of peak-to-peak). */ double R = 0.075; /* Min ripple radius. */ double W = 1.5*M_PI/M_LN2; /* Frequency scale factor. */ double R2 = R*R; double D2 = x*x + y*y + R*R; double t = log(D2/R2); double dtdx = 2*x/D2; double dtdy = 2*y/D2; double F = A*cos(W*t); double dFdt = -A*W*sin(W*t); (*z) = F; if (dz != NULL) { dz->c[0] = dFdt*dtdx; dz->c[1] = dFdt*dtdy; } } void function_babel(r2_t p, double RI, double RO, double N, double EF, double *z, r2_t *dz) { /* Tower of Babel: */ double x = p.c[0], y = p.c[1]; double GZ = 0.1; /* Height of ground floor. */ double TZ = 0.80; /* Height of top platform. */ double DR = (RO - RI)/(N+1); /* Width of ramp. */ double DZ = (TZ - GZ)/N; /* Height difference between adjacent turns. */ double EP = EF*DR; /* Width of shoulder. */ /* Convert {x,y} to polar {r,t}: */ double r2 = x*x + y*y; double r = sqrt(r2); /* Distance from center. */ double t = atan2(y, x)/(2*M_PI); /* Argument of {(x,y)} in turns. */ /* Compute height {mz} of ramp, ignoring shoulders and Z-clipping: */ double rpt = RI + (1 - t)*DR; /* Radius of platform at this {t}. */ double s = r - rpt; /* Radial distance from platform. */ int k = floor(s/DR); /* Index of ramp turn (0 = nearest to center). */ double ph = fmax(0,s - k*DR); /* Distance from inner edge of ramp, ignoring shoulder. */ double mz = TZ - DZ*(1 + k - t); /* Height without shoulders or clipping. */ if (mz >= TZ) { /* Central platform: */ (*z) = TZ; if (dz != NULL) { dz->c[0] = dz->c[1] = 0.00; } } else if ((mz <= GZ) && (ph > EP)) { /* The ramp here is below ground and we are not on the shoulder, so we are on the ground: */ (*z) = GZ; if (dz != NULL) { dz->c[0] = dz->c[1] = 0.00; } } else if (mz <= GZ - DZ) { /* The ramp here is buried so deeply that there is not even the shoulder, so we are on the ground: */ (*z) = GZ; if (dz != NULL) { dz->c[0] = dz->c[1] = 0.00; } } else { /* We are on the ramp or shoulder (including the outermost shoulder): */ double drdx = x/r; double drdy = y/r; double dtdx = -y/r2/(2*M_PI); double dtdy = +x/r2/(2*M_PI); double dmdx = DZ*dtdx; double dmdy = DZ*dtdy; /* First compute {z,dz} ignoring the shoulder fill: */ if (mz < GZ) { /* Ramp is buried so we would be on the ground: */ (*z) = GZ; if (dz != NULL) { dz->c[0] = dz->c[1] = 0.00; } } else { /* We would be on the ramp: */ (*z) = mz; if (dz != NULL) { dz->c[0] = dmdx; dz->c[1] = dmdy; } } if (ph < EP) { /* We are on the shoulder region, add the shoulder fill: */ /* Get the relative position {w} across the shoulder, in [0_1]: */ double w = ph/EP; /* Compute the the derivatives of {w}: */ double dsdx = drdx + DR*dtdx; double dsdy = drdy + DR*dtdy; double dwdx = dsdx/EP; double dwdy = dsdy/EP; /* We need a cubic Hermite shoulder {h(w)} of unit width and height: */ double h = (2*w + 1)*(w - 1)*(w - 1); double dhdw = 6*w*(w - 1); double dhdx = dhdw*dwdx; double dhdy = dhdw*dwdy; /* Now get the height {a} of the shoulder: */ double a, dadx, dady; if (mz + DZ > TZ) { /* Innermost shoulder: */ a = TZ - mz; dadx = -dmdx; dady = -dmdy; } else if (mz < GZ) { /* Outermost shoulder. */ a = mz + DZ - GZ; dadx = dmdx; dady = dmdy; } else { /* Full-height shoulder: */ a = DZ; dadx = dady = 0; } /* Now add the shoulder of height {a} to {mz}: */ (*z) += a*h; if (dz != NULL) { dz->c[0] += dadx*h + a*dhdx; dz->c[1] += dady*h + a*dhdy; } } } } void function_09(r2_t p, double *z, r2_t *dz) { /* Tower of Babel with soulders: */ double N = 3.0; /* Number of full turns. */ double RI = 0.05; /* Inner radius of top platform (excl. shoulder). */ double RO = 0.90; /* Outer radius of tower (excl. shoulder). */ double EF = 1.0/3.0; /* Relative width of soft shoulder. */ function_babel(p, RI, RO, N, EF, z, dz); } void function_10(r2_t p, double *z, r2_t *dz) { /* Diverging parabolic ramps, with discontinuity: */ double x = p.c[0], y = p.c[1]; double Cx = 1.0, Cy = 1.5; double S = Cx*x + Cy*y; double T = -Cy*x + Cx*y; double A = 0.10; if (S < 0) { /* Flat region: */ (*z) = 0; if (dz != NULL) { dz->c[0] = 0; dz->c[1] = 0; } } else if (T > 0) { /* Rising parabolic ramp: */ (*z) = + A*S*S; if (dz != NULL) { dz->c[0] = + 2*A*S*Cx; dz->c[1] = + 2*A*S*Cy; } } else { /* Descending parabolic ramp: */ (*z) = - A*S*S; if (dz != NULL) { dz->c[0] = - 2*A*S*Cx; dz->c[1] = - 2*A*S*Cy; } } } void function_wave(r2_t *p, r2_t *f, double phase, double *z, r2_t *dz) { double t = 2*M_PI*r2_dot(p, f) + phase; double fm = r2_norm(f); double A = 1.5/fmax(0.5, 2.0*M_PI*fm); (*z) = A*sin(t); if (dz != NULL) { double ct = A*2*M_PI*cos(t); dz->c[0] = ct*f->c[0]; dz->c[1] = ct*f->c[1]; } } void function_11(r2_t p, double *z, r2_t *dz) { /* A low-frequency sinusoid wave: */ r2_t f = (r2_t){{ 0.25, 0.50 }}; function_wave(&p, &f, 0.0, z, dz); } void function_12(r2_t p, double *z, r2_t *dz) { /* A low-frequency co-sinusoid wave: */ r2_t f = (r2_t){{ 0.25, 0.50 }}; function_wave(&p, &f, M_PI/2, z, dz); } void function_13(r2_t p, double *z, r2_t *dz) { /* A medium-frequency sinusoid wave: */ r2_t f = (r2_t){{ 3.00, 2.50 }}; function_wave(&p, &f, 0.0, z, dz); } void function_14(r2_t p, double *z, r2_t *dz) { /* A medium-frequency co-sinusoid wave: */ r2_t f = (r2_t){{ 3.00, 2.50 }}; function_wave(&p, &f, 0.0, z, dz); } void function_15(r2_t p, double *z, r2_t *dz) { /* A medium-frequency sinusoid wave: */ r2_t f = (r2_t){{ 5.00, 7.00 }}; function_wave(&p, &f, 0.0, z, dz); } void function_16(r2_t p, double *z, r2_t *dz) { /* A medium-frequency co-sinusoid wave: */ r2_t f = (r2_t){{ 5.00, 7.00 }}; function_wave(&p, &f, 0.0, z, dz); } void function_21(r2_t p, double *z, r2_t *dz) { /* First, low frequency sinusoid component*/ double z_low; r2_t dz_low; r2_t f_low = (r2_t){{ 0.25, 0.50 }}; function_wave(&p, &f_low, 0.0, &z_low, &dz_low); /* Second the high frequency co-sinuisoid component*/ double z_high; r2_t dz_high; r2_t f_high = (r2_t){{ 7.00, -5.00 }}; function_wave(&p, &f_high, M_PI/2, &z_high, &dz_high); *z = z_low + z_high; if(dz != NULL) r2_add(&dz_low,&dz_high,dz); } void function_22(r2_t p, double *z, r2_t *dz) { /* A low-frequency co-sinusoid wave: */ r2_t f = (r2_t){{ 1.1, 0.0 }}; function_wave(&p, &f, 0.0, z, dz); } void function_23(r2_t p, double *z, r2_t *dz) { /* A low-frequency co-sinusoid wave: */ r2_t f = (r2_t){{ 1.1, 1.1 }}; function_wave(&p, &f, 0.0, z, dz); } void function_17(r2_t p, double *z, r2_t *dz) { /* Tower of Babel without soulders: */ double N = 2.25; /* Number of full turns. */ double RI = 0.05; /* Inner radius of top platform. */ double RO = 0.90; /* Outer radius of tower. */ double EF = 0.0; /* Relative width of soft shoulder. */ function_babel(p, RI, RO, N, EF, z, dz); } void function_18(r2_t p, double *z, r2_t *dz) { /* Tilted bicubic ramp: */ double x = p.c[0], y = p.c[1]; /* Tilt {theta}, half-side {A} and height {H} of ramp: */ double theta = M_PI/7; double Ct = cos(theta); double St = sin(theta); double A = 0.75/(fabs(Ct) + fabs(St)); double H = 0.80; /* Remap {x,y} to {u,v} tilted by theta: */ double u = (+ Ct*x + St*y)/A; double v = (- St*x + Ct*y)/A; double Vmax = 0.75; double Umax = 1.50; if ((u < -1.00) || (u > +Umax) || (fabs(v) > Vmax)) { /* Ground floor: */ (*z) = 0; if (dz != NULL) { dz->c[0] = 0; dz->c[1] = 0; } } else if (u > +1.00) { /* Upper floor: */ (*z) = H; if (dz != NULL) { dz->c[0] = 0; dz->c[1] = 0; } } else { /* Ramp: */ (*z) = 0.25*H*(2 + 3*u - u*u*u); if (dz != NULL) { double dzdu = 0.75*H*(1 - u*u)/A; dz->c[0] = dzdu*Ct; dz->c[1] = dzdu*St; } } } void function_19(r2_t p, double *z, r2_t *dz) { /* Buried sphere with undefined regions and scattered outliers: */ /* Compute the buried sphere without noise: */ function_05(p, z, dz); if (dz == NULL) { return; } /* Now add a large value (20.0) to the slopes at random places: */ bool_t bad = FALSE; /* For now. */ double x = p.c[0], y = p.c[1]; /* Make it bad outside the dome: */ double R = 0.95; /* Slightly bigger than sphere radius. */ double D2 = x*x + y*y; if (D2 > R*R) { bad = TRUE; } /* Make a pizza slice bad: */ if ((y > 0) && (y < -x/3)) { bad = TRUE; } /* Drill some holes in various places: */ int M = 2; double xf = M*x - floor(M*x); double yf = M*y - floor(M*y); double xc = 0.5 + 0.2*sin(4.615*(x+y)); double yc = 0.5 + 0.2*sin(4.634*(2*x-3*y)); double hr = 0.1 + 0.05*sin(22.01*x); double hd = hypot(xf-xc, yf-yc); if (hd <= hr) { bad = TRUE; } /* Add 1.0 to the slope where we want it to be bad: */ double Big = 20; /* Hopefully bigger than any actual slope. */ if (bad) { dz->c[0] += Big; dz->c[1] += Big; } } void function_20(r2_t p, double *z, r2_t *dz) { /* Buried annulus connected by straight ramps: */ double x = p.c[0], y = p.c[1]; auto void fix_ramp(double r, double z0, double z1, double t, double w, double d, double s, double *z, r2_t *dz); /* Defines the {z,dz} values inside a ramp of width {w} and total length {d}. The ramp connects a region at level {z0} (bottom) to a region at level {z1} (top), separated by a circle with radius {r}. The ramp is located at angle {t} from the X-axis. The ramp consists of a flat region with length {s} and a slanted plane with length {d-s}. The ramp is directed outwards if {d>0}, inwards if {d<0}. The sign of {s} must be the same as that of {d}. */ void fix_ramp(double r, double z0, double z1, double t, double w, double d, double s, double *z, r2_t *dz) { /* Rotate {x,y} so that the ramp is at +X: */ double ct = cos(t), st = sin(t); double xr = + ct*x + st*y; double yr = - st*x + ct*y; if ((fabs(yr) <= w/2) && (xr > 0)) { /* Compute X of beginning of ramp: */ double xbeg = (d > 0 ? sqrt(r*r - w*w/4) - 0.00001*r : 1.00001*r); /* Compute the X position relative to the begin & end of ramp: */ double u = (xr - xbeg)/d; if ((u >= 0) && (u < 1.0)) { double ubot = s/d; if (u < ubot) { (*z) = z0; if (dz != NULL) { (*dz) = (r2_t){{ 0, 0 }}; } } else { double v = (u - ubot)/(1.0 - ubot); (*z) = (1-v)*z0 + v*z1; if (dz != NULL) { /* Compute gradient relative to {xr,yr}: */ double dzdxr = (z1 - z0)/(d - s); /* Rotate gradient to original position: */ dz->c[0] = dzdxr*ct; dz->c[1] = dzdxr*st; } } } } } double zA = 0.10; /* Background level. */ double zB = 0.00; /* Annulus level. */ double zC = 0.05; /* Inner platform level. */ /* Initialize {*dz} with a plane: */ if (dz != NULL) { dz->c[0] = dz->c[1] = 0.00; } /* Determine the basic floor level: */ double rA = 0.75; /* Outer annulus radius. */ double rB = 0.50; /* Outer annulus radius. */ double r2 = x*x + y*y; if (r2 > rA*rA) { (*z) = zA; } else if (r2 < rB*rB) { (*z) = zC; } else { (*z) = zB; } /* Excavate outer ramp: */ double tA = 0.80; /* Position in radians. */ double wA = 0.20; /* Width. */ double dA = +0.30; /* Depth. */ double sA = +0.05; /* Threshold depth. */ fix_ramp(rA, zB, zA, tA, wA, dA, sA, z, dz); /* Excavate inner ramp: */ double tB = 2.50; /* Position in radians. */ double wB = 0.16; /* Width. */ double dB = -0.24; /* Depth. */ double sB = -0.04; /* Threshold depth. */ fix_ramp(rB, zB, zC, tB, wB, dB, sB, z, dz); } options_t *parse_options(int 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 *)malloc(sizeof(options_t)); argparser_skip_parsed(pp); argparser_get_keyword(pp, "-function"); o->ZFunc = argparser_get_next_int(pp, MIN_ZFUNC, MAX_ZFUNC); argparser_get_keyword(pp, "-size"); o->NX = argparser_get_next_int(pp, 1, 4095); o->NY = argparser_get_next_int(pp, 1, 4095); o->smoothZ = parse_sampling_options(pp, "-smoothZ"); o->smoothG = parse_sampling_options(pp, "-smoothG"); o->numGrad = argparser_keyword_present(pp, "-numGrad"); argparser_get_keyword(pp, "-maxGDiff"); o->maxGDiff = argparser_get_next_double(pp, 0.0, +DBL_MAX); argparser_get_keyword(pp, "-noiseG"); o->sigmaG = argparser_get_next_double(pp, 0.0, +DBL_MAX); argparser_get_keyword(pp, "-noiseW"); o->sigmaW = argparser_get_next_double(pp, 0.0, +DBL_MAX); argparser_get_keyword(pp, "-outPrefix"); o->outPrefix = argparser_get_next(pp); argparser_finish(pp); return o; } sampling_opts_t parse_sampling_options(argparser_t *pp, char *key) { sampling_opts_t smooth; if (argparser_keyword_present(pp, key)) { smooth.L = argparser_get_next_int(pp, 0, 2); smooth.N = argparser_get_next_int(pp, 1, (smooth.L == 0 ? 1 : 65)); } else { smooth.L = 0; smooth.N = 1; } return smooth; } /* COPYRIGHT, AUTHORSHIP, AND WARRANTY NOTICE: ** ** Copyright � 2005 by the State University of Campinas (UNICAMP). ** ** Created on 2005-08-15 by Jorge Stolfi, IC-UNICAMP. ** ** Permission to use, copy, modify, and redistribute this software and ** its documentation for any purpose and without fee is hereby ** granted, provided that: (1) the copyright notice at the top of this ** file and this copyright, authorship, and warranty notice is retained ** in all derived source files and documentation; (2) no executable ** code derived from this file is published or distributed without the ** corresponding source code; and (3) these same rights are granted to ** any recipient of such code, under the same conditions. ** This software is provided "as is", WITHOUT ANY EXPLICIT OR IMPLICIT ** WARRANTIES, not even the implied warranties of merchantibility and ** fitness for a particular purpose. END OF NOTICE. */