Patent Publication Number: US-RE43206-E

Title: Apparatus and method for providing panoramic images

Description:
FIELD OF THE INVENTION 
     The present invention relates to an apparatus, algorithm, and method for stitching different pieces of images of a scene into a panoramic environment map. 
     BACKGROUND 
     The most common way to electronically represent the real world is with image data. Unlike traditional graph-based systems, there are systems which use panoramic images to construct a virtual world. The major advantage of a system which uses panoramic images is that very vivid and photo-realistic rendering results can be obtained even when using PCs. In addition, the cost of constructing the virtual world is independent of scene complexity. In such systems, panoramic images are stitched together into a panoramic map from several individual images which are acquired by rotating a camera horizontally or vertically. This panoramic map can be used in different applications such as movie special effects, the creation of virtual reality, or games. A typical problem is how to stitch the different pieces of a scene into a larger picture or map. One approach to address this problem is to manually establish correspondences between images to solve unknown parameters of their relative transformation. Because manual methods are tedious for large applications, automatic schemes are preferably used for generating a seamless panoramic image from different pieces of images. 
     One proposed approach uses a nonlinear minimization algorithm for automatically stitching panoramic images by minimizing the discrepancy in intensities between images. This approach has the advantage of not requiring easily identifiable features. However, this technique does not guarantee finding the global minimum if the selection of starting points is not proper. Further because the optimization process is time-consuming, the approach is inefficient. In this invention, the domain of images under consideration is panoramic images. 
     SUMMARY OF THE INVENTION 
     The invention allows users to generate panoramic images from a sequence of images acquired by a camera rotated about its optical center. In general, the invention combines feature extraction, correlation, and relaxation techniques to get a number of reliable and robust matching pairs used to derive registration parameters. Based on the obtained registration parameters, different pieces of consecutive images can-be stitched together to obtain a seamless panoramic image. 
     In a first aspect, a method of merging a pair of images to form a seamless panoramic image includes the following steps. A set of feature points along the edges of the images is extracted, each feature point defining an edge orientation. A set of registration parameters is obtained by determining an initial set of feature points from a first one of the images which matches a set of feature points of a second one of the images. A seamless panoramic image is rendered using the first and second images with the set of registration parameters. 
     The invention provides a feature-based approach for automatically stitching panoramic images acquired by a rotated camera and obtaining a set of matching pairs from a set of feature points for registration. Since the feature points are extracted along the edges, each feature point specifies an edge orientation. Because the orientation difference between two panoramic images is relatively small, the difference between edge orientations of two feature points is also small if they are good matches. Based on this assumption, edge information of feature points can be used to eliminate in advance many false matches by checking their orientation difference. Moreover, many unnecessary calculations involving cross-correlation can be screened in advance, thereby significantly reducing the search time needed for obtaining correct matching pairs. After checking, by calculating the value of correlation of the remaining matching pairs, a set of possible matches can be selected with a predefined threshold. The set of possible matches are further verified through a relaxation scheme by calculating the quality of their matches. Once all of the correct matching pairs are found, they are then used to derive registration parameters. In this invention, an iterative scheme is applied to increase the reliability in providing matching results. Since only three iteration or fewer are needed and only a few feature points are involved in the matching pairs, the whole procedure can be accomplished very efficiently. Also, as discussed above, because the orientation difference of two feature points is checked in advance (before matching). Many calculations involving cross-correlation are not required and the efficiency of stitching is significantly improved. Compared with conventional algorithms, the proposed scheme offers improved efficiency and reliability for stitching images. 
     Embodiments of this aspect of the invention may include one or more of the following features. In one embodiment, a set of feature points are first extracted through wavelet formations. Among other advantages, the invention uses wavelet transforms to obtain a number of feature points with edge orientations. Such edge information can speed up the entire registration process by eliminating many impossible matches in advance and avoiding many unnecessary calculations of correlation. The method determines a number of reliable and robust matching pairs through relaxation. The method also measures the quality of a matching pair, imposes angle consistency constraint for improving the robustness of registration, and uses a voting concept to get the desired solution from the set of final matching results. 
     In other embodiments, the method forms the final panoramic image with the help of the registration results. In particular, the method adjusts the intensity differences between consecutive input images and blends the intensities of adjacent images to obtain a seamless panoramic image. The final panoramic images can then be used to build a virtual world. 
     Still other embodiments may include one or more of the following features: 
     For example, the method selects a number of feature points through wavelet transforms. Each feature point is associated with an edge orientation so that the speed of the registration process is increased. 
     The method uses an angle constraint to construct a set of matching pairs, which are used to obtain reliable matching results through a relaxation and a voting technique. The set of matching results are then used to form the final seamless panoramic image. 
     Constructing an initial set of matching pairs for registration includes comparing the edge orientation differences of feature points in one image and its corresponding feature points in another, calculating the values of correlation of each possible matching pair, and thresholding them with a predefined threshold. 
     Getting reliable matching results through relaxation and a voting technique includes calculating the quality of a matching pair, imposing angle consistency constraint to filter out impossible matching pairs, updating matching results through relaxation, and using the voting technique to obtain the reliable registration parameters. In addition, it refines the final registration results by using the correlation technique with a proper starting point. 
     Forming the final panoramic images includes dynamically adjusting and properly blending the intensity differences between adjacent images. 
     In another aspect, the invention features a system for merging pairs of images to form a panoramic image. The system includes an imaging device which, in operation, acquires a series of images, a storage for storing a series of images, a memory which stores computer code, and at least one processor which executes computer code to extract a set of feature points along the edges of the images, each feature point defining an edge orientation and to obtain a set of registration parameters by determining an initial set of feature points from a first one of the images which matches a set of feature points of a second one of the images, and to render a seamless panoramic image using the first and second images with the set of registration parameters. 
     Other advantages and features of the invention will become apparent from the following description, including the claims and the drawings. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  illustrates an approach for generating a virtual panoramic world. 
         FIG. 2  shows the geometry for cylindrical mapping. 
         FIG. 3  is a flowchart illustrating the steps used to extract feature points from images using wavelet transform. 
         FIG. 4  is a flowchart illustrating the steps used to obtain an initial set of matches using edge orientations and correlation values. 
         FIG. 5  is a flowchart illustrating the relaxation procedure used to obtain a set of reliable matching pairs. 
         FIG. 6  illustrates adjusting intensities of I a  and I b  such that their intensity differences are small. 
         FIG. 7  illustrates an example for a blending technique. 
         FIG. 8  shows an architecture for implementing the present invention. 
         FIG. 9  shows a series of input images for stitching. 
         FIG. 10  shows the stitching result obtained by the proposed stitcher. 
         FIG. 11  shows images used to demonstrate the ghost-like effect. 
         FIG. 12  shows how the ghostlike effect is removed: (a) with a wider blending width and (b) with a narrower blending width. 
     
    
    
     DESCRIPTION 
     Referring first to the flow diagram of  FIG. 1 , a method of the general steps of generating a seamless panoramic image from different pieces of an image is shown. 
     Input Image Warping 
     In general, it is difficult to seamlessly stitch two adjacent images together to form a panoramic image due to perspective distortion introduced by a camera. To remove the effects of this distortion, these images are preferably reprojected onto a simple geometry, e.g., a cube, a cylinder, or a sphere. In many applications, a cylindrical geometry is preferable since its associated geometrical transformation is simple. In this example, the cylindrical geometry is used.  FIG. 2  illustrates the relationship between the cylindrical surface ▭pqrs and film plane ▭PQRS of the camera. Note that the plane ▭PQRS is tangent to the cylindrical surface ▭pqrs. Let O denote the optical center, C the center of image plane, f the focal length and r the radius of the cylinder. Further, assume that P(x, y) is a pixel in the image plane and p(u, v) is its corresponding pixel in the cylindrical map. Using  FIG. 2 , the coordinates of p (u, v) can be obtained as follows: 
                     u   =       r∠C   ⁢           ⁢   OB     =       r   ⁢           ⁢     tan     -   1       ⁢       CB   _       OC   _         =     r   ⁢           ⁢     tan     -   1       ⁢     x   f             ,           (   1   )               
