Method of Mosaicing Ultrasonic Volumes for Visual Simulation

The present invention relates to image processing for medical imaging and, more particularly, to a method for registering images, especially ultrasound images, of a patient, acquired from several viewpoints.

The problem addressed by the present invention is that of mosaicing a large volume from a set of small volumes acquired from a 3D ultrasound device. The formation of large volumes is necessary for an ultrasound visual simulator, such as the one described in U.S. patent application Ser. No. 08/316,841. Trainees using the simulator are able to learn how to identify and diagnose a wide range of medical cases, including rare pathologies, by operating a simulated ultrasound device on a mannequin, without any need for actual patients. The principle of the simulator is simple but effective. In an off-line preprocess, a volume buffer is generated from real ultrasound images, and then the simulated images are generated by slicing the volume buffer on-line. Such images can be generated very rapidly, including postprocessing enhancements, and can produce images which are, in most cases, indistinguishable from real ultrasound images.

The ultrasonic volume buffer has to be big enough, and has to represent a large portion of the mannequin, to permit the unbounded practice of a real life diagnosis. Contemporary ultrasound devices do not provide the capability of obtaining the entire volume in a single acquisition. This implies that the volume buffer has to be reconstructed from several sub-volumes obtained from different viewpoints. The registration of mono-modal datasets has been extensively investigated in medical applications where atlas data are used. However, ultrasonic datasets are far more problematic than other medical modalities, such as CT or MRI, because-the ultrasound values are significantly more noisy, blurred, deformed, and have view dependent variations which are described below.

The basic technique of a mosaicing operation is to align and register two given volumes with a significant overlap into a single volume which smoothly combines the information from both. The type of registration technique that can be appropriately applied is directly dependent on the type of variation between the two volumes. Thus, to design a registration method for ultrasonic volumes, it is necessary to know the type of variation exhibited by ultrasonic volumes.

The typical size of an ultrasound image generated by common ultrasonic devices is limited to a maximum width of 12-15 cm. The acquisition of a volume is therefore reconstructed from a series of two-dimensional slices. There are two main methods of collecting the series of slices: a free-hand collection and a mechanical collection. In a free-hand collection the location and orientation of the slice is tracked by a six-degree-of-freedom device. The slices are stored in the volume and the gaps between the slices are filled by interpolation. This method is impractical for large volumes because of the long acquisition time and the inaccuracy caused by the numerous unsampled gaps.

A better approach is to attach the transducer probe to a mechanical motor that sweeps the slice along some type of trajectory (e.g., fan, rotation). In particular, with a parallel sweep, a series of parallel slices are obtained with no inter-slice gaps. The desired image resolution may be obtained by sweeping the slice sufficiently slowly. Such parallel dense slices provide small volumes of good quality. A series of such volumes needs to be collected and assembled to form a large volume.

The registration of two volumes requires the detection of the changes between the two images and the design of a transformation that deforms the images in order to remove or reduce the variations between them. The source variations can be classified into the following three types.

1. Directional Variations

These variations are due to changes in the view point. They cause a misalignment that can be simply corrected by a rigid transformation. However, as we noted above, the acquisition of the same volume from a different view point causes other effects that are not compensated for by spatial transformation. For example, shadows are cast with a strong correlation with the probe viewing direction.

2. Volumetric Variations

These are caused by the characteristics of ultrasonic technology, for example, depth gain compensation and gain distortions, and the inherently noisy and blurred ultrasound signal. These effects are difficult to model and to remove. One can attempt to reduce them by tuning the acquisition parameters.

3. Geometric Variations

Geometric deformations are caused by the movements of the body of the patient during acquisition. Some movements are forced by the acquisition device, because the ultrasound probe must have good contact with the body. The human body is soft and not flat, and it is rather difficult to maintain contact without causing the patient's muscles to contract.

Another unavoidable deformation is caused by breathing and other natural behavior of the patient. Periodic deformation (like that of the heart) can be overcome by "gating". In gating, the acquisition is synchronized with the period of the deformation, and the slices are acquired in phase with the deformation, using equipment similar to E.C.G. that monitors heart activity.

