Patent Publication Number: US-6219462-B1

Title: Method and apparatus for performing global image alignment using any local match measure

Description:
This application claims benefit of U.S. provisional patent application Ser. No. 60/046,069, filed May 9, 1997 and incorporated herein by reference. 
    
    
     The present invention generally relates to image processing systems and, more particularly, to a method and apparatus for aligning images within a image processing system. 
     BACKGROUND OF THE DISCLOSURE 
     In many image processing systems it is necessary that images be aligned with one another to perform image merging or analysis. The phrase image processing, as used herein, is intended to encompass the processing of all forms of images including temporally unrelated images as well as images (frames) of a video signal, i.e., a sequence of temporally related images. Image alignment in an image processing system is necessary to create mosaics of multiple images, perform some forms of image compression, perform motion estimation and/or tracking and the like. Alignment (also known as registration) of images begins with determining a displacement field that represents the offset between the images and then warping one image to the other to remove or minimize the offset. The images may be taken from the same sensor or from two entirely different sensors, possibly of different modality (e.g., an infrared sensor and a visible sensor). Often, the displacement field that defines the offset between the images can be described as a global parametric transformation between the two images, such as an affine, quadratic, or a projective transformation. Many techniques have been developed for the parametric transformation of a pair of images. 
     Most flow-based techniques divide the registration process into two steps: first a flow-field is estimated, then, using regression, the global parametric transformation which best describes the flow field is found. However, often the local flow estimates are noisy and unreliable, resulting in poor registration accuracy and a lack of robustness. 
     To overcome this problem, direct gradient-based techniques have been developed. These techniques estimate the global transformation parameters by directly using local image intensity information without first computing a local flow-field. They achieve high registration accuracy since they avoid the noisy step of computing a local flow-estimation. However, these techniques assume that the intensity values of corresponding pixels in the two images are the same, which is known as the “brightness constancy assumption”. As a result, the applicability of direct gradient-based techniques is limited to situations when the images to be registered are substantially similar in appearance. Consequently, the direct gradient-based techniques cannot handle large changes in illumination and contrast between the images. Because the images need to be substantially similar to be registered using direct gradient-based techniques, images produced by sensors having different modality and/or containing a substantial range of motion cannot be accurately registered by direct gradient-based techniques. 
     Therefore, a need exists in the art for a method and apparatus that aligns images having substantial illumination differences between the images and/or a substantial amount of motion and/or other image differences that would otherwise make registration difficult. 
     SUMMARY OF THE INVENTION 
     The disadvantages of the prior art are overcome by the present invention of a method, which when given any local match measure, applies global estimation directly to the local match measure data, without first performing an intermediate step of local flow estimation. Any local match measure can be used as part of the inventive method, such as correlation, normalized-correlation, squared or absolute brightness difference, statistical measures such as mutual information, and the like. Global estimation constrains the analysis of the local match measure, thereby avoiding noisy local motion estimates, while still providing an accurate result for image registration. 
     In one embodiment of the invention, the inventive generalized global alignment method is used with a normalized-correlation match-measure to result in a global correlation-based alignment method which combines the robustness and accuracy of global alignment with the broad applicability of the normalized-correlation match measure. The inventive method overcomes many of the limitations of existing gradient-based and flow-based techniques. In particular, the novel method can handle large appearance differences between the images and large image motion within the image scene. Also, in contrast to the flow-based methods, the invention accurately registers imagery with sparse texture or feature content. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The teachings of the present invention can be readily understood by considering the following detailed description in conjunction with the accompanying drawings, in which: 
     FIG. 1 depicts a block diagram of an image processing system of the present invention; 
     FIG. 2 depicts a flow diagram representing the method of the present invention; and 
     FIG. 3 depicts a flow diagram representing a second embodiment of the method of the present invention. 
    
    
     To facilitate understanding, identical reference numerals have been used, where possible, to designate identical elements that are common to the figures. 
     DETAILED DESCRIPTION 
     FIG. 1 depicts a block diagram of an image processing system  100  comprising at least one image sensor  104  and an image processor  106 . The present invention can be used to process images produced by multiple sensors (e.g., sensor  104  and  104 A) and is especially useful in processing images that are produced by sensors having different modalities, e.g., an infrared wavelength sensor  104 A and a visible wavelength sensor  104 . The sensors  104  and  104 A are intended to provide imagery of a three dimensional scene  102 . 
     To process these images, the image processor  106  comprises a central processing unit (CPU)  108 , a memory device  110 , conventional CPU support circuits  112  and input/output (I/O) peripherals  114 . The CPU  108  is a general purpose computer that, when executing specific routines that are recalled form memory  110 , becomes a specific purpose computer. The CPU can be any high-performance processor such as a PENTIUM II processor manufactured by Intel Corporation or a POWER PC processor manufactured by Motorola Inc. The memory  110  can be random access memory (RAM), read only memory (ROM), a hard disk drive, a floppy disk drive or any combination thereof. The support circuits  112  include various conventional circuits such as frame grabber circuits, analog-to-digital (A/D) circuits, clock circuits, cache, power supplies, and the like. The I/O peripherals  114  generally include a keyboard, a mouse, and a display, but may also include a video tape recorder, a video disk player, and the like. The images that are processed by the image processor  106  may not be sourced by the sensors  104 , but may also be sourced from pre-recorded material such as would be provided by a video tape recorder or other image storage device. 
     The present invention is a routine  116  that, when executed by the CPU  108 , provides a method, which when given any local match measure, applies global estimation directly to the local match measure data, without first performing an intermediate step of local flow estimation. Any local match measure can be used as part of the inventive method, such as correlation, normalized-correlation, squared or absolute brightness difference, statistical measures such as mutual information, and the like. Global estimation constrains the analysis of the local match measure, thereby avoiding noisy local motion estimates, while still providing an accurate result for image registration. 
     More specifically, when a pair of images representing a three dimensional scene  102  are to be aligned and merged, i.e., the images are to be registered, the induced image motion between two images depends on the cameras&#39; internal parameters such as zoom and external parameters such as actual camera motion, as well as parallax motion within the 3D scene structure. In many practical situations, the induced motion field between the image pair can be modeled in terms of one or a small number of parametric transformations. Such transformations are well known in the art as represented in J. R. Bergen et al., “Hierarchical Model-based Motion Estimation,” European Conference on Computer Vision, pp. 237-252, Santa Margarita Ligure, May 1992. The present invention expands upon this parametric transformation model by using a two-dimensional parametric transformation technique, although the approach generalizes to other classes of models as well. 
     Specifically, when the motion field is a linear function of a few unknown parameters {p i }, then the motion vector {right arrow over (u)}(x,y)=(u(x,y),v(x,y)) T  can be expressed as: 
     
       
           {right arrow over (u)} ( x,y;{right arrow over (p)} )= X ( x,y )· {right arrow over (p)}   (1) 
       
     
     where X(x,y) is a matrix which depends only on the pixel coordinates (x,y), and {right arrow over (p)}=(p 1 , . . . , p n ) T  is the parameter vector. For example, for an affine transformation:          [           u        (     x   ,     y   ;     p   →         )                 v        (     x   ,     y   ;     p   →         )             ]     =     [             p   1     +       p   2        x     +       p   3        y                   p   4     +       p   5        x     +       p   6        y             ]                     
     therefore, in this case: {right arrow over (p)}=(p 1 ,p 2 ,p 3 ,p 4 ,p 5 ,p 6 ) T  and        X   =     [         1       x       y       0       0       0           0       0       0       1       x       y         ]                     
     and for a quadratic transformation:          [           u        (     x   ,     y   ;     p   →         )                 v        (     x   ,     y   ;     p   →         )             ]     =     [             p   1     +       p   2        x     +       p   3        y     +       p   7          x   2       +       p   8        xy                   p   4     +       p   5        x     +       p   6        y     +       p   7        xy     +       p   8          x   2               ]                     
     therefore: {right arrow over (p)}=(P 1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 , P 8 ) T  and        X   =     [         1       x       y       0       0       0         x   2         xy           0       0       0       1       x       y       xy         y   2           ]                     
     Given two images f and g, the local match measure m(u,v;x,y) computes the similarity (or dissimilarity) between the pixels f(x,y) and g(x+u, y+v). Some of the common match measures used in existing registration techniques include: (i) correlation, (ii) normalized-correlation, which is generally preferred over regular correlation, since it is invariant to local changes in mean and contrast, (iii) squared brightness differences, which is applicable under the brightness constancy assumption, (iv) sum of squared brightness differences (SSD), and (v) mutual information which measures the statistical correlation between the two signals. 
     The present invention performs parametric registration of two images, f and g, by determining a parametric transformation {right arrow over (p)} that maximizes (or minimizes) a global match-measure M({right arrow over (p)}):                M        (     p   →     )       =       ∑     x   ,   y            m        (       u        (     x   ,     y   ;     p   →         )       ,       v        (     x   ,     y   ;     p   →         )       ;   x     ,   y     )                 (   2   )                         
     where (u(x,y;{right arrow over (p)})),v(x,y;{right arrow over (p)})) is the motion field described by the parametric transformation {right arrow over (p)}, and m(u,v;x,y) measures the local match between pixel (x,y) in f and pixel (x+u, y+v) in g. For example, when a correlation match-measure is used, then:            m        (     u   ,     v   ;   x     ,   y     )            =   def            ∑       (     i   ,   j     )        εW              f        (       x   +   i     ,     y   +   j       )            g        (       x   +   u   +   i     ,     y   +   v   +   j       )             ;                   
     when a squared-brightness-difference measure is used, then:          m        (     u   ,     v   ;   x     ,   y     )            =   def            (       f        (   x   )       -     g        (       x   +   u     ,     y   +   v       )         )     2                     
     The invention determines a global transformation {right arrow over (p)} that maximizes (or minimizes) the sum of local match-measures m. The routine of the present invention estimates {right arrow over (p)} via regression directly on surfaces of the local match-measure m, without first computing locally to specific image displacements {(u(x,y),v(x,y))}. 
     The inventive routine, first defines a match-measure surface: Let S (x,y) (u,v) denote a local match-measure surface corresponding to pixel (x,y) in f. For any shift (u,v) of g relative to f, S (x,y)  is defined as:            S     (     x   ,   y     )            (     u   ,   v     )            =   def          m        (     u   ,     v   ;   x     ,   y     )                       
     For example, when the local match-measure m(u,v;x,y) is derived using a correlation technique, then S (x,y)  is the correlation surface of pixel (x,y). Equation (2) can therefore be rewritten in terms of the collection of surfaces {S (x,y) }:                      M        (     p   →     )       =                  ∑     (     x   ,   y     )              S     (     x   ,   y     )            (       u        (     x   ,     y   ;     p   →         )       ,     v        (     x   ,     y   ;     p   →         )         )                     =                  ∑     (     x   ,   y     )              S     (     x   ,   y     )       (       u   →          (     x   ,     y   ;     p   →         )                         (   3   )                         
     For compactness of notation, {right arrow over (u)} denotes the two-dimensional vector (u,v). To solve for {right arrow over (p)} that maximizes M({right arrow over (p)}), the routine uses Newton&#39;s method, which is an iterative method. Let {right arrow over (p)} o  denote the parametric transformation that was computed in the previous iteration step. A second order Taylor expansion of M({right arrow over (p)}) around {right arrow over (p)} o  yields: 
     
       
           M ( {right arrow over (p)} )= M ( {right arrow over (p)}   0 +(∇ {right arrow over (p)}   M ( {right arrow over (p)}   0 )) T {right arrow over (δ)} p +{right arrow over (δ)} p   T   H   M ( {right arrow over (p)}   0 ){right arrow over (δ)} p   (4) 
       
     
     where, {right arrow over (δ)} p ={right arrow over (p)}−{right arrow over (p)} o  is the unknown refinement step of {right arrow over (p)} 0  that is to be solved for, ∇ p M denotes the gradient of M (i.e., first derivatives of M), and H M  denotes the Hessian of M (i.e., second derivatives of M):            ▽   p          M        (     p   →     )              =   def            (         ∂   M       ∂     p   1         ,       ∂   M       ∂     p   2         ,   …   ,       ∂   M       ∂     p   n           )     T                     
     In order to find the parametric transformation {right arrow over (p)} that maximizes M({right arrow over (p)}), the right-hand side of Equation (4) is differentiated with respect to {right arrow over (δ)} p  and set equal to zero: 0=∇ p M({right arrow over (p)} 0 )+H M ({right arrow over (p)} 0 ){right arrow over (δ)} p . Therefore: 
     
       
         {right arrow over (δ)} p =−( H   M ( {right arrow over (p)}   0 )) −1 ·∇ p   M ( {right arrow over (p)}   0 )  (5) 
       
     
     where {right arrow over (δ)} p * denotes {right arrow over (δ)} p  which maximizes the quadratic approximation of M({right arrow over (p)}) around {right arrow over (p)} o  (Equation (4). Hence, given {right arrow over (p)} o  from the previous iteration, the system estimates a refinement step {right arrow over (δ)} p * to obtain a better estimate {right arrow over (δ)} p ={right arrow over (p)}−{right arrow over (δ)} p *. The new {right arrow over (p)} is then used as an initial estimate {right arrow over (p)} o  for the next iteration. 
     Expressing the right-hand-side of Equation 5 in terms of the measurable quantities {S (x,y) }                 ▽   p          M        (     p   →     )         =         ∑     x   ,   y            (       ▽   p          S        (     u   →     )         )       =       ∑     x   ,   y            (         X   T     ·     ▽   u            S        (     u   →     )         )                 (   6   )                   H   M          (     p   →     )       =       ∑     x   ,   y            (       X   T     ·       H   S          (     u   →     )       ·   X     )                                         
     When X is the matrix defined in Equation (1) (i.e., {right arrow over (u)}(x,y;{right arrow over (p)})=X(x,y)·{right arrow over (p)}); ∇ u S({right arrow over (u)}) is the gradient of S (x,y) {right arrow over (u)}(x,y;{right arrow over (p)}), and H S  is the Hessian of S (x,y) {right arrow over (u)}(x,y;{right arrow over (p)}), then Equation 6 becomes:            ▽   u          S        (     u   →     )              =   def            (         ∂     S     (     x   ,   y     )           ∂   u       ,       ∂     S     (     x   ,   y     )           ∂   v         )       T                              H   S          (     u   →     )            =   def          [               ∂   2          S     (     x   ,   y     )           ∂     u   2                   ∂   2          S     (     x   ,   y     )             ∂   u          ∂   v                       ∂   2          S     (     x   ,   y     )             ∂   u          ∂   v                   ∂   2          S     (     x   ,   y     )           ∂     v   2               ]                     
     Substituting Equation (6) into Equation (5) provides an expression for the refinement step {right arrow over (δ)} p * in terms of the local match-measure surfaces {S (x,y) }.                    δ   →     p     *     =       -       (       ∑     (     x   ,   y     )              X   T     ·       H   S          (       u   →     0     )       ·   X       )       -   1         ·     (       ∑     (     x   ,   y     )                X   T     ·     ▽   u            S        (       u   →     0     )           )               (   7   )                         
     where {right arrow over (u)} 0 ={right arrow over (u)}(x,y;{right arrow over (p)} 0 ) is the displacement induced at pixel (x,y) by the estimated parametric transformation {right arrow over (p)} o  at the previous iteration. The process represented by Equation (7) functions as follows: at each iteration, for each pixel (x,y), the local quadratic approximation of the pixel&#39;s match-measure surface S (x,y) ({right arrow over (u)}) is computed around its previously estimated displacement {right arrow over (u)} 0 . These local quadratic approximations are then used in Equation (7) to solve for the global parametric refinement {right arrow over (δ)} p * 
     To account for large misalignments between the two images, the system performs multi-resolution coarse-to-fine estimation, e.g., as in J. R. Bergen et al., “Hierarchical Model-based Motion Estimation,” European Conference on Computer Vision, pp. 237-252, Santa Margarita Ligure, May 1992. To facilitate a coarse-to-fine alignment process, a Gaussian (or a Laplacian) pyramid is constructed for each of the two images. Although Gaussian and/or Laplacian pyramids are generally used, other forms of image decomposition may be appropriate such as the use of wavelets and the like. As such, the regressive process described above is used at each level of the pyramids to accurately process the images. 
     FIG. 2 depicts a flow diagram of a routine  116  representing the method of the present invention as executed by the CPU of FIG.  1 . The routine begins at step  200  and proceeds to step  202  where the pyramid processing is performed to decompose a pair of images f and g into a plurality of pyramid levels. Images f and g are supplied by one or more sensors, as described above, and digitally stored in memory as arrays of pixel values. 
     In the following description, f l  and g l  denote the images at resolution level l in the pyramids of images f and g, respectively. At step  204 , the routine selects the coarsest level of the pyramids of each image. Starting at the selected coarsest resolution level with {right arrow over (p)} o  initially set to 0, the following steps are performed at each resolution level of the pyramids. 
     At step  206 , for each pixel (x,y) at f l  compute a local match measure surface around {right arrow over (u)} o  (i.e., around the displacement estimated at the previous iteration: {right arrow over (u)} 0 ={right arrow over (u)}(x,y;{right arrow over (p)} 0 )=X{right arrow over (p)} 0 ). In practice, the match-measure surface is estimated only for displacements {right arrow over (u)}of g l  within a radius d around {right arrow over (u)} 0 , i.e.: 
     
       
           S   l   (x,y) ( {right arrow over (u)} )= m ( u,v;x,y ),∀ {right arrow over (u)} =( u,v ) s.t.∥{right arrow over (u)}−{right arrow over (v)}   0    ∥≦d   
       
     
     where the distance d of possible {right arrow over (u)} from {right arrow over (u)} 0  is determined by the size of the masks used for discretely estimating the first and second order derivatives of S (x,y) ({right arrow over (u)}) at {right arrow over (u)}. Although other masking functions can be used, in the present embodiment of the invention, Beaudet&#39;s masks were used to estimate the first and second order derivatives of the surfaces. Beaudet&#39;s masks are described in detail in Paul R. Beaudet, “Rotationally Invariant Image Operators,” International Conference on Pattern Recognition, pp. 579-583 (1978). Satisfactory results have been found when using 3×3 masks (i.e., d=1) or with 5×5 masks (i.e., d=2). 
     At step  208 , the routine performs the regression step of Equation. (7) to determine the parametric refinement {right arrow over (δ)} p *. 
     At step  210 , the parametric refinement is used to update {right arrow over (p)} o :{right arrow over (p)} o :={right arrow over (p)} o +{right arrow over (δ)} p *. To refine the computational accuracy, the routine returns to step  206 . The loop represented by steps  206 ,  208  and  210  is repeated (step  218 ) until a predefined parametric transform accuracy is achieved for a given pyramid level. The predefined parametric transform accuracy varies from application to application and generally depends up the image resolution, the processing power of the computer performing the calculations, the complexity of the images being aligned and the like. After repeating the process loop of steps  206 ,  208 , and  210  for a few iterations (typically four), the parameters {right arrow over (p)} are propagated to the next resolution level, and the process is repeated at that resolution level. 
     The routine then queries, at step  212  whether another pyramid level is to be processed. If the query is affirmatively answered, the routine selects, at step  214 , another level of the pyramid, e.g., the next finest level, and returns to step  208 . Otherwise, the routine ends at step  216 . The process is generally stopped when the iterative process at the highest resolution level is completed. 
     In practice, to improve performance, an image warping step  220  can be added before each iteration. An image warping process that would be useful in this step is disclosed in the Bergen et al. paper cited above. To implement image warping, the image g (an inspection image) is warped towards the image f (the reference image) according to the current estimated parametric transformation {right arrow over (p)} o . After warping the images, {right arrow over (p)} o  is again set to 0, and {right arrow over (δ)} p * is estimated between the two warped images. Warping compensates for the spatial distortions between the two images (e.g., scale difference, rotations, etc.), and hence improves the quality of the local match-measures which are generally based on a window around a pixel, such as in correlation. 
     To further condition the regression step of Equation (7) and make the process more robust, only pixels (x,y), where a quadratic approximation of S (x,y) ({right arrow over (u)}) around {right arrow over (u)} 0  is concave, are used in the regression process, i.e., only concave surfaces are selected for use in the regression and all others are ignored. If M({right arrow over (p)}) is minimized rather than maximized with respect to {right arrow over (p)}, then only convex (rather than concave) pixels are used. All other pixels are ignored. This automatic outlier rejection mechanism provides the routine with a strong locking property. The locking property occurs because of the well-known phenomenon that motion in the image that is close to the current global motion estimate influences the subsequent motion estimates more effectively than does motion estimates that are very far from the current estimate, i.e., the iterative technique “locks on” to the motion that is close to the estimate. In the present invention, this locking property is augmented by the active rejection of local match surfaces whose shape is inconsistent with the hypothesis that there is a local minimum of the error (or maximum of the match) in the neighborhood of the current estimate. 
     This selective use of the match surfaces can be extended such that the system selects only those local match surfaces that meet a specific criteria to be used in the regression process. As such, once all the match surfaces are computed in step  206  of FIG. 2, that step may optionally prune the number of match surfaces used in the regression process using criteria such as the gradient (slope) of the surface, the surface peak, the surface shape, and the like. Such selection will remove surfaces that are, for example, flat and prone to be noisy. 
     In a specific embodiment of the invention, the generalized process described above can be used with a normalized-correlation function as the local match-measure. The normalized-correlation function is represented as:          ∑       (       f        (       x   +   i     ,     y   +   j       )       -       f   _     W       )     ·     (       g        (       x   +   u   +   i     ,     y   +   v   +   j       )       -       g   _     W       )                 (       f        (       x   +   i     ,     y   +   j       )       -       f   _     W       )     2       ·         (       g        (       x   +   u   +   i     ,     y   +   v   +   j       )       -       g   _     W       )     2                         
     where the summations are performed for (i,j)εW, and {overscore (fw)} and {overscore (gw)} denote the mean brightness value within corresponding windows around pixel (x,y) in f and pixel (x+u,y+v) in g, respectively. Normalized correlation is widely applicable, since it is invariant to local changes in mean and contrast, i.e., when the images f and g are linearly related within the window W. Since the windows W are usually chosen to be small, global alignment based on local normalized correlation is invariant to a variety of non-linear global brightness variations between the two images f and g. The global correlation-based alignment process therefore combines the robustness and accuracy of global estimation with the broad applicability and large capture range of normalized correlation. 
     In the present embodiment of the invention, the normalized-correlation is estimated using 9×9 correlation windows applied to Laplacian pyramid images. The quadratic approximations are based on 3×3 Beaudet&#39;s masks (i.e., d=1). The global correlation-based alignment process of the present invention has been applied to many image pairs, consisting of a variety of types of scenes and motions including pairs of images containing actual image motion, external camera motion, and internal camera motion, e.g., camera zoom. The process of the present invention registered the images with superior results over the prior art flow field- and gradient-based registration techniques. 
     Due to the robust nature of the present invention, it can be used to register imagery produced by sensors of varying types, i.e., sensors having differing modalities such as an image from a visible camera and an image from an infrared camera. In an alternative embodiment of the invention, the global correlation-based alignment method presented above is applied to registering multi-sensor images. The adaptation of the present invention to multi-sensor images is based on a number of observations regarding multi-sensor imagery. First, when the images f and g are imaged by sensors of different modalities, the relationship between the brightness values of corresponding pixels is complex and unknown. Contrast reversal may occur in some parts of the images, while not in others. Furthermore, visual features present in one sensor image may not appear in the other image, and vice versa. Moreover, multiple brightness values in f may map to a single brightness value in g, and vice versa. In other words, the two multi-sensor images are usually not globally correlated, and often not even statistically correlated. Locally, however, (i.e., within small image patches), it is realistic to assume that the two images are correlated, up to a sign reversal. 
     Since a multi-resolution (coarse-to-fine) search is used in the present invention to handle large misalignments between the two images, it is important that the local match measure m(u,v;x,y) be applicable at all resolution levels (i.e., applicable to low-pass or band-pass filtered images). When the two images are obtained by the same sensor, or by two different cameras of same modality, then typically the corresponding signals at all resolution levels are correlated. In a multi-sensor image pair, however, the signals are correlated primarily in high resolution levels, where the details that correspond to the physical structure in the scene are captured. Low-pass filtered multi-sensor images are usually not correlated, as low resolution features depend heavily on the photometric and physical imaging properties of the sensors, which can differ substantially between the two images. To apply coarse-to-fine search to multi-sensor images, it is therefore necessary to project high-resolution scene information into low-resolution levels of the pyramid. 
     Consequently, these observations require the foregoing method of FIG. 2 to be slightly altered to produce a multi-sensor registration process. FIG. 3 depicts a block diagram of the multi-sensor registration routine  300 . In routine  300 , the images are pre-processed prior to being subjected to the global alignment routine  116  of FIG.  2 . Specifically, the routine  300  begins at step  302  and proceeds to step  304  to accomplish image preprocessing. The images are pre-filtered to assure that when the pyramid is constructed, high-resolution information will in fact be projected into low resolution levels. In particular, in step  306 , each of the two multi-sensor images is first high-pass filtered using a Laplacian filter, then, at step  306 , the filtered images are squared. The squaring process assures that high resolution features appear in low resolution levels. Processing to achieve high resolution feature shifting is disclosed in P. J. Burt, “Smart Sensing With a Pyramid Vision Machine, Proceeding of the IEEE, 76:1006-1015 (1988). Since the Laplacian-energy image is invariant to contrast reversal in the original image, a local match measure based on normalized-correlation is adequate for this kind of imagery. 
     At step  310 , the global alignment routine  116  of FIG. 2 is applied with a normalized-correlation local match measure to the Laplacian-energy images (the preprocessed images) rather than to the original images f and g. 
     Further increase in the sensitivity of the local match measure can be achieved by using directional-derivative energy images instead of the Laplacian-energy images. The global correlation-based registration was extended to apply simultaneously to a collection of n multi-sensor pairs of directional-derivative-energy images. The extension is performed by estimating a single parametric transformation {right arrow over (p)}, which simultaneously maximizes the local normalized-correlations of n image pairs: 
     
       
           M ( {right arrow over (p)} )=Σ (x,y) Σ i=1   n   S   i   (x,y) ( u )( x,y;{right arrow over (p)} ), v ( x,y;{right arrow over (p)} ))  (8) 
       
     
     where S i   (x,y)  is the match-measure surface estimated for pixel (x,y) in the i-th image pair. In particular, this can be applied to n corresponding pairs of directional-derivative-energy images. 
     In practice, the invention has be used with four (n=4) directional derivatives (horizontal, vertical, and the two diagonals). Applying the global correlation-based alignment simultaneously to the 4 corresponding pairs provides the global alignment with increased robustness, which allows for successful multi-sensor registration even in very challenging situations, i.e., when there are significant differences in image content between the two sensor-images. Apart from having significantly different appearance, there can be many non-common features in the multi-sensor image pair, which can lead to false matches. These are overcome by the automatic outlier mechanism of the present invention. 
     Although various embodiments which incorporate the teachings of the present invention have been shown and described in detail herein, those skilled in the art can readily devise many other varied embodiments that still incorporate these teachings.