/** @file find_register.c * Registration of a pair of overlapping images. * Given a pair of images that are known to overlap, we want to * register them exactly, so that within the overlapping area each * pixel in one image represents the same physical location as a pixel * in the other. We handle only the simplest case, where the two * images already have identical spatial scales, are undistorted * relative to some common coordinate plane, and are not rotated * relative to each other. There's even some a priori information * about the the relative displacement of one image relative to the * other, so we simply search all possible relative positions in the * neighborhood of this initial suggestion and accept the one that * yields the best value of some goodness-of-fit function. * * @author University of Arizona Digital Image Analysis Lab * @date 2003 * @version $Id: find_register.c,v 1.4 2004/09/07 23:35:08 mmunro Exp $ */ /* This file is part of TREES. TREES is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. TREES is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with TREES; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. */ #if HAVE_CONFIG_H #include #endif /* HAVE_CONFIG_H */ #include #include #include "proto.h" #include #include #if WITH_DMALLOC #include #endif /* WITH_DMALLOC */ /** * Represents a range of indicies, e.g., within a particular image dimension * (bands, lines, pixels, etc.). */ typedef struct dim_range { uint32_t off; /**< offset, the lower bound of the range. */ uint32_t len; /**< the length of the range. */ } Dim_range_t; typedef struct dim_range *dim_rangep_t; /** * Sums of pixel values within a rectangular region of interest (ROI) * (limited to a single band and interval). */ typedef struct bandroi { const IMAGE *img; /**< the image within which the ROI is located. */ uint32_t interval; /**< interval index, currently unused. */ uint32_t band; /**< band index (the ROI is within this band). */ Dim_range_t lin; /**< the range of lines covered by the ROI. */ Dim_range_t pix; /**< the range of pixels covered by the ROI. */ long double sum; /**< sum of pixel values within the ROI. */ long double sumsq; /**< sum of squared pixel values within the ROI. */ unsigned long int n; /**< total number of pixels == lin.len * pix.len. */ } Bandroi_t; typedef struct bandroi *bandroip_t; /** * Represents the overlap interval between two images. * The interval is on a single image dimension (pixels, lines, etc.); * an offset from the image origin and a length (less than the * image size in that dimension) are the basis of this representation. */ typedef struct overlap { uint32_t soff; /**< coordinate offset in the first image. */ uint32_t toff; /**< coordinate offset in the second image. */ uint32_t len; /**< length of the interval in both images. */ } Overlap_t; typedef struct overlap *overlapp_t; /** * Forward or backward displacement both horizontal and vertical. */ typedef struct xydisplace { int dx; /**< positive or negative horizontal displacement. */ int dy; /**< positive or negative vertical displacement. */ } Xydisplace_t; typedef struct xydisplace *xydisplacep_t; /** * The goodness of fit metric value from matching images, with their positions. */ typedef struct mpos { double t; /**< goodness of fit value at a particular position. */ Xydisplace_t pos; /**< the relative position. */ } Mpos_t; typedef struct mpos *mposp_t; /** * A fixed-size ordered queue (a heap structure) of goodness of fit values. */ typedef struct mheap { size_t n_max; /**< maximum number of entries in the queue. */ size_t n; /**< current number of entries in the queue. */ mposp_t q; /**< array representation of a complete binary tree (heap). */ } Mheap_t; typedef struct mheap * mheapp_t; /** * Registration parameters for a particular downsampled image resolution. */ typedef struct regparam { uint32_t x_radius; /**< horizontal search radius. */ uint32_t y_radius; /**< vertical search radius. */ } Regparam_t; typedef struct regparam * regparamp_t; /** * A collection of settings for registering images at various resolutions. * The parameters in @c plevel[0] are for the highest resolution, * those in @c plevel[nlevels-1] for the lowest. */ typedef struct regparam_pyramid { const size_t nlevels; /**< the total number of resolution levels. */ const uint32_t ratio; /**< resolution ratio between levels in the pyramid. */ const regparamp_t plevel; /**< array of level search parameters. */ } Regparam_pyramid_t; typedef struct regparam_pyramid * regparam_pyramidp_t; /** * A pyramid of images, ordered by resolution. * The image in @c ilevel[0] is at the highest resolution, those in * @c ilevel[1] to @c ilevel[nlevels-1] are downsampled to progressively * lower resolutions. */ typedef struct regimage_pyramid { size_t nlevels; /**< the total number of resolution levels. */ uint32_t ratio; /**< resolution ratio between levels in the pyramid. */ IMAGE * * ilevel; /**< an array of pointers to the SADIE images. */ } Regimage_pyramid_t; typedef struct regimage_pyramid * regimage_pyramidp_t; /* /\** */ /* * Default horizontal and vertical search radii for various resolutions. */ /* *\/ */ /* static Regparam_t default_radii[] = */ /* { */ /* { 1, 1 }, */ /* { 1, 1 }, */ /* { 2, 2 }, */ /* { 11, 11 } */ /* }; */ /* /\** */ /* * Default registration settings (for a modest four-level image pyramid). */ /* *\/ */ /* static Regparam_pyramid_t default_param = */ /* { */ /* 4, 2, default_radii */ /* }; */ /** * Default horizontal and vertical search radii for various resolutions. */ static Regparam_t default_radii[] = { { 1, 1 }, { 1, 1 }, { 2, 2 }, { 11, 11 } }; /** * Default registration settings (for a modest four-level image pyramid). */ static Regparam_pyramid_t default_param = { 4, 2, default_radii }; /** * Marker value guaranteed to be less than any goodness of fit value. */ static const double MARKER_T = -1e36; /** * Arbitrary lower bound for goodness of fit values > @a MARKER_T. */ static const double MIN_T = -1e35; /** * Error codes for registration failures. */ typedef enum { REGERR_OK = 0, /**< no error (zero, to make @c if @c (!err) tests work). */ REGERR_EMPTY, /**< unexpected null data structure. */ REGERR_MISMATCH, /**< image pyramid depth or resolution mismatch. */ REGERR_BADPARAM, /**< corrupt search parameters. */ REGERR_NOLAP, /**< no overlap between the images. */ REGERR_TOOSMALL /**< the overlap area is too small. */ } regerr_t; /** * Make a fixed-size ordered queue to hold goodness of fit values and * image positions. Use an array representation of a complete binary * tree for the queue: @c q[i] represents a node of the tree having * child nodes at @c q[2*i] and @c q[2*i+1], and a parent node at * @c q[i/2]; all nodes @c q[k] for @c k>n don't exist. Maintain the * heap property for the goodness of fit values, @a t: each node has * @a t >= the @a t value of its parent. There's a constant dummy * node at @c q[0] to avoid having to treat the root node as a special * case. I've not managed to dig the original reference out of Knuth * vol. 3, so this is based on Chapter 11 of Sedgewick, R. * @e Algorithms Addison-Wesley: Reading MA: 1983. There's also * something at * http://www.nist.gov/dads/HTML/heap.html * @param size the maximum number of entries the queue can hold. * @return a pointer to the newly-created queue (initially empty). * @see dispose_mheap() @see fix_mheap() @see min_mheap() @see add_mheap() */ static mheapp_t make_mheap (size_t size) { mheapp_t themheap; themheap = malloc (sizeof (Mheap_t)); if (themheap) { themheap->n_max = size; themheap->n = 0; themheap->q = calloc (size + 1, sizeof (Mpos_t)); if (themheap->q) { themheap->q[0].t = MIN_T; themheap->q[0].pos.dx = 0; themheap->q[0].pos.dy = 0; } else { free (themheap); themheap = NULL; } } return themheap; } /** * Clean up a fixed-size ordered queue. * @param themheap a pointer to the queue data structure to destroy. * @see make_mheap() */ static void dispose_mheap (mheapp_t themheap) { if (themheap) { free (themheap->q); free (themheap); } } /** * Restore the heap property to a heap with a modified leaf node. * @param themheap a pointer to the heap to be fixed. * @param i index of the node that may have invalidated the heap property. * @see make_heap() */ static void leaf_fix_mheap (mheapp_t themheap, size_t i) { Mpos_t p; if ((themheap == NULL) || (themheap->q == NULL) || (i > themheap->n_max)) { fprintf (stderr, "leaf_fix_mheap: invalid data structure or index.\n"); exit (-2); } p = themheap->q[i]; /* suspect leaf node */ while (themheap->q[i / 2].t > p.t) { themheap->q[i] = themheap->q[i / 2]; i = i / 2; } themheap->q[i] = p; } /** * Restore the heap property to a heap with a modified root node. * @param themheap a pointer to the heap to be fixed. * @see make_heap() */ static void root_fix_mheap (mheapp_t themheap) { Mpos_t p; size_t i; size_t j; if ((themheap == NULL) || (themheap->q == NULL)) { fprintf (stderr, "root_fix_mheap: invalid data structure or index.\n"); exit (-2); } i = 1; p = themheap->q[i]; /* suspect root node */ while (i <= themheap->n / 2) { j = i * 2; /* left child node */ if ((j < themheap->n) && (themheap->q[j].t > themheap->q[j + 1].t)) j++; /* right child node, if it's smaller */ if (p.t <= themheap->q[j].t) break; /* heap property satisfied */ themheap->q[i] = themheap->q[j]; /* move smaller child towards root */ i = j; } themheap->q[i] = p; } /** * The minimum t-value currently in the queue. * If the queue exists & isn't empty, the root node holds the minimum t-value. * @param themheap pointer to the queue. * @return the double minimum t-value, or a value lower than any possible t. */ static inline double min_mheap (const mheapp_t themheap) { return (themheap && themheap->n) ? themheap->q[1].t : MIN_T; } /** * Whether or not the queue has reached its maximum (fixed) capacity. * @param thedheap pointer to the queue. * @return non-zero iff the queue is valid and is full to capacity. */ static inline int full_mheap (const mheapp_t themheap) { return (themheap && themheap->n_max && (themheap->n >= themheap->n_max)); } /** * Add a candidate node to a heap data structure. * @param themheap the heap to modify. * @param cand the candidate node to add. */ static void add_mheap (mheapp_t themheap, const mposp_t cand) { size_t i; if ((themheap == NULL) ||(themheap->q == NULL) || (cand == NULL)) { fprintf (stderr, "add_mheap: invalid data structure or index.\n"); exit (-2); } if ((min_mheap (themheap) >= cand->t) && full_mheap (themheap)) return; if (themheap->n_max == 1) { /* frequently encountered degenerate case */ themheap->n = 1; themheap->q[1] = *cand; return; } /* Do a crude linear search for the candidate coordinates. */ i = themheap->n; while ((i != 0) && ((cand->pos.dx != themheap->q[i].pos.dx) || (cand->pos.dy != themheap->q[i].pos.dy))) i--; if (i == 0) { /* The candidate node position isn't yet in the heap. */ if (!full_mheap (themheap)) { themheap->n++; themheap->q[themheap->n] = *cand; leaf_fix_mheap (themheap, themheap->n); } else { themheap->q[1] = *cand; root_fix_mheap (themheap); } } else { /* Found the same displacement already within the heap. */ if (cand->t > themheap->q[i].t) { /* Candidate is better: update the old node. */ themheap->q[i].t = cand->t; leaf_fix_mheap (themheap, i); } } } /** * Create a pyramid of progressively downsampled images. * Assumes the images at different levels in the pyramid are related to * each other by a fixed downsampling ratio. * @param full_img a pointer to the full-resolution (level 0) SADIE image. * @param levs the number of levels in the pyramid. * @param ratio the ratio between resolutions in adjacent pyramid levels. * @return the newly-created pyramid, or NULL on failure. */ regimage_pyramidp_t make_regpyramid (const IMAGE * full_img, size_t levs, uint32_t ratio) { regimage_pyramidp_t thepyramid; int i; thepyramid = (regimage_pyramidp_t) malloc (sizeof (Regimage_pyramid_t)); if (thepyramid) { thepyramid->nlevels = levs; thepyramid->ratio = ratio; thepyramid->ilevel = (IMAGE **) calloc (levs, sizeof (IMAGE *)); if (thepyramid->ilevel) { thepyramid->ilevel[0] = (IMAGE *)full_img; i = 0; while (thepyramid->ilevel[i] && (i < (levs - 1))) { i++; RESAMPL (thepyramid->ilevel[i - 1], ratio, ratio, ratio, ratio, &thepyramid->ilevel[i]); } if (thepyramid->ilevel[i] == NULL) { /* Clean up after an allocation failure. */ while (--i > 0) RELMEM (thepyramid->ilevel[i]); /* But be sure to leave the full image, level[0], intact. */ free (thepyramid->ilevel); free (thepyramid); thepyramid = NULL; } } } return thepyramid; } /** * Clean up an image pyramid. * @param thepyramid a pointer to the pyramid structure to destroy. */ static void dispose_regpyramid (regimage_pyramidp_t thepyramid) { int i; if (thepyramid) { if (thepyramid->ilevel) { for (i = (thepyramid->nlevels - 1); i > 0; i--) RELMEM (thepyramid->ilevel[i]); /* But be sure to leave the full image, level[0], intact. */ free (thepyramid->ilevel); } free (thepyramid); } } /** * Initialize a single-band rectangular region of interest from the image data. * Offsets from the image origin to one corner of the rectangle and the * lengths of its sides define the region of interest, which is clipped * if it extends beyond the image boundaries. This computes the total number * of pixels within the rectangle, the sum of the pixel values, and the sum * of the squared values. * @param theroi the single-band region of interest to initialize. * @param theimg the SADIE image to which the ROI applies. * @param ival interval within which the band lies (presently unused). * @param band the ROI lies within this band. * @param loff offset from the image origin to the first line of the ROI. * @param nl the requested number of lines in the region. * @param poff offset to the first ROI pixel in each line. * @param np the requested number of ROI pixels in each line. */ static void initialize_bandroi (bandroip_t theroi, const IMAGE * theimg, uint32_t ival, uint32_t band, uint32_t loff, uint32_t nl, uint32_t poff, uint32_t np) { int i, j; long double px; if (theroi == NULL) return; theroi->interval = ival; theroi->band = band; theroi->lin.off = loff; theroi->pix.off = poff; theroi->lin.len = theroi->pix.len = 0; theroi->sum = theroi->sumsq = 0; if ((theimg != NULL) && (band < theimg->nbnd) && (loff < theimg->nlin) && (poff < theimg->npix)) { theroi->lin.len = ((loff + nl) <= theimg->nlin) ? nl : (theimg->nlin - loff); theroi->pix.len = ((poff + np) <= theimg->npix) ? np : (theimg->npix - poff); for (i = loff; i < (loff + theroi->lin.len); i++) for (j = poff; j < (poff + theroi->pix.len); j++) { px = theimg->data[band][i][j]; theroi->sum += px; theroi->sumsq += px * px; } } theroi->img = theimg; theroi->n = theroi->lin.len * theroi->pix.len; } /** * Copy the details of one single-band region of interest to another. * Equivalent to a bitwise copy of one ROI stucture to another, * but with some sanity checks. * @param dest the destination ROI (details overwritten). * @param src the source ROI, containing valid ranges and sums to copy. * @return the destination ROI. */ bandroip_t copy_bandroi (bandroip_t dest, const bandroip_t src) { if (src == NULL) { fprintf (stderr, "copy_bandroi: null ROI structure.\n"); exit (-2); } if (src->img == NULL) { fprintf (stderr, "copy_bandroi: ROI without an associated image.\n"); exit (-2); } if ((src->band >= src->img->nbnd) || (src->lin.off >= src->img->nlin) || ((src->lin.off + src->lin.len) > src->img->nlin) || (src->pix.off >= src->img->npix) || ((src->pix.off + src->pix.len) > src->img->npix)) { fprintf (stderr, "copy_bandroi: ROI outwith image bounds.\n"); exit (-2); } if (src->n != ((long int) src->pix.len * (long int) src->lin.len)) { fprintf (stderr, "copy_bandroi: ROI total number of pixels incorrect.\n"); exit (-2); } if (src->sumsq < 0) { fprintf (stderr, "copy_bandroi: ROI sum of value squares corrupt.\n"); exit (-2); } *dest = *src; return dest; } /** * Compute the sum of the products of the pixel values in two identical ROIs. * @param s the first single-band rectangular ROI. * @param t the second single-band rectangular ROI. * @return the sum of the products of corresponding pixel values. */ long double product_bandroi (bandroip_t s, bandroip_t t) { long double ps = 0; size_t i, j; if ((s->lin.len != t->lin.len) || (s->pix.len != t->pix.len)) { fprintf (stderr, "ROI mismatch %d, %d versus %d, %d\n", s->lin.len, s->pix.len, t->lin.len, t->pix.len); exit (-2); } for (i = 0; i < s->lin.len; i++) for (j = 0; j < s->pix.len; j++) { ps += (s->img->data[s->band][s->lin.off + i][s->pix.off + j] * t->img->data[t->band][t->lin.off + i][t->pix.off + j]); } return ps; } /** * Determine the ovelap interval in a particular dimension (lines or pixels). * @param off the displacement of the second image relative to the first. * @param ssize the first image length in this dimension. * @param tsize the second image length. * @param lap set to the newly-computed overlap interval. */ void clip_overlap (long int off, uint32_t ssize, uint32_t tsize, overlapp_t lap) { if (off >= 0) { lap->soff = (off < (long int) ssize) ? (uint32_t) off : ssize; lap->toff = 0; lap->len = (tsize < (ssize - lap->soff)) ? tsize : (ssize - lap->soff); } else { lap->soff = 0; lap->toff = ((-off) < (long int) tsize) ? (uint32_t) (-off) : tsize; lap->len = (ssize < (tsize - lap->toff)) ? ssize : (tsize - lap->toff); } } /** * Shift the index range covered by a rectangular ROI in a single dimension. * Given a previously established range of indices for an ROI in a * single dimension (e.g., pixels or lines), and a relative * displacement for this range, determine the new range for the ROI as * a whole, and two secondary ranges: the range whose values must be * subtracted to update the ROI, and the range whose values must be * added. The actual additions and subtractions happen elsewhere. * @param delta the relative displacement of the range. * @param limit clip the ranges to fit in the interval [0..@a limit). * @param old the previously established ROI range to be shifted. * @param all the complete ROI range after shifing. * @param sub the subrange whose values must be subtracted in the update. * @param add the subrange whose values must be added in the update. */ void shift_ranges (long int delta, uint32_t limit, dim_rangep_t old, dim_rangep_t all, dim_rangep_t sub, dim_rangep_t add) { long int new_off; new_off = old->off + delta; if (new_off >= limit) { all->off = limit; all->len = 0; } else if (new_off < 0) { all->off = 0; all->len = ((old->len + new_off) > 0) ? (uint32_t) (old->len + new_off) : 0; } else { all->off = (uint32_t) new_off; all->len = ((new_off + old->len) > limit) ? (limit - all->off) : old->len; } if (delta < 0) { sub->off = all->off + all->len; add->off = all->off; add->len = old->off - all->off; } else { sub->off = old->off; add->off = old->off + old->len; add->len = (all->off + all->len) - (old->off + old->len); } sub->len = (uint32_t) (labs (delta)); } /** * Displace a rectangular region of interest by the specified number of lines. * Given a previously establised single-band rectangular ROI within an image, * displace it by a certain number of lines, clipping the revised ROI at * the edge of the image. If the old and new ROIs overlap by a substantial * amount (arbitrarily set to at least 2/3 of the area), derive the new * ROI from the old by adding and subtracting values rather than * re-computing it from scratch. * @param theroi the region of interest to displace. * @param delta the relative displacement (plus or minus number of lines). */ static void lineshift_bandroi (bandroip_t theroi, long int delta) { Dim_range_t all, sub, add; uint32_t ladd, lsub, ipx, jpx; long double add_px, sub_px; if (theroi == NULL) return; shift_ranges (delta, theroi->img->nlin, &theroi->lin, &all, &sub, &add); /* It's simpler to re-initialize if more than 1/3 of the area will change. */ if ((3 * labs (all.off - theroi->lin.off)) > theroi->lin.len) initialize_bandroi (theroi, theroi->img, theroi->interval, theroi->band, all.off, all.len, theroi->pix.off, theroi->pix.len); else { for (lsub = sub.off; lsub < (sub.off + sub.len); lsub++) { for (ipx = theroi->pix.off; ipx < theroi->pix.off + theroi->pix.len; ipx++) { sub_px = theroi->img->data[theroi->band][lsub][ipx]; theroi->sum -= sub_px; theroi->sumsq -= sub_px * sub_px; } theroi->n -= theroi->pix.len; } for (ladd = add.off; ladd < (add.off + add.len); ladd++) { for (jpx = theroi->pix.off; jpx < theroi->pix.off + theroi->pix.len; jpx++) { add_px = theroi->img->data[theroi->band][ladd][jpx]; theroi->sum += add_px; theroi->sumsq += add_px * add_px; } theroi->n += theroi->pix.len; } theroi->lin = all; } } /** * Displace a rectangular region of interest by the specified number of pixels. * Given a previously establised single-band rectangular ROI within an * image, displace it by a certain number of pixels within each line * (keeping the range of lines constant), and clip the revised ROI at the * edge of the image. If the old and new ROIs overlap by a * substantial amount (arbitrarily set to at least 2/3 of the area), * derive the new ROI from the old by adding and subtracting values * rather than re-computing it from scratch. * @param theroi the region of interest to displace. * @param delta the relative displacement (plus or minus number of pixels). */ static void pixelshift_bandroi (bandroip_t theroi, long int delta) { Dim_range_t all, sub, add; uint32_t line, iadd, isub; long double add_px, sub_px; int line_change; if (theroi == NULL) return; shift_ranges (delta, theroi->img->npix, &theroi->pix, &all, &sub, &add); /* It's simpler to re-initialize if more than 1/3 of the area will change. */ if ((3 * labs (all.off - theroi->pix.off)) > theroi->pix.len) initialize_bandroi (theroi, theroi->img, theroi->interval, theroi->band, theroi->lin.off, theroi->lin.len, all.off, all.len); else { line_change = (int) (add.len - sub.len); for (line = theroi->lin.off; line < theroi->lin.off + theroi->lin.len; line++) { for (isub = sub.off; isub < (sub.off + sub.len); isub++) { sub_px = theroi->img->data[theroi->band][line][isub]; theroi->sum -= sub_px; theroi->sumsq -= sub_px * sub_px; } for (iadd = add.off; iadd < (add.off + add.len); iadd++) { add_px = theroi->img->data[theroi->band][line][iadd]; theroi->sum += add_px; theroi->sumsq += add_px * add_px; } theroi->n += line_change; } theroi->pix = all; } } /** * Compute the classic product-moment correlation coefficient, @a r, * then use Fisher's z-transformation to compute a @a t-value from this * (to get a goodness of fit index that's independent of the degrees * of freedom, which in this case are determined by the sizes of the * regions of interest). Use one of the standard textbook formulae, * computing a covariance term from all the values in both search and * target ROIs, and correcting this by the previously computed * variance terms for each ROI. Update the queue of optimum * @a t-values and their locations if necessary. * @param search a particular search region of interest. * @param target the one invariant target region of interest. * @param ss_target the sum of squares for the target ROI. * @param i the current search pattern column (pixel) index. * @param j the current search row (line) index. * @param apriori estimate of the relative overlap of the two images. * @param matches priority queue of previously determined @a r-value maxima. */ void maximize_r (const bandroip_t search, const bandroip_t target, long double ss_target, int i, int j, const xydisplacep_t apriori, mheapp_t matches) { long double ss_search; long double sum_products; long double r; long double t_temp; Mpos_t newmatch; ss_search = (search->sumsq - ((search->sum * search->sum) / search->n)); sum_products = (product_bandroi (search, target) - ((search->sum * target->sum) / search->n)); r = sum_products / sqrtl (ss_search * ss_target); t_temp = atanhl (r) * sqrtl (search->n - 3); newmatch.t = (t_temp < MIN_T) ? MIN_T : (double) t_temp; if ((!full_mheap (matches)) || (min_mheap (matches) < newmatch.t)) { newmatch.pos.dx = i + apriori->dx; newmatch.pos.dy = j + apriori->dy; add_mheap (matches, &newmatch); } } /** * Find the maximum goodness of fit values in a set of image overlaps. * Given a pair of images at a suggested overlap position, examine all the * other overlap positions close to this, and identify the positions that * maximize the goodness of fit index, @a t. Note that the search ROI * is not the same as the search pattern: the former is a substantial part of * one of the images, identical in size to the target ROI; the latter * is a (generally much smaller) set of coordinate pairs, each pair defining * a potential overlap position for the two images. The search order, * outwards from the center of the pattern, reduces the maximum length * of the chains of updating operations needed by the algorithm. Use any * previously recorded optima in the heap structure that @a matches points to, * replacing them or leaving them in place as necessary. * @param search_image the image within which to move a search ROI. * @param target_image contains a fixed ROI to match with the search ROI. * @param prm points to the registration search parameter structure. * @param apriori points to the structure defining the a priori position. * @param xlap defines the horizontal overlap interval in the images. * @param ylap defines the vertical overlap interval in the images. * @param matches a heap structure recording the optimal match positions. */ void find_optima (const IMAGE * search_img, const IMAGE * target_img, const regparamp_t prm, const xydisplacep_t apriori, const overlapp_t xlap, const overlapp_t ylap, mheapp_t matches) { Bandroi_t target; Bandroi_t search_centre, search_line, search_pixel; long double ss_target; int x_size, y_size; int i_lo, i_hi, j_lo, j_hi, k_lo, k_hi; x_size = 1 + 2 * prm->x_radius; y_size = 1 + 2 * prm->y_radius; initialize_bandroi (&target, target_img, 0, 0, (xlap->toff + prm->x_radius), (ylap->len - y_size), (ylap->toff + prm->y_radius), (xlap->len - x_size)); ss_target = target.sumsq - (target.sum * target.sum) / target.n; initialize_bandroi (&search_centre, search_img, 0, 0, (ylap->soff + prm->y_radius), (ylap->len - y_size), (xlap->soff + prm->x_radius), (xlap->len - x_size)); copy_bandroi (&search_line, &search_centre); for (j_hi = prm->y_radius; j_hi < y_size; j_hi++) { copy_bandroi (&search_pixel, &search_line); for (i_hi = prm->x_radius; i_hi < x_size; i_hi++) { maximize_r (&search_pixel, &target, ss_target, i_hi, j_hi, apriori, matches); if (i_hi < (x_size - 1)) pixelshift_bandroi (&search_pixel, +1); } copy_bandroi (&search_pixel, &search_line); for (i_lo = (prm->x_radius - 1); i_lo >= 0; i_lo--) { pixelshift_bandroi (&search_pixel, -1); maximize_r (&search_pixel, &target, ss_target, i_lo, j_hi, apriori, matches); } if (j_hi < (y_size - 1)) lineshift_bandroi (&search_line, +1); } copy_bandroi (&search_line, &search_centre); for (j_lo = (prm->y_radius - 1); j_lo >= 0; j_lo--) { lineshift_bandroi (&search_line, -1); copy_bandroi (&search_pixel, &search_line); for (k_hi = prm->x_radius; k_hi < x_size; k_hi++) { maximize_r (&search_pixel, &target, ss_target, k_hi, j_lo, apriori, matches); if (k_hi < (x_size - 1)) pixelshift_bandroi (&search_pixel, +1); } copy_bandroi (&search_pixel, &search_line); for (k_lo = (prm->x_radius - 1); k_lo >= 0; k_lo--) { pixelshift_bandroi (&search_pixel, -1); maximize_r (&search_pixel, &target, ss_target, k_lo, j_lo, apriori, matches); } } } /** * Determine registration offsets between a pair of images. * Initialize the structures that define a particular overlap position, * and update the priority queue of optimum goodness of fit indices (and * their associated positions) in @a matches. * @param search_img pointer to a SADIE image, the first input image. * @param target_img pointer to a SADIE image, the second input image. * @param prm points to the registration search parameter structure. * @param est_off Xydisplace_t estimate of the x- and y-offset of the images. * @param matches pointer to a heap data structure (queue of max. t-values). * @return a non-zero error code on failure. */ static regerr_t multi_register (const IMAGE * search_img, const IMAGE * target_img, const regparamp_t prm, const xydisplacep_t est_off, mheapp_t matches) { Overlap_t xlap, ylap; Xydisplace_t apriori; if (((int64_t)est_off->dx > (int64_t)search_img->npix) || ((int64_t)est_off->dy > (int64_t)search_img->nlin)) return REGERR_NOLAP; clip_overlap (est_off->dx, search_img->npix, target_img->npix, &xlap); clip_overlap (est_off->dy, search_img->nlin, target_img->nlin, &ylap); if ((xlap.len < (1 + 2 * prm->x_radius)) || (ylap.len < (1 + 2 * prm->y_radius))) return REGERR_TOOSMALL; apriori.dx = est_off->dx - prm->x_radius; apriori.dy = est_off->dy - prm->y_radius; find_optima (search_img, target_img, prm, &apriori, &xlap, &ylap, matches); return REGERR_OK; } /** * Find a single optimum overlap position between two image pyramids. * Given two image pyramids and a set of search parameters (radii of * the search pattern in both horizontal and vertical directions), * refine an estimate of the position of the target images relative to the * search images. Rank the ovelap positions by a goodness of fit * index for the images, and discard all but a fixed number * of best positions by comparing the coarsest resolution images at all * possible positions in a search pattern for that resolution, centered on * the initial estimate of the overlap. Use these estimates as the centers * of searches in successively finer resolutions down the image pyramid, * culminating in a search restricted to a single optimum overlap position * at the finest resolution. * @param search_pyr points to the search image pyramid. * @param target_pyr points to the target image pyramid. * @param param_pyr points to search parameters for each pyramid level. * @param est_off an initial estimate of the relative image overlap. * @param actual_off set to the refined estimate of the image overlap. * @return zero on success, or an error code on failure. */ static regerr_t pyramid_register (const regimage_pyramidp_t search_pyr, const regimage_pyramidp_t target_pyr, const regparam_pyramidp_t param_pyr, const xydisplacep_t est_off, xydisplacep_t actual_off) { size_t qsize; int scale; int i_l, j, k, s; mheapp_t sought, found; Xydisplace_t scaled_est_off; regerr_t err; if (!search_pyr || !target_pyr || !param_pyr || !est_off || !actual_off) return REGERR_EMPTY; if ((search_pyr->nlevels != target_pyr->nlevels) || (target_pyr->nlevels != param_pyr->nlevels) || (search_pyr->ratio != target_pyr->ratio) || (target_pyr->ratio != param_pyr->ratio)) return REGERR_MISMATCH; i_l = param_pyr->nlevels - 1; if (i_l < 0) return REGERR_BADPARAM; /* We will handle the case where i_l == 0 (a degenerate 1-level pyramid). */ found = NULL; qsize = ((i_l) ? ((size_t) (rint (sqrt ((param_pyr->plevel[i_l].x_radius * param_pyr->plevel[i_l].x_radius) + (param_pyr->plevel[i_l].y_radius * param_pyr->plevel[i_l].y_radius))))) : 1); sought = make_mheap (qsize); for (s = i_l, scale = 1; s > 0; s--) scale *= param_pyr->ratio; /* Search at the coarsest resolution of the pyramid first. */ scaled_est_off.dx = (int) rint ((double) est_off->dx / (double) scale); scaled_est_off.dy = (int) rint ((double) est_off->dy / (double) scale); err = multi_register (search_pyr->ilevel[i_l], target_pyr->ilevel[i_l], ¶m_pyr->plevel[i_l], &scaled_est_off, sought); /* Refine the search results at successively finer resolutions. */ while (!err && (--i_l > 0)) { dispose_mheap (found); found = sought; sought = make_mheap (qsize); for (j = 1; !err && (j <= found->n); j++) { scaled_est_off.dx = param_pyr->ratio * found->q[j].pos.dx; scaled_est_off.dy = param_pyr->ratio * found->q[j].pos.dy; err = multi_register (search_pyr->ilevel[i_l], target_pyr->ilevel[i_l], ¶m_pyr->plevel[i_l], &scaled_est_off, sought); } } /* For the final search, find a single optimum at the finest resolution. */ if (!err && (i_l == 0)) { dispose_mheap (found); found = sought; sought = make_mheap (1); for (k = 1; !err && (k <= found->n); k++) { scaled_est_off.dx = param_pyr->ratio * found->q[k].pos.dx; scaled_est_off.dy = param_pyr->ratio * found->q[k].pos.dy; err = multi_register (search_pyr->ilevel[0], target_pyr->ilevel[0], ¶m_pyr->plevel[0], &scaled_est_off, sought); } } *actual_off = (err) ? *est_off : sought->q[1].pos; dispose_mheap (found); dispose_mheap (sought); return err; } /** * Determine the best registration offsets of two overlapping images. * Starting with some a priori values, search various possible * horizontal and vertical offsets for the combination that give the * best goodness of fit index between the ovelapping areas of the images. * Determine the offset and bias correction once the accurate offsets * are available. * @param img1 the first input image. * @param img2 the second input image. * @param est_x_off the estimated x-offset of the second image. * @param est_y_off the estimated y-offset of the second image. * @param x_off the actual x-offset of the second image. * @param y_off the actual y-offset of the second image. * @param bias_adj the bias adjustment for the second image. * @param gain_adj the gain adjustment for the second image. */ void FIND_REGISTER (IMAGE * img1, IMAGE * img2, int est_x_off, int est_y_off, int *x_off, int *y_off, PIXEL * bias_adj, PIXEL * gain_adj) { char msg[SLEN]; regimage_pyramidp_t search_pyr, target_pyr; Xydisplace_t est_off, actual_off; struct timeval start, find_first_reg, findgain, end; gettimeofday (&start, NULL); if (!CHECKIMG (img1)) { MESSAGE ('E', " Can't identify first input image."); return; } if (!CHECKIMG (img2)) { MESSAGE ('E', " Can't identify second input image."); return; } if (NAMES) { MESSAGE ('I', ""); MESSAGE ('I', "FIND_REGISTER"); MESSAGE ('I', ""); sprintf (msg, " First input image: %s", img1->text); MESSAGE ('I', msg); sprintf (msg, " Second input image: %s", img2->text); MESSAGE ('I', msg); MESSAGE ('I', ""); MESSAGE ('I', ""); sprintf (msg, " Estimated x-offset: %d", est_x_off); MESSAGE ('I', msg); sprintf (msg, " Estimated y-offset: %d", est_y_off); MESSAGE ('I', msg); } if (((int64_t)est_x_off > (int64_t)img1->npix) || ((int64_t)est_y_off > (int64_t)img1->nlin)) { MESSAGE ('E', " Images cannot be registered because they do not overlap!"); return; } search_pyr = make_regpyramid (img1, default_param.nlevels, default_param.ratio); if (search_pyr) { target_pyr = make_regpyramid (img2, default_param.nlevels, default_param.ratio); if (target_pyr == NULL) { dispose_regpyramid (search_pyr); search_pyr = NULL; } } if (search_pyr == NULL) { MESSAGE ('E', " Unable to create the image pyramids to register."); return; } gettimeofday (&find_first_reg, NULL); est_off.dx = est_x_off; est_off.dy = est_y_off; switch (pyramid_register (search_pyr, target_pyr, &default_param, &est_off, &actual_off)) { case REGERR_OK: *x_off = actual_off.dx; *y_off = actual_off.dy; printf ("Registration results: (%d,%d) ==> (%d,%d)\n", est_x_off, est_y_off, *x_off, *y_off); fprintf (stdout, "deltamax = (%d,%d)\n", *x_off - est_x_off, *y_off - est_y_off); gettimeofday (&findgain, NULL); FINDGAINADJ (img1, img2, *x_off, *y_off, gain_adj, bias_adj); gettimeofday (&end, NULL); if (NAMES) { MESSAGE ('I', ""); sprintf (msg, " Computed x-offset: %d", *x_off); MESSAGE ('I', msg); sprintf (msg, " Computed y-offset: %d", *y_off); MESSAGE ('I', msg); MESSAGE ('I', ""); sprintf (msg, "Time to compute match = %ld ms.", delay (find_first_reg, findgain)); MESSAGE ('I', msg); MESSAGE ('I', msg); sprintf (msg, "Time to compute gain/bias registration = %ld ms.", delay (findgain, end)); MESSAGE ('I', msg); sprintf (msg, "Total time required = %ld ms.", delay (start, end)); MESSAGE ('I', msg); MESSAGE ('I', " ..............."); } break; case REGERR_EMPTY: MESSAGE ('E', " Unexpected empty data structure (supposedly impossible)."); break; case REGERR_MISMATCH: MESSAGE ('E', " Mismatch in image pyramid settings."); break; case REGERR_BADPARAM: MESSAGE ('E', " Corrupt search parameters."); break; case REGERR_NOLAP: MESSAGE ('E', " The images do not overlap."); break; case REGERR_TOOSMALL: MESSAGE ('E', " The image overlap area is too small."); break; } dispose_regpyramid (search_pyr); dispose_regpyramid (target_pyr); }