RELATED APPLICATIONS 
     This application is a 35 U.S.C. 371 national stage filing of International Application No. PCT/EP2005/013488 filed 15 Dec. 2005, which claims priority to Great Britain Patent Application No. 0428175.4 filed on 23 Dec. 2004 in Great Britain. The contents of the aforementioned applications are hereby incorporated by reference. 
     FIELD OF THE INVENTION 
     The present invention relates to image processing methods. The methods provided by the invention find particular use in the processing of images obtained from medical investigative tools such as CT, MRI, PET, portal images, and others. 
     BACKGROUND ART 
     One of the most commonly used image processing methods is that of co-registering pictures to each other, i.e. transforming one image to the reference system of another image so that they can be compared directly or even merged. A common example is for a CT image taken on one occasion to be compared with an MRI image taken on another occasion. The process of registration of one image with the other makes it possible for the user directly to compare the images correctly, or enables a fusion of the images whereby both images are combined to a single image. 
     Both automatic and manual methods of doing so exist at present. 
     SUMMARY OF THE INVENTION 
     Hitherto, methods for the registration of images have assumed that there are only two images. Where more than two images are to be registered, these have been co-registered in appropriate pairs, using the same techniques. 
     The present invention is based on a realisation that this assumption is often incorrect. In the medical field (for example) there are often many images of a patient, derived at different times via different investigative apparatus. The registration of these multiple images provides an opportunity to improve the accuracy of registration of chosen pairs of images within that multiplicity. 
     In general, the registration process involves the identification of (in general terms) ‘n’ transformation steps, such as a translation in both the x and y directions, a rotation, a scaling, and other degrees of freedom, all of which when applied together and in an appropriate order will transform one image into the other. Clearly, ‘n’ will vary with the intended accuracy of the transformation, in terms of the number and type of steps that are to be carried out, and with the number of dimensions of each of the images. A pair of two-dimensional images that are translated in two directions, rotated, and scaled in two directions will thus give rise to a five-dimensional vector. In general, the factors can be considered to form a single n-dimensional vector quantity, representing the overall transformation that is required. This vector will be referred to herein as a ‘registration vector’. 
     This vector approach allows the registration process to be conceptualised in the familiar terms of vector relationships, from which an understanding of the invention can be derived. Assuming a simple case of three images (A, B and C), these will define three points in the relevant n-dimensional space and the registration vectors will represent the displacements necessary to move from one point to another. Thus, to reach B from A, it is possible to use the single transformation that links the two points, or to move first via the transformation A to C followed by the transformation C to B. 
     In general, having obtained registration vectors for each of the three paths linking the three images, and thus both the single path vector A-B and the composite vector A-C-B, these two vectors will not be precisely equal as theory demands. This difference between the two possible paths from A to B represents the errors present in both; each will be an imperfect approximation to a theoretically perfect but unknown vector that represents the actual (ideal) transformation. However, given that it is to be expected that random errors in the vectors will outweigh systematic errors, it will be apparent that an average of the two path vectors will more closely approximate to the actual transformation. 
     Thus, the present invention is based on the realisation that by using the remaining images to obtain a choice of paths between the two selected images, the true transformation between the two can be determined with greater accuracy by averaging those paths. This is, at first sight, counterintuitive. It would be expected that the contribution of the indirect paths would result in a reduction in overall accuracy, given that they rely on the summation of a plurality of imperfect approximations. However, this is not in fact the case, since the third image (and any subsequent images) offer new information not present in the two selected images, which can be used to improve the accuracy of the overall method. 
     It is to be preferred that, when averaging the paths, greater weight is given to paths whose accuracy is known or reasonably believed to be greatest. It may be possible to place an a priori expectation of accuracy on certain paths, since they may involve the registration of images that have greater fidelity, or which have registration marks included in the image, for example. It is also possible to estimate the registration accuracy based on analysis of the image detail. 
     It will be appreciated that the method is readily applicable to more than three images, with the greater number of possible paths allowing a wider average to be taken. 
     Iteration of the process may be possible, where the available computation power is available. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       An embodiment of the present invention will now be described by way of example, with reference to the accompanying figures in which; 
         FIG. 1  shows a schematic set of three images and the transformation required between them; 
         FIG. 2  shows a representation in two dimensions of a hypothetical path set consisting of three images; 
         FIG. 3  shows the averaging process and the consequent reduction in error; and 
         FIG. 4  shows a representation in two dimensions of a hypothetical path set consisting of four images. 
     
    
    
     DETAILED DESCRIPTION OF THE EMBODIMENTS 
     The challenge lies in obtaining the most accurate co-registration as possible. Many contexts require an error of less than 1 mm for the resulting output to be clinically useful. In non-ideal situations where the images have a poor resolution, poor contrast or contain artefacts, or when difficult combinations such as PET to CT are presented, known methods can face difficulties in obtaining the necessary level of accuracy. 
     A further problem with existing automatic methods of registration is that there is no established quality measure for the predicted level of accuracy. In manual registration, where the user places a number of landmarks on both images, it is possible to calculate an estimation of the registration error. No corresponding possibility for automatic methods exists. The user must therefore rely on a purely visual inspection of the result. 
     All present registration techniques, both automatic and manual, use a target image and a reference image. The aim with the co-registration process is to transform the target image on to the reference image. Often, however, the user has access to more than two images. An example might be a frameless PET image taken prior to treatment, which is to be co-registered with a CT image and an MRI image taken, with a frame and indicator box, just prior to treatment. The co-registration of the CT and MR-image is obtained directly by identifying the fiducials of the indicator box. The target image (the PET image) is thereafter typically registered on either the CT image or the MRI image acting as a reference image. 
     This presents an opportunity to improve the process, compared to prior art methods, by using the information that the third image contains. 
     Referring to  FIG. 1 , this shows three sample images by way of schematic representation. The first image  10 , designated ‘A’ contains no registration information but simply includes anatomical detail  12  of the patient. This is typical of a PET image, for example. The second image  14 , designated ‘B’, contains the same anatomical detail  16  but in a different reference frame. In simple terms, this means that the images are essentially of the same general area of the patient, but may be subject to (for example) translation, rotation and scaling so that the two are not directly comparable. Image B may, for example, represent a CT image and is thus shown with reference marks  18 . 
     Finally, the third image  20 , designated ‘C’, again contains the same anatomical detail  22  but is set in a third reference frame. This might, for example, be an MRI image. Images of the same reference markers  18  are visible, in different image locations reflecting the differences between the second and third reference frames. 
     Thus, the co-registration between A and B can be calculated directly or be obtained by the merging of co-registrations A-C and C-B. By calculating the weighted mean value (in the manner described below) between the co-registration A-B and the merged co-registration A-C-B, the expected error can be lowered by a factor 
     
       
         
           
             1 
             - 
             
               
                 
                   
                     k 
                     
                       a 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       c 
                     
                     2 
                   
                   + 
                   
                     k 
                     cb 
                     2 
                   
                 
               
               
                 
                   1 
                   + 
                   
                     k 
                     
                       a 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       c 
                     
                     2 
                   
                   + 
                   
                     k 
                     cb 
                     2 
                   
                 
               
             
           
         
       
     
     as explained below, where k ac  is the standard deviation for the error in the A-C coregistration divided with the error in the A-B coregistration and where k cb  is the standard deviation of the error in the C-B coregistration divided with the error in the A-B coregistration. 
     If we suppose that all registrations are just as good, i.e. k ac =k cb =1, the expected error is lowered by 18%. Of course, in this example the coregistration C-B is expected to have half the error of the A-B and A-C coregistrations, since it is calculated using the indicator box, thereby giving a better result (particularly if the images are of poor quality). In this case, the expected error is lowered by more than 25%. 
     It is noteworthy that the expected error is always lower, even if the coregistrations of A-C and B-C have an expected error greater than that of coregistration of A-B. The method still works well with more than three images, and then gives even higher gains. 
     The method can be summarized as follows. 
     Let each image be a node in a graph. A calculated co-registration between two images forms an edge between corresponding nodes. Each image is thus registered on each other image. A coregistration between two selected images is obtained as the weighted mean value of all paths between the images/nodes, where the weights are calculated from an estimation of the variance of each registration, and is chosen so as to minimize the variance in the result. The deviation between the different paths can be used as a measure of the quality of the registration. If landmark-based registration is used, then any of the the usual error estimates could be used to estimate the variance of each registration. Otherwise, an estimation using the images&#39; voxel size could be used. The expected error of the resulting registration is always lower compared with the direct registration. 
       FIG. 2  shows such a graph. Image A is denoted at  100 , Image B at  102  and Image C at  104 . The direct registration vector from A to B is shown as edge  106 . However, it will be seen that it is also possible to reach B via image C, i.e. from  100  via edge  108  to image C at  104 , then via edge  110  to reach image B at  102 . 
     As shown in  FIG. 3 , an enlarged part of  FIG. 2 , image B lies somewhere in the region denoted  102 ′. The end points of registration vectors  106  and  110  do not co-incide exactly; there is a gap  112  between them. An average point  114  along the line connecting the end points of registration vectors  106  and  110  will give a better indication of the end point of the “true” registration vector  116 . As shown, the average point  114  is midway between the two end points, but the weighting applied may affect this. In addition, the true registration vector  116  is shown between the measured vectors  106 ,  110  based on the assumption that systematic errors are outweighed by random errors. 
       FIG. 4  shows an example involving four images A, B, C and D. In this case, there are multiple paths from A to B; 
     
       
         
               
               
               
             
           
               
                   
                   
               
               
                   
                 Images 
                 Route 
               
               
                   
                   
               
             
             
               
                   
                 A-B 
                 200: 208: 202 
               
               
                   
                 A-C-B 
                 200: 210: 204: 218: 202 
               
               
                   
                 A-D-B 
                 200: 216: 206: 214: 202 
               
               
                   
                 A-C-D-B 
                 200: 210: 204: 212: 206: 214: 202 
               
               
                   
                 A-D-C-B 
                 200: 216: 206: −212: 204: 218: 202 
               
               
                   
                   
               
             
          
         
       
     
     It may be preferable to use some or all of these additional paths. The greater number of paths will offer more composite registration vectors, but will increase the computational load. A balance may have to be struck in the light of the available computational effort, the required accuracy, and the time available. 
     Calculation of the Weighted Mean Value 
     Let a, b and c be three images, and suppose that a should be registered to b. Let T ab  be the registration between a and b, and define T ac  and T cb  accordingly. Let T acb  be the registration between a and b that is obtained by composing T ac  T cb . The registration T between a and b is given by the weighted average
 
 T=αT   ab +(1−α) T   acb , αε(0,1)
 
     Let X ab , X ac  and X cb  be the stochastic error in T ab , T ac  and T cb . The stochastic error X for T is thus (since the errors are additive)
 
 X=αX   ab +(1−α)( X   bc   +X   cb )
 
     It can be assumed that X ab , X ac  and X cb  are independent and have expectation values of 0. The error X is thus minimized when the variance V[ X ] is minimized according to 
     
       
         
           
             
               
                 
                   min 
                   ⁢ 
                   
                       
                   
                 
                 
                   α 
                   ∈ 
                   
                     ( 
                     
                       0 
                       , 
                       1 
                     
                     ) 
                   
                 
               
               ⁢ 
               
                 V 
                 ⁡ 
                 
                   [ 
                   X 
                   ] 
                 
               
             
             = 
             
               
                 
                   α 
                   2 
                 
                 ⁢ 
                 
                   σ 
                   ab 
                   2 
                 
               
               + 
               
                 
                   
                     ( 
                     
                       1 
                       - 
                       α 
                     
                     ) 
                   
                   2 
                 
                 ⁢ 
                 
                   ( 
                   
                     
                       σ 
                       
                         a 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         c 
                       
                       2 
                     
                     + 
                     
                       σ 
                       cb 
                       2 
                     
                   
                   ) 
                 
               
             
           
         
       
     
     where σ ab   2 =V[X ab ], σ ac   2 =V[X ac ] and σ cb   2 =V[X cb ]. These can be estimated from (for example) the voxel sizes of the images, or by using one of the existing error estimates for manual co-registrations. Let σ ac =k ac σ ab  and σ cb =k cb σ ab , then 
     
       
         
           
             
               
                 
                   min 
                   ⁢ 
                   
                       
                   
                 
                 
                   α 
                   ∈ 
                   
                     ( 
                     
                       0 
                       , 
                       1 
                     
                     ) 
                   
                 
               
               ⁢ 
               
                 V 
                 ⁡ 
                 
                   [ 
                   X 
                   ] 
                 
               
             
             = 
             
               
                 σ 
                 ab 
                 2 
               
               ⁢ 
               
                 ⌊ 
                 
                   
                     α 
                     2 
                   
                   + 
                   
                     
                       
                         ( 
                         
                           1 
                           - 
                           α 
                         
                         ) 
                       
                       2 
                     
                     ⁢ 
                     
                       ( 
                       
                         
                           k 
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             c 
                           
                           2 
                         
                         + 
                         
                           k 
                           cb 
                           2 
                         
                       
                       ) 
                     
                   
                 
                 ⌋ 
               
             
           
         
       
     
     which assumes its minimum for 
     
       
         
           
             
               
                 
                   
                     α 
                     = 
                     
                       
                         
                           k 
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             c 
                           
                           2 
                         
                         + 
                         
                           k 
                           cb 
                           2 
                         
                       
                       
                         1 
                         + 
                         
                           k 
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             c 
                           
                           2 
                         
                         + 
                         
                           k 
                           cb 
                           2 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     This 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     gives 
                   
                 
               
             
             
               
                 
                   
                     V 
                     ⁡ 
                     
                       [ 
                       X 
                       ] 
                     
                   
                   = 
                   
                     
                       σ 
                       ab 
                       2 
                     
                     ⁢ 
                     
                       
                         
                           k 
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             c 
                           
                           2 
                         
                         + 
                         
                           k 
                           cb 
                           2 
                         
                       
                       
                         1 
                         + 
                         
                           k 
                           
                             a 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             c 
                           
                           2 
                         
                         + 
                         
                           k 
                           cb 
                           2 
                         
                       
                     
                   
                 
               
             
           
         
       
     
     The standard deviation D[X]=√{square root over (V[X])} is thus lower than D[X ab ] with a factor 
     
       
         
           
             
               1 
               - 
               
                 
                   
                     
                       k 
                       
                         a 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         c 
                       
                       2 
                     
                     + 
                     
                       k 
                       cb 
                       2 
                     
                   
                 
                 
                   
                     1 
                     + 
                     
                       k 
                       
                         a 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         c 
                       
                       2 
                     
                     + 
                     
                       k 
                       cb 
                       2 
                     
                   
                 
               
             
             &gt; 
             0 
           
         
       
     
     Given that this factor is greater than zero even for k ac =k cb =1, it follows that the use of additional images to provide alternate paths is a useful means for improving the registration accuracy of the selected pair of images. 
     It will of course be understood that many variations may be made to the above-described embodiment without departing from the scope of the present invention.
Methods for the registration of images typically assume that there are only two images, an assumption that is not always valid. By using the remaining images to obtain a choice of paths between two selected images, the transformation between the two can be determined with greater accuracy by averaging those paths. When averaging the paths, greater weight can be given to paths whose accuracy is known or reasonably believed to be greatest. Iteration of the process may be possible, where the available computation power is available.