/* Last edited on 2007-10-28 19:49:04 by stolfi */ /* INTERNAL PROTOTYPES */ /* IMPLEMENTATIONS */ void foo(void) { /* Check image channel counts: */ demand(nc > 0, "invalid channel count"); for (i = 0; i < n; i++) { /* Check channel count of {img[i]}: */ demand(nc <= img[i]->sz[0], "invalid channel count"); } for (k = 0; k < 2; k++) { for (i = 0; i < n; i++) { int npi = img[i]->sz[k+1]; /* Number of pixels along axis {k}. */ demand(p[i].c[k] - dmax >= 0 - fudge, "point too close to low edge"); demand(p[i].c[k] + ns.c[k]*step.c[k] <= npi + fudge, "point too close to high edge"); } double dmax = ns.c[k]*step.c[k]; /* Max displacement along axis {k}. */ demand(nr.c[k] >= 0, "window radius must be non-negative"); } /* Compute weights of window points: */ double w0[nr.c[0]]; /* {w1[r]} is weight for sample point {ħr} from window center. */ double w1[nr.c[1]]; /* {w1[r]} is weight for sample point {ħr} from window center. */ for (k = 0; k < 2; k++) { int r; for (r = 0; r <= nr.c[k]; r++) { (k == 0 ? w0 : w1)[r] = fial_sample_weight(r, nr.c[k]); } } } void fial_generate_all_alignments(int nims, int dmax, double step, int nal, alignment_t *al) { int sc, sr; int *dc = (int *)malloc(nims*sizeof(int)); int *dr = (int *)malloc(nims*sizeof(int)); int k; int carry; /* boolean */ int nfound = 0; k = nims; sc = 0; sr = 0; carry = 0; while(1) { /* Here {sc} is the sum of {dc[k..nims-1]}, ditto for {sr,dr}. */ /* Generate the smallest solution with current {dc,dr[k..nfound-1]}: */ int j; while (k > 0) { /* Move down and set digit to minimum value: */ k--; dc[k] = imax(-dmax, -(k*dmax + sc)); dr[k] = imax(-dmax, -(k*dmax + sr)); sc += dc[k]; sr += dr[k]; } /* Store this solution: */ if (nfound >= nal) { pnm_error("num of alignments underestimated"); } if ((sc != 0) || (sr != 0)) { pnm_error("displacements don't add to zero"); } { alignment_t a = alloc_alignment(nims); fprintf(stderr, "\n"); for (j = 0; j < nims; j++) { a.dcol[j] = step*dc[j]; a.drow[j] = step*dr[j]; a.mism[j] = (interval_t){0.00, 0.25}; fial_debug_displacement(" ", j, a.dcol[j], a.drow[j], a.mism[j], "\n"); } a.tvar = (interval_t){0.00, 0.25}; al[nfound] = a; nfound++; } /* Find first position {k} that can be incremented, and increment it: */ carry = 1; while (carry) { /* Now {sc,sr} is sum of {dc,dr[k..nims-1]} */ if (k >= nims) { /* Exhausted all possibilities, stop: */ if (nfound != nal) { pnm_error("number of alignments doesn't check"); } free(dc); free(dr); return; } sc -= dc[k]; sr -= dr[k]; if (dr[k] < imin(dmax, k*dmax - sr)) { /* Bump {dr[k]}, reset carry: */ dr[k]++; carry = 0; } else if (dc[k] < imin(dmax, k*dmax - sc)) { /* Bump {dc[k]}, reset {dr[k]}, reset carry: */ dc[k]++; dr[k] = imax(-dmax, -(k*dmax+sr)); carry = 0; } else { /* This position can't be incremented, carry to the next one: */ k++; } /* Now {sc,sr} is sum of {dc,dr[k+1-carry..nims-1]} */ } sc += dc[k]; sr += dr[k]; } } void fial_adjust_alignment ( int nims, image_rec_t *iml, double dmax, double epsilon, double radius, alignment_t *a ) { if (nims >= 2) { int arad = ceil(3*radius); /* Weights are negligible beyond this. */ int nwt = 2*arad+1; Image *avg = allocate_image(nwt, nwt, fial_max_chns(iml)); double *gwt = fial_gaussian_distr(nwt, radius); double s1 = ((double)(nims-1))/((double)nims); double corr = s1*s1; /* Variance correction factor. */ image_rec_t *impi; double step; int i; /* Loop until convergence: */ i = 0; impi = iml; step = 0.5*(dmax + epsilon); while(1) { int krowi, kcoli; interval_t prev_best; /* Adjust displacements so that they add to zero: */ fial_recenter_displacements(nims, al->dcol, al->drow); /* Compute the average of all images excluding image {i}: */ pnm_message("averaging images except %d...", i); fial_average_images(nims, iml, i, avg, al->dcol, al->drow); if (avg->cols < 20) { debug_image("avg", "", avg); } /* Compare image {i} against average, in all possible displacements: */ pnm_message("computing current mismatch of image %d...", i); al->mism[i] = fial_compute_mismatch(impi->im, al->dcol[i], al->drow[i], avg, corr, nwt, gwt); fial_debug_displacement("INIT", i, al->dcol[i], al->drow[i], al->mism[i], "\n"); prev_best = al->mism[i]; for (krowi = -1; krowi <= +1; krowi++) { double drowi = al->drow[i] + step*krowi; for (kcoli = -1; kcoli <= +1; kcoli++) { double dcoli = al->dcol[i] + step*kcoli; interval_t s = fial_compute_mismatch(impi->im, dcoli, drowi, avg, corr, nwt, gwt); if (fial_compare_mismatches(&s, &(al->mism[i])) < 0) { fial_debug_displacement("good", i, dcoli, drowi, s, "\n"); al->mism[i] = s; al->dcol[i] = dcoli; al->drow[i] = drowi; } } } if (fial_compare_mismatches(&(al->mism[i]), &prev_best) < 0) { fial_debug_displacement("BEST", i, al->dcol[i], al->drow[i], al->mism[i], "\n"); } impi = impi->next; i++; if (i >= nims) { i = 0; impi = iml; if (step <= 1.00001*epsilon) { fial_recompute_pixel_variance(nims, iml, al, avg, nwt, gwt); return; } step /= 2; } } } } void write_optimal_positions(int nims, image_rec_t *iml, double *dcol, double *drow) { int i; image_rec_t *impi; for (i = 0, impi = iml; i < nims; i++, impi = impi->next) { double tcol = dcol[i] + impi->ccol; double trow = drow[i] + impi->crow; printf("%03d %8.3f %8.3f\n", i, tcol, trow); } fflush(stdout); } void fial_recenter_displacements(int nims, double *dcol, double *drow) { double tcol = 0, trow = 0; int i; for (i = 0; i < nims; i++) { tcol += dcol[i]; trow += drow[i]; } tcol /= (double)nims; trow /= (double)nims; for (i = 0; i < nims; i++) { dcol[i] -= tcol; drow[i] -= trow; } } void fial_align_var ( int n, float_image_t *img[], /* Images to align. */ int nc, /* Number of channels to consider. */ i2_t ns, /* Number of adjustment steps along each axis. */ r2_t step, /* Precision required along each axis. */ r2_t p[], /* (IN/OUT) Corresponding points in each image. */ i2_t nr, /* Half-extent of comparison window along each axis. */ double *QP /* (OUT) Discrepancy measure for best alignment. */ ) { bool_t debug = TRUE; double fudge = 1.0e-7; /* Fudge term to compensate roundoff errors. */ int i, k; /* Check image channel counts: */ demand(nc > 0, "invalid channel count"); for (i = 0; i < n; i++) { /* Check channel count of {img[i]}: */ demand(nc <= img[i]->sz[0], "invalid channel count"); } /* Check search grid parameters: */ for (k = 0; k < 2; k++) { demand(ns.c[k] >= 0, "search step count must be non-negative"); if (ns.c[k] > 0) { demand(step.c[k] > 0, "search {step} must be positive"); } double dmax = ns.c[k]*step.c[k]; /* Max displacement along axis {k}. */ for (i = 0; i < n; i++) { int npi = img[i]->sz[k+1]; /* Number of pixels along axis {k}. */ demand(p[i].c[k] - dmax >= 0 - fudge, "point too close to low edge"); demand(p[i].c[k] + ns.c[k]*step.c[k] <= npi + fudge, "point too close to high edge"); } demand(nr.c[k] >= 0, "window radius must be non-negative"); } /* Compute weights of window points: */ double w0[nr.c[0]]; /* {w1[r]} is weight for sample point {ħr} from window center. */ double w1[nr.c[1]]; /* {w1[r]} is weight for sample point {ħr} from window center. */ for (k = 0; k < 2; k++) { int r; for (r = 0; r <= nr.c[k]; r++) { (k == 0 ? w0 : w1)[r] = fial_sample_weight(r, nr.c[k]); } } i2_t d[n]; /* Adjustment step counts to apply to {p[0..n-1]}. */ i2_t dbest[n]; /* Best step count tuple found so far. */ double Qbest = INF; /* Mismatch of best alignment. */ r2_t q[n]; /* Work: Displaced points {q[i] = p[i] + d[i]}. */ /* Initialize the adjustment step counts {d[0..n-1]}: */ for (i = 0; i < n; i++) { d[i] = (i2_t){{0,0}}; } /* Enumeration loop: */ while (1) { /* Evaluate {Q} for the adjustment step count tuple {d[0..n-1]}: */ fial_displace_points(n, p, d, q, step); double Q = fial_mismatch_var(n, img, q, nc, nr, w0, w1); if (debug) { fial_debug_displacement("trial", n, d, q, Q, "\n"); } if (Q < Qbest) { Qbest = Q; for (i = 0; i < n; i++) { dbest[i] = d[i]; } } /* Advance to the next tuple: */ if (fial_next_tuple(n, d, ns)) { /* Enumeration finished, return the best solution: */ fial_displace_points(n, p, dbest, p, step); (*QP) = Qbest; return; } } } double fial_sample_weight(int r, int rmax); /* Returns the relative weight of a sample that is displaced by distance {r} from the center sample, assuming that the maximum {r} within the comparison window is {rmax}. */ alignment_t *fial_alloc_alignment_set(int ncan); void fial_generate_all_alignments(int nims, int dmax, double step, int nal, alignment_t *al); /* Generate all relative alignments for {nims} images whose entries range in {[-dmax .. +dmax]} times {step}. Expects {nal} to be the exact number of such alignments. Returns the alignments in {al[0..*nal-1]}. */ image_rec_t *fial_reduce_images(image_rec_t *iml); /* Reduces all imagaes in {iml} by a linear factor of 2; returns a list of the results. Also computes the coordinates of the centers of the reduced images. */ int fial_num_rel_shifts(int nims, int dmax, int sum); /* Computes the number of integer vectors in {[-dmax..+dmax]^nims} which add to {sum}. */ void fial_recenter_displacements(int nims, double *dcol, double *drow); /* Shifts all displacements so that their sum is zero. */ void fial_adjust_alignment ( int nims, image_rec_t *iml, double dmax, double epsilon, double radius, alignment_t *a ); /* Adjusts the shifts {(al->dcol[i],al->drow[i])} for each image {i}, so as to minimize the average pixel variance {al->tvar}. The maximum adjustment is {dmax} and the nominal precision required is {epsilon}, in each direction. The pixel variances are weighted with a Gaussian of characteristic radius {radius}. A displacement of zero means to align the center of the image at the origin. The displacements are adjusted so as to add to zero. */ void fial_average_images(int nims, image_rec_t *iml, int iex, Image *avg, double *dcol, double *drow); /* Computes average of all images, excluding image {iex}, shifting each one so that its center displaced by {dcol[j],drow[j]} matches the center of {avg}. */ interval_t fial_compute_mismatch ( Image *tim, double dcol, double drow, Image *sim, double scale, int nwt, double *gwt ); /* Computes the mismatch between the image {tim}, shifted by {(dcol,drow)}, against the square image {sim}. The vector {wt} is a unidimensional weight distribution. Assumes that the length of {wt} is the same as the number of rows and columns of {sim}. The mismatch is the sum of the squared pixel differences, each multiplied by the {scale} factor and the bidimensional weight distribution {wt[col]*wt[row]}. */ void fial_recompute_pixel_variance(int nims, image_rec_t *iml, alignment_t *al, Image *avg, int nwt, double *wt); /* Recomputes the mismatches {al->mism[i]} and the total pixel variance {al->tvar}, from scratch. */ int fial_max_chns(image_rec_t *iml); /* Returns the maximum {chns} among all images in {iml}. */ void fial_sort_and_prune_alignments(int nims, int *ncanP, alignment_t *can, double epsilon, int max); double *fial_gaussian_distr(int width, double sigma); /* Gaussian distribution {g[i] = A*exp(-((i+0.5)-m)^2/sigma^2)} where {m = width/2} and {A} is such that sum is 1. */ void fial_debug_displacement(char *label, int i, double col, double row, interval_t s, char *tail); double fial_compare_cands(int nims, alignment_t *al, alignment_t *bl); /* Returns the maximum difference between the shifts of {al} and {bl}. */ int fial_compare_mismatches(interval_t *a, interval_t *b); /* Returns +1 if {a} looks smaller than {b}, -1 if {a} looks larger than {b}, and 0 if they look equal or incomparable. */ /* IMPLEMENTATIONS */ void fial_optimal_shifts ( int nims, float_image_t *iml, double dxmax, double dymax, double epsilon. double rx, double ry, double *dx, double *dy ); void fial_multiscale_optimal_shifts ( int nims, float_image_t *iml, double dxmax, double dymax, double epsilon, double rx, double ry, int mult, alignment_t *opt, int *noptP ) { double guess_step = 0.5; /* Initial guesses are spaced every half-pixel: */ int idmax = (int)(floor(dmax/guess_step)); double allshifts = (double)fial_num_rel_shifts(nims, idmax, 0); /* Watch oflow... */ double allsols = allshifts*allshifts; /* Total possible alignments. */ alignment_t *can = NULL; int ncan; int j, k; fprintf(stderr, "enter compute_multiscale_optimal_shifts\n"); fprintf(stderr, " allsols = %g\n", allsols); if (mult >= allsols) { /* The caller asked for all possible solutions. */ int nall = (int)(floor(allsols+0.5)); fprintf(stderr, " computing all %d solutions...\n", nall); can = fial_alloc_alignment_set(nall); fial_generate_all_alignments(nims, idmax, guess_step, nall, can); ncan = nall; } else { /* Generate {4*mult} solutions for coarser scale, then refine and select. */ image_rec_t *rdl = fial_reduce_images(iml); /* Compute 4*mult solutions for coarse images */ can = fial_alloc_alignment_set(4*mult); fprintf(stderr, " recursion to smaller scale...\n"); compute_multiscale_optimal_shifts(nims, rdl, 4*mult, dmax/2, guess_step, radius/2, can, &ncan); for (k = 0; k < ncan; k++) { /* Map displacements to current scale: */ alignment_t *al = &(can[k]); image_rec_t *rdli, *imli; int i; for (i=0, imli=iml, rdli=rdl; inext, rdli=rdli->next) { al->dcol[i] = 2*al->dcol[i] + (rdli->ccol - imli->ccol); al->drow[i] = 2*al->drow[i] + (rdli->crow - imli->crow); al->mism[i] = (interval_t){0.00, 0.25}; } } } /* Refine each sol by moving at most {epsilon} from current pos: */ fprintf(stderr, " local refinement...\n"); for (k = 0; k < ncan; k++) { fial_adjust_alignment(nims, iml, guess_step, epsilon, radius, &(can[k])); } /* Select {mult} best solutions: */ fprintf(stderr, " sorting and pruning...\n"); fial_sort_and_prune_alignments(nims, &ncan, can, epsilon, mult); (*noptP) = (ncan < mult ? ncan : mult); for (j = 0; j < (*noptP); j++) { opt[j] = can[j]; } for (j = (*noptP); j < ncan; j++) { opt[j].dcol = NULL; opt[j].drow = NULL; free(can[j].dcol); free(can[j].drow); free(can[j].mism); } free(can); fprintf(stderr, "exit compute_multiscale_optimal_shifts\n"); } int fial_num_rel_shifts(int nims, int dmax, int sum) { fprintf(stderr, "(%d)", nims); if (abs(sum) > nims*dmax) { return 0; } else if (nims <= 1) { return 1; } else if (nims == 2) { return 2*dmax + 1; } else if (nims == 3) { return 3*dmax*(dmax + 1) + 1; } else if (nims == 4) { return 8*dmax*(dmax*(2*dmax + 3) + 1)/3 + 1; } else if (nims == 5) { /* There must be a better way... */ int i, s = 0; for (i=-dmax; i <= dmax; i++) { s += fial_num_rel_shifts(nims-1,dmax,sum-i); } return s; } else { pnm_error("too many images"); return 0; } } image_rec_t *fial_reduce_images(image_rec_t *iml) { fprintf(stderr, "r"); if (iml == NULL) { return NULL; } else { Image *im = iml->im; double rcol = iml->ccol; double rrow = iml->crow; Image *rd = reduce_image(im, &rcol, &rrow); return cons_image(rd, rcol, rrow, fial_reduce_images(iml->next)); } } void fial_sort_and_prune_alignments(int nims, int *ncanP, alignment_t *can, double epsilon, int max) { /* !!! must eliminate duplicates !!! */ int ncan = *ncanP; int nok, i; /* Select the best {max} alignments: */ nok = 1; for (i = 1; i < ncan; i++) { int j; alignment_t t = can[i]; j = nok-1; while ((j >= 0) && (fial_compare_cands(nims, &(can[j]), &t) > epsilon/2)) { j--; } if (j < 0) { j = nok; while ((j > 0) && (fial_compare_mismatches(&(can[j-1].tvar), &(t.tvar)) > 0)) { if (j < max) { can[j] = can[j-1]; } j--; } can[j] = t; if (nok < max) { nok++; } } } (*ncanP) = nok; } double fial_compare_cands(int nims, alignment_t *al, alignment_t *bl) { int i; double dmax = 0; for (i = 0; i < nims; i++) { double dc = fabs(al->dcol[i] - bl->dcol[i]); double dr = fabs(al->drow[i] - bl->drow[i]); if (dc > dmax) { dmax = dc; } if (dr > dmax) { dmax = dr; } } return dmax; } int fial_compare_mismatches(interval_t *a, interval_t *b) { double amin = 0.75*a->lo + 0.25*a->hi; double bmin = 0.75*b->lo + 0.25*b->hi; return cmp_double(&amin, &bmin); } alignment_t *fial_alloc_alignment_set(int ncan) { int i; alignment_t *als = (alignment_t *)malloc(ncan*sizeof(alignment_t)); if (als == NULL) { pnm_error("out of memory for alignments"); } for (i = 0; i < ncan; i++) { alignment_t *alsi = &(als[i]); alsi->dcol = NULL; alsi->drow = NULL; alsi->mism = NULL; alsi->tvar = (interval_t){0.00, 0.25}; } return als; } alignment_t alloc_alignment(int nims) { alignment_t al; int i; al.dcol = (double *)malloc(nims*sizeof(double)); al.drow = (double *)malloc(nims*sizeof(double)); al.mism = (interval_t *)malloc(nims*sizeof(interval_t)); if ((al.dcol == NULL) || (al.drow == NULL) || (al.mism == NULL)) { pnm_error("out of memory for displacement vectors"); } for (i=0; i < nims; i++) { al.dcol[i] = 0; al.drow[i] = 0; al.mism[i] = (interval_t){0.00, 0.25}; } al.tvar = (interval_t){0.00, 0.25}; return al; } int imin(int x, int y) { return (x < y ? x : y); } int imax(int x, int y) { return (x > y ? x : y); } void fial_generate_all_alignments(int nims, int dmax, double step, int nal, alignment_t *al) { int sc, sr; int *dc = (int *)malloc(nims*sizeof(int)); int *dr = (int *)malloc(nims*sizeof(int)); int k; int carry; /* boolean */ int nfound = 0; k = nims; sc = 0; sr = 0; carry = 0; while(1) { /* Here {sc} is the sum of {dc[k..nims-1]}, ditto for {sr,dr}. */ /* Generate the smallest solution with current {dc,dr[k..nfound-1]}: */ int j; while (k > 0) { /* Move down and set digit to minimum value: */ k--; dc[k] = imax(-dmax, -(k*dmax + sc)); dr[k] = imax(-dmax, -(k*dmax + sr)); sc += dc[k]; sr += dr[k]; } /* Store this solution: */ if (nfound >= nal) { pnm_error("num of alignments underestimated"); } if ((sc != 0) || (sr != 0)) { pnm_error("displacements don't add to zero"); } { alignment_t a = alloc_alignment(nims); fprintf(stderr, "\n"); for (j = 0; j < nims; j++) { a.dcol[j] = step*dc[j]; a.drow[j] = step*dr[j]; a.mism[j] = (interval_t){0.00, 0.25}; fial_debug_displacement(" ", j, a.dcol[j], a.drow[j], a.mism[j], "\n"); } a.tvar = (interval_t){0.00, 0.25}; al[nfound] = a; nfound++; } /* Find first position {k} that can be incremented, and increment it: */ carry = 1; while (carry) { /* Now {sc,sr} is sum of {dc,dr[k..nims-1]} */ if (k >= nims) { /* Exhausted all possibilities, stop: */ if (nfound != nal) { pnm_error("number of alignments doesn't check"); } free(dc); free(dr); return; } sc -= dc[k]; sr -= dr[k]; if (dr[k] < imin(dmax, k*dmax - sr)) { /* Bump {dr[k]}, reset carry: */ dr[k]++; carry = 0; } else if (dc[k] < imin(dmax, k*dmax - sc)) { /* Bump {dc[k]}, reset {dr[k]}, reset carry: */ dc[k]++; dr[k] = imax(-dmax, -(k*dmax+sr)); carry = 0; } else { /* This position can't be incremented, carry to the next one: */ k++; } /* Now {sc,sr} is sum of {dc,dr[k+1-carry..nims-1]} */ } sc += dc[k]; sr += dr[k]; } } void fial_adjust_alignment ( int nims, image_rec_t *iml, double dmax, double epsilon, double radius, alignment_t *a ) { if (nims >= 2) { int arad = ceil(3*radius); /* Weights are negligible beyond this. */ int nwt = 2*arad+1; Image *avg = allocate_image(nwt, nwt, fial_max_chns(iml)); double *gwt = fial_gaussian_distr(nwt, radius); double s1 = ((double)(nims-1))/((double)nims); double corr = s1*s1; /* Variance correction factor. */ image_rec_t *impi; double step; int i; /* Loop until convergence: */ i = 0; impi = iml; step = 0.5*(dmax + epsilon); while(1) { int krowi, kcoli; interval_t prev_best; /* Adjust displacements so that they add to zero: */ fial_recenter_displacements(nims, al->dcol, al->drow); /* Compute the average of all images excluding image {i}: */ pnm_message("averaging images except %d...", i); fial_average_images(nims, iml, i, avg, al->dcol, al->drow); if (avg->cols < 20) { debug_image("avg", "", avg); } /* Compare image {i} against average, in all possible displacements: */ pnm_message("computing current mismatch of image %d...", i); al->mism[i] = fial_compute_mismatch(impi->im, al->dcol[i], al->drow[i], avg, corr, nwt, gwt); fial_debug_displacement("INIT", i, al->dcol[i], al->drow[i], al->mism[i], "\n"); prev_best = al->mism[i]; for (krowi = -1; krowi <= +1; krowi++) { double drowi = al->drow[i] + step*krowi; for (kcoli = -1; kcoli <= +1; kcoli++) { double dcoli = al->dcol[i] + step*kcoli; interval_t s = fial_compute_mismatch(impi->im, dcoli, drowi, avg, corr, nwt, gwt); if (fial_compare_mismatches(&s, &(al->mism[i])) < 0) { fial_debug_displacement("good", i, dcoli, drowi, s, "\n"); al->mism[i] = s; al->dcol[i] = dcoli; al->drow[i] = drowi; } } } if (fial_compare_mismatches(&(al->mism[i]), &prev_best) < 0) { fial_debug_displacement("BEST", i, al->dcol[i], al->drow[i], al->mism[i], "\n"); } impi = impi->next; i++; if (i >= nims) { i = 0; impi = iml; if (step <= 1.00001*epsilon) { fial_recompute_pixel_variance(nims, iml, al, avg, nwt, gwt); return; } step /= 2; } } } } void fial_recompute_pixel_variance(int nims, image_rec_t *iml, alignment_t *al, Image *avg, int nwt, double *wt) { int i = 0; image_rec_t *impi = iml; double qt = 1.0/((double)nims); al->tvar = (interval_t){0, 0}; fial_average_images(nims, iml, nims, avg, al->dcol, al->drow); for (i = 0, impi=iml; i < nims; i++, impi = impi->next) { interval_t s = fial_compute_mismatch(impi->im, al->dcol[i], al->drow[i], avg, 1.0, nwt, wt); al->mism[i] = s; al->tvar.lo += qt*s.lo; al->tvar.hi += qt*s.hi; } } void fial_recenter_displacements(int nims, double *dcol, double *drow) { double tcol = 0, trow = 0; int i; for (i = 0; i < nims; i++) { tcol += dcol[i]; trow += drow[i]; } tcol /= (double)nims; trow /= (double)nims; for (i = 0; i < nims; i++) { dcol[i] -= tcol; drow[i] -= trow; } } double *fial_gaussian_distr(int width, double sigma) { double m = ((double)width)/2.0; double *g = (double *)malloc(width*sizeof(double)); double s = 0; int i; if (g == NULL) { pnm_error("out of memory for gaussian weight"); } for (i=0; i < width; i++) { double x = (i + 0.5 - m)/sigma; double y = exp(-x*x); g[i] = y; s += y; } for (i=0; i < width; i++) { g[i] /= s; } return g; } void write_optimal_positions(int nims, image_rec_t *iml, double *dcol, double *drow) { int i; image_rec_t *impi; for (i = 0, impi = iml; i < nims; i++, impi = impi->next) { double tcol = dcol[i] + impi->ccol; double trow = drow[i] + impi->crow; printf("%03d %8.3f %8.3f\n", i, tcol, trow); } fflush(stdout); } interval_t fial_compute_mismatch ( Image *tim, double dcol, double drow, Image *sim, double scale, int nwt, double *gwt ) { int srow, scol; double wt2; interval_t d2; interval_t td2 = (interval_t){0.0, 0.0}; /* Total mismatch. */ if ((sim->rows != nwt) || (sim->cols != nwt)) { pnm_error("ref image has wrong size"); } for (srow = 0; srow < sim->rows; srow++) { for (scol = 0; scol < sim->cols; scol++) { int sindex = (srow * sim->cols + scol)*sim->chns; interval_t *sv = &(sim->el[sindex]); interval_t tv[MAX_CHANNELS]; double tcol = scol + ((double)(tim->cols - sim->cols))/2.0 - dcol; double trow = srow + ((double)(tim->rows - sim->rows))/2.0 - drow; /* debug = ((scol % 56) == 23) & ((srow % 56) == 23); */ eval_at_point(tim, tcol, trow, &(tv[0])); if (debug) { debug_floatized_pixel("imi", tcol, trow, sim->chns, &(tv[0]), "\n"); } if (tim->chns != sim->chns) { convert_pixel(tim->chns, &(tv[0]), sim->chns); } subtract_pixel(sim->chns, &(sv[0]), &(tv[0])); if (debug) { debug_floatized_pixel("dif", scol, srow, sim->chns, &(tv[0]), "\n\n"); } d2 = pixel_norm(sim->chns, &(tv[0])); wt2 = scale*wt[srow]*wt[scol]; td2.lo += wt2*d2.lo; td2.hi += wt2*d2.hi; } } return td2; } void fial_average_images(int nims, image_rec_t *iml, int iex, Image *avg, double *dcol, double *drow) { int row, col; int navg = ((iex >= 0) && (iex < nims) ? nims-1 : nims); double qt = 1.0/((double)navg); double tqt; for (row = 0; row < avg->rows; row++) { for (col = 0; col < avg->cols; col++) { int aindex = (row * avg->cols + col)*avg->chns; interval_t *fva = &(avg->el[aindex]); int j, k; image_rec_t *impj; for (k = 0; k < avg->chns; k++) { fva[k] = (interval_t){ 0.0, 0.0 }; } for (j = 0, impj = iml; j < nims; j++, impj = impj->next) { if (j != iex) { Image *imj = impj->im; interval_t fvj[MAX_CHANNELS]; double colj = col + ((double)(imj->cols - avg->cols))/2.0 - dcol[j]; double rowj = row + ((double)(imj->rows - avg->rows))/2.0 - drow[j]; eval_at_point(imj, colj, rowj, &(fvj[0])); if (imj->chns != avg->chns) { convert_pixel(imj->chns, &(fvj[0]), avg->chns); } accum_pixel(avg->chns, &(fvj[0]), qt, fva, &tqt); } } } } } int fial_max_chns(image_rec_t *iml) { int chns = 0; while (iml != NULL) { Image *im = iml->im; if (im->chns > chns) { chns = im->chns; } iml = iml->next; } return chns; } double fial_sample_weight(int r, int rmax) { if ((rmax == 0) || (r == 0)) { return 1.0; } else { /* Hann window: */ double z = r/(1.0 + rmax); return 0.5*(1 + cos(M_PI*z)); } }