and
 
                   v   =       r   ⁢           ⁢       PB   _       OB   _         =     r   ⁢           ⁢       x         x   2     +     f   2           .                 (   2   )               
Moreover, since the radius r is equal to f, Equations (1) and (2) can be rewritten as follows:
 
u=f tan −1 x/f
 
and
 
                   v   =     f   ⁢           ⁢       x         x   2     +     f   2           .               (   4   )               
Based on Eqs. (3) and (4), input images are provided and then (step 10) warped into a cylindrical map for further registration to construct a complete panoramic scene (step 11).
 
     Feature Extraction Using Wavelet Transform 
     Referring again to  FIG. 1  and  FIG. 3 , once the input images have been warped, useful features required for registration are extracted (step 12). This process advantageously uses a wavelet transform for feature extraction (step 31). In particular, all of the feature points are detected along the edges. First, let S(x, y) be a 2-D smoothing function. Two wavelets, ψ 1 (x, y) and ψ 2 (x, y) are the partial derivatives of the smoothing function S(x,y) in the x and y directions, respectively, where 
                 ψ   1     ⁡     (     x   ,   y     )       =       ∂           ⁢     S   ⁡     (     x   ,   y     )           ∂   x             
and
 
                 ψ   2     ⁡     (     x   ,   y     )       =         ∂           ⁢     S   ⁡     (     x   ,   y     )           ∂   y       .           
