/* Last edited on 2024-12-07 01:33:09 by stolfi */ scene->bg = (frgb_t){{ 0.0, 0.0, 0.0 }}; scene->fg = (frgb_t){{ 1.0, 1.0, 1.0 }}; double cbd[2]; /* Splitting coordinates. */ /* A tree node defines two planes {P[0],P[1]} perpendicular to the given {axis} at coordinates {cbd[0]} and {cbd[1]}. Objects in subtree [child[0]} lie entirely in the low half-space of the plane {P[0]}, meaning that every point of their surface has coordinate {axis} less than or equal to {cbd[0]}. Likewise all objects in subtree {child[1]} lie entirely in the high half-space of {P[1]}. */ /* Note that one may have {cbd[0] > cbd[1]}, meaning that the two halfspaces intersect. If {child[0]} is {NULL} then {cbd[0]} is {-INF}, and if {child[1]} is {NULL} then {cbd[1]} is {+INF}. If both are null then {cbd[0],cbd[1]} are irrelevant, and are set to {NAN,NAN}. */ interval_t bbox[3]; /* Bounding box of all objects. */ /* The {bbox} field is a tight bounding box of all objects {objs[0..NO-1]}, excluing the floor. If {NO} is zero, {bbox} is an empty box. */ /* and their bounding box {bbox} will be set to empty. */ scene->bbox[j] = (interval_t){{ +1, -1 }}; /* Empty interval. */ tr->cbd[ic] = dir*INF; if (ko_min > ko_max) { tr->child[ic] = NULL; } else { for (int32_t ko = ko_min; ko <= ko_max; ko++) { multifok_scene_object_t *ok = &(objs[k0]); double ck = ok->ctr.c[tr->axis] - dir*ok->rad; if (dir*(ck - tr->cbd[ic]) < 0) { tr->cbd[ic] = ck; } } tr->cbd[ic] -= dir*fudge; /* Shift bounding plane to allow for roundoff. */ interval_t bboxi2[2]; bboxi2[0] = box[0]; bboxi2[1] = box[1]; bboxi2[tr->axis].end[1-ic] = tr->cbd[ic]; } } /* Consistency check of {box}: */ for (int8_t ax = 0; ax < 2; ax++) { assert(! interval_is_empty(&(box[ax]))); assert(tr->obj->ctr.c[ax] - tr->obj->rad >= box[ax].end[0] - fudge); assert(tr->obj->ctr.c[ax] + tr->obj->rad <= box[ax].end[1] + fudge); } /* contained in two sub-boxes {B[0],B[1]} of {B = box[0]×box[1]}, cut from {B} by halfspaces delimited by planes perpendicular to the {axis} */ void multifok_scene_get_object_bbox(multifok_scene_object_t *obj, interval_t bbox[]) { /* Get the preliminary range of positions {[blo[0..2] _ bho[0..2]} for the object: */ double blo[3], bhi[3]; for (int32_t j = 0; j < 3; j++) { double cctr = obj->ctr.c[j]; double crad = (obj->flat ? 0.0 : obj->rad); double cmin = cctr - crad; double cmax = cctr + crad; bbox[j] = (interval_t){{ cmin, cmax }} } } /* Elimination of {ws} from dot products. */ double ws[] /* with weights {ws[0..NS-1]} ws[i]* */ double ws[], double wk = ws[k]; double wk = ws[k]; double *multifok_basis_sample_weights(int32_t NW) { /* Create a unidimensional weight table: */ double u[NW]; /* Unidimensional weights. */ bool_t normalize = TRUE; wt_table_fill_binomial(NW, u, normalize); /* Now fill the bidimensional table: */ int32_t NS = multifok_basis_num_samples(NW); double *ws = notnull(malloc(NS*sizeof(double)), "no mem"); for (int32_t iy = 0; iy < NW; iy++) { for (int32_t ix = 0; ix < NW; ix++) { double wxy = u[iy]*u[ix]; ws[iy*NW + ix] = wxy; } } return ws; } void multifok_basis_normalize_window_samples ( int32_t NW, double x[], double ws[], double noise, double *avg_P, double *dev_P ) { int32_t NS = NW*NW; double sum_wx = 0; double sum_w = 0; for (int32_t ks = 0; ks < NS; ks++) { sum_wx += ws[ks]*x[ks]; sum_w += ws[ks]; } assert(sum_w > 0); double avg = sum_wx/sum_w; double sum_wd2 = 0; for (int32_t ks = 0; ks < NS; ks++) { double d = x[ks] - avg; sum_wd2 += ws[ks]*d*d; } double dev = sqrt(sum_wd2/sum_w); noise = fmax(1.0e-200, noise); /* To avoid division of zero by zero. */ double mag = hypot(dev, noise); for (int32_t ks = 0; ks < NS; ks++) { x[ks] = (x[ks] - avg) / mag; } (*avg_P) = avg; (*dev_P) = dev; } double multifok_basis_prod(int32_t NW, double x[], double y[], double ws[]) { int32_t NS = multifok_basis_num_samples(NW); double prod = 0.0; for (int32_t k = 0; k < NS; k++) { double xk = x[k]; double yk = y[k]; double wk = ws[k]; prod += wk*xk*yk; } return prod; } double multifok_basis_dist_sqr(int32_t NW, double x[], double y[], double ws[]) { int32_t NS = multifok_basis_num_samples(NW); double d2 = 0.0; for (int32_t k = 0; k < NS; k++) { double dk = x[k] - y[k]; double wk = ws[k]; d2 += wk*dk*dk; } return d2; } /* and a set of {sample weight} {ws[0..NS-1]} */ void multifok_basis_normalize_prod_range(int32_t NW, double ws[], int32_t NB, double *bas[]) { int32_t NS = NW*NW; for (int32_t kb = 0; kb < NB; kb++) { /* Determne the range {[sMin _ sMax]} of the dot product of {bas[kb]} and: */ double sMax = 0.0; /* Sum of positive basis coordinates. */ double sMin = 0.0; /* Sum of negative basis coordinates. */ for (int32_t ks = 0; ks < NS; ks++) { double wpk = ws[ks]*bas[kb][ks]; if (wpk > 0.0) { sMax += wpk; } if (wpk < 0.0) { sMin += wpk; } } assert((sMin <= 0) && (sMax >= 0)); /* Rescale basis samples to ensure the dot product is in the range {[-1 _ +1]}: */ double scale = fmax(-sMin, sMax); assert(scale > 0.0); for (int32_t ks = 0; ks < NS; ks++) { bas[kb][ks] /= scale; } } } /* Assumes that {coeff[0..NB-1]} are the coefficients of the window samples in some local operator basis. Then, if {NT} is positive, the procedure forms {NT} term values {term[0..NT-1]} from the coeefs {coeff[0..NB-1]} as specified in the vector {tix[0..NT-1]}. Namely, let {tix[k]} be {(ib,jb)}. If {ib} and {jb} are non-negative, then {term[k]} is the quadratic term {coeff[ib]*coeff[jb]}. If only {ib} is non-negative, then {term[k]} is {coeff[ib]}. If both are negative, then {term[k]} is 1. Finally, the procedure However, if {squared} is true, assumes that the combination of the terms is an estimate of the squared sharpness, so that the returned {score} is the square root of that combination. In any case the returned value is clipped to the range {[0 _ 1]}. If {NT} is zero, the parameters {tix} and {wt} are ignored, and the sharpness score is set to zero. The samples are destroyed in the process. */ else { /* Scale each element so that the dot product with any window in {[0_1]^NS} is always in {[-1 _ +1]}: */ multifok_basis_normalize_prod_range(NW, NB, bas); } /* Old junk */ int32_t NT = 0; string_vec_t tname = string_vec_new(50); double_vec_t wt = double_vec_new(50); while (TRUE) { fget_skip_spaces(rd); int32_t c = fgetc(rd); if (c == EOF) { break; } else if (c == '#') { do { c = fgetc(rd); } while ((c != '\n') && (c != EOF)); if (c == EOF) { break; } continue; } else if (c == '\n') { continue; } else { ungetc(c, rd); string_vec_expand(&tname,NT); double_vec_expand(&wt,NT); tname.e[NT] = fget_string(rd); wt.e[NT] = fget_double(rd); if (verbose) { fprintf(stderr, " %3d %+16.8f * %s\n", kt, wt.e[kt], tname.e[kt]); } NT++; fget_comment_or_eol(rd, '#'); } } string_vec_trim(&tname,NT); double_vec_trim(&wt,NT); fclose(rd); /* Extract product index table {prix} from term names: */ int32_t NP; multifok_term_prod_t *prix = NULL; multifok_term_indices_from_names(NB, belName, NT, tname.e, &NP, &prix, verbose); if (verbose) { fprintf(stderr, " * %-12s", tnk); fprintf(stderr, " = %3d %3d\n", prix[kt].jb1, prix[kt].jb2); } } int32_t multifok_basis_fill_LAPL(int32_t NW, double *bas[], char *belName[]); /* Requires {NW=3}. Sets {*NB_P} to {NB=5} and {bas[0..NB-1][0..NS-1]} to the sample coefficients used to compute the derivative operators "FX", "FY", "FXX", "FXY", and "FYY". */ int32_t multifok_basis_fill_CUBE(int32_t NW, double *bas[], char *belName[]); /* Requires {NW=3}. Sets {*NB_P} to {NB=7} and {bas[0..NB-1][0..NS-1]} to the sample coefficients used to compute the derivative operators "FX", "FY", "FXX", "FXY", "FYY", "C3", and "S3" (3-saddles). */ int32_t multifok_basis_fill_LAPL(int32_t NW, double *bas[], char *belName[]) { demand(NW == 3, "{LAPL} basis requires {NW=3}"); int32_t NB = 5; /* Actual basis size. */ assert(NB <= NW*NW); int32_t kb = 0; setbas(bas,belName, kb, "FX", -1.0, 00.0, +1.0, -1.0, 00.0, +1.0, -1.0, 00.0, +1.0f ); kb++; /* d/dX. */ setbas(bas,belName, kb, "FY", -1.0, -1.0, -1.0, 00.0, 00.0, 00.0, +1.0, +1.0, +1.0f ); kb++; /* d/dY. */ setbas(bas,belName, kb, "FXX",-1.0, +1.0, -1.0, -1.0, +1.0, -1.0, -1.0, +1.0, -1.0f ); kb++; /* d^2/dX^2. */ setbas(bas,belName, kb, "FYY",-1.0, -1.0, -1.0, +1.0, +1.0, +1.0, -1.0, -1.0, -1.0f ); kb++; /* d^2/dY^2. */ setbas(bas,belName, kb, "FXY",-1.0, 00.0, +1.0, 00.0, 00.0, 00.0, +1.0, 00.0, -1.0f ); kb++; /* d^2/dXdY. */ assert(kb == NB); return NB; } int32_t multifok_basis_fill_CUBE(int32_t NW, double *bas[], char *belName[]) { demand(NW == 3, "{CUBE} basis requires {NW=3}"); int32_t NB = 7; assert(NB <= NW*NW); int32_t kb = 0; setbas(bas,belName, kb, "FX", -1.0, 00.0, +1.0, -1.0, 00.0, +1.0, -1.0, 00.0, +1.0f ); kb++; /* d/dX. */ setbas(bas,belName, kb, "FY", -1.0, -1.0, -1.0, 00.0, 00.0, 00.0, +1.0, +1.0, +1.0f ); kb++; /* d/dY. */ setbas(bas,belName, kb, "FXX",-1.0, +1.0, -1.0, -1.0, +1.0, -1.0, -1.0, +1.0, -1.0f ); kb++; /* d^2/dX^2. */ setbas(bas,belName, kb, "FYY",-1.0, -1.0, -1.0, +1.0, +1.0, +1.0, -1.0, -1.0, -1.0f ); kb++; /* d^2/dY^2. */ setbas(bas,belName, kb, "FXY",-1.0, 00.0, +1.0, 00.0, 00.0, 00.0, +1.0, 00.0, -1.0f ); kb++; /* d^2/dXdY. */ setbas(bas,belName, kb, "S3", -1.0, +1.0, -1.0, 00.0, 00.0, 00.0, +1.0, -1.0, +1.0f ); kb++; /* Saddle3. */ setbas(bas,belName, kb, "C3", +1.0, 00.0, -1.0, -1.0, 00.0, +1.0, +1.0, 00.0, -1.0f ); kb++; /* Cosaddle3. */ assert(kb == NB); return NB; } typedef struct multifok_test_object_t { r3_t ctr; /* Center of disk/sphere. */ double rad; /* Radius of disk/sphere. */ bool_t flat; /* True if disk, false if sphere. */ frgb_t bg; /* Background color of object pattern. */ frgb_t fg; /* Foreground color of object pattern. */ } multifok_test_object_t; /* Specifies a sphere or horizontal disk with center {ctr}, radius {rad}, and color {bg[]}. */ typedef struct multifok_test_scene_t { interval_t box[3]; int32_t ND; multifok_test_object_t *obj; bool_t back_tilt; frgb_t bg; frgb_t fg; } multifok_test_scene_t; /* A test scene to generate test images from. It has {ND} objects entirely contained in the given {box}. The lowest {Z} of the {box} must be non-negative. It also has a backplane with colors ranging between {bg} and {fg}. If {back_tilt} is false, the backplane is horizontal at {Z=0}. If {back_tilt} is true, the backplane contains the bottom left and top right edges of the {box} .*/ typedef void multifok_test_pattern_t(double x, double y, double z, int32_t iobj, int32_t NC, float fs[]); /* Type of a function that computes pixel samples {fs[0..NC-1]} as a function of an arbitrary /object index/ {iobj} and {(x,y,z)}. */ #define multifok_test_SCENE_DEPTH 10.0 /* Total depth of simulated scene. */ multifok_test_scene_t *multifok_test_scene_throw ( interval_t box[], bool_t back_tilt, double rMin, double rMax, double minSep, bool_t verbose ); /* Generates a test scene with several objects picked at random inside the box {box[0..2]}. The backplane will be horizontal at {Z=0} if {back_tilt} is false, or tilted across the {box} if {back_tilt} is true (see the {back_tilt} attribute of {multifok_test_scene_t}). In the second case there will be no objects, and the parameters {rMin,rMax,minSep} are ignored. The backplane's {bg} and {fg} colors are chosen at random. If {back_tilt} is false, the number of objects is chosen by the procedure. The objects will be generated with {multifok_test_object_throw(box,rMin,rMax)}. If {minSep} is negative, the objects may overlap in {X} and {Y} and may extend partially outside the box in those directions. If {minSep} is zero or positive, the {X} and {Y} projections of the objects will be disjoint and separated by at least {minSep} pixels, and will be at least {minsep} away from the vertical walls. In any case the {box} must be strictly wider than {2*rMax} in all three axes. Also {box[2]} must be contained in {[0_ZMAX]} where {ZMAX=multifok_test_SCENE_DEPTH}. */ multifok_test_object_t multifok_test_object_throw(interval_t box[], double margin, double rMin, double rMax); /* Generates a random sphere or disk entirely contained in the given 3D domain {box[0..2]}. The object will have a random center and radius, and random but contrasting {bg} and {fg} colors. The radius {rad} of the object will in principle be chosen randomly in {[rMin _ rMax]}. The {Z} interval {box[2]} has somewhat different meaning that {box[0]} and {box[1]}. For disks, the center's {Z} will be chosen at random in {box[2]}. For balls the center's {Z} will be chosen so that the ball can fit entirely inside the {box[2]} interval; which may require slao reducing the radius {rad} to half the width of that interval {Y}. If the given {margin} is non-negative, the center's {X} and {Y} coordinates will be chosen to fit entirely within the {boxx[0]} and {box[1]} intervals, at least {margin} away from all four sides. In this case too the object's radius may have to be reduced to satisfy this constraint. If {margin} is negative, the center's {X} and {Y} will be chosen within the {box[0]} and {box[1]} intervals, and the radius will not be reduced, even if the object ends up extending partially ouside those ranges. In any case, the {box} must be wide enough to make the placement possible. */ double pt_rbl = (isfinite(pt_shr) ? 1/pt_shr : 0.0); /* Blur radius at hit point. */ pt_rbl = hypot(0.5*rPix, pt_rbl); pt_shr = rPix/pt_rbl; float vMin = -1.0e-38f; /* To avoid {NAN} values if {img} is all zeros. */ float vMax = +1.0e-38f; /* To avoid {NAN} values if {img} is all zeros. */ float_image_update_sample_range(img, 0, &vMin, &vMax); /* Ignores infinities. */ float vR = (float)fmax(fabs(vMin), fabs(vMax)); vMin = -vR; vMax = +vR; /* ---------------------------------------------------------------------- */ void multifok_score_read_term_weights_and_names ( FILE *rd, int32_t NB, char *belName[], int32_t *NP_P, multifok_term_prod_t **prix_P, int32_t *NT_P, char ***termName_P, bool_t verbose ); /* Reads from {rd} a a certain number {NT} of terms -- quadratic local operators -- and the weights of each term in a sharpness score. Each line must have fields "{kt} {wt[kt]} {termName[kt]}" where {kt} is the term index in {0..NT-1}, {wt[kt]} is a numeric weight, and {termName[kt]} is the term's name (formula). The latter must be a sum of products of pairs of basis element names, such as "FX*FX+FY*FY". See {multifok_term_indices_from_names} for the valid format of {namek} and semantics. The procedure analyzes the formula {termName[kt]}, determining the number {NP} of coeff pair products that appear in them and building a table {prix[0..NP-1} as required by {multifok_term_values_from_basis_coeffs}. The values of {NP,prix,NT,wt,termName} are returned in {*NP_P,*prix_P,*NT_P,*wt_P,*termName_P}. Comments starting with "#" and blank lines are ignored. If {verbose} is true, also prints the data to stderr. */ void multifok_score_read_term_weights_names_indices ( FILE *rd, int32_t NB, char *belName[], int32_t *NP_P, multifok_term_prod_t **prix_P, int32_t *NT_P, double **wt_P, char ***termName_P, bool_t verbose ) { int32_t NT; multifok_term_read_set_and_weights(rd, &NT, wt_P, termName_P); int32_t NP; string_vec_t termName = string_vec_new(50); double_vec_t wt = double_vec_new(50); while (TRUE) { bool_t ok = fget_test_comment_or_eol(rd, '#'); if (ok) { continue; } if (fget_test_eof(rd)) { break; } /* There is something there: */ char *tnk = fget_string(rd); double wtk = fget_double(rd); fget_comment_or_eol(rd, '#'); /* Save in tables: */ int32_t kt = NT; if (verbose) { fprintf(stderr, "%3d %+16.12f %s", kt, wtk, tnk); } string_vec_expand(&termName,NT); double_vec_expand(&wt,NT); termName.e[kt] = tnk; wt.e[kt] = wtk; NT++; } string_vec_trim(&termName,NT); double_vec_trim(&wt,NT); /* Convert term names to indices: */ int32_t NP; multifok_term_prod_t *prix; multifok_term_indices_from_names(NB, belName, NT, termName.e, &NP, &prix, verbose); if (verbose) { fprintf(stderr, "product index table:\n"); for (int32_t ip = 0; ip < NP; ip++) { multifok_term_prod_t *pri = &(prix[ip]); fprintf(stderr, "%3d %3d %3d %3d %s\n", ip, pri->jb1, pri->jb2, pri->kt, pri->pname); } } (*NT_P) = NT; (*termName_P) = termName.e; (*prix_P) = prix; (*wt_P) = wt.e; } void multifok_score_write_term_weights_and_names(FILE *wr, int32_t NT, double wt[], char *termName[]) { for (int32_t kt = 0; kt < NT; kt++) { fprintf(wr, "%3d %+12.4f %s\n", kt, wt[kt], termName[kt]); } fflush(wr); } void multifok_test_write_term_names(char *outPrefix, int32_t NT, char *termName[]); /* Writes the term names {termName[0..NT-1]} to a file "{outPrefix}-tnames.txt", one per line. */ /* ---------------------------------------------------------------------- */ /* Scene size in user units: */ double WX_scene = 2*interval_rad(&(scene->dom[0])); double WY_scene = 2*interval_rad(&(scene->dom[1])); /* Scale from image pixel {XY} to scene {XY}: */ double scale_img = fmax(((double)NX_img)/WX_scene, ((double)NY_img)/WY_scene); r3_t ctr_scene = (r3_t){{ interval_mid(&(scene->dom[0])), interval_mid(&(scene->dom[1])), interval_mid(&(scene->dom[2])) }}; r3_t pix_pos = multifok_test_pixel_to_scene((double)ix, (double)iy, NX_img, NY_img, &ctr_scene, scale_img); if (ws != NULL) { /* Apodize basis elements: */ for (int32_t kb = 0; kb < NB; kb++) { double *bask = bas[kb]; for (int32_t ks = 0; ks < NS; ks++) { bask[ks] *= ws[ks]; } } } multifok_term_set_t *multifok_score_read_term_weights_names_get_indices ( FILE *rd, int32_t NB, char *belName[], int32_t *NT_P, double **wt_P, char ***termName_P, int32_t *NP_P, multifok_term_prod_t **prix_P, bool_t verbose ) { multifok_term_read_set_and_weights(rd, NT_P, wt_P, termName_P, verbose); multifok_term_indices_from_names(NB, belName, *NT_P, *termName_P, NP_P, prix_P, verbose); } multifok_term_set_t *multifok_term_read_product_table ( FILE *rd, multifok_basis_t *basis, int32_t *NP_P, multifok_term_prod_t **prix_P, int32_t *NT_P, char ***termName_P, bool_t verbose ) { int32_t NB = basis->NB; int32_t NP_max = NB*(NB+1)/2 + NB + 1; /* Number of coeffs, canonical coeff pairs, and the unit term. */ multifok_term_prod_t *prix = talloc(NP_max, multifok_term_prod_t); for (int32_t ip = 0; ip < NP_max; ip++) { prix[ip] = (multifok_term_prod_t){ -1,-1, -1, NULL }; } int32_t NJJ = NB*NB; /* Number of possible basis element pairs {jb1,jb2} in any order. */ bool_t jjseen[NJJ]; /* {jjseen[jb1*NB+jb2]} is true if the pair {jb1,jb2} already occurred. */ for (int32_t jjb = 0; jjb < NJJ; jjb++) { jjseen[jjb] = FALSE; } int32_t NT_max = NP_max; int32_t nr_kt[NT_max]; /* Number of products assigned to each term. */ for (int32_t kt = 0; kt < NT_max; kt++) { nr_kt[kt] = -1; } int32_t NT = 0; /* Number of terms seen in file. */ int32_t NP = 0; /* Number of products read from file. */ auto void read_data_line(int32_t *jb1_P, int32_t *jb2_P, int32_t *kt_P, char **pr_P); /* Reads a non-blank data line and returns the data in {*jb1_P,*jb2_P,*kt_P,*pr_P}. */ if (verbose) { fprintf(stderr, "reading index triples...\n"); } while (TRUE) { bool_t ok = fget_test_comment_or_eol(rd, '#', NULL); if (ok) { continue; } if (fget_test_eof(rd)) { break; } /* There is something there: */ int32_t jb1, jb2, kt; char *pr; read_data_line(&jb1, &jb2, &kt, &pr); if (kt == NT) { /* Should be a new term: */ nr_kt[kt] = 1; NT++; } else { /* Previously started term: */ assert(nr_kt[kt] > 0); nr_kt[kt]++; } prix[NP] = (multifok_term_prod_t){ jb1, jb2, kt, pr }; if (verbose) { fprintf(stderr, " %3d %3d %3d %3d %s\n", NP, jb1, jb2, kt, pr); } NP++; } if (NP < NP_max) { prix = (multifok_term_prod_t*)notnull(realloc(prix, NP*sizeof(multifok_term_prod_t)), "no mem"); } if (verbose) { fprintf(stderr, "building term names...\n"); } char **termName = (char **)notnull(malloc(NT*sizeof(char*)), "no mem"); for (int32_t kt = 0; kt < NT; kt++) { termName[kt] = NULL; } for (int32_t ip = 0; ip < NP; ip++) { char *pr = prix[ip].pname; int32_t kt = prix[ip].kt; char *otm = termName[kt]; char *tm = NULL; if (otm == NULL) { asprintf(&tm, "%s", pr); } else { asprintf(&tm, "%s+%s", otm, pr); free(otm); } termName[kt] = tm; } if (verbose) { for (int32_t kt = 0; kt < NT; kt++) { fprintf(stderr, " %3d %s\n", kt, termName[kt]); } } (*NP_P) = NP; (*prix_P) = prix; (*NT_P) = NT; (*termName_P) = termName; return; /* INTERNAL IMPLEMENTATIONS */ void read_data_line(int32_t *jb1_P, int32_t *jb2_P, int32_t *kt_P, char **pr_P) { int32_t ip = fget_int32(rd); demand(ip == NP, "unexpected product index"); int32_t jb1 = fget_int32(rd); demand((jb1 >= 0) && (jb1 < NB), "invalid basis coeff index {jb1}"); int32_t jb2 = fget_int32(rd); demand((jb2 >= 0) && (jb2 < NB), "invalid basis coeff index {jb2}"); demand(jb1 <= jb2, "non-canonical product {jb1>jb2}"); int32_t jjb = jb1*NB + jb2; demand(! jjseen[jjb], "repeated produc {jb1,jb2}"); jjseen[jjb] = TRUE; int32_t kt = fget_int32(rd); demand(kt > 0, "invalid {kt}"); demand(kt <= NT, "{kt} skipped term"); char *pr_read = fget_string(rd); /* Validate {pr_read} with {jb1,jb2}: */ char *pr = NULL; asprintf(&pr, "%s*%s", belName[jb1], belName[jb2]); demand(strcmp(pr_read, pr) == 0, "product name does not match {jb1,jb2}"); free(pr_read); fget_comment_or_eol(rd, '#', NULL); demand(kt <= NT, "{kt} skipped term"); (*jb1_P) = jb1; (*jb2_P) = jb2; (*kt_P) = kt; (*pr_P) = pr; } } /* PARSED TERM TABLE I/O */ void multifok_term_read_index_table ( FILE *rd, int32_t NB, char *belName[], int32_t *NP_P, multifok_term_prod_t **prix_P, int32_t *NT_P, char ***termName_P, bool_t verbose ); /* Reads from {rd} a table describing how a certain number {NT} of terms of degree 2 are to be composed from {NP} products of pairs of {NB} basis coefficients. The file has one line for each coeff product to be evaluated. Each line has the format "{ip} {jb1} {jb2} {kt} {pname}" where {ip} is the product index in {0..NP-1}, {jb1} and {jb2} are the indices of the two basis coeffs, with {0 <= jb1 <= jb2 < NB}; {kt} is the index of the term to which this product is part; and {pname} is the textual product "{el1}*{el2}", e.g. "Xom*Xpp"; where {el1} must be {belName[jb1]} and {el2} must be {belName[jb2]}. The pairs {jb1,jb2} must be all distinct so the maximum value of {NP} is {NB*(NB+1)/2}. The quadruples {(jb1, jb2, kt, pname)} are saved in a table {prix[0..NP-1]}. The number of terms {NT}, which is between 1 and {NP}, is inferred from the indices {kt}. The procedure also assembles a formula {termName[kt]} for each term, denoting the sum of all pairs of products assigned to that term, e.g. "Xom*Xop+Xmo*Xpo". It assumes that {belName[jb]} is the name of basis element {jb}. The values of {NP,prix,NT,termName} are returned in {*NP_P,*prix_P,*NT_P,*termName_P}. Comments starting with "#" and blank lines are ignored. If {verbose} is true, also prints the data to stderr. */ void multifok_test_read_term_weights_names_get_indices ( char *fname, int32_t NB, char *belName[], int32_t *NT_P, double **wt_P, char ***termName_P, int32_t *NP_P, multifok_term_prod_t **prix_P, bool_t verbose ) { FILE *rd = open_read(fname, TRUE); multifok_score_read_term_weights_names_get_indices ( rd, NB, belName, NT_P, wt_P, termName_P, NP_P, prix_P, verbose ); fclose(rd); } void multifok_test_read_term_index_table ( char *fname, int32_t NB, char *belName[], int32_t *NP_P, multifok_term_prod_t **prix_P, int32_t *NT_P, char ***termName_P, bool_t verbose ) { FILE *rd = open_read(fname, TRUE); multifok_term_read_index_table(rd, NB, belName, NP_P, prix_P, NT_P, termName_P, verbose); fclose(rd); } void multifok_test_read_term_index_table ( char *fname, int32_t NB, char *belName[], int32_t *NP_P, multifok_term_prod_t **prix_P, int32_t *NT_P, char ***termName_P, bool_t verbose ); /* Same as {multifok_term_read_index_table}, but reads from a file {fname} instead of a {FILE} descriptor. */ float_image_t *multifok_image_grayscale_read(char *prefix, char *code, float vMin, float vMax) { char *fname = NULL; asprintf(&fname, "%s-%s.pgm", prefix, code); bool_t isMask = FALSE; double gama = multifok_image_gamma; double bias = multifok_image_bias; bool_t yup = TRUE; bool_t warn = TRUE; bool_t verbose = FALSE; float_image_t *img = float_image_read_pnm_named(fname, isMask, gama, bias, yup, warn, verbose); int32_t NC = (int32_t)img->sz[0]; demand(NC == 1, "grayscale image should have 1 channels"); if ((vMin != 0.0) || (vMax != 1.0)) { float_image_rescale_samples(img, 0, 0.0,1.0, vMin,vMax); } free(fname); return img; } void multifok_image_grayscale_write(float_image_t *img, char *prefix, char *code, float vMin, float vMax) { int32_t NC = (int32_t)img->sz[0]; demand(NC == 1, "grayscale image should have 1 channel"); float_image_t *cpy = float_image_copy(img); if ((vMin != 0.0) || (vMax != 1.0)) { float_image_rescale_samples(cpy, 0, vMin, vMax, 0.0, 1.0); } char *fname = NULL; asprintf(&fname, "%s-%s.pgm", prefix, code); bool_t isMask = FALSE; /* Assume uniform distr. of pixel values in encoding/decoding. */ double gamma = multifok_image_gamma; double bias = multifok_image_bias; bool_t yup = TRUE; bool_t warn = TRUE; bool_t verbose = TRUE; float_image_write_pnm_named(fname, cpy, isMask, gamma, bias, yup, warn, verbose); float_image_free(cpy); free(fname); } /* Conversion from image to scene coordinates: */ double scale = 2.0*fmin(interval_rad(&(dom[0]))/NX, interval_rad(&(dom[1])/NY)); r2_t ctr_scene = (r2_t){{ interval_mid(&(dom[0])), interval_mid(&(dom[1]) }}; r2_t ctr_img = (r2_t){{ 0.5*NX, 0.5*NY }}; That is, {pRay} should be {(xRay,yRay,zRay)} where {xRay = ctr[0] + s*(pSmp.c[0] - NX/2)}, {yRay = ctr[1] + s*(pSmp.c[1], assumes that the image plane is the horizontal plane in the scene at {Z} coordinate {zFoc}, and the image domain {[0 _ NX] x [0 _ NY]} on the scene plane {Z=0} fits snugly in the scene coordinates rectangle {dom[0] x dom[1]} and is concentric with it. And it returns in {*blur@@_P} the nominal blurring indicator {blur@@(R)} of the surface's image at that point, defined as the square of {2*(zHit-zFoc)/zDep}. In particular, if {zDep} is {+INF} (vertical rays), or {zHit(R)} is equal to {zFoc}, then {blur@@(R)} is {0.0}. double dz = zHit_ray - zFoc; double rdz = ((zDep == +INF) || (dz == 0) ? 0.0 : dz/zDep); double blur@@_ray = 1.0 - rdz*rdz;/* Sharpness at hit point. */ (*blur@@_P) = blur@@_ray; } void multifok_scene_estimate_pixel_zAve_zDev ( int32_t ix, int32_t iy, multifok_scene_t *scene, int32_t NX, int32_t NY, double *zAve_P, double *zDev_P ) { ?? r3_t d_scene = (r3_t){{ 0.0, 0.0, -1.0 }}; /* Ray direction in scene coords. */ double zFoc = ZMAX + 1.0; int32_t NR = 9; /* Number of rays to trace. */ double zHit[NR]; /* {Z} surface coordinates at the subpixel rays. */ int32_t kz = 0; for (int32_t ex = -1; ex <= +1; ex++) { for (int32_t ey = -1; ey <= +1; ey++) { /* Compute a point {pRay%%} equal to center of pixel {(ix,iy)} perturbed by {0.5*(ex,ey)}: */ r3_t p3_img = (r3_t){{ ix + 0.5*(ex + 1), iy + 0.5*(ey + 1), 0.0 }}; r3_t p3_scene = multifok_scene_coords_from_image_coords(&p3_img, NX, NY, scene->dom, zFoc); /* Ray-trace the scene from {p} straight down: */ multifok_scene_object_t *oHit_ray; /* Object hit, or {NULL} if floor. */ r3_t pHit_ray; /* Hit point. */ multifok_scene_ray_trace(scene, &p_scene, &d_scene, FALSE, &oHit_ray, &pHit_ray); assert(kz < NR); zHit[kz] = pHit_ray.c[2]; /* {Z} of hit points. */ kz++; } } assert(kz == NR); double sum_zHit = 0.0; for (int32_t kz = 0; kz < NR; kz++) { sum_zHit += zHit[kz]; } double zAve = sum_zHit/NR; double sum_dz2 = 1.0e-200; for (int32_t kz = 0; kz < NR; kz++) { double dz = zHit[kz] - zAve; sum_dz2 += dz*dz; } double zDev = sqrt(sum_dz2/NR); (*zAve_P) = zAve; (*zDev_P) = zDev; } void multifok_scene_estimate_pixel_zAve_zDev ( int32_t ix, int32_t iy, multifok_scene_t *scene, int32_t NX, int32_t NY, double *zAve_P, double *zDev_P ); /* Compute a quick estimate of the average {zAve} and deviation {zDev} of {Z} of scene at image pixel {ix,iy}. Returns them in {*zAve_P} and {*zDev_P}. */ /* Round up {HR} up so that {2*HR+1} is an (odd) multiple of {2*HS+1}: */ int32_t KR = (int32_t)ceil(((double)2*HR + 1)/((double)2*HS+1)); KR = 2*(KR/2) + 1; /* Round up to odd. */ int32_t HR_new = (KR*(2*HR+1) - 1)/2; if (verbose) { fprintf(stderr, "rounding {HR} up from %d to %d\n", HR, HR_new); } HR = HR_new; void multifok_scene_object_throw_colors(frgb_t *tone, frgb_t *bg, frgb_t *fg) { double wmax = +INF; for (int32_t ic = 0; ic < 3; ic++) { double vi = tone->c[ic]; bg->c[ic] = (float)(0.25*vi); wmax = fmin(wmax, 1.0 - bg->c[ic]); } for (int32_t ic = 0; ic < 3; ic++) { fg->c[ic] = (float)(bg->c[ic] + wmax); assert(fg->c[ic] <= 1.00001); fg->c[ic] = (float)fmin(fg->c[ic], 1.0); } } /* ---------------------------------------------------------------------- */ /* multifok_sampling.c */ oid multifok_sampling_choose_ray_tilts_and_weights ( uint32_t KR, uint32_t NS, uint32_t *NR_P, r2_t **tRay_P, double **wRay_P, bool_t verbose ) { uint32_t HR; if (KR == 0) { if (verbose) { fprintf(stderr, "using one vertical ray for all sampling points\n"); } HR = 0; } else { demand(KR >= 0, "invalid {KR}"); demand(NS >= 1, "invalid {NS}"); demand((NS % 2) == 1, "invalid {NS}"); /* Factor {NS} into {SS^2*MS} where {MS} is squarefree: */ uint32_t SS = 1, MS = NS; uint32_t d = 3; while (TRUE) { uint32_t d2 = d*d; if (MS < d2) { break; } else if ((MS % d2) == 0) { MS /= d2; SS *= d; } else { d += 2; } } if (verbose) { fprintf(stderr, "factored NS = %d into *%d^2*%d\n", NS, SS, MS); } assert(SS*SS*MS == NS); assert((SS*MS % 2) == 1); /* Now choose {KR_new} so that {KR_new*MS} is an odd square and {KR_new >= KR}: */ uint32_t TR = (uint32_t)ceil((sqrt(((double)KR)/((double)MS)) - 1.0)/2); uint32_t LR_loc = (2*TR+1)*MS*SS; uint32_t NR_loc = LR_loc*LR_loc; uint32_t KR_new = NR_loc/NS; HR = (LR_loc - 1)/2; if (verbose) { if (KR_new != KR) { fprintf(stderr, "rays per sampling point {KR} rounded up from %d to %d\n", KR, KR_new); } } KR = KR_new; assert(KR >= 1); } uint32_t LR = 2*HR + 1; uint32_t NR = LR*LR; fprintf(stderr, "generating %d × %d = %d total rays (HR = %d)\n", LR, LR, NR, HR); if (NR != 1) { assert((NR % NS) == 0); } uint32_t NR_gen; r2_t *tRay; double *wRay; multifok_sampling_generate_2D_hann_samples(HR, 1.0, &NR_gen, &tRay, &wRay); assert(NR_gen == NR); if (NR > 1) { /* Normalize the tilts to unit RMS radius: */ double sum_w = 0; double sum_w_r2 = 0; for (uint32_t ir = 0; ir < NR; ir++) { sum_w += wRay[ir]; double r2 = r2_norm_sqr(&(tRay[ir])); sum_w_r2 += wRay[ir]*r2; } assert(sum_w >= 1.0); double r_avg = sqrt(sum_w_r2/sum_w); assert(r_avg > 0.0); for (uint32_t ir = 0; ir < NR; ir++) { tRay[ir].c[0] /= r_avg; tRay[ir].c[1] /= r_avg; } } ... } /* multifok_sampling_generate_quadrant_hann_sampling_points_and_weights */ /* Decide the weight cutoff {wCut}: */ double wMin = wHan[0]; /* Smallest 1D weight. */ assert(wHan[NW-1] == wMin); double wSec = wHan[HV-1]; /* Second-largest weight. */ assert(wHan[HV+1] == wSec); double wCut = 0.99*wMin*wSec; if (debug) { fprintf(stderr, "discarding sampling vectors with weight less than %10.8f\n", wCut); } /* ---------------------------------------------------------------------- */ uint32_t kq = 0; /* First-quadrant vectors generated and accepted so far. */ for (uint32_t ix = 1; ix <= +HV; ix++) { double vx = R*((double)ix)/(HV+1); for (uint32_t iy = 0; iy <= +HV; iy++) { assert(kq < NVQ_max); double vy = R*((double)iy)/(HV+1); double wk = wHan[ix+HV]*wHan[iy+HV]; if (wk >= wCut) { assert(kq < NVQ_max); vq[kq] = (r2_t){{ vx, vy }}; wq[kq] = wk; kq++; } } } uint32_t NVQ = kq; /* Number of sample points actually retained in the first quadrant. */ if (debug) { fprintf(stderr, "retained %d sampling vectors in first quadrant\n", NVQ); }