# Patent

## Title
Method and device for joining a plurality of individual digital images to form a total image

## Abstract

In a device and a corresponding method for joining a plurality of individual digital images to form a total image, a plurality of features is determined in a first individual image by means of a selection unit using a feature-based algorithm and then tracked in a second individual image by means of a tracking unit. A transformation matrix, with which the individual images are joined in an output unit to form the total image, is calculated from the determined feature correspondences in a transformation unit. The individual images can be joined in real time and with a high degree of accuracy by means of the feature-based algorithm in combination with a robust algorithm to calculate the transformation matrix.

## Background
CROSS REFERENCE TO RELATED APPLICATIONS
This application is a United States National Phase application of International Application PCT/EP2010/062251 and claims the benefit of priority under 35 U.S.C. §119 of German Patent Application Serial No. 10 2009 039 251.3 filed Aug. 28, 2009, the entire contents of which are incorporated herein by reference.
FIELD OF THE INVENTION
The invention relates to a method and a device for joining a plurality of individual digital images to form a total image.
BACKGROUND OF THE INVENTION
Endoscope systems are used to inspect cavities in the industrial and medical sector that are difficult to access. The restricted image field of endoscope systems makes the orientation and the coordination more difficult when inspecting the cavities, so the entire inspection process becomes significantly more complex and time-consuming. In order to extend the viewing field of endoscope systems, methods have been developed which produce a total image from a sequence of endoscopically recorded, individual digital images. The joining is also called stitching or mosaicing. The viewing field of the endoscope system is artificially extended by the joined total image, which is also called an image mosaic.
A real time-capable method for joining endoscopically recorded individual images to form an image mosaic, which is based on a calculation of the optical flow, is known from the specialist article “Real-Time Image Mosaic for Endoscopic Video Sequences” by Konen et al., published in “Bildverarbeitung für die Medizin”, 2007. A vector field is called an optical flow and gives the movement direction and speed between two respective individual images for each pixel. It is assumed here that the pixel values between the individual images do not change, but merely displace. The drawback in this method is that calculation of the optical flow is laborious.
SUMMARY OF THE INVENTION
The invention is based on an object of providing a method for joining a plurality of individual digit images to form a total image, which is fast, in particular real time-capable and accurate.
This object is achieved by a method for joining a plurality of individual digital images to form a total image, comprising the steps:

    
    
        providing a first individual digital image and a second individual digital image overlapping the latter,
        determining a plurality of features of the first individual image, wherein the features are each allocated first feature coordinates,
        determining the features in the second individual image, wherein the features are each allocated second feature coordinates,
        determining a transformation matrix having at least six degrees of freedom between the individual images, in that
        
            a first part quantity of the features is selected,
            a test transformation matrix is calculated from the first and second feature coordinates of this part quantity,
            a quality value is determined for the test transformation matrix with a second part quantity of the features, and
            the transformation matrix is determined proceeding from the test transformation matrix, when the quality value fulfils a quality criterion, and
        
        
        joining the individual images to form a total image by means of the transformation matrix.
    
    


The method according to the invention is based on a feature-based algorithm, with which a plurality of features are determined in the first individual image and these features are pna again determined or followed in the second individual image. The following of the features determined in the first individual image is also called tracking. The KLT-algorithm, the Harris-Corner detector, the monotonicity operator, the SIFT algorithm (Scale Invariant Feature Transformation) or the SURF algorithm (Speeded Up Robust Features) can be used, for example, as the feature-based algorithm.
A transformation matrix, which is used to join the individual images to form the total image, is calculated from the determined and tracked features. The transformation matrix has at least six degrees of freedom and is, for example, an affine transformation matrix with precisely six degrees of freedom or a projective transformation matrix with precisely eight degrees of freedom. The transformation matrix, which moves the two individual images into one another, is calculated from the quantity of feature correspondences determined, in other words, the first and associated second feature coordinates to the features M1 to Mn. The total quantity of feature correspondences is called

A={(p1a,p1b),(p2a,p2b), . . . ,(pna,pnb)}  (1)

pna are the first feature coordinates of the feature Mn and
pnb are the associated second feature coordinates of the feature Mn.
If the basic transformation is the perspective transformation, there applies

pb=Hb,apb  (2)

wherein

    Hb,a is a transformation matrix with eight degrees of freedom, which connects the points pa=(xa, ya, 1)T and pb=(xb, yb, 1)T to one another.


Writing out the elements of the matrix H produces




  
    
      
        
          (
          
            
              
                
                  x
                  a
                
              
            
            
              
                
                  y
                  a
                
              
            
            
              
                1
              
            
          
          )
        
        =
        
          
            (
            
              
                
                  
                    h
                    11
                  
                
                
                  
                    h
                    12
                  
                
                
                  
                    h
                    13
                  
                
              
              
                
                  
                    h
                    21
                  
                
                
                  
                    h
                    22
                  
                
                
                  
                    h
                    23
                  
                
              
              
                
                  
                    h
                    31
                  
                
                
                  
                    h
                    32
                  
                
                
                  
                    h
                    33
                  
                
              
            
            )
          
          ⁢
          
            (
            
              
                
                  
                    x
                    b
                  
                
              
              
                
                  
                    y
                    b
                  
                
              
              
                
                  1
                
              
            
            )
          
        
      
    
    
      
        (
        3
        )
      
    
  