Let
 
               ψ     2   j     1     =       1     4   j       ⁢       ψ   1     ⁡     (       x     2   j       ,     y     2   j         )               
and
 
                 ψ     2   j     2     =       1     4   j       ⁢       ψ   2     ⁡     (       x     2   j       ,     y     2   j         )           ,         
At each scale 2 j , the 2-D wavelet transform of a function I(x,y) in L 2 (R 2 ) can be decomposed into two independent directions as follows:
 
W 2     j     1 I(x, y)=I*ψ 2     j     1 (x, y) and W 2     j     2 I(x, y)=I*ψ 2     j     2 (x, y).
 
Basically, these two components are equivalent to the gradients of I(x,y) smoothed by S(x,y) at scale 2 j  in the x and y directions. At a specific scale s=2 i , the modulus of the gradient vector of f(x,y) can be calculated as follows:
 
M 2     j   I(x, y)=√{square root over (|W 1   2     j   I(x, y)| 2 +|W 2   2     j   I(x,y)| 2 )}{square root over (|W 1   2     j   I(x, y)| 2 +|W 2   2     j   I(x,y)| 2 )}.
 
If the local maximum of M 2     j   I(x, y) are located and thresholded with a preset value, then the edge points of I(x,y) at scale 2 j  can be detected. Since we are interested in specific feature points for scene registration, additional constraints have to be introduced. In general, noise is the main cause of false detection of edge points. In order to suppress the effect of noise, a criterion called edge correlation is introduced (step  32 ):
 
                   R   n     ⁡     (     j   ,   x   ,   y     )       =         ∏     n   -   1         i   =   0       ⁢       M     2     j   +   1         ⁢     I   ⁡     (     x   ,   y     )             ,         
where n is a positive integer indicating the number of scales involved in the multiplication, and j represents the initial scale for edge correlation. R n  reveals a peak whenever a true edge exists and is suppressed by the multiplication process if a point at location (x,y) is not a true edge. Thus, using the relationship for R n (j, x, y) above, the noise in an image can be suppressed while the true edges can be retained. In one embodiment of the invention, the number of scales for multiplication is chosen to be 2. In order to conserve the energy level, R n (j, x, y) should be normalized as follows:
 
                     R   _     n     ⁡     (     j   ,   x   ,   y     )       =         R   n     ⁡     (     j   ,   x   ,   y     )       ⁢         MP   ⁡     (   j   )           RP   n     ⁡     (   j   )               ,         
where
 
               MP   ⁡     (   j   )       =       ∑     x   ,   y       ⁢              M     2   j       ⁢           ⁢     I   ⁡     (     x   ,   y     )              2             