In most cases, the registration of ultrasound images cannot be based on image features. In ultrasonic data, it would be very difficult for even a human expert, let alone an automatic procedure, to identify useful features, because 3D features can appear significantly different in two different images. Therefore, a technique for registering ultrasound images must be based directly on image pixel values, and not on the interpretation of those values. A reasonable strategy would be to base the registration technique on some "majority". That is, the success of the technique depends on the average success over the whole image, even though part of the image may not contain any useful information. Such techniques are classified as automatic direct registration techniques. One such technique, for two dimensional images, was published by Barber (D. C. Barber (1992), Registration of low resolution medical images, Phys. Med. Biol., vol. 37 no. 7, pp. 1485-1498, which is incorporated by reference for all purposes as if fully set forth herein). Barber's technique is not directly suitable for use with ultrasound data sets, however, because of the unique characteristics of ultrasound data sets noted above.

There is thus a widely recognized need for, and it would be highly desirable to have, a method for automatic direct registration of ultrasound images.

According to the present invention there is provided a method for warping an input target image to match an input reference image, both images being at an original level of resolution, both images being composed of a plurality of pixels, the method comprising the steps of: (a) convolving the input target image with a low pass filter, thereby obtaining a filtered target image; (b) convolving the input reference image with said low pass filter, thereby obtaining a filtered reference image; (c) subsampling said filtered target image, thereby obtaining a subsampled target image; (d) subsampling said filtered reference image, thereby obtaining a subsampled reference image; (e) determining an affine transformation that maximizes a normalized cross-correlation between said subsampled target image and said subsampled reference image; and (f) scaling said affine transformation to apply to the original level of resolution, thereby obtaining a scaled factor of said affine transformation.

According to the present invention there is provided a method for warping an input target image to match an input reference image, both images being of a dimensionality D, both images being at an original level of resolution, both images being composed of a plurality of pixels, the method comprising the steps of: (a) subdividing the input target image into a plurality of target subimages; (b) subdividing the input reference image into a plurality of reference subimages, there being a one-to-one correspondence between said reference subimages and said target subimages; (c) for each of said target subimages, determining a local transformation that warps said target subimage to match said corresponding reference subimage.

According to the present invention there is provided a method for blending a three-dimensional target image with a three-dimensional reference image, the target image including a plurality of points, the method comprising the steps of: (a) determining at least one seam plane shared by the target image and the reference image; (b) for each of said at least one seam plane: (i) extracting an input target slice from the target image along said seam plane, (ii) extracting an input reference slice from the reference image along said seam plane, (iii) subdividing the input target slice into a plurality of target subslices, (iv) subdividing the input reference slice into a plurality of reference subslices, there being a one-to-one correspondence between said reference subslices and said target subslices, (v) for each of said target subslices, determining a local transformation that warps said target subimage to match said corresponding reference subimage, (vi) recursively performing both said subdividing and said determining of said local transformation at a plurality of levels of localization, said recursion further including the steps of: (A) for each of said target subslices, composing said local transformation with a prior transformation, thereby obtaining a local composed transformation, and (B) forming a weighted sum of said local composed transformations, thereby obtaining a global transformation; (c) applying a weighted sum of said at least one global transformation to at least one of the points.

The method of the present invention is built on an affine mapping of a target image into a reference image. Because this affine mapping is based on the gradients of the images, it is referred to herein as the "gradient technique". The naive gradient technique is extended herein by two improvements, referred to herein as "multiresolution and localization". The derivations herein are 2D derivations: their extension to 3D and higher dimensions will be obvious to one skilled in the art.

The gradient technique assumes that the misalignment of the two images is not too large. That is, the starting position of the two images is "good enough", where the exact meaning of "good enough" is discussed below. There are other techniques which can register two images with arbitrary partial overlapping. Such techniques usually perform a coarse alignment, possibly by rigid transformation. The present invention includes a partial matching method which applies rigid transformation to bring the two images to a satisfactory starting position for the fine and elastic parts of the registration.

The naive gradient technique finds an affine operator T that transforms the coordinates of the target image into the coordinates of the reference image in a way that maximizes the overlap of the two images. This technique is based on a Taylor expansion of the exact mapping, under the assumption that there is a majority of target pixels which are close enough to the corresponding positions in the reference image. The Taylor expansion requires, as its input, gradients of the images These gradients are approximated numerically. Whenever the displacement is large, however, the Taylor expansion is no longer valid. Barber suggested a technique to converge to the solution by iterations. The series of transformations are computed, where at each step k the solution for T.sup.k is found over an image warped by the previous transformation T.sup.k-1. The transformations must include some filtering and not cause unnecessary loss of information. To minimize applying many transformations, Barber suggested composing the transformations. Because Barber uses bilinear transformations, which are not closed under composition, he approximates his composition by another bilinear transformation.