and the coordinates are calculated as





  
    
      
        
          
            x
            a
          
          =
          
            
              
                
                  h
                  11
                
                ⁢
                
                  x
                  b
                
              
              +
              
                
                  h
                  12
                
                ⁢
                
                  y
                  b
                
              
              +
              
                h
                13
              
            
            
              
                
                  h
                  31
                
                ⁢
                
                  x
                  b
                
              
              +
              
                
                  h
                  32
                
                ⁢
                
                  y
                  b
                
              
              +
              
                h
                33
              
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        4
        )
      
    
  
  
    
      
        
          y
          a
        
        =
        
          
            
              
                h
                21
              
              ⁢
              
                x
                b
              
            
            +
            
              
                h
                22
              
              ⁢
              
                y
                b
              
            
            +
            
              h
              23
            
          
          
            
              
                h
                31
              
              ⁢
              
                x
                b
              
            
            +
            
              
                h
                32
              
              ⁢
              
                y
                b
              
            
            +
            
              h
              33
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




Simple transformation of these equations leads to

h11xb+h12yb+h13−h31xbxa−h32ybxa−h33xa=0  (6)

and

h21xb+h22yb+h23−h31xbya−h32ybya−h33ya=0  (7)

Each feature correspondence provides two homogeneous equations of this type in order to determine the nine parameters h11 . . . h33. However, the matrix Hb,a is only clearly defined to a desired scaling factor, which is why it has eight degrees of freedom. If the homogeneous equations are divided by the element h33, eight unknowns remain, as this element is cancelled. For this reason, a quantity of four feature correspondences is sufficient to determine the transformation matrix between two individual images.
The calculation of the transformation matrix can be carried out by means of a singular value partition. From the vectors

h=(h11,h12,h13,h21,h22,h23,h31,h32,h33)T  (8)

αi=(xib,yib,1,0,0,0,−xib,xia,−yibxia,−xia)T  (9)

βi=(0,0,0xib,yib,1,−xib,yia,−yibyia,−yia)T  (10)

the homogeneous equations can be written as

αiTh=0  (11)

βiTh=0.  (12)

From the n feature correspondences wherein n≧4, the system matrix Q for the singular value partition can now be constructed




  
    
      
        Q
        =
        
          (
          
            
              
                
                  α
                  11
                
              
              
                
                  α
                  12
                
              
              
                ...
              
              
                
                  α
                  19
                
              
            
            
              
                
                  α
                  21
                
              
              
                
                  α
                  22
                
              
              
                ...
              
              
                
                  α
                  29
                
              
            
            
              
                ⋮
              
              
                ⋮
              
              
                ⋮
              
              
                ⋮
              
            
            
              
                
                  α
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    1
                  
                
              
              
                
                  α
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    2
                  
                
              
              
                ...
              
              
                
                  α
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    9
                  
                
              
            
            
              
                
                  β
                  11
                
              
              
                
                  β
                  12
                
              
              
                ⋮
              
              
                
                  β
                  19
                
              
            
            
              
                
                  β
                  21
                
              
              
                
                  β
                  22
                
              
              
                ⋮
              
              
                
                  β
                  29
                
              
            
            
              
                ⋮
              
              
                ⋮
              
              
                ⋮
              
              
                ⋮
              
            
            
              
                
                  β
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    1
                  
                
              
              
                
                  β
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    2
                  
                
              
              
                ⋮
              
              
                
                  β
                  
                    n
                    ⁢
                    
                      
                    
                    ⁢
                    9
                  
                
              
            
          
          )
        
      
    
    
      
        (
        13
        )
      
    
  




and the homogeneous equations from equation (11) can be written as

Qh=0  (14)

The solution to this equation system is the zero space of Q. The solution is unambiguous, if four or more feature correspondences, which are not co-linear, are contained in the quantity. If more than four such correspondences exist, the solution with the smallest error squares can also be calculated with the aid of the singular value partition. The singular value partition factorizes Q in three matrices Q=UΣVT. The diagonal matrix Σ contains the decreasingly sorted singular values σ1 . . . σ9. In the case of a perfect solution there applies σ9=0 and σi≠0 for 1≦i≦8. In the case of an overdetermined equation system, there applies σ9≈0. In both cases, the corresponding last column vector of V is the solution vector. This vector contains nine elements, which produce the looked for transformation matrix if this is filled line-wise with the elements.
For the affine transformation model, the transformation matrix can be determined in a corresponding manner. As this transformation matrix only has six degrees of freedom, a quantity of at least three non-co-linear feature correspondences is sufficient here. There applies σ8=0 and σ9=0 and the corresponding right singular vectors are the eighth or ninth unit vectors. In the case of a perfect solution, there applies σ7=0 and in the case of a solution with the smallest error squares σ7≈0. The solution is the seventh column vector from V here.
As the feature-based algorithm, because of its simplicity, in particular in endoscopically recorded individual images, because of their poor image quality, leads to incorrect allocations in the feature correspondences, in other words to incorrectly tracked features, a robust calculation of the transformation matrix is necessary. For this purpose, the so-called RANSAC algorithm (RANdom SAmple Consensus) is used, with which incorrectly allocated feature correspondences remain disregarded in the calculation of the transformation matrix. Firstly, a test transformation matrix is calculated using a first part quantity, for example a minimum part quantity of required feature correspondences. A quality value, the so-called support, is now determined for this test transformation matrix. This means that it is determined whether the test transformation matrix for a second part quantity of the features, for example for the remaining features, satisfies a predetermined quality criterion. For example, it is determined what number of features of the second part quantity satisfies the quality criterion. The features that satisfy the quality criterion are designated so-called inliers; the other features are correspondingly designated outliers. For a test transformation matrix, a feature correspondence is, for example, classified as an inlier if there applies

∥pia−hpip∥≦t  (15)

wherein ∥ . . . ∥ designates the Euclidean distance.