and
 
                 RP   n     ⁡     (   j   )       =       ∑     x   ,   y       ⁢                R   n     ⁡     (     j   ,   x   ,   y     )            2     .             
During the feature point selection process, an edge point is recognized as a candidate if its corresponding normalized edge correlation R 2 (j, x, y) is larger than its corresponding modulus value. Basically, the above mentioned process is equivalent to detecting an edge point with the strongest edge response in a local area. The three conditions which will be used to judge whether a point P(x,y) is a feature point or not are as follows:
         Condition 1: P(x,y) must be an edge point of the image I(x,y). This means that P(x,y) is a local maxima of M 2     j   I(x, y), and M 2     j   I(x, y)&gt;a threshold (step  34 ),   Condition 2:  R   2 (1, x, y)&gt;M 2     j   I(x, y) (step  36 ),   Condition 3:       

                   M     2   j       ⁢     I   ⁡     (     x   ,   y     )         =       max       (     x   ,   y     )     ∈     N   p         ⁢     {       M     2   j       ⁢     I   ⁡     (     x   ,   y     )         }         ,         
where N p  is the neighborhood of P(x,y) (step  34 ).
 
     Obtaining Initial Matches Using Edge Orientations and Cross Correlation 
     Referring to  FIGS. 3 and 4 , as part of step  12  of  FIG. 1 , a set of feature points is extracted using wavelet transform (step  31 ). Let FP I     a   ={p i =(p x   i , p y   i )} and FP I     b   ={q i =(q x   i , q y   i )} represent two sets of feature points extracted from two partially overlapping images I a  and I b  respectively. Assume that N I     a    and N I     b    represent the number of elements in the set of points FP I     a    and FP I     b   , respectively. Let. A(u) be the angle of an edge point u. A standard method for estimating the orientation of an edge-based feature point u at scale 2 j  can be expressed as follows:
 
Arg(W 2     j     2 I(x, y)+W 2     j     2 I(x, y)).
 
However, the above representation can be very sensitive to noise. Therefore, an edge tracking technique plus a line-fitting model is used to solve the noise problem (step  41 ).
 
     Let P be a feature point and N e  be its neighborhood. Since P is an edge point, there should exist an edge line passing through it. By considering P as a bridge point, an edge line passing through P can be determined by searching in all the directions from P. All the edge points on the edge are then used as candidates for determining the orientation of the edge line. During the searching process, the edge connection constraint and the direction consistency constraint are used to restrict the searching domain. The edge connection constraint means that if N e  contains another edge l but l does not pass P, all edge points in l will not be included in estimating the orientation of P. In certain cases, however, there will exist more than one edge line passing through P. In these cases, the first line detected is adopted to estimate the orientation. Let l 1  denote this line. The direction consistency constraint means all the edge points along other edge lines whose orientations are inconsistent with l 1  are not included to estimate the orientation of P. In this way, a set of edge points can be selected and then used to estimate the orientation of P using a line-fitting model. In other embodiments, other edge tracking technique can also be applied to provide a better estimation. 
     After the orientation is estimated, all feature points u will associate with an edge orientation A(u). For a feature point p i  in the set of points FP I     a    and q j  in the set of points FP I     b   , the orientation difference between them is calculated as follows (step  42 ):
 
θ ij =A(q j )−A(p i ).
 
In fact, if p i  and q i  provide a good match, the value of θ ij  will be small since the orientation of image I a  is similar to that of I b . Assuming this is the case, step  14  ( FIG. 1 ) is performed to eliminate impossible false matches in advance, avoiding many unnecessary correlation calculations.
 
     For a feature point p i  in FP I     a    and a feature point q j  in FP I     b   , if they form a good match (step  44 ), the following condition will be satisfied:
 
|A(p i )−A(q i )|&lt;10°.   Condition 1:
 