The solution based on iterations assumes that at each step the majority of pixels contribute a displacement towards the correct solution. However, in many cases in practice, the initial displacement is too large. The present invention includes a multiresolution method to better handle large misalignments. The transformation is computed on a resolution pyramid and the results from the low resolution transformation are used to guide the computation of finer levels. A resolution pyramid consists of the original image and a number of copies of it at lower resolutions. At lower resolutions, adjacent pixels and local gradients represent large distances of the original image. Thus, the low resolution transformation is scaled up according to image resolution. Because all the transformations at all levels are affine transformations, their composition is also an affine transformation. Thus, the final image is a transformation of the original image, and the multiresolution method does not cause any information loss.

Once the final transformation is computed it is applied to the target image. Then, the transformed image is partitioned into four sub-images and the registration process is repeated for each of the sub-images to yield four local transformations which refine the registration locally. This process is recursively repeated for each sub-image and the registration is further refined. The aggregation of the independent local affine transformations forms a global elastic transformation with a high degree of freedom.

The above registration method is static in the sense that all the pixels of the image contribute to the solution. Many pixels may not contain any useful information, or their local environment may not have a matching counterpart in the other image. This situation often occurs in ultrasound images because of hidden parts or noise. Therefore, the present invention includes a rejection mechanism to cull away pixels that are estimated to be irrelevant to the solution. Such pairs of pixels are not assigned to the linear system of equations. This mechanism significantly improves the registration of partially occluded scenes.

As noted above, the present invention relies heavily on the Taylor expansion and the correctness of the gradients. Thus, the images should be continuous and differential, which is not necessarily the case in many discrete images. However, the registration method need not be affected by sharp local variations. Local variations may be softened and blurred by a filter as a preprocess to the gradient calculation. The filter also removes part of the noise from the images and smoothes them. This improves the exactness of the Taylor expansion and hence the quality of the results. However, after the transformation is computed over the filtered images, it is applied to the original images during the mosaicing stage.

The invention is herein described, by way of example only, with reference to the accompanying drawings, wherein

FIG. 1 is a diagram of the data flow of multiresolution registration;

FIG. 2A is a reference image;

FIG. 2B is a target image to be warped to match the reference image of FIG. 2A;

FIG. 2C is the result of warping the target image of FIG. 2B without using multiresolution;

FIG. 2D is the result of warping the target image of FIG. 2B using multiresolution;

FIG. 2E is the image of FIG. 2B after one multiresolution step;

FIG. 2F is the image of FIG. 2B after two multiresolution steps;

FIG. 3A is a reference image;

FIG. 3B is a target image to be warped to match the reference image of FIG. 3A;

FIG. 3C is the result of warping the target image of FIG. 3B using a global transformation, including multiresolution;

FIG. 3D is the result of warping the target image of FIG. 3B using one level of localization;

FIG. 3E is the result of warping the target image of FIG. 3B, using two levels of localization;

FIG. 4 shows the data flow of the localization process;

FIG. 5 shows schematically the overlap of a target volume and a reference volume, including two seam planes;

FIG. 6 shows the volume mosaicing pipeline.

The present invention is of a method for registration of medical images, particularly ultrasound images. Specifically, the present invention can be used to create an internally consistent numerical database of ultrasound images from a collection of disparate, overlapping digital ultrasound subimages.

The principles and operation of medical image registration according to the present invention may be better understood with reference to the drawings and the accompanying description.

The Gradient Method