This process is repeated until a test transformation matrix with an adequately high support is found or a maximum number of iterations has been run through. The final transformation matrix is then, for example, the best test transformation matrix. Alternatively, the transformation matrix can be calculated from the inliers present for the best test transformation matrix.
The feature-based algorithm therefore allows endoscopic image mosaics or an endoscopic total image of individual images to be set up, an outlier recognition for incorrectly followed or tracked features being provided by the RANSAC algorithm. The method according to the invention therefore ensures, due to the feature-based algorithm, a fast joining of the individual images and, because of the RANSAC algorithm, an accurate joining of the individual images, even if they only have a low image quality, as is the case in endoscopically recorded individual images. The method according to the invention is, in particular, real time-capable, real time-capability being taken to mean at least ten joined individual images per second when using a convention commercial PC.
A method, in which the feature coordinates are determined iteratively in different resolution stages of the individual images, allows a determination of the feature coordinates in the case of relatively large displacements between the individual images. Proceeding from a lowest resolution stage, a feature is localized and tracked, for example with sub-pixel accuracy. The result is projected into the next resolution stage, where the feature is again determined and tracked. The displacements in pixels are therefore limited to a certain amount per resolution stage.
A method, in which the features of the first individual image are determined, in that in each case a pixel window is formed from a plurality of adjacent pixels, a matrix is formed from the gradients for these pixels, the eigenvalues are calculated from the matrix, the pixel window is characterized as the feature, if the eigenvalues exceed a threshold value, and the feature is allocated first feature coordinates, and in which the features in the second individual image are determined, in that in each case an error vector is formed for the feature as a function of the gradients and the differences in the pixel values of the individual images in the pixel window, a displacement vector is calculated from the error vector and the matrix, and the feature is allocated second feature coordinates, provides a good compromise between rapidity and accuracy. In this feature-based algorithm, which is also called the KLT algorithm, the determination and tracking of the features are optimally matched to one another, as a pixel window or feature window is only characterized as a feature if the linear equation system to be solved for the tracking of the feature is stably solvable. For this purpose, a matrix is formed from the gradients to the pixels of the pixel window in the first individual image and the eigenvalues thereof calculated. The matrix has the rank 2. If the two eigenvalues exceed a predefined threshold, the pixel window is characterized as a feature, to which first feature coordinates are allocated. As the features are also called feature points, the feature coordinates are also correspondingly called point coordinates. These feature coordinates may be the pixel coordinates of one of the pixels located in the pixel window. The matrix is then used to track or determine the feature in the second individual image. For this purpose, a displacement vector is calculated from the matrix and an error vector by solving a linear equation system. The displacement vector describes the displacement of the pixel window or feature in the second individual image and applies to all the pixels of the pixel window.
The feature-based algorithm is described in detail in the specialist article “Detection and Tracking of Point Features” by Carlo Tomasi and Takeo Kanade (Technical Report CMU-CS-91-132, Carnegie Mellon University, 1991), to which reference is made here and the entire contents of which are incorporated by reference. The algorithm is also known as the KLT (Kanade-Lukas-Tomasi) tracking algorithm or as the KLT tracker or KLT algorithm. The KLT algorithm is based on the linear equation system

G·d=e  (16)

wherein

    G is a symmetrical matrix with rank 2,
    d is a displacement vector (dx, dy) and
    e is an error vector.


For two respective consecutive individual images I (hereinafter also Ia) and J (hereinafter also Ib), the matrix G can be calculated by determining the gradient g from the first individual image as follows:




  
    
      
        G
        =
        
          
            ∫
            W
          
          ⁢
          
            
              gg
              T
            
            ⁢
            w
            ⁢
            
              d
              A
            
          
        
      
    
    
      
        (
        17
        )
      
    
  




wherein

    g is a gradient vector to a pixel of the pixel window W,
    w is a weighting function for weighting the pixels of the pixel window W and
    dA=dx·dy is the integration differential of the double integral.


The two-dimensional error vector e is calculated from the calculated gradients g and the differences of the pixel values of the individual images in the pixel window W, as follows:




  
    
      
        e
        =
        
          
            ∫
            W
          
          ⁢
          
            
              
                (
                
                  
                    I
                    ⁡
                    
                      (
                      x
                      )
                    
                  
                  -
                  
                    J
                    ⁡
                    
                      (
                      x
                      )
                    
                  
                
                )
              
              ·
              gw
            
            ⁢
            
              d
              A
            
          
        
      
    
    
      
        (
        18
        )
      
    
  




wherein

    I(x) is the pixel value of the first individual image in the pixel x,
    J(x) is the pixel value of the second individual image in the pixel x,
    g is the gradient vector to a pixel of the pixel window W,
    w is a weighting function for weighting the pixels of the pixel window W and
    dA=dx·dy is the integration differential of the double integral.


The double integrals in the equations (16) to (18) are to be replaced by double sums in the x- and y-direction of the image plane for discrete pixels.
When calculating the eigenvalues, two small eigenvalues indicate homogeneous image areas, while a large and a small eigenvalue stand for unidirectional features, such as, for example, edges. Two large eigenvalue indicate corners or high-frequency textures, which are features that are good to track. As the larger eigenvalue is upwardly limited by the maximum pixel values or intensity values, it is sufficient for the smaller of the two eigenvalues to be above a specific threshold. Consequently, during the feature selection, a pixel window is accepted if for the eigenvalues λ1 and λ2 thereof, there applies that

min(λ1,λ2)>λt  (19)

wherein

    λt is a predefined threshold.