Adding this criterion will significantly speed up the search time. On the other hand, if p i  and q i  form a good match, the similarity degree between p i  and q i  should be larger. A cross-correlation which can be used to measure the similarity degree between p i  and q j  (step  46 ) and is defined as follows (step  32 ):
 
     
       
         
           
             
               
                 
                   
                     
                       C 
                       ⁡ 
                       
                         ( 
                         
                           
                             p 
                             i 
                           
                           ; 
                           
                             q 
                             j 
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                         1 
                         
                           
                             σ 
                             i 
                           
                           ⁢ 
                           
                             
                               
                                 σ 
                                 j 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     2 
                                     ⁢ 
                                     M 
                                   
                                   + 
                                   1 
                                 
                                 ) 
                               
                             
                             2 
                           
                         
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           
                             ∑ 
                             
                               x 
                               , 
                               
                                 y 
                                 = 
                                 M 
                               
                             
                           
                           
                             x 
                             , 
                             
                               y 
                               = 
                               
                                 - 
                                 M 
                               
                             
                           
                         
                         ⁢ 
                         
                           
                             [ 
                             
                               
                                 
                                   I 
                                   a 
                                 
                                 ( 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   
                                     x 
                                     + 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       p 
                                       x 
                                       i 
                                     
                                   
                                   , 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     y 
                                     + 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       p 
                                       y 
                                       i 
                                     
                                   
                                 
                                 ) 
                               
                               - 
                               
                                   
                               
                               ⁢ 
                               
                                 u 
                                 i 
                               
                             
                             ] 
                           
                           [ 
                           
                               
                           
                           ⁢ 
                           
                             
                               
                                 I 
                                 b 
                               
                               ( 
                               
                                   
                               
                               ⁢ 
                               
                                 
                                   x 
                                   + 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     q 
                                     x 
                                     i 
                                   
                                 
                                 , 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   y 
                                   + 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     q 
                                     y 
                                     i 
                                   
                                 
                               
                               ⁢ 
                               
                                   
                               
                               ) 
                             
                             ⁢ 
                             
                                 
                             
                             - 
                             
                                 
                             
                             ⁢ 
                             
                               u 
                               j 
                             
                           
                           ⁢ 
                           
                               
                           
                           ] 
                         
                       
                     
                   
                   ⁢ 
                   
                       
                   
                   , 
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     where u i  and u i u j  are the local means of p i  and q j , respectively, σ i  and σ i  are the local variances of p i  and q j , respectively; and (2M+1) 2  represents the area of matching window. Based on this correlation measure, a pair {p i     q   j } is qualified as a possible matching pair if the following conditions are satisfied (step  33 ): 
                   Condition  2:       ⁢       C       I   a     ⁢     I   b         ⁡     (       p   i     ;     q   j       )         =       Max       q   k     ∈     FP     I   b           ⁢     C   ⁡     (       p   i     ;     q   k       )           ,     
     ⁢         Condition  3:       ⁢       C       I   a     ⁢     I   b         ⁡     (       p   i     ;     q   j       )         =       Max       p   k     ∈     FP     I   a           ⁢     C   ⁡     (       p   k     ;     q   i       )           ,         
and
 
C i     a     I     b   (p i ;q j )&gt;T e , where T c =0.65 (step  35 ).   Condition 4:
 
