/* Last edited on 2023-03-27 04:39:23 by stolfi */ #include void check_satisfaction(int m, int n, double M[], int p, double X_ref[]); /* Checks whether the {(n-p) × p} matrix {X} is still a solution of the system {A X = B} where {A} is the first {n-p} columns of {M}, and {B} is the last {p} columns. */ void check_satisfaction(int m, int n, double M[], int p, double X[]) { fprintf(stderr, " checking system's solution...\n"); int q = n - p; int i, j, k; for (k = 0; k < p; k++) { double tol = 1.0e-6; for (i = 0; i < m; i++) { double s = 0.0; for (j = 0; j < q; j++) { s += M[i*n + j]*X[j*p + k]; } s -= M[i*n + q + k]; if (fabs(s) > tol) { fprintf(stderr, "(A X - B)[%d,%d] = %22.16e tol = %22.16e\n", i, k, s, tol); demand(FALSE, "** system is no longer satisfied"); } } } } void minn_constr_throw_ortho_vector ( int32_t n, int32_t p, double C[], int32_t k, double H[] ); /* Assumes that {C} is an {n} by {p} matrix, with orthonormal columns; that {H} is a {d=n-p} by {n} matrix; that {k} is less than {d}; and that the first {k} rows of {H} are orthonormal vectors orthogonal to the columns of {C}. Stores into row {k} of {H} a another unit {n}-vector {h} that is orthogonal to all those rows and all columns of {C}. */ void minn_constr_throw_ball_vector(int32_t n, double rmin, double rmax, double a[]) { demand ((0 <= rmin) && (rmin <= 0.80*rmax), "bad norm range"); while (TRUE) { /* Throw a random alignment vector in the cube of side {2*rmax} centered at the origin, compute its squared norm: */ double sum2 = 0.0; /* Sum of squares of all coordinates. */ for (int32_t i = 0; i < n; i++) { for (int32_t j = 0; j < 2; j++) { double pij = rmax*(2*drandom() - 1.0); a[i].c[j] = pij; sum2 += pij*pij; } } /* Reject if the norm is outside the requested range:: */ if ((sum2 <= rmax*rmax) && (sum2 >= rmin*rmin)) { return; } } } void minn_constr_throw_ortho_vector(int32_t n, double rad[], int32_t k, double H[]) { bool_t debug = TRUE; if (debug) { fprintf(stderr, " computing H[%d]\n", k); } double *hk = &(H[n*k]); int32_t nit = 0; /* Counts iterations, for safety. */ while (TRUE) { nit++; assert(nit < 1000); /* Should never happen if {k < d} */ /* Generate a vector {hk} uniformly distributed in the unit ball that is not too short: */ double rmin = 0.5; /* To reduce the effect of roundoff noise on normalization. */ minn_constr_throw_ball_vector(n, rmin, 1.0, hk); /* Project {hk} perpendicular to coordinates where {rad} is zero: */ for (int32_t i = 0; i < n; i++) { for (int32_t j = 0; j < 2; j++) { double rij = rad[i].c[j]; demand (rij >= 0, "invalid {rad}"); if (rij == 0.0) { hk[i].c[j] = 0.0; } } } /* Project {hk} perpendicular to the all-ones vectors, preserving conformity: */ for (int32_t j = 0; j < 2; j++) { double sum = 0.0; int32_t nv = 0; for (int32_t i = 0; i < n; i++) { double rij = rad[i].c[j]; if (rij != 0.0) { sum += hk[i].c[j]; nv++; } } if (nv > 0) { double avg = sum/nv; for (int32_t i = 0; i < n; i++) { double rij = rad[i].c[j]; if (rij != 0.0) { hk[i].c[j] -= avg; } } } } /* Project {hk} perpendicular to the previous delta vectors. */ /* Since the previous vectors are conformal to {rad} and balanced, the projection preserves these properties: */ for (int32_t r = 0; r < k; r++) { double *hr = &(H[n*r]); double sdot = minn_constr_dot(n, hk, hr); for (int32_t i = 0; i < n; i++) { for (int32_t j = 0; j < 2; j++) { hk[i].c[j] -= sdot*hr[i].c[j]; } } } /* Check if norm is still large enough: */ double sum2 = minn_constr_norm_sqr(n, hk); if (sum2 >= rmin*rmin) { /* Normalize and return this vector: */ double norm = sqrt(sum2); for (int32_t i = 0; i < n; i++) { for (int32_t j = 0; j < 2; j++) { hk[i].c[j] /= norm; } } return; } /* Remaining vector {hk} was too short. Try again. */ } } void minn_constr_print_search_basis_vector(FILE *wr, int32_t n, int32_t ix, double Uk[], double uraddk) { int32_t nz = 0; for (int32_t i = 0; i < n; i++) { fprintf(wr, " %s", name); if (ix >= 0) { fprintf(wr, "[%d]", ix); } fprintf(wr, "[%d] = (", i); for (int32_t j = 0; j < 2; j++) { double pij = a[i].c[j]; fprintf(wr, " %12.8f", pij); if (pij == 0) { nz++; } } fprintf(wr, " )\n"); } if (ctvars) { int32_t nc = 2*n; i2_t nv = minn_constr_count_variable_coords(n, a); int32_t d = minn_constr_count_degrees_of_freedom (n, a); fprintf(wr, " nc = %d nz = %d nv = %d+%d = %d d = %d\n", nc, nz, nv.c[0], nv.c[1], nv.c[0]+nv.c[1], d); } } /* Choose two main direction vectors: */ demand(d >= 2, "not enough degrees of freedom to plot"); int32_t k0 = 0, k1 = 1; double *U0 = &(U[k0*n]); double urad0 = urad[k0]; double *U1 = &(U[k1*n]); double urad1 = urad[k1]; int32_t nd = minn_count_degrees_of_freedom(n, arad); /* Save the initial guess: */ double ctr[n]; /* Center of ellipsoid (aved initial guess). */ for (int32_t i = 0; i < n; i++) { ctr[i] = v[i]; } /* Compute the axes and radii of the search ellipsoid {\RF}: */ double U[nd*n]; /* Basis of conformal balanced delta vectors. */ double urad[nd]; minn_compute_search_ellipsoid (n, arad, nd, U, urad); /* !!! Convert {U} to the packed lattice basis !!! */ minn_enum_grid(n, F2, ctr, arad, nd, U, urad, tol, v, F2valP); void minn_enum_grid ( int32_t n, /* Number of images to align. */ minn_goal_t *F2, /* Function that evaluates the goal between the images. */ double ctr[], /* (IN/OUT) Corresponding points in each image. */ double arad[], /* Max delta vector coordinates for each image. */ int32_t nd, /* Degrees of freedom (dimension of {\RF}). */ double U[], /* Adjstment vectors parallet to axes of {\RF}. */ double urad[], /* Radii of {\RF} along those directions. */ double step, /* Grid step. */ double v[], /* (OUT) Optimum alignment vector. */ double *F2valP /* (OUT) Value of {F2} for the optimum alignment vector. */ ) { bool_t debug = FALSE; /* Compute the number of samples along each main axis of {\RF}: */ int32_t n[n]; /* Number of sample points on each side of 0 along each axis of {\RF}. */ double fudge = 1.0e-10; /* Fudge to ensure we include the edges of the grid. */ for (int32_t k = 0; k < n; k++) { n[k] = (int32_t)floor(urad[k]/step + fudge); if (debug) { fprintf(stderr, "plot axis %d length = %.8f steps = %d\n", k, urad[k], n[k]); double *uk = &(U[k*n]); minn_print_vector(stderr, n, "u", k, uk, FALSE); } } /* Enumerate all integer tuples {t[0..n-1]} where {t[k]} ranges in {-n[k]..+n[k]}: */ (*F2valP) = +INF; /* Minimum goal found so far. */ int32_t t[n]; /* Enumeration variables. */ for (int32_t k = 0; k < n; k++) { t[k] = -n[k]; } double psmp[n]; /* Sampling point. */ int32_t knext = -1; /* Next tuple elem to be incremented is {t[knext]}. */ while (knext < n) { if (debug) { fprintf(stderr, " t = ( "); for (int32_t k = 0; k < n; k++) { fprintf(stderr, " %d", t[k]); } fprintf(stderr, " )\n"); } /* !!! Should avoid generating tuples outside {\RF} !!! */ /* Compute the {urad}-relative squared norm {sum2} of the displacement: */ double sum2 = 0; /* Squared norm of {(v-ctr)/arad}. */ for (int32_t k = 0; k < nd; k++) { double ek = t[k]*step/urad[k]; sum2 += ek*ek; } if (sum2 <= 1.0 + 1.0e-8) { /* Sample point will be inside {\RF}. */ /* Build the sample point {psmp}: */ for (int32_t i = 0; i < n; i++) { for (int32_t j = 0; j < 2; j++) { psmp[i].c[j] = ctr[i].c[j]; double rij = arad[i].c[j]; if (rij != 0) { double dij = 0; for (int32_t k = 0; k < nd; k++) { double *uk = &(U[k*n]); dij += t[k]*step*uk[i].c[j]; } psmp[i].c[j] += dij; } } } /* Evaluate the goal function at {psmp}: */ double F2val = F2(n, psmp); if (F2val < (*F2valP)) { /* Update the current optimum: */ for (int32_t i = 0; i < n; i++) { v[i] = psmp[i]; } (*F2valP) = F2val; } } /* Get the next tuple {t[0..nd-1]}: */ if ((knext < 0) || (t[knext] >= n[knext])) { /* Can't increment this, try next one: */ knext++; } else { /* Increment {t[knext]} and clear all previous {t[k]}: */ t[knext]++; while (knext > 0) { knext--; t[knext] = -n[knext]; } } } } double minn_max_rel_diff (int32_t n, double a[], double b[], double arad[]); /* Given two {n}-vectors {a,b}, returns the maximum differences between the coordinates, relative to the radius vector {arad[0..n-1]}. That is, returns the maximum of {fabs(a[i] - b[i])/arad[i]} for all {i} in {0..n-1}. If {a} or {b} is {NULL}, assumes a vector of {n} zeros. Thus the search box {\RB} consists of all vectors {v} such that {minn_max_rel_diff(n,v,NULL,arad) <= 1}. is all vectors {a} such that . */ double minn_rel_dist_sqr (int32_t n, double a[], double b[], double arad[]); /* Given two {n}-vectors {a,b}, returns the total squared coordinate differences between them, relative to the radius vector {arad[0..n-1]}. That is, returns the sum of {((a[i] - b[i])/arad[i])^2} for all {i} in {0..n-1}. If {a} or {b} is {NULL}, assumes a vector of {n} zeros. Thus the search ellipsoid {\RE} consists of all vectors {a} such that {minn_rel_dist_sqr(n,a,NULL,arad) <= 1}. */ double arad[], /* Radius of search domain along each axis. */ /* INTERNAL PROTOTYPES */ void minn_enum_grid ( int32_t ni, /* Number of images to align. */ minn_mismatch_t *F2, /* Function that evaluates the mismatch between the images. */ r2_t ctr[], /* (IN/OUT) Corresponding points in each image. */ r2_t arad[], /* Max delta vector coordinates for each image. */ int32_t nd, /* Degrees of freedom (dimension of {\RF}). */ r2_t U[], /* Adjstment vectors parallet to axes of {\RF}. */ double urad[], /* Radii of {\RF} along those directions. */ double step, /* Grid step. */ r2_t p[], /* (OUT) Optimum alignment vector. */ double *F2valP /* (OUT) Value of {F2} for the optimum alignment vector. */ ); /* Enumerates every point {psmp} that is {ctr} plus a linear combination of the lattice vectors {U[k]} for {k} in {0..nv-1} with coefficients that are {step} times integers in the range {-tmax..+tmax}; where {U[k][i] = U[k*ni + i]} for {i} in {0..ni-1}. The {arad} parameter is used only to skip coordinates where {\RE} has size zero. Then evaluates {F2(ni, psmp)} at every alignment vector {psmp} that is inside the ellipsoid with center {ctr} and radius vector {arad}. Sets {*F2valP} to the minimum of those values, and sets {p[0..ni-1]} to the alignment vector {psmp} that realizes the minimum. */ /* IMPLEMENTATIONS */ double fudge = 1.0e-10; /* Fudge to ensure we include the edges of the grid. */ for (int32_t k = 0; k < n; k++) { n[k] = (int32_t)floor(urad[k]/tol + fudge); if (debug) { fprintf(stderr, "axis %d domain radius = %.8f steps = %d\n", k, arad[k], n[k]); } } double minn_constr_rel_norm_sqr (int32_t n, double a[], double b[], double arad[]); /* Given two {n}-vectors {a,b}, returns the total squared coordinate differences between them, relative to the radius vector {arad}. That is, returns the sum of {((a[i] - b[i])/arad[i])^2} for all {i} in {0..n-1}. If {b} is {NULL}, assumes a vector of {n} zeros. Thus the basic ellipsoid {\CE} consists of all vectors {a} such that {minn_constr_rel_norm_sqr(n,a,NULL,arad) <= 1}. */ void minn_constr_print_search_basis_vector(FILE *wr, int32_t n, int32_t k, double Uk[], double uradk) { fprintf(wr, "direction U[%d,*] = ", k); rn_gen_print(wr, n, Uk, "%12.7f", "[", " ", "]\n"); fprintf(wr, "radius urad[%d] = %.8f\n", k, uradk); } /* ---------------------------------------------------------------------- */ /* minn_quad.c */ /* INTERNAL PROTOTYPES */ void minn_quad_grid ( int32_t ni, /* Number of images to align. */ minn_goal_t *F2, /* Function to be minimized. */ r2_t ctr[], /* (IN/OUT) Corresponding points in each image. */ r2_t arad[], /* Max delta vector coordinates for each image. */ int32_t nd, /* Degrees of freedom (dimension of {\RF}). */ r2_t U[], /* Adjstment vectors parallet to axes of {\RF}. */ double urad[], /* Radii of {\RF} along those directions. */ double step, /* Grid step. */ r2_t p[], /* (OUT) Optimum alignment vector. */ double *F2val_P /* (OUT) Value of {F2} for the optimum alignment vector. */ ); /* Enumerates every point {psmp} that is {ctr} plus a linear combination of the lattice vectors {U[k]} for {k} in {0..nv-1} with coefficients that are {step} times integers in the range {-tmax..+tmax}; where {U[k][i] = U[k*ni + i]} for {i} in {0..ni-1}. The {arad} parameter is used only to skip coordinates where {\RE} has size zero. Then evaluates {F2(ni, psmp)} at every alignment vector {psmp} that is inside the ellipsoid with center {ctr} and radius vector {arad}. Sets {*F2val_P} to the minimum of those values, and sets {p[0..ni-1]} to the alignment vector {psmp} that realizes the minimum. */ /* IMPLEMENTATIONS */ /* Compute the number of samples along each main axis of {\RF}: */ int32_t n[n]; /* Number of sample points on each side of 0 along each axis of {\RF}. */ double fudge = 1.0e-10; /* Fudge to ensure we include the edges of the grid. */ for (int32_t k = 0; k < n; k++) { n[k] = (int32_t)floor(urad[k]/tol + fudge); if (debug) { fprintf(stderr, "axis %d domain radius = %.8f steps = %d\n", k, arad[k], n[k]); } } /* Enumerate all integer tuples {t[0..n-1]} where {t[k]} ranges in {-n[k]..+n[k]}: */ (*F2val_P) = +INF; /* Minimum goal found so far. */ int32_t t[n]; /* Enumeration variables. */ for (int32_t k = 0; k < n; k++) { t[k] = -n[k]; } double u[n]; /* Sampling vector. */ int32_t knext = -1; /* Next tuple elem to be incremented is {t[knext]}. */ while (knext < n) { if (debug) { fprintf(stderr, " t = ( "); for (int32_t k = 0; k < n; k++) { fprintf(stderr, " %d", t[k]); } fprintf(stderr, " )\n"); } /* !!! Should avoid generating tuples outside {\RF} !!! */ /* Build the sample vector {u}: */ for (int32_t i = 0; i < n; i++) { u[i] = ((double)t[i])*tol; } /* Determine whether the point is inside the domain: */ bool_t inside = (box ? TRUE : minn_rel_dist_sqr(n, u, NULL, arad) <= 1 + fudge) if (inside) { /* Evaluate the goal function at {v}: */ double F2val = F2(n, u); if (F2val < (*F2val_P)) { /* Update the current optimum: */ rn_copy(n, u, v); (*F2val_P) = F2val; } } /* Get the next tuple {t[0..nd-1]}: */ if ((knext < 0) || (t[knext] >= n[knext])) { /* Can't increment this, try next one: */ knext++; } else { /* Increment {t[knext]} and clear all previous {t[k]}: */ t[knext]++; while (knext > 0) { knext--; t[knext] = -n[knext]; } } } double image_eval(int32_t i, double vOpt[], double xsc, double ysc) { /* Returns the mother image {img}, displaced by {sopt[j]}, shrunk by {iscale.c[j]} along axis {j}, then evaluated at {(xsc,ysc)}. */ double scale = (double){{ pow(2.0, iscale.c[0]), pow(2.0, iscale.c[1]) }}; double dxsc = sopt[i].c[0]/scale.c[0]; double dysc = sopt[i].c[1]/scale.c[1]; double psc = (double){{ xsc - dxsc, ysc - dysc }}; return tmnn_tools_eval_mother_image(&psc, &scale, mom_NF, mom_frq, mom_phi, mom_amv); } double F2_imgmis(int32_t n, double v[], double vOpt[]) { double F2v = tmnn_compute_test_goal_func ( n, &image_eval, v, iscale, wsize, wtx, wty ); return F2v; } void tmnn_choose_initial_guess(double d, double urad[], double sini[]); /* Chooses and stores into {sini[0..d-1]} the coeeficients of the initial guess, in the search basis {U}; assuming that the radii if the search ellipsoid {\RF} are {urad[0..d-1]}. */ /* Choose base radii and constraints: */ double arad[n]; /* Domain radius along each axis. */ tmnn_choose_radii(n, FALSE, arad, it); /* No zeros. */ tmnn_choose_radii(n, arad); double A[q*n]; tmnn_choose_constraints(n, q, A); tmnn_test_minn_subspace(n, d, U, urad, it, verbose); bool_t rad_zeros = TRUE; /* Some radii may be zero. */ tmnn_choose_radii(n, TRUE, arad); /* Allow zeros. */ } int32_t v; double *C = NULL; minn_normalize_constraints(n, arad, q, A, verbose, &v, &C); tmnn_check_constraint_normalization(n, arad, q, A, v, C); /* Compute the search ellipsoid: */ int32_t d = n - v; double U[d*n]; double urad[d]; minn_compute_search_ellipsoid(n, arad, v, C, d, U, urad); tmnn_check_search_ellipsoid(n, arad, v, C, d, U, urad); /* The mother image: */ /* Data for constituent waves: */ int32_t mom_NF = 30; /* Number of waves. */ double mom_amv[mom_NF]; /* Raw amplitude. */ double mom_phi[mom_NF]; /* Initial wave phase (fraction of cycle). */ double mom_frq[mom_NF]; /* Wave frequency in X and Y (cycles per pixel). */ /* Set auxiliary parameters {sini[],sopt[],arad[],astv[]}: */ double sopt[n]; /* Actual optimum alignment. */ double sini[n]; /* Initial alignment. */ double arad[n]; /* Search radius for each coordinate. */ double astv[n]; /* Search step or precision for each coordinate. */ tmnn_choose_sini_arad_astp(n, adj_rad, adj_stp, sini, arad, astp); /* Choose the goal function {F2_raw} for optimization according to {fname}. */ if (strcmp(fname, "indiff") == 0) { F2_raw = &F2_indiff; /* For this function, the optimum {sopt} is the initial point {sini}: */ for (int32_t i = 0; i < n; i++) { sopt[i] = sini[i]; } compare_to_sopt = bias; /* If {bias} is true, the solution should be {sopt}. */ } else if (strcmp(fname, "optdst") == 0) { F2_raw = &F2_optdst; /* Choose an optimum {sopt} near the initial point {sini}: */ tmnn_choose_optimum(n, sini, arad, sopt); compare_to_sopt = ! bias; /* If {bias} is false, the solution should be {sopt}. */ } else if (strcmp(fname, "imgmis") == 0) { F2_raw = &F2_imgmis; /* Initialize the parameters of the mother image: */ tmnn_tools_initialize_mother_image(mom_NF, mom_frq, mom_phi, mom_amp, TRUE); /* Choose the optimum displacements {sopt} near the initial point {sini}: */ tmnn_choose_optimum(n, sini, arad, sopt); compare_to_sopt = ! bias; /* If {bias} is false, the solution should be {sopt}. */ /* Write the images out for debuging: */ for (int32_t isc = -1; isc < 3; isc++) { i2_t ipscale = ( isc < 0 ? iscale : (i2_t){{ isc, isc }} ); /* Plotting scale. */ for (int32_t i = 0; i < n; i++) { tmnn_write_test_image(image_eval, i, ipscale, cmp_rad, sini[i], "ini", arad[i]); tmnn_write_test_image(image_eval, i, ipscale, cmp_rad, sopt[i], "opt", arad[i]); } } } else { demand(FALSE, "invalid goal function name"); } /* Print raw function value and bias term for initial point: */ double F2ini = F2_raw(n, sini, iscale); double b2ini = (bias ? bias_term(sini) : 0); tmnn_debug_points(0, "initial guess ", n, "sini", sini, NULL, NULL, arad, astp, F2ini, b2ini); /* Print raw function value and bias term for optimum point: */ double F2opt = F2_raw(n, sopt, iscale); double b2opt = (bias ? bias_term(sini) : 0); tmnn_debug_points(0, "actual optimum", n, "sopt", sopt, sini, NULL, arad, astp, F2opt, b2opt); /* Plot the goal function in the neighborhood of the initial guess and optimum point: */ for (int32_t isc = -1; isc < 3; isc++) { i2_t ipscale = ( isc < 0 ? iscale : (i2_t){{ isc, isc }} ); /* Plotting scale. */ fprintf(stderr, "plotting goal function at scale (%d,%d)...\n", ipscale.c[0], ipscale.c[1]); nF2 = 0; tmnn_plot_goal(F2_full, n, ipscale, sini, "ini", arad); fprintf(stderr, "%d calls to {F2_full} for plotting.\n", nF2); nF2 = 0; tmnn_plot_goal(F2_full, n, ipscale, sopt, "opt", arad); fprintf(stderr, "%d calls to {F2_full} for plotting.\n", nF2); } /* If search step is zero, just plot, skip optimization: */ if ((adj_stp.c[0] == 0.0) && (adj_stp.c[1] == 0.0)) { fprintf(stderr, "no optimization: search step is zero\n"); return; } /* Call the optimizer: */ fprintf(stderr, "optimizing"); debug_points = TRUE; double psol[n]; /* Computed alignment vector. */ double F2sol; /* Goal function at {psol} including bias. */ nF2 = 0; for (int32_t i = 0; i < n; i++) { psol[i] = sini[i]; } if (monoscale) { if (quadopt) { fprintf(stderr, " (single_scale_quadopt)...\n"); r2_opt_single_scale_quadopt(n, iscale, F2_full, arad, astp, psol, &F2sol, debug_opt); } else { fprintf(stderr, " (single_scale_enum)...\n"); r2_opt_single_scale_enum(n, iscale, F2_full, arad, astp, psol, &F2sol, debug_opt); } } else { fprintf(stderr, " (multi_scale)...\n"); r2_opt_multi_scale(n, F2_full, quadopt, arad, astp, psol, &F2sol, debug_opt); } fprintf(stderr, "done optimizing.\n"); fprintf(stderr, "%d calls to {F2_full}.\n", nF2); /* Print raw function value and bias term for computed optimum: */ double F2_raw_sol = F2_raw(n, psol, iscale); double bias_sol = (bias ? bias_term(psol) : 0); tmnn_debug_points(0, "computed optimum", n, "psol", psol, sini, sopt, arad, astp, F2_raw_sol, bias_sol); assert(F2sol == F2_raw_sol + bias_sol); /* Print again raw function value and bias term for optimum point: */ tmnn_debug_points(0, "actual optimum", n, "sopt", sopt, sini, NULL, arad, astp, F2opt, b2opt); bool_t ok = tmnn_check_ellipsoid(n, psol, F2sol, sini, (compare_to_sopt ? sopt : NULL), arad, astp); if(! ok) { exit(1); } return; /* INTERNAL IMPLEMENTATIONS */ } int32_t nF2 = 0; /* Counts calls to {nF2}. */ auto double F2(int32_t n, double v[]); /* Goal function. */ double vSol[]; /* Computed minimum. */ dpuble F2Sol = NAN; minn_uniform(n, F2, dMax, box, atol, meth, vSol, &F2Sol); tmnn_check_solution(n, F2, vSol, F2Sol, vOpt, atol, NAN); } void tmnn_make_weight_table(int32_t hw, int32_t *nwp, double **wtp) { demand (hw >= 0, "invalid sampling radius"); int32_t nw = 2*hw + 1; double *wt = notnull(malloc(nw*sizeof(double)), "no mem"); bool_t norm = TRUE; wt_table_fill_hann(nw, wt, norm); (*nwp) = nw; (*wtp) = wt; } void tmnn_choose_sini_arad_astp ( int32_t n, double adj_rad, double adj_stp, double sini[], double arad[], double astp[] ) { /* Coordinates {kfix + r*mfix} for integer {r} will be fixed: */ int32_t kfix = 1; int32_t mfix = 7; /* Define the parameters: */ int32_t k = 0; /* Counts coordinates. */ for (int32_t i = 0; i < n; i++) { /* Choose the initial guess {sini}: */ sini[i] = (double){{ 50.0, 50.0 }}; for (int32_t j = 0; j < 2; j++) { if (((k - kfix) % mfix == 0) || (adj_rad.c[j] == 0)) { /* Do not optimize this coordinate: */ arad[i].c[j] = 0.0; astp[i].c[j] = 0.0; } else { /* Choose the search radius {arad[i].c[j]} by wriggling {adj_rad.c[j]}: */ arad[i].c[j] = (0.85 + 0.15*sin(k))*adj_rad.c[j]; if (adj_stp.c[j] == 0.0) { /* Plot but do not optimize this coordinate: */ astp[i].c[j] = 0.0; } else { /* Choose the enum step {astp[i].c[j]} by wriggling {adj_stp.c[j]}: */ astp[i].c[j] = (0.80 + 0.20*sin(k+0.3))*adj_stp.c[j]; } } } } tmnn_debug_params(0, "rad", n, arad); tmnn_debug_params(0, "stp", n, astp); } void tmnn_choose_optimum(int32_t n, double sini[], double arad[], double sopt[]) { /* Set {sopt[0..n-1]} to {sini[0..n-1]}, balanced. */ for (int32_t j = 0; j < 2; j++) { double da[n]; /* Displacements {sopt[i]-sini[i]} along axis {j}. */ /* Set {da[0..n-1]} to random displcements; compute their sum: */ double sumdvar = 0; /* Sum of all {da[0..n-1]} excluding fixed ones. */ int32_t nvar = 0; /* Number of variable coordinates. */ for (int32_t i = 0; i < n; i++) { double ra = arad[i].c[j]; /* Search radius for {v[i]} along axis {j}. */ if (ra == 0.0) { /* Coordinate is fixed: */ da[i] = 0.0; } else { /* Coordinate is variable: */ assert(ra > 0.0); da[i] = (2*drandom()-1)*ra; sumdvar += da[i]; nvar ++; } } /* Rebalance {da} to zero sum, add to {sini} to get {sopt}: */ double avg = sumdvar/nvar; for (int32_t i = 0; i < n; i++) { double ra = arad[i].c[j]; if (ra > 0.0) { da[i] = da[i] - avg; } sopt[i].c[j] = sini[i].c[j] + da[i]; } } } double tmnn_compute_test_goal_func ( int32_t n, /* Number of images being compared. */ tmnn_image_eval_proc_t *eval, /* Image evaluator. */ double v[], /* Unscaled sampling grid center for each image. */ double vOpt[], /* Image shrink scale in each axis. */ i2_t wsize, /* Width of comparison window in each axis (must be odd). */ double wtx[], /* Weight table for X displacement. */ double wty[] /* Weight table for Y displacement. */ ) { /* Get half-widths of sampling window: */ demand((wsize.c[0] % 2) == 1, "window width must be odd"); demand((wsize.c[1] % 2) == 1, "window height must be odd"); int32_t hwx = (wsize.c[0]-1)/2; int32_t hwy = (wsize.c[1]-1)/2; double scale = (double){{ pow(2.0, iscale.c[0]), pow(2.0, iscale.c[1]) }}; /* Enumerate the sampling points, and get the variances of the values: */ double sum_wtxy_var = 0; double sum_wtxy = 0; for (int32_t jy = -hwy; jy <= +hwy; jy++) { for (int32_t jx = -hwx; jx <= +hwx; jx++) { /* Sample the images at this grid sampling point: */ double val[n]; for (int32_t i = 0; i < n; i++) { /* Get samples of image {i}: */ double *pi = &(v[i]); double xsc = pi->c[0]/scale.c[0] + jx; double ysc = pi->c[1]/scale.c[1] + jy; val[i] = eval(i, iscale, xsc, ysc); } /* Compute mean and variance, and sampling point weight {wtxy}: */ double avg, var; tmnn_compute_avg_var(n, val, &avg, &var); double wtxy = wtx[jx+hwx]*wty[jy+hwy]; sum_wtxy_var += wtxy*var; sum_wtxy += wtxy; } } return (sum_wtxy == 0 ? 0.0 : sum_wtxy_var / sum_wtxy); } bool_t tmnn_check_ellipsoid ( int32_t n, double psol[], double F2sol, double sini[], double sopt[], double arad[], double astp[] ) { auto void compare_to_point(char *tag, double pref[]); /* Prints the RMS absolute difference {psol[i] - pref[i]}, and the RMS difference {psol[i] - pref[i]} relative to {arad[i]}, over {i} in {0..n-1}. Ignores coordinate {j} of {psol[i],pref[i]} iff {arad[i].c[j]} is zero. */ int32_t ind = 0; /* Indentation */ fprintf(stderr, "%*schecking the solution\n", ind, ""); bool_t valid = TRUE; if (sopt != NULL) { compare_to_point("optimum", sopt); } if (sini != NULL) { compare_to_point("initial", sini); } fprintf(stderr, "\n"); return valid; void compare_to_point(char *tag, double pref[]) { int32_t nv = 0; /* Number of nonzero coordinates in {arad[0..n-1]}. */ double sum_d2 = 0; /* Total abs squared disp between {psol} and {pref}. */ double sum_e2 = 0; /* Total rel squared disp between {psol} and {pref}. */ for (int32_t i = 0; i < n; i++) { double *q = &(psol[i]); double *o = &(pref[i]); double *r = &(arad[i]); double d, e; r2_sub(q, o, &d); for (int32_t j = 0; j < 2; j++) { if (r->c[j] == 0) { d.c[j] = 0; e.c[j] = 0; } else { e.c[j] = d.c[j]/r->c[j]; nv++; } } sum_d2 += r2_norm_sqr(&d); sum_e2 += r2_norm_sqr(&e); } double rms_d2 = (nv == 0 ? 0.0 : sqrt(sum_d2/nv)); double rms_e2 = (nv == 0 ? 0.0 : sqrt(sum_e2/nv)); /* Check whether displacement from initial point add to (0,0): */ fprintf(stderr, "%*s Abs and rel RMS displacement from %s = %10.6f %10.6f\n", ind, "", tag, rms_d2, rms_e2); } } void tmnn_compute_avg_var(int32_t nz, double z[], double *avgP, double *varP) { /* Compute the average {avg}: */ double sum_z = 0; /* Sum of {z[ki]} */ for (int32_t i = 0; i < nz; i++) { sum_z += z[i]; } double avg = sum_z/nz; /* Compute the sample variance: */ double sum_du2 = 0; for (int32_t i = 0; i < nz; i++) { double dui = z[i] - avg; sum_du2 += dui*dui; } double var = sum_du2/nz; (*avgP) = avg; (*varP) = var; } void minn_throw_arad(int32_t n, double zfrac, double rmin, double rmax, double rad[]) { for (int32_t j = 0; j < 2; j++) { /* Decide how many coordinates {rad[0..n-1].c[j]} will be set to zero. */ int32_t nzgoal = (int32_t)floor(zfrac.c[j]*n + 0.5); if (n >= 2) { /* Cook {nzgoal} to be {0} or {n} iff {zfrac} is 0.0 or 1.0, resp.: */ if ((nzgoal == 0) && (zfrac.c[j] > 0.0)) { nzgoal = 1; } if ((nzgoal == n) && (zfrac.c[j] < 1.0)) { nzgoal = n-1; } } assert(nzgoal <= n); int32_t nc = n; /* Number of coordinates remaining. */ int32_t nz = nzgoal; /* Number of coordinates still to be set to zero. */ for (int32_t i = 0; i < n; i++) { assert(nz <= nc); double rij; if ((nz == nc) || (nc*drandom() < nz)) { rij = 0.0; nz--; } else { rij = rmin + (rmax - rmin)*drandom(); } rad[i].c[j] = rij; nc--; } } } bool_t tmnn_check_ellipsoid ( int32_t n, double arad[], int32_t v, double C[], int32_t d, double U[], double urad[] ); /* Check correctness of the search basis {U} and the search ellipsoid radii {urad[0..d-1]}, assuming that the base ellipsoid radii are {arad[0..n-1]} and the (orthonormal) constraint basis is {C} (a {v} by {n} matrix). */ double tmnn_compute_test_goal_func ( int32_t n, /* Dimension of base space. */ double v[], /* A vector in base space. */ double vopt[] /* Optimum vector in base space. */ ); /* Computes a function of the vector {v} that has a quadratic minimum in the neighnborhood of vector {vopt}.*/ /* ---------------------------------------------------------------------- */ /* 0200_test_minn/Makefile */ # Test runs # Format: ${FUNC}:${BIAS}:${MONO}:${QUAD}:${NI}:${RADX}:${RADY}:${STPX}:${STPY}:${SCLX}:${SCLY} TESTS := \ imgmis/0/1/1/02/02/02/00.25/00.25/0/0 PLOT_TEST := \ imgmis/0/1/0/02/02/02/00.00/00.00/0/0 ALL_TESTS := \ optdst/0/1/0/02/04/04/01.00/01.00/0/0 \ optdst/0/1/1/02/04/04/01.00/01.00/0/0 \ indiff/1/1/1/02/04/04/01.00/01.00/0/0 \ imgmis/0/1/0/02/02/02/00.25/00.25/0/0 \ imgmis/0/1/0/02/20/20/01.00/01.00/4/4 TEST_ARGS := ${subst /, ,${TEST}} FUNC := ${word 1, ${TEST_ARGS}} BIAS := ${word 2, ${TEST_ARGS}} MONO := ${word 3, ${TEST_ARGS}} QUAD := ${word 4, ${TEST_ARGS}} NI := ${word 5, ${TEST_ARGS}} RADX := ${word 6, ${TEST_ARGS}} RADY := ${word 7, ${TEST_ARGS}} STPX := ${word 8, ${TEST_ARGS}} STPY := ${word 9, ${TEST_ARGS}} SCLX := ${word 10, ${TEST_ARGS}} SCLY := ${word 11, ${TEST_ARGS}} -${PROG} \ ${FUNC} ${BIAS} ${MONO} ${QUAD} ${NI} \ ${RADX} ${RADY} ${STPX} ${STPY} ${SCLX} ${SCLY}