Given two images, a target image I and a reference image J, the goal is to find a transformation T such that L(T(I,J) is maximal, where L is the normalized cross-correlation operator ##EQU1## and I and J are the average values of the target and reference images, respectively. In the development it is assumed that the images are smooth differentiable functions. In the algorithm, a discrete version of the derivatives involved is used.

The transformation T can be defined as a polynomial which maps the x,y coordinates of I to the u,v coordinates of J:

In the preferred embodiments of the present invention, an affine transformation is used:

but other degrees can be used as well. For example, Barber used a bilinear transformation.

To find the coefficients of the affine transformation, a linear system of equations is constructed from a large set of pairs of points p.sub.i .epsilon.I and q.sub.i .epsilon.J. Those point are implicitly described by their displacement, which can be expressed in terms of the transformation.

Let T(x,y)=(u,v). If T warps I onto J, it is assumed that I(xy)=J(u,v) holds for the majority of the image pixels. Define dx=u-x and dy=v-y. Then I(x,y)=J(x+dx, y+dy). Assuming the displacement vector (dx,dy) is small enough to be approximated by the linear term of the Taylor expansion, then

However, because

it follows that

Assuming that Equation (6) holds for all the pixels in the image, the system of equations can be written in matrix form. Denote by p.sub.i and q.sub.i the values of the target and reference pixels, respectively. Let Q be the vector {q.sub.i } and A the vector of unknowns:

where a.sub.00 compensates for the amplitude scale of the two images:

The system of equations in Equation (6) can be written in a matrix form: ##EQU2## then, ideally,

The system of equations is overdetermined, because there are many more pixels than unknown coefficients in the vector A. Using a least squares technique, A is obtained by a pseudo-inverse method:

where D.sup.t is the transpose of matrix D. Because the images are discrete, the partial derivatives are approximated by forward and backward differences:

In the following it will be shown how it is possible to exploit multiresolution to improve the registration between two images. Multiresolution improves the performance of the method of the present invention for images whose initial misalignment is rather large. Because the basic method is based on the Taylor expansion, the displacement values must be small enough for the gradient values to be relevant. At lower resolutions, adjacent pixels and local gradients represent large distances of the original image. A displacement computed on a low resolution image corresponds to a larger displacement on the highest resolution of the original image. These larger displacements may yield transformations that compensate for larger misalignments. However, those are only approximate transformations, because they are based on coarse representations of the original images.

A resolution pyramid consists of the original image and several copies of it at lower resolutions. The multiresolution registration is computed on the resolution pyramid, and the results from the low resolution transformation are used to guide the computation of the finer levels. The low resolution displacements and the corresponding transformations are scaled up according to image resolution.

Before the multiresolution method is elaborated, it first is shown how the k-th factor low resolution images are subsampled. Let I.sub.k denote the k-th factor low resolution image of I, a convolution of I with a low-pass filter having a kernel H: ##EQU3## A simple box filter H=1/k.sup.2 yields satisfactory results.

Homogeneous coordinates are used to represent points. A 2D point is represented by the triple (x,y,w), whose corresponding Cartesian coordinates are (x/w,y/w). Accordingly, an affine transformation such as Equation (3) may be represented by a matrix: ##EQU4##

Given an affine transformation B, its f-th scaled factor B.sup.f is defined as: ##EQU5## and it is easy to verify that if Bp=q, where p=(p.sub.x,p.sub.y,1).sup.t and q=(q.sub.x,q.sub.y,1).sup.t, then B.sup.f p.sup.f =q.sup.f, where p.sup.f =(f * px, f * py, 1).sup.t and qf=(f * qx, f * qy, 1).sup.t. Thus, if B.sub.k is a transformation between two low resolution images I.sub.k and J.sub.k, then B.sub.k.sup.k is an approximation of the transformation between images I and J. That is, B.sub.k.sup.k .apprxeq.B and dB.sub.k .apprxeq.dB/k, where dB.sub.k denotes the displacement values

The quality of the approximation depends on the size of the factor k because of the inevitable loss of information in low resolution images.

The multiresolution registration process starts by selecting a scale factor k.sub.0 >1, where k.sub.0 is not too large, so that the scaled down images contain enough information to calculate a meaningful, although coarse, approximation of the transformation between the two images. Denote by I and J the target and reference full scale (input) images, respectively. Let T.sub.k.sbsb.0 be the transformation between I.sub.k.sbsb.0 and J.sub.k.sbsb.0. T.sub.k.sbsb.0.sup.k.sbsp.0 is the first coarse approximation for the transformation T between I and J. Denote I.sup.k.sbsp.0 =T.sub.k.sbsp.0.sup.k.sbsp.0 (I). I.sup.k.sbsp.0 is the result of applying the first approximation to the input target image. Now the approximation is refined by applying the method to I.sub.k.sbsb.1.sup.k.sbsp.0 and J.sub.k.sbsb.1, where k.sub.1 <k.sub.0, to yield the transformation T.sub.k.sbsb.1. The transformations T.sub.k.sbsb.0.sup.k.sbsp.0 and T.sub.k.sbsb.1.sup.k.sbsp.1 are composed to yield a refined affine mapping M.sup.k.sbsp.1 between the original images. This process continues iteratively, such that at each step T.sub.k.sbsb.j is calculated and used to refine M.sup.k.sbsp.j :

Note that T.sub.k.sbsb.j is calculated between I.sub.k.sbsb.j.sup.k.sbsp.j.sup.-1 and J.sub.k.sbsb.j, where:

Because the affine transformations are closed under composition, M.sup.k.sbsp.j is an affine transformation, and I.sup.k.sbsp.j is a direct mapping of the input image, so that the iterative process does not accumulate errors. The reduction factors k.sub.i can be selected arbitrarily. Preferably , however, the binary reduction, k.sub.i =k.sub.i-1 /2, i>0, is used. This continues until the original resolution is restored. FIG. 1 is a diagram of the data flow of multiresolution registration.

The multiresolution method improves the performance of the registration in terms of the initial misalignment of the source and target images. FIGS. 2A through 2F present an example of the use of multiresolution to improve alignment. FIG. 2A is a reference image. FIG. 2B is a target image to be warped to match the reference image of FIG. 2A. Note that FIG. 2B has a small square in its upper left hand corner for reference. FIG. 2C shows the results of warping FIG. 2B using the gradient method without multiresolution. The normalized cross-correlation between FIGS. 2A and 2C is only 0.82. FIG. 2D shows the results of warping FIG. 2B using multiresolution. The normalized cross-correlation between FIGS. 2A and 2D is 0.93. This multiresolution registration was done in three steps. FIG. 2E shows the image of FIG. 2B after one multiresolution step, with a normalized cross-correlation of 0.81 with FIG. 2A. FIG. 2F shows the image of FIG. 2B after two multiresolution steps, with a normalized cross-correlation of 0.91.

The registration method described above is global in the sense that it considers all the pixels of the source and target images. The solution of the equation system is based on the least squares criterion, in which all the pixel values contribute to the solution. Some of the pixels are irrelevant to the solution, as the corresponding pixels in the other image bear no resemblance to them, or may not even exist. The "vote" of these pixels not only cannot contribute, but may even add noise to the equation system. The present invention includes a rejection process that culls away pixels that are estimated to not contribute to the correct solution. Specifically, a fast local similarity measure is used: ##EQU6## where I.sub.s and I.sub.t are the average values in some local neighborhood around I.sub.s (x,y) and I.sub.t (x,y), respectively. A pixel at (x,y) whose similarity value E(xy) is smaller than a given threshold is rejected. Determining an appropriate local neighborhood and an appropriate threshold is straightforward for one ordinarily skilled in the art.

In most real life cases it is quite rare to find a global polynomial that transforms one image to another. In most cases the transformation characteristic is more general and nonuniform, with different variations over the image. It therefore is necessary to construct the transformation by an aggregation of many local transformations, which define a general elastic transformation. The present invention finds local transformations using the above global registration method over small windows and then interpolates the local transformations to define the global transformation. The localization permits using transformations of low degree, whose combination can define a global elastic transformation of high degree. The present invention preferably uses local affine transformations, taking advantage of the fact that affine transformations are closed under composition.

However, there is a problem of working directly on small areas. A small window in one image is not the transformation of the corresponding window in the other image. This means that it is necessary to deal with partial matching. Moreover, it may occur that images of some windows are mismatched while the entire image matches. Thus, an independent registration between the windows most likely will fail. This suggests the necessity of defining a progressive process in which the registration is defined by a series of refined transformations. In the following, an hierarchical localization and a (weighted) interpolation are introduced. The process is illustrated by the example of FIGS. 3A through 3E.

The localization process starts with the definition of a global transformation between the initial two images, the reference image of FIG. 3A and the target image of FIG. 3B. The target image of FIG. 3B is mapped to an intermediate image, FIG. 3C, by this global transformation . The intermediate image of FIG. 3C and the reference image of FIG. 3A are each partitioned into four windows. Then the corresponding windows are registered. This process continues recursively in each window and is iterated down to a predefined level of localization.

At localization level k, the images are partitioned into 2.sup.k by 2.sup.k windows. Let L.sub.i,j.sup.k denote the local transformation centered at the (i,j)-window of the k-th level. First, the global transformation T.sup.0 is computed between the input images I and J. The transformed image I.sup.1 =T.sup.0 (I) and the reference image J each are partitioned into four images: I.sub.00.sup.1, I.sub.01.sup.1, I.sub.10.sup.1, I.sub.11.sup.1 and J.sub.00.sup.1, J.sub.01.sup.1, J.sub.10.sup.1, J.sub.11.sup.1, respectively. The four local transformations L.sub.i,j.sup.1 are computed between the corresponding four pairs of images. The local transformations are then composed with the global transformation to yield T.sub.i,j.sup.1, the local composed transformations:

A global transformation M.sup.1 is defined as the weighted average of these four transformations: ##EQU7## where the w.sub.i,j are weights defined below.

FIG. 4 shows the data flow of the localization process. In general, at level k, each local transformation L.sub.i,j.sup.k is composed with the upper level transformation

and the global transformation M.sup.k is defined by ##EQU8## The transformed image I.sup.k-1 =M.sup.k (I) and the reference image J are then partitioned again and the process is repeated.

The local composed transformations are combined to form a global elastic mapping by an inverse distance weighted interpolation: ##EQU9## where w.sub.i,j satisfies: 1. w.sub.i,j (x.sub.ij,y.sub.ij)=1, ##EQU10## w.sub.i,j (x,y).gtoreq.0 and (x.sub.ij,y.sub.ij) is the sub-image center. Shepard (D. Shepard (1968), A two-dimensional interpolation function for irregularly spaced data, Proc. 23rd Natl. Conf. Of the ACM, ACM Press, N.Y., pp. 517-524) suggested the following function: ##EQU11## for any (x,y).noteq.(x.sub.ij,y.sub.ij) and 1 otherwise. The Euclidean distance between points a and b is denoted by d(a,b), and the value of the exponent u determines the smoothness of the transformation. Preferably, u=2 to emphasize the local influence.

Note that the global transformation M.sup.k is neither an affine transformation nor a linear transformation, but a general elastic transformation. The computation of each local transformation L.sub.i,j.sup.k is computed with the multiresolution and culling techniques.

The example in FIGS. 3A through 3E shows the importance of localization in cases where it is necessary to register images between which no affine transformation exists, for example, in cases where the deformation is local to a small part of one of the images. The global warped image of FIG. 3C has a normalized cross-correlation with the image of FIG. 3A of only 0.17. FIG. 3D shows the results of using four local registrations (one level of localization): the normalized cross-correlation between the images of FIGS. 3A and 3D is 0.64. FIG. 3E shows the results of using sixteen local registrations (two levels of localization): the normalized cross-correlation between the images of FIGS. 3A and 3E is 0.81.

Although the preferred mode of localization is by halving the linear dimension of the subwindow at every localization level, as described above, the scope of the present invention includes all ways of subdividing the reference and target images into nonoverlapping subimages, for example, dividing the linear dimensions of the subwindow by an integer greater than 2 at every localization level.

As was mentioned above, the sampled ultrasonic volumes are deformed during acquisition. Because unavoidable deformation is inherently a 3D deformation, the warp transformation that can align the deformed volume must be a true 3D warp. Otherwise, it would be like trying to register two slices lying on different planes.

The extension of the gradient method to 3D is straightforward. All that is required is to treat the third dimension analogously to the other two dimensions. Define a 3D affine transformation: ##EQU12##

Assume, as in the 2D case, that I and J are smooth differentiable functions from R.sup.3 to R. What is needed is an affine transformation T(x,y,z)=(s,t,w) such that I(x,y,z)=J(s,t,w). Define the displacements dx=s-x, dy=t-y, and dz=w-z, so:

If the displacement vector (dx,dy,dz) is small enough to be approximated by the linear term of a Taylor expansion, then

Because

it follows that

Because it is assumed that Equation (28) holds for most of the voxels, there are obtained, as in the 2D case, a large number of equations which are solved by the pseudo inverse method, similar to the 2D case (Equation (9)). The multiresolution, culling, and localization techniques are extended and applied in a similar way. In particular, localization is done using 3D nonoverlapping subimages, preferably by partitioning the images initially into eight subimages, and, in each succeeding level of localization, partitioning each of the subimages of the previous level into eight subimages.

This formalism is easily extended by one skilled in the art to images of dimensionality higher than 3, and the scope of the present invention includes the registration of images of dimensionality higher than 3. In the localization step for images of dimensionality D, if the subimages are obtained by dividing the linear dimensions of the parent images by an integer N at every localization level, then each parent image is divided into N.sub.2D subimages at every localization level.

The gradient method defines an elastic transformation which warps the target volume to agree with the geometry of the reference volume. The method relies on the assumption that the two volumes are relatively close. Given two volumes, it is necessary to recover their relative positions. With 2D images this can be done manually, but it is almost impossible to have the same success in 3D because of the lack of proper 3D user interfaces.

A 3D rigid transformation is a function of six variables, three for translation and three for orientation. Unlike the elastic gradient method, here an indirect method, i.e., a search method, is used. Because the two volumes are of the same type, a cross-correlation can be used as a similarity measure. The search for a global maximum is carried out in a six dimensional domain

The problem of finding a global maximum for a multidimensional function is well known. This is a classic problem of optimization in which the main difficulty is to avoid converging to some local maximum. There are many solutions for non linear systems, such as Newton or Taylor techniques, which assume an initial proximity to the global solution, because they do not distinguish between global and local extrema. The preferred method in the present invention is the well established technique of "simulating annealing", which can elude local extrema. The parameters of the technique require data-dependent tuning, but this can be done easily by one skilled in the art. A simulating annealing method can also be applied to search for a non-rigid transformation. However, because search methods are sensitive to the number of parameters, the search time is impractical.

After the coarse rigid alignment has been applied, the location of the target volume in terms of the reference volume coordinate system is known, as shown in FIG. 5. This defines a clipping box contained within the overlapping portion of the two volumes. The target volume is warped by the 3D gradient method based on the context of two subvolumes, one from each volume, bounded by the clipping box. These two steps, i.e., rigid alignment and the 3D elastic warp, are the first two stages along the mosaicing pipeline illustrated in FIG. 6. Equation (22) is assumed valid throughout the target volume, not just within the clipping box. The rest of the pipeline stitches the warped volume over the reference volume.

The warped target volume has to be composed over the reference volume in a seamless way. In most cases, adjacent volumes have different intensities, mainly because of directional variations. Some parts in one volume bear no resemblance to the corresponding part of the other volume or often do not exist. This requires the application of a 2D registration of the two slices at the seam face to further refine the misalignments at the connecting planes, and to blend (cross-dissolve) the two volumes to blur away the seam faces.

The determination of the seam planes is not as simple as it may seem. The seam plane should pierce the clipping volume but not necessarily in the middle. Because we apply a 2D registration along this plane, it should contain significant information for a valid registration. A simple heuristic may be applied to look for the slice parallel to the volume face which has the most meaningful pixels, where a meaningful pixel is defined as a pixel whose close neighborhood has a variation above some predefined threshold. Note, that there may be more than one seam plane. In many cases there are three seam planes, and if the target volume is set deep into the reference volume there may be up to six seam planes.

Once the seam planes are determined, two slices are extracted, one from each volume for each seam plane. The images of the slices are registered by the 2D version of the elastic multiresolution gradient method described above. The 2D registration defines a set of transformations M.sub.i (Equation (21)) one for each seam plane. Then, the target volume is deformed by a weighted blend of the 2D transformations defined as follows. Let p be an arbitrary point in the target volume and p.sub.i the closest point (the orthogonal projection) on the i-th seam plane. Let w.sub.i be an inverse power of the Euclidean distance between p and p.sub.i : ##EQU13## where the exponent u controls the distance decay. Then ##EQU14## The blended registration function M improves the geometric registration in the proximity of the seam planes and reduces the distortions away from them. The gray levels of the two images are blended by a similar weighted sum to give the final gray level at point p.

While the invention has been described with respect to a limited number of embodiments, it will be appreciated that many variations, modifications and other applications of the invention may be made.