Condition 2 means that given a feature point p i  it is desired to find a point q i εFP I     b    such that a the value of C I     a     I     b   (p i ;q j ) is maximized for all points q k εFP I     b   . Condition 3 means that given a feature points q i , it is desired to find a point p i εFP I     a    such that the value of C I     a     I     b   (p i ;q j ) is maximized. If only Condition 2 is used, it is possible that several points p i  match with single point q j . Conversely, if only Condition 3 is used, several points q j  possibly will match with single point p i . As for Condition 3, it forces the value of C I     a     I     b    of a matching pair to be larger than a threshold. In a preferred implementation, the orientation constraint will be checked first. If the constraint is not satisfied, it is not necessary to check Conditions 2, 3, and 4. In this way, only a few pairs are needed to calculate the cross-correlation measure C I     a     I     b   , which is considered a time bottleneck of the whole process.
 
     Eliminating False Matches through Relaxation 
     Referring to  FIGS. 1 and 5 , a set of reliable matching pairs is obtained through relaxation. Let MP I     a     I     b   ={p i   q i } i=1, 2 . . .  be the set of matching pairs which satisfy Conditions 1, 2, 3, and 4 above, where p i  is a point in image I a (x, y) and q i  is another point in image I b (x, y). Let Ne I     a   (p k ) and Ne I     b   (q k ) denote the neighbors of p k  and q k  within an area of radius R, respectively. Assume that NP p     i     q     j   ={n k   1   n k   2 } k=1, 2 . . .  is the set of matching pairs, where n k   1 εNe I     a   (p i ), n k   2 εNe I     b   (q j ), and all elements of NP p     i     q     j    belong to MP I     a     I     b   . The proposed method is based on a concept that if {p i   q i } and {p j   q j } provide a pair of good matches, the distance between p i  and p j  should be similar to the one between q i  and q j . Based on this assumption, we can measure the quality of a matching pair {p i   q i } according to how many matches {n j   1   n j   2 } in NP P     i     q   j  whose distance d(p i , n j   1 ) is similar to the distance d(q i , n j   2 ), where d(u i , u j )=|u i −u j ∥, the Euclidean distance between two points u i  and u j . With this concept, the measure of the quality of a match {p i   q i } is defined as follows: 
                         G       I   a     ⁢     I   b         ⁡     (   i   )       =       ∑       {       n   k   i     ⇔     n   k   2       }     ∈     NP   Pi4i         ⁢       r   ⁡     (     i   ,   k     )         1   +     dist   ⁡     (     i   ,   k     )               ,           (   6   )               