A method, in which the calculation of the respective displacement vector takes place iteratively, in that the individual images are displaced relative to one another by a first displacement vector, and in that a second displacement vector is calculated by means of the displaced individual images, increases the accuracy when calculating the displacement vectors. In the solution of the linear equation system (16), the displacement vector d(i) will have a residual error, as the Taylor approximation through the linear term, in other words through the gradient vector g, is only a precise solution when the pixel values linearly depend on x and y. This is not the case above all in pixel areas with high frequencies, which are well suited for following or tracking. The approximation solution for the displacement vector d(i) may, however, be used as the starting point for further iteration steps, in which the equation system (1) is again solved and a new displacement vector d(i+1) is calculated. In order to check whether this is necessary, the respective pixel window is firstly displaced by the specific displacement vector in one of the individual images. As this position generally lies between discrete pixel coordinates, new pixel values for the pixel window are calculated by means of bilinear interpolation and the error, in other words the sum of the quadratic differences of the pixel values or image intensities, is then determined therefrom. If this error exceeds a predefined tolerance value, further iteration steps are carried out, proceeding from the current displacement vector d, until the error falls below the tolerance value or a maximum number of iteration steps has been reached. Typically, less than five iteration steps are required to precisely localize the feature looked for. As, after each iteration step, the pixel values or intensity values are redetermined by linear interpolation, the displacement vector d can be calculated with sub-pixel accuracy.
A method, in which the calculation of the respective displacement vector takes place in different resolution stages of the individual images, in that a first displacement vector is calculated in a first, low resolution stage of the individual images, in that the individual images are displaced relative to one another in a second, higher resolution stage by the first displacement vector, and in that a second displacement vector is calculated by means of the displaced individual images in the second resolution stage, allows a calculation of the displacement vectors in the case of relatively large displacements between the individual images. The feature-based algorithm, because of the linear approach on which it is based, is only in a position to determine comparatively small displacements. Larger displacements can also be determined by a pyramid-based approach, in which the individual images are present in a plurality of resolution stages j. Proceeding from a lowest resolution stage, a feature is localized, for example with sub-pixel accuracy, and the position calculated by the displacement vector d(j) is projected to the next resolution stage, where it is used as the starting point for the further calculation of the displacement vector d(j+1). Thus, the algorithm approaches an exact solution for a displacement vector d with each resolution stage, wherein the displacements to be determined per resolution stage are adequately small. The algorithm can, for example, be controlled with a predeterminable search radius r, which states how large the displacements may be, which are to be recognized by the algorithm. Using the search radius, the algorithm automatically calculates the number of resolution stages and the factor of sub-samplings between two resolution stages. At the same time, a limit range is defined in the x- and y-direction, which is ignored by the algorithm in the further calculation. This is because edge regions are produced, which are undefined, by the Gaussian filter, which is applied to every resolution stage.
A method, in which a pixel window is only characterized as feature if a minimum distance from the already determined features is exceeded, ensures that the features to be tracked are distributed as homogeneously as possible over the individual images. As a result, displacements in the most varied directions can easily be tracked.
A method, in which the transformation matrix has precisely eight degrees of freedom and characterizes a projective transformation, ensures a high image quality of the total image. The projective transformation is also called a perspective transformation or homography. By means of the projective transformation, the free movement of an imaging sensor over a flat object can be modeled, or any scene, in which the movement of the imaging sensor is limited to the rotation about its own axis.
A method, in which the feature coordinates are normalized before the calculation of the transformation matrix, increases the numerical stability in the calculation of the transformation matrix.
A method, in which the normalization takes place in such a way that the feature coordinates belonging to an individual image have an average distance of √{square root over (2)} from a focal point, ensures a high numerical stability in connection with a high image quality of the total image. Firstly, the focal point is in each case calculated from the features of the individual images j. The focal point C=(Cx,Cy) is calculated as




  
    
      
        c
        =
        
          
            1
            n
          
          ⁢
          
            
              ∑
              
                i
                =
                1
              
              n
            
            ⁢
            
              p
              i
              j
            
          
        
      
    
    
      
        (
        20
        )
      
    
  




Furthermore, for each image j, the mean distance d from the focal point can be calculated:




  
    
      
        d
        =
        
          
            1
            n
          
          ⁢
          
            
              ∑
              
                i
                =
                1
              
              n
            
            ⁢
            
              
              
                
                  p
                  i
                  j
                
                -
                c
              
              
            
          
        
      
    
    
      
        (
        21
        )
      
    
  




The aim is to displace the coordinate origin of the calculated features from each individual image into the focal point and to scale the features in such a way that the distance from the origin is on average √2. This takes place by multiplication by the normalization matrices Nj, which are calculated as follows:




  
    
      
        
          N
          j
        
        =
        
          (
          
            
              
                
                  
                    2
                  
                  d
                
              
              
                0
              
              
                
                  
                    -
                    
                      
                        2
                        ⁢
                        
                          c
                          x
                        
                      
                    
                  
                  d
                
              
            
            
              
                0
              
              
                
                  
                    2
                  
                  d
                
              
              
                
                  
                    -
                    
                      
                        2
                        ⁢
                        
                          c
                          y
                        
                      
                    
                  
                  d
                
              
            
            
              
                0
              
              
                0
              
              
                1
              
            
            
              
                
                  
                
              
              
                
                  
                
              
              
                
                  
                
              
            
          
          )
        
      
    
    
      
        (
        22
        )
      
    
  




Using the normalized feature coordinates pij, the normalized transformation matrix H is now calculated so that there applies

 pia≅ H pib  (23)

The non-normalized transformation matrix H can now be calculated from H. There applies:

H=Na−1 HNb.  (24)

A method, in which the calculation of the transformation matrix takes place in such a way that for a repetition count a first part quantity of the features is in each case selected, a test transformation matrix belonging to this part quantity is calculated, the number of features of a second part quantity is determined for the respective test transformation matrix as the quality value and these features are imaged by the test transformation matrix with a defined accuracy, the test transformation matrix with the highest number is selected, and the transformation matrix is calculated from the associated features, which is imageable with the defined accuracy, ensures a robust and accurate calculation of the transformation matrix. An optimal accuracy of the transformation matrix is achieved by the concluding calculation of the transformation matrix from all the inliers of the best test transformation matrix.
A method, in which the individual images are joined in a planar projection face, ensures a simple and efficient joining of the individual images to form the total image, which is significant, in particular for the real time capability.
A method, in which the individual image which is later in terms of time during joining is superimposed on the individual image that is earlier in terms of time, ensures simple visualization of the total image in real time. The total image is a hybrid view of a static and a dynamic fraction. The static fraction is individual images going back in time, which are superimposed by the current individual image, as a dynamic fraction. The current individual image is preferably in each case displayed in the monitor center of a monitor. The image mosaic is therefore always centered on the most current view.
A method, in which, for joining, the x-coordinate and the y-coordinate of the pixels of the displaced individual image are rounded to the nearest integral value, ensures a rapid and calculation-optimized joining of the individual images. This is significant, in particular, for joining in real time.
A method, in which the individual images are recorded by means of an endoscope system, ensures a rapid and precise joining to form a total image, even if the individual images have poor image quality. Individual images recorded by an endoscope system are generally distorted by the wide angle lenses used. In addition, the image quality is impaired by the movement of the imaging sensor and by highlight effects.
A method, in which the individual images are rectified before the determination of the features, improves the image quality of the total image. Imaging systems with lenses, in particular wide angle lenses, cause distortions in the individual images. This is also the case in endoscope systems. The rectification takes place by means of an inverse distortion model. With regard to the rectification, reference is made to the specialist article “Precise Radial Un-distortion of Images” by John Mallon and Paul F. Whelan, published in “Proceedings of 17th International Conference on Pattern Recognition”, 2004, the entire contents of which are incorporated herein by reference.
A method, in which a highlight mask is set up for each individual image, in that pixel values are detected, which exceed a highlight threshold value and no features are permitted for the calculation of the transformation matrix in the highlight masks, avoids the determination and tracking of incorrect features. Endoscopically recorded individual images regularly have a large fraction of highlight as a result of the smooth and/or damp surfaces to be inspected. Highlight occurs where the light emitted by the endoscope system impinges perpendicularly on the object surface and is reflected from there. As these highlights have high pixel values or intensity values, the strongest contrasts can often be found at the limits of the highlights, so that the feature-based algorithm finds the supposedly best features here. However, these features do not belong to the scene, but are produced during the recording of the individual images and are therefore accordingly dependent on the movement of the imaging sensor with regard to their movement. As the highlight only changes slightly from individual image to individual image, the associated features can be tracked well, the RANSAC algorithm also not being in a position to eliminate them entirely. This ultimately leads to a poor transformation matrix. As highlights clearly stand out with regard to their pixel values from the remainder of the individual image, these can be recognized relatively easily by a threshold test and subsequent erosion. For each individual image, all the pixel values are compared with a highlight threshold. The pixels, the pixel values of which exceed the highlight threshold, are allocated to a highlight mask. The highlight mask thus produced is expanded by erosion. This corresponds to a morphological dilation with a round structure element. Using the binary highlight mask, features in the region thereof are excluded. For this purpose, no features are determined in the region of the highlight mask, or features determined therein are not tracked, or not used to calculate the transformation matrix.
The invention is also based on the object of providing a device for joining a plurality of individual digital images to form a total image, which is fast, in particular real time-capable, and precise.
This object is achieved by a device for joining a plurality of digital images to form a total image with

    
    
        an input unit for providing a first individual digital image and a second individual digital image overlapping the latter,
        a selection unit to determine a plurality of features of the first individual image, which is configured in such a way that the features are each allocated first feature coordinates,
        a tracking unit to determine the features in the second individual image, which is configured in such a way that the features are each allocated second feature coordinates,
        a transformation unit to determine a transformation matrix having at least six degrees of freedom between the individual images, which is configured in such a way that
        
            a first part quantity of the features is selected,
            a test transformation matrix is calculated from the first and second feature coordinates of this part quantity,
            a quality value is determined for the test transformation matrix with a second part quantity of the features, and
            the transformation matrix is determined proceeding from the test transformation matrix if the quality value satisfies a quality criterion, and
        
        
        an output unit to join the individual images to form a total image by means of the transformation matrix.
    
    


The advantages of the device according to the invention correspond to the advantages already described of the method according to the invention. In particular, the device can also be developed according to the method of the invention.
A device, in which a rectifying unit is provided to rectify the individual images, increases the precision in joining the individual images to form an overall image, in particular if the individual images are recorded by means of an endoscope system.
A device, in which a highlight unit is provided to determine highlight masks of the individual images, increases the precision when joining the individual images, if these contain highlights.
An endoscope system with an imaging sensor arranged on a shaft to record individual digital images and a device for joining a plurality of individual digital images to form a total image according to the invention, has a virtually extended viewing field. The endoscope system can be used in the technology for analyzing and checking components and for the documentation of faulty locations. Moreover, the endoscope system can be used in medicine for investigating hollow organs and for documenting lesions and in minimally invasive surgery.
The various features of novelty which characterize the invention are pointed out with particularity in the claims annexed to and forming a part of this disclosure. For a better understanding of the invention, its operating advantages and specific objects attained by its uses, reference is made to the accompanying drawings and descriptive matter in which preferred embodiments of the invention are illustrated.



BRIEF DESCRIPTION OF THE DRAWINGS
In the drawings:
FIG. 1 is a schematic view of an endoscope system with a device for joining individual digital images to form a total image;
FIG. 2 is a schematic view of two consecutively recorded individual images;
FIG. 3 is a schematic view of a total image joined from the individual images in FIG. 2;
FIG. 4 is a schematic flowchart of the method for joining the individual images to form the total image; and
FIG. 5 is a schematic flowchart of the feature-based algorithm for determining and tracking features in the individual images.



DESCRIPTION OF THE PREFERRED EMBODIMENTS
An endoscope system 1 has a shaft 2, on which a handle 4 is arranged on a first end 3 and an imaging system 6 is arranged on a second end 5. The shaft 2 is flexible. Alternatively, the latter may also be rigid.
The optical system 6 comprises a radiation source 7, a wide angle lens 8 and an imaging camera sensor 9. The radiation source 7 and the camera sensor 9 can be controlled by the handle 4. The imaging sensor 9 is, for example, configured as a miniaturized video sensor in Tip-Chip technology.
The individual images Ia, Ib recorded by the sensor 9 are transmitted by a signal line 10 from the sensor 9 to a device 11, which is used to join the individual digital images Ia, Ib to form a total image G. The signal line 10 is configured, for example, as an electric image transmission line.
The device 11 is configured as a personal computer, on which an input unit 12, a rectifying unit 13, a highlight unit 14, a selection unit 15, a tracking unit 16, a transformation unit 17 and an output unit 18 are implemented. To store data, the device additionally has a storage unit 19. To visualize the total image G, the device 11 is connected to a monitor 20.
Using the camera sensor 9, during a time interval Δt, a first individual image Ia and a second individual image Ib that is later with respect to time are initially recorded. The camera sensor 9 was displaced between the recording of the individual images Ia and Ib, so that the scene located in front of the camera sensor 9 was recorded from different positions.
The individual image Ia is firstly read in step S1 into the input unit 12 and then rectified in the rectifying unit 13 in step S2 by means of an inverse distortion model. The rectified individual image Ia is stored in the storage unit 19 in step S3. The rectified individual image Ia is then passed to the highlight unit 14, which produces, in step S4, a highlight mask La, which is stored in the storage unit 19 in step S5. To produce the highlight mask La, the pixel values pa(x,y) are compared with a highlight threshold t1 and those pixels are allocated to the highlight mask La, the pixel values pa(x,y) of which exceed the highlight threshold value t1. This pixel range or these pixel ranges is or are then extended by morphological dilation with a round structure element and the highlight mask La is thus formed.
The individual image Ia is then fed in step S6 to the selection unit 15 and the tracking unit 16, in which the feature-based algorithm, for example the KLT algorithm, is implemented. The precise mode of functioning of the KLT algorithm will be described in more detail below.
The selection unit 15 in step S7 outputs n features M1 to Mn, which were determined in the individual image Ia. The features M1 to Mn are characterized by associated feature coordinates p1a to p1a, which are stored in the storage unit 19. The individual image Ia is fed in step S8 to the output unit 18, with which this is shown in step S9 in a planar projection face. No joining takes place yet for the first individual image Ia, as no further individual image is present. The origin of the first individual image Ia can be re-determined in step S9. For this purpose, an origin transformation matrix Ha,0 is read out from the storage unit 19 in step S10. In step S11, a return is then made to step S1, where the following second individual image Ib is read in step S1. The steps S2 to S7 are then repeated for the individual image Ib.
The features M1 to Mn determined in the first individual image Ia are then determined in the second individual image Ib. This process is also called tracking or following. In step S12, the features M1 to Mn and parameters for the KLT algorithm, such as, for example, a search radius r, are read from the storage unit 19 into the tracking unit 16. The mode of functioning of the tracking unit 16 will be described below in detail.
The result of the KLT algorithm, which is output by the tracking unit 16, is feature correspondences, which for each of the features M1 to Mn consist of the first feature coordinates p1a to pna and associated second feature coordinates p1b to pnb. The feature correspondences are accordingly (p1a, p1b) to (pna, pnb). The feature correspondences are formed in step S13, the highlight masks La and Lb entering in step S14, so that no features M are permitted, which lie within one of the highlight masks La or Lb.
In step S15, the transformation matrix Hb,a is calculated in the transformation unit 17. The calculation takes place in accordance with equations (1) to (24) by means of the RANSAC algorithm, for which purpose, in step S16, parameters for the RANSAC algorithm, such as, for example, the size of the first part quantity T1 and the second part quantity T2 of features Mi and the number N of iterations are read in. For subsequent individual images, a feature history can additionally be read in step S16. The feature history contains all the coordinates and values of a feature over the time and may be useful when features become lost. Using the feature history, additional feature correspondences between past individual images and the current individual image can be generated in order to stabilize the total image or the production process thereof. The feature coordinates of the lost feature are further transformed by means of the transformation matrix or over a plurality of consecutive individual images by means of a plurality of transformation matrices. As a result, these features can be continued in the list of feature correspondences.
To calculate the transformation matrix Hb,a, the feature coordinates are firstly normalized and the transformation matrix Hb,a is then calculated by means of the RANSAC algorithm. For this purpose, a first part quantity T1 is selected from the total quantity of features M1 to Mn. In the case of a projective transformation matrix, eight unknowns are to be determined, so at least four feature coordinates have to be selected in order to clearly calculate the transformation matrix Hb,a. The RANSAC algorithm carries out N-iteration steps, a first part quantity T1 being randomly selected in each iteration step and a test transformation matrix T being calculated. The calculated test transformation matrix T is then checked with a second part quantity T2 of features Mi, the part quantities T1 and T2 not overlapping. The number NI of the so-called inliers is used as the quality value for the respective test transformation matrix T. A feature Mi of the second part quantity T2 is an inlier when it can be imaged within a predefined accuracy by the calculated test transformation model. The Euclidean distance is calculated for this purpose with the test transformation matrix T. The test transformation matrix T with the highest number NI of inliers is used as the starting point to calculate the transformation matrix Hb,a with all the features Mi characterized as inliers of the second part quantity T2.
If no transformation matrix is found, the individual image Ib is discarded in step S17 and a further individual image is read in step S11. If no further individual image is present, the method ends in step S18. If a transformation matrix Hb, a is found, this is passed in step S19 to the output unit 18, in which the individual images Ia and Ib are joined to form the total image G. The total image G is visualized by the monitor 20, the current individual image Ib preferably being superimposed on the preceding individual image(s) Ia. The most current individual image Ib is preferably shown centrally on the monitor 20.
The transformation matrix Hb, a is stored in the storage unit 19 in step S20. If no further individual image is present, the method ends in step S18.
The mode of functioning of the selection unit 15 and the tracking unit 16 will be described in more detail below with reference to FIG. 5. In step S6, the individual image Ia is firstly read in. The individual image Ia is then smoothed in step S61 with a smoothing filter and then a resolution pyramid is set up in step S62 from a plurality of resolution stages j. The number of resolution stages j is calculated for example, from a predetermined search radius r.
Features M1 to Mn are then determined or extracted in accordance with equations (16) to (19) from the individual image Ia. The features M1 to Mn are calculated all the resolution stages j, so that they can also be tracked in the case of large displacements in the subsequent individual image Ib. The displacement vector is firstly determined in the lowest resolution stage and then projected into the next higher resolution stage. The displacement in pixels can therefore be limited in each resolution stage to a desired amount.
A pixel window W with a predefined width b in the x-direction and a predefined height h in the y-direction is displaced over the individual image Ia in step S63, wherein in each or a plurality of positions of the pixel window W, a matrix G is calculated in accordance with equation (17). The matrix G is produced from the sum of the products ggT over all the pixels of the pixel window W, g being the gradient vector for a specific pixel of the pixel window W. The pixel window W of the feature M1 with a gradient vector g is drawn by way of example in FIG. 2. In addition, a weighting function w can be taken into account. If this is not desired, w=1 can be selected. For each pixel window W, the eigenvalues λ1 and λ2 of the matrix G are calculated and a decision is made with the aid of equation (19) whether the pixel window W characterizes a feature M. In addition, a feature M being newly added has to maintain a minimum distance ΔM from the already determined features M. In step S64, storage of the determined features M1 to Mn is initiated, so that the latter are transmitted to the storage unit 19 in step S7. A return is then made in step S65 to step S61, where the second individual image Ib is read in and smoothed and a resolution pyramid is set up for this in step S62.
For the following or tracking of the features M1 to Mn in the tracking unit 16, these are selected one after the other in step S66. Firstly, for the first feature M1, the individual images Ia and Ib in the lowest resolution stage j=1 are selected in S67. The equation system (16) is then solved for this resolution stage. The solution of the equation system (16) is the displacement vector d for the first feature M1. The matrix G is already known for step S68. The error vector e is calculated according to equation (18) in step S68, the pixel window W being assumed in a first iteration step i=1 in the second individual image Ib at a location corresponding to the first individual image Ia. The displacement vector d(i, j) determined in the first iteration step still has errors. For a second iteration step i=2, the pixel window P in the second individual image Ib is assumed at a location displaced by the displacement vector d determined in the first iteration step. The equation system (16) is solved again with a newly calculated error vector e, the second displacement vector d determined in the second iteration step being an improved solution of the equation system (16). The described iteration step S69 is repeated until the solution for the displacement vector d converges and the feature M1 in the second individual image Ib is determined in step S70. The solving of the equation system (16) can take place by means of the Newton-Raphson algorithm. If no converging solution for the displacement vector d is found after a maximum number of iterations, the associated feature M1 is discarded in step S71 and a return is made in step S72 to step S66, where the next feature M2 is selected.
If the feature M1 is determined in step S70, a return is made in step S73 to step S67, where the next highest resolution stage j=2 of the individual images Ia and Ib is selected. The steps S68, S69 and S70 are then repeated for this resolution stage taking into account the displacement vector d determined in the preceding resolution stage. If this loop has been run through for all the resolution stages j, a check is made in step S74 whether further features M to be determined are present. If this is the case, step S72 is repeated. If all the features M were determined in the individual image Ib, lost features M are replaced in step S75. A feature M can be lost if the latter is no longer contained in the individual image Ib because of the displacement of the sensor 9, or it was discarded as unusable in step S71. In step S75, lost features M are replaced in such a way that the number of features M is again n. The replacement of the features M in step S75 takes place in accordance with step S63. A storage process is initiated for the features M to be newly added in step S67. If further individual images are present, step S65 is repeated. If no further individual images I are present, the KLT algorithm is left in step S76 and a transfer is made to step S13.
To join the individual images Ia and Ib, the x- and y-coordinates of the pixels of the displaced individual image Ib are rounded to the nearest integral value. This method is also called the “nearest neighbor” method. The following individual images can be calculated by multiplication of the current transformation matrix Hi,i-1 by the preceding transformation matrix Hi-1,0 as follows:

Hi,0=Hi-1,0·Hi,i-1  (25)

Using the method according to the invention and the corresponding device 11, 6 to 25 images per second can be achieved on a PC of the type Intel Core 2 6420 with 2.13 GHz and 2 GB RAM.
While specific embodiments of the invention have been shown and described in detail to illustrate the application of the principles of the invention, it will be understood that the invention may be embodied otherwise without departing from such principles.

## Claims

1. A method for joining a plurality of individual digital images to form a total image, the method comprising the steps:
providing a first individual digital image and a second individual digital image overlapping the first individual digital image;
determining a plurality of features of the first individual image, wherein the features are each allocated first feature coordinates;
determining the features in the second individual image, wherein the features are each allocated second feature coordinates;
normalizing the feature coordinates before the calculation of a transformation matrix, wherein the normalization takes place in such a way that the feature coordinates belonging to an individual image have an average distance of √{square root over (2)} from a focal point;
determining the transformation matrix having at least six degrees of freedom between the individual images, wherein:
a first part quantity of the features is selected;
a test transformation matrix is calculated from the first feature coordinates and the second feature coordinates of this part quantity;
a quality value is determined for the test transformation matrix with a second part quantity of the features; and
the transformation matrix is determined proceeding from the test transformation matrix, when the quality value fulfils a quality criterion; and

joining the individual images to form a total image by means of the transformation matrix.

2. A method according to claim 1, wherein the feature coordinates are determined iteratively in different resolution stages of the individual images.

3. A method according to claim 1, wherein the features of the first individual image are determined, wherein in each case:
a pixel window is formed from a plurality of adjacent pixels;
a matrix is formed from the gradients for the plurality of adjacent pixels;
the eigenvalues are calculated from the matrix;
the pixel window is characterized as the feature, if the eigenvalues exceed a threshold value; and
the feature is allocated first feature coordinates; and
the features in the second individual image are determined, wherein in each case:
an error vector is formed for the feature as a function of the gradients and the differences in the pixel values of the individual images in the pixel window;
a displacement vector is calculated from the error vector and the matrix; and
the feature is allocated second feature coordinates.

4. A method according to claim 3, wherein the calculation of the respective displacement vector takes place iteratively, wherein:
the individual images are displaced relative to one another by a first displacement vector; and
a second displacement vector is calculated by means of the displaced individual images.

5. A method according to claim 3, wherein the calculation of the respective displacement vector takes place in different resolution stages of the individual images, wherein:
a first displacement vector is calculated in a first, low resolution stage of the individual images;
the individual images are displaced relative to one another in a second, higher resolution stage by the first displacement vector; and
a second displacement vector is calculated by means of the displaced individual images in the second resolution stage.

6. A method according to claim 1, wherein a pixel window is only characterized as feature if a minimum distance from the already determined features is exceeded.

7. A method according to claim 1, wherein the transformation matrix has precisely eight degrees of freedom and characterizes a projective transformation.

8. A method according to claim 1, wherein the calculation of the transformation matrix takes place in such a way that for a repetition count:
a first part quantity of the features is in each case selected;
a test transformation matrix belonging to said first part quantity is calculated;
the number of features of a second part quantity is determined for the respective test transformation matrix as the quality value and the features of said second part quantity are imaged by the test transformation matrix with a defined accuracy;
the test transformation matrix with the highest number is selected; and
the transformation matrix is calculated from the associated features, which is imageable with the defined accuracy.

9. A method according to claim 1, wherein the individual images are joined in a planar projection face.

10. A method according to claim 1, wherein the individual image which is later in terms of time during joining is superimposed on the individual image that is earlier in terms of time.

11. A method according to claim 1, wherein, for joining, the x-coordinate and the y-coordinate of the pixels of the displaced individual image are rounded to the nearest integral value.

12. A method according to claim 1, wherein the individual images are recorded by means of an endoscope system.

13. A method according to claim 1, wherein the individual images are rectified before the determination of the features.

14. A method according to claim 1, wherein a highlight mask is set up for each individual image, wherein pixel values are detected, which exceed a highlight threshold value and no features are permitted for the calculation of the transformation matrix in the highlight masks.

15. A method according to claim 1, wherein the transformation matrix is calculated via a RANSAC algorithm after the feature coordinates are normalized.

16. A device for joining a plurality of digital images to form a total image, the device comprising:
an input unit for providing a first individual digital image and a second individual digital image overlapping the first individual digital image;
a selection unit to determine a plurality of features of the first individual image, which is configured in such a way that the features are each allocated first feature coordinates;
a tracking unit to determine the features in the second individual image, which is configured in such a way that the features are each allocated second feature coordinates;
a transformation unit to determine a transformation matrix having at least six degrees of freedom between the individual images, which is configured in such a way that:
the feature coordinates are normalized before the calculation of the transformation matrix, wherein the normalization takes place in such a way that the feature coordinates belonging to an individual image have an average distance of √{square root over (2)} from a focal point;
a first part quantity of the features is selected;
a test transformation matrix is calculated from the first feature coordinates and the second feature coordinates of the first part quantity;
a quality value is determined for the test transformation matrix with a second part quantity of the features; and
the transformation matrix is determined proceeding from the test transformation matrix if the quality value satisfies a quality criterion; and

an output unit to join the individual images to form a total image by means of the transformation matrix.

17. A device according to claim 16, wherein a rectifying unit is provided to rectify the individual images, the transformation matrix being calculated via a RANSAC algorithm after the feature coordinates are normalized.

18. A device according to claim 16, wherein a highlight unit is provided to determine highlight masks of the individual images.

19. An endoscope system, comprising:
an imaging sensor arranged on a shaft to record individual digital images;
a device for joining a plurality of individual digital images to form a total image with an input unit for providing a first individual digital image and a second individual digital image overlapping the first individual digital image;
a selection unit to determine a plurality of features of the first individual image, which is configured in such a way that the features are each allocated first feature coordinates;
a tracking unit to determine the features in the second individual image, which is configured in such a way that the features are each allocated second feature coordinates;
a transformation unit to determine a transformation matrix having at least six degrees of freedom between the individual images, which is configured such that:
the feature coordinates are normalized before the calculation of the transformation matrix, wherein the normalization takes place in such a way that the feature coordinates belonging to an individual image have an average distance of √{square root over (2)} from a focal point;
a first part quantity of the features is selected;
a test transformation matrix is calculated from the first feature coordinates and the second feature coordinates of the first part quantity;
a quality value is determined for the test transformation matrix with a second part quantity of the features; and
the transformation matrix is determined proceeding from the test transformation matrix if the quality value satisfies a quality criterion; and

an output unit to join the individual images to form a total image by means of the transformation matrix.

20. An endoscope system according to claim 19, wherein the transformation matrix is calculated via a RANSAC algorithm after the feature coordinates are normalized.