where
 
                 dist   ⁡     (     i   ,   k     )       =       [       d   ⁡     (       p   i     ,     n   k   1       )       +     d   ⁡     (       q   i     ,     n   k   2       )         ]     /   2       ,       r   ⁡     (     i   ,   k     )       =     {           ⅇ       -     u   ⁡     (     i   ,   k     )         /     T   1                 if   ⁢           ⁢     μ   ⁡     (     i   ,   k     )         &lt;     T   2                 0   ,         otherwise                   
with the two predefined thresholds
 
                 T   1     ⁢           ⁢   and   ⁢           ⁢     T   2       ,       and   ⁢           ⁢     u   ⁡     (     i   ,   k     )         =                    d   ⁢       (p     i       ,     r   k   1       )     -     d   ⁡     (       q   i     ,     r   k   2       )                dist   ⁡     (     i   ,   k     )         .             
The contribution of a pair {n k   1   n k   2 } in NP p     i     q     i    monotonically decreases based on the value of dist(i,k). Further, if the value of u(i,k) is larger than the threshold T 2 , the contribution of {n k   1   n k   2 } is set to zero.
 
     Referring to  FIG. 5 , after the quality of match of each pair {p i   q i } in MP I     a     ,I     b    (step  51 ), relative quality value of each pair is obtained to G I     a     I     b    (i) to further eliminate false matches (step  51 ). Now, based on the quality value of each candidate match, a relaxation technique is used to eliminate false candidates for further registration. If we define the energy function as follow: 
             F   =       ∑              p   i     ⇔     q   i            |     ∈     MP       I   a     ⁢     I   b               ⁢       G       I   a     ⁢       I   b     ⁡     (   i   )           .             
then the relaxation procedure can be formulated as follow:
 
                                         Iterate {             -Compute the quality for each candidate match             -Choose the best possible candidates for minimizing F             according to the quality value G I     a     I     b    (i).           } until F converges.                    
There are several strategies for updating the matching candidates. In one application an update strategy, referred to here as “some-looser-take-nothing” is used to update the matching candidates (step  52 ). First, according to the quality value of G I     a     I     b    (i), we sort each element of MP I     a     ,I     b    is sorted in increasing order (step  52 ). Then, a predetermined percentage Q% of matches in are eliminated as the impossible matches. Thus, the remaining (100−Q%) of matches are selected as potential matches for further relaxation (step  54 ). In our implementation, Q is set to 25. Three iterations are generally sufficient for achieving relaxation (step  53 ).
 
     On the other hand, in order to make the matching results more reliable, the method includes an “angle consistency” constraint (step  53 ) within the first iteration of the relaxation process to further eliminate impossible matching pairs. That is, if {p i   q i } and {p j   q j } are well matched, the angle between {right arrow over (p i p j )} and {right arrow over (q i q j )} must be close to zero. In this case, during this first iteration, for each element {p i   q i } in MP I     a     ,I     b   , a counter CA i  is used to record the number of matches {p k   q k } in MP I     a     ,I     b    where the angle between {right arrow over (p i p k )} and {right arrow over (q i q k )} is less than a predefined threshold θ. According to the value of CA i , the elements of MP I     a     ,I     b    are sorted in increasing order. The first Q% of matches in MP I     a     ,I     b    are considered to be impossible matches. 
     Obtaining the Desired Offset from Matches 
     After applying the relaxation process, a set of reliable matches is obtained. Referring to  FIG. 1 , the method uses a “voting” concept to derive a desired offset from the set of reliable matching pairs. For example, assume that this set is CP I     a     ,I     b   ={u i   v i } i=1, 2 . . . Ne , where N e  is the total number of elements in CP I     a     ,I     b   . In general, the 2-D point sets {u i } and {v i } satisfy the following relation:
 
v i =u i +T, for i=1,2,3 . . . , N e ,
 
where T is the desired solution. However, in real cases, different pairs {u i   v i } will lead to different offsets T i . Therefore, a voting technique is used to measure the quality of different solutions T i . Let S(i) denote a counter which records the number of solutions T k  consistent with T i . Two solutions T i  and T k  are said to be consistent if the distance between T i  and T k  is less than a predefined threshold. Since there are N e  elements in CP I     a     ,I     b   , the total number of consistency tests will be N e (N e −1)/2. After applying the consistency test, the offset  T  associated with the maximum value of S(i) is chosen as the desired solution.
 
     Note that due to noise and image quality, the positions of feature points will not be precisely located and the accuracy of T is affected. In practical implementation, the y-component  T y    of  T  is generally affected more then  T x   . Therefore, the following method is proposed to refine and correct the solution of  T y    for the y-component. Let {u v} be the matching pair of CP I     a     ,I     b    having the highest quality value. Given point u and the offset  T , the refined offset can be found by searching the local neighborhood of the point (u+  T ) in another image I h  using the correlation technique. However, if there is little texture information within the local neighborhood of the point u, the above approach will not necessarily provide a satisfactory solution. Let 
                 g   x     ⁡     (   k   )       =         ∑   5       i   =   1       ⁢              I   a     ⁡     (         u   x     +   i     ,   k     )       -       I   b     ⁡     (         u   x     -   i     ,   k     )                      
be the horizontal gradient of the point (u x ,k) in I h , where u x  is the x-coordinate of u. Instead of using the starting point u directly, we use another starting point ū for refining the desired offset  T  by searching the point ū whose a horizontal gradient g x  is largest along a column of pixels with the same x-coordinate u x . Based on the starting point ū, the final offset can then be accurately obtained with the correlation technique.
 
     Rendering the Final Panoramic Image 
     In general, when stitching two adjacent images, discontinuities of intensity exist between their common areas. Therefore, step  17  of  FIG. 1  is used to eliminate such intensity discontinuities. The scheme used in step  17  can be divided into two stages. The first stage is used to adjust the intensities of two adjacent images such that their intensities are similar. The second stage is used to blend their image intensities according to a distance measure such that the final composite image appears smooth. Assume I a  and I b  are two adjacent images with their widths in w a  and w b , respectively. Let ΔI be the average intensity difference between the overlapping area 
     of I a  and I b , that is, 
                     ΔI   =       1        A          ⁢       ∑     i   ∈   A       ⁢     (         I   b     ⁡     (     q   ⁡     (   i   )       )       -       I   a     ⁡     (     p   ⁡     (   i   )       )         )           ,           (   7   )               
where A is the overlapping area of I a  and I b , |A| is the number of pixels in A, p(i) is a pixel in I a , and q(i) is its corresponding pixel in I b .
 
     In particular, referring to  FIGS. 6A and 6B , the gap of average intensity between I a  and I b  is about ΔI. According to ΔI, w a , and w b , the first stage is used to adjust the intensities of I a  and I b  as follows: 
                         I   a     ⁡     (     p   ⁡     (     x   ,   y     )       )       =         I   a     ⁡     (     p   ⁡     (     x   ,   y     )       )       +     x   ⁢           ⁢       Δ   ⁢           ⁢   I       2   ⁢           ⁢     w   1               ,           (   8   )               
and
 
     
       
         
           
             
               
                 
                   
                     
                       I 
                       b 
                     
                     ⁡ 
                     
                       ( 
                       
                         q 
                         ⁡ 
                         
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         I 
                         b 
                       
                       ⁡ 
                       
                         ( 
                         
                           q 
                           ⁡ 
                           
                             ( 
                             
                               x 
                               , 
                               y 
                             
                             ) 
                           
                         
                         ) 
                       
                     
                     + 
                     
                       
                         ( 
                         
                           x 
                           - 
                           
                             w 
                             2 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           ΔI 
                           
                             2 
                             ⁢ 
                             
                               w 
                               2 
                             
                           
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     After this adjusting step, the intensities of I a  and I b  in  FIG. 6(a)  will be gradually changed to approach the intensity line  EF  shown in  FIG. 6(b) , thereby bringing the intensities between I a  and I b  closer. As will be described below, in order to further smooth the intensity discontinuity between I a  and I b , a blending technique is then applied. 
     The second stage uses a ray-casting method to blend different pixel intensities together. Referring to  FIG. 7 , p i  is a pixel in I a , q i  is its corresponding pixel in I h , and I a  and I h  are two boundary lines in I a  and I h , respectively. With p i  and q i , the intensity of the corresponding pixel r i  in the composite image I can be obtained as follows: 
                     I   ⁡     (     r   i     )       =           d   b   t     ⁢       I   a     ⁡     (     p   i     )         +       d   a   t     ⁢       I   b     ⁡     (     q   i     )               d   a   t     +     d   b   t                 (   10   )               
Where d a  is the distance between p i  and I a , d h  the distance between q i  and I h , and t is an adjustable parameter. Using Equation (10), the intensities in I a  are gradually changed to approach the intensities of pixels in I h  such that the final composite image I looks very smooth. In fact, if the blending area is chosen too large, a “ghostlike” effect will occur, particularly when moving objects in the common overlapping area exist between I a  and I b . However, since the intensities of I a  and I b  have been adjusted, the blending width can be chosen small such that the so-called ghostlike effect is significantly reduced. In one preferred embodiment, the blending width is chosen as one-third of the original width of the overlapping area.
 
     Architecture for Implementation 
     Referring to  FIG. 8 , one preferred architecture for implementing the real-time stitcher apparatus, described above, includes input devices  60  (e.g., digital cameras or scanners) to acquire a series of panoramic images. Then, the panoramic images are stored into external storage  62  such as hard disks for further processing or being directly provided to one or more microprocessors  64  for stitching. The microprocessors  64  perform the stitching including warping, feature extraction, edge orientation estimation, correlation calculation, relaxation, rendering and blending the final panoramic images, etc. All the temporary data are stored in the system RAM memory to speed up stitching. Finally, the stitching result is sent to the display engine for displaying. In many applications this architecture can be implemented using a general personal computer. 
     Performance of the Invention 
     Referring to  FIG. 9 , to analyze the performance of the real-time stitcher apparatus, described above, a series of original panoramic images  9 a- 9 e captured by a rotated camera were provided. 
     Referring to  FIG. 10 , the stitched panoramic image obtained by our proposed method is shown. In addition, if moving objects in two adjacent images exist, the proposed method can significantly minimize ghost-like effects.  FIG. 11  shows two adjacent images with moving objects. The posture of the man in (a) is clearly different from the one in (b). Besides, large intensity differences exist between them. In order to smooth such large intensity differences, the blending width should be chosen to be large. However, with a large area, the ghost-like effect will appear in the final composite image (see  FIG. 12(a) ). Contrary to the previous technique, with our proposed blending technique, since the intensities of input images have been adjusted before blending, the blending width can be chosen smaller. Therefore, the ghostlike effect is significantly lessened, the quality of (b) is clearly improved over the quality shown in  FIG. 12(a) .