Patent Application
20060034500
The present invention relates to a method of computing the transformation for transforming two images, in particular medical MR- or CT-images of a patient, one into the other. Moreover, the present invention relates to a corresponding device and to a computer program for implementing said method.
In numerous medical and non-medical applications of imaging methods a problem is encountered in that two images formed of the same object have to be analyzed as regards which elements in the images correspond to one another and how these elements have shifted and/or have become deformed from one image to the other. Such comparisons of two images are intended notably for the analysis of flexible objects and shapes. In the context of the automatic analysis of pairs of images, that mathematical transformation is computed which transforms the two images one into the other. Such a transformation may also be referred to as a motion and deformation field, since it indicates how every point of the first image has moved in the other image or how a surface element or volume element of the first image has become deformed in the second image. In the case of medical applications, the local distribution of the deformation or motion in an image can be directly visualized so as to support the diagnosis of, for example the growth of a tumor. Furthermore, the deformation can be used as the basis for cardiac applications, for a comparison of images formed before and after a treatment, and for the compensation of motions of a patient. From the deformation field, local elastic properties can be deduced. These elastic properties can reflect pathologies, for instance rigid tumors in soft environment.
Various methods have been proposed for the automatic computation of a motion and deformation field which transforms two different images of the same object one into the other. Many algorithms attempt to establish a correspondence between parting edges (J. Declerck, G. Subsol, J.-P. Thirion, N. Ayache, Automatic retrieval of anatomical structures in 3d medical images, Lecture Notes in Computer Science 950 (1995), p. 153) or boundaries (D. Davatzikos, J. L. Prince, R. N. Bryan, Image Registration Based on Boundary Mapping, IEEE Trans. Medical Imaging 15 (1996), p. 112) that have been extracted in both images. According to these methods, however, only a part of the information contained in the images is used. Therefore, the methods yield transformation parameters only for the lines of the selected characteristics and the transformation must subsequently be expanded to cover the entire volume.
An alternative approach consists in elastic registration that is based on gray values and utilizes the entire image contents. The calculation of a motion and deformation field means that a respective set of transformation parameters must be assigned to each sub-volume of the image (or even to each voxel). The calculation, therefore, is actually an optimization method for determining the field that produces maximum similarity between the images. A motion and deformation field thus represents a very large number of degrees of freedom. Overall optimization schemes for the elastic registration change the overall motion and deformation field in every step of the computation and hence also change a large number of parameters. A large number of local optima is to be expected during these methods and, therefore, in many cases the optimization method becomes stuck in one of the local optima instead of reaching the desired overall optimum.
In another method, that is, the so-called block matching, the image is subdivided into small sections which are separately assigned (P. J. Kostelec, J. B. Weaver, D. M. Healy, Jr., Multiresolution Elastic Image Registration, Med. Phys. 25 (1998), p. 1593). The method, illustrated on the basis of 2D images and proposed for 3D images, commences with a rigid registration of the overall image. Only a shift (translation and rotation of the image) takes place during a rigid transformation, but not a deformation, that is, a location-dependent expansion or compression of lengths. After the rigid registration, the image is successively subdivided into ever smaller sections that are then rigidly assigned again. Each section utilizes the registration parameters found in the “parent block” as initial estimated values for the own registration. Even though this method yields acceptable results for 2D images with comparatively small deformations, in 3D images not even an approximately rigid registration of image sections that are larger than approximately 1 cm is possible in the presence of large deformations. The quality of the initial estimated values, therefore, is very likely inadequate for successful registrations in the last steps of the algorithm where a given overlap is required between small corresponding blocks in the source image and in the target image.
Particularly in the case of global optimisation schemes, a large number of parameters, e.g. the positions of some 104 spline control points for typical 3D volumes, have to be varied. This results in a large number of local optima of the similarity measure and in a large number of function evaluations per optimisation step. In order to achieve tolerable computation times and to increase robustness, multi-resolution approaches are used. However, re-sampling the image at courser resolutions results in blurring of tissue boundaries and might result in unrealistic deformation fields.
In case of non-rigid registration algorithms an optimisation of the deformation field for each pixel or voxel of the reference image, e.g. about 134-106 voxels for a 5123 image, would take an unacceptable amount of time, so that the calculation of the deformation field is usually based on a set of sub-volumes, as proposed for block-matching, on control points of a continuous function, e.g. splines, or on the parameters of a global polynomial transformation. The continuity of these functions can successfully describe deformations where the spatial variation of local transformation parameters is smooth. Examples are deformations of the female breast. However, if anatomical entities move relative to each other, e.g. the liver moving relative to the ribs due to respiration, the real local transformation parameters are not continuous at the liver boundary. Current methods interpolate the dense deformation field at an image position from all control points close to this position even if the control points are located in different organs.
U.S. 2002/0054699 A1 discloses a method of computing the transformation which transforms two different images of an object one into the other while describing the motion or deformation of the object. In accordance with the method first the local transformation parameters are computed for sub-regions. By choosing the sub-regions so as to be sufficiently small, a rigid transformation can be used for this purpose. Starting from at least one first predetermined sub-region, the local transformation parameters of the neighboring sub-regions are successively computed, the starting values on which each computation is based being the already determined local transformation parameters of neighboring sub-regions. Using the local transformation parameters thus determined for sub-regions, the overall transformation can subsequently be computed. This template propagation procedure implicitly assumes a continuous variation of local transformation parameters, an assumption that is not justified in cases where anatomical entity is moved relative to each other so that discontinuities in the displacement field turn up.
Considering the foregoing it is an object of the present invention to provide a method of computing the transformation for transforming two images, in particular of the same object one into the other which avoids the above described problems, in particular avoids unwanted interpolation artifacts for points located close to tissue boundaries, increases computational efficiency and accuracy and reduces the required time for registration of the images.
This object is achieved according to the present invention by a method as claimed in claim 1 comprising the steps of:
A corresponding device according to the present invention is defined in claim 9. Furthermore, the present invention relates to a computer program for implementing the method as claimed in claim 10. Preferred embodiments of the invention are defined in the dependent claims.
The method according to the invention is intended to compute the mathematical transformation that transforms two different images one into the other. The images may notably be medical images which have been acquired, for example by means of an X-ray computed tomography (CT) apparatus or by means of a magnetic resonance (MR) tomography apparatus. In this context a transformation between two images is to be understood to mean a function which assigns the points of one image to the points of the other image while leaving the neighbor relations of the points unchanged. Thus, a transformation of this kind is a continuous function which is also referred to as a motion and deformation field, because it describes the motion of the points of the images from one image to the other in relation to the deformation of surface elements or volume elements. The function is preferably objective, so that it assigns each point of one image reversibly unambiguously to a point of the other image.
The present invention is based on the idea to add a new processing step in the context of non-rigid registration of a 2D or 3D target image to a reference image. This additional processing step shall be referred to as control point clustering and is performed as explained in the following.
Given a set of point pairs pi,ref and pi,trans (i=1, . . . , N) in the reference and target image respectively, for each of the points in the reference image, a set of transformation parameters ti is first derived by global optimization steps or by a block matching method. Rather than treating each point on its own, neighboring points are grouped together in clusters Cj (j=1, . . . , m) such that all points of a cluster have similar transformation parameters:
Pi,refεCj:N(ti, Tj)<s
where Tj denotes the transformation parameters assigned to the cluster Cj, the deviation is measured by a norm N and the maximum deviation from Cj is given by threshold value s. For example, Tj can be assigned to the average transformation parameters of all cluster points and N can be the length of the difference vector ti-TJ.
Not all points pi,ref need to belong to a cluster; one point can belong to more than one cluster. This is required to treat irregular motion patterns, e.g. in the lumen of MR cardiac images, and to avoid binary cluster borders which would be heavily affected by noise and inaccuracies in the values of transformation parameters.
In known interpolation methods the transformation parameters of a voxel at position rref is influenced by the properties of all control points pi,ref in a certain area around r where the contribution of each control point depends only on the distance |r-pi|. This often gives plausible values for positions at the centre of organs. However, the treatment of points closed to tissue boundaries is more problematic as the interpolation for these positions is based on control points belonging to different anatomical structures. This problem appears for the case of respiratory motion compensation where inaccuracies in the liver area turn up because control points positioned on the ribs are used to calculate the deformation field within the liver. This leads to inaccuracies because during inspiration the ribs follow the expansion of the chest while the motion of the liver is dominated by the movement of the diaphragm. The proposed method addresses this difficulty by imposing an additional cluster-based classification step during interpolation as defined in steps d1-d3:
For control points r belonging to more than one cluster, the transformation parameters are derived in two steps. First, for each cluster containing the control point, transformation parameters at the position r are deduced. This results in k sets of intermediate transformation parameters. In a second step, these intermediate transformation parameters are combined for determining the real transformation parameters (step d3).
Applying this method during global optimisation leads to a dynamically updated subdivision of the image into regions characterized by similar or substantially identical transformation properties. Applied after optimisation, the clustering step avoids unwanted interpolation artifacts for control points located close to tissue boundaries.
According to a preferred embodiment the transformation parameters in step d3 are determined by a combination and weighting of the intermediate transformation parameters. Further, in order to assign a control point to one or more clusters it can be checked, as proposed according to a further preferred embodiment, whether the control point lies within the convex hull of the cluster under consideration. Another possibility is to use the distance between the control point and the closest point of the cluster as a criterion.
According to another embodiment the transformation parameters are determined by a selection of one of the intermediate transformation parameters based on a similarity measure of the image information belonging to the control points under consideration. The similarity measure therein indicates, preferably based on image information, to which image object or to which cluster a particular control point should belong.
According to another aspect of the invention the clustering step is repeated for optimisation of the clusters in case a control point has been assigned to a particular cluster in previous step d. Thus, the clustering can be optimised.
The proposed clustering can be incorporated into a template propagation method as, for instance, described in P. Rösch, T. Netsch, M. Quist, G. P. Penney, D. L. G. Hill and J. Weese, “Robust 3D Deformation Field Estimation by Template Propagation”, vol. 1935, pages 521-530, MICCAI, Springer, 2000. In this case, clustering is performed dynamically during the propagation process. Rather than propagating starting estimates to neighbouring templates, the hypotheses that a set of templates forms a cluster is tested on the basis of local transformation parameters, and the clustering information is taken into account during propagation to avoid that starting values derived from templates belonging to one anatomical entity are used for another anatomical structure. An advantage of this method is that in addition to a set of correspondences a segmentation of the image in clusters is performed. This segmentation information can be refined afterwards.
According to another preferred embodiment control points can be interactively assigned to clusters by a user. Particularly in cases where the automatic clustering as inaccurate this will be advantageous and will improve accuracy of the method. In cases where it is not clear to which a control point belongs the method can be adapted to ask the user to make the decision or to ignore the particular control point.
From a conceptual point of view, the proposed invention and the advantages obtained by performing a clustering during non-rigid registration can be seen as an attempt to support image segmentation by the results of elastic registration.
The invention also relates to a computer program for computing the transformation that transforms two digitized images of an object, preferably two medical images acquired by means of computed tomography or magnetic resonance tomography, one into the other. The computer program is characterized in that it carries out a computation in conformity with one of the methods described above.
The invention also relates to a data carrier for a computer program on which a computer program of the kind set forth is stored. The data carrier may notably be a magnetic data carrier (disc, magnetic tape, hard disc), an optical data carrier (CD), a semiconductor memory (RAM, ROM . . . ) or the like. The data carrier may notably form part of a computer in which the computer program stored on the data carrier is executed.
Finally, the invention relates to a device for computing the transformation which transforms two digitized images of an object, preferably acquired by means of computed tomography or magnetic resonance tomography, one into the other. The device comprises a central processing unit and at least one memory unit with which the central processing unit is connected and to which it has access for reading and writing of data and commands. The memory unit may especially store the images to be transformed as well as a computer program to be executed by the central processing unit. The memory unit may notably be a magnetic data carrier (disc, magnetic tape, hard disc), an optical data carrier (CD), a semiconductor memory (RAM, ROM . . . ) or the like. The program that is stored in memory and that controls the central processing unit is adapted to calculate the transformation on the central processing unit by a method as it was explained above, i.e. the central processing unit executes the steps a) to d) as explained above. Moreover, the above-mentioned improvements of the method may be implemented in the computer program.
The invention will now be explained in more detail with reference to the drawings in which
As shown in
In a subsequent step, which is illustrated in
Thereafter a clustering is performed as shown in
Subsequently, as illustrated in
In case of control point 6, which is part of cluster C2, transformation parameters of neighboring control points which are part of the same cluster C2 are used for interpolation, i.e. transformation parameters t3 and t4 of control points 3, 4 are used for interpolation of transformation parameters t6 describing the motion of control point 6 in the reference image 10 into control point 6′ in the target image 10′.
For control points which belong to more than one cluster, additional criteria have to be considered. The situation is shown in
Under the assumption that control point 7 belongs to cluster C, the position 7b′ would thus be found, and the transformation parameters t7b would be selected or determined as transformation parameters while a conventional interpolation method would result in position 7a′ and corresponding transformation parameters t7a. Thus, the above described problem at tissue boundaries can be solved.
A particular application of the proposed method is to support the optimization procedure for deformation field estimation if global optimization schemes are used. In that optimization procedure in a first step a set of control points as described above is initialized. Thereafter an optimization step is performed which leads to a modified control point distribution yielding a larger numerical value of the applied similarity measure, e.g. mutual information or local correlation. Thereafter a clustering of the corresponding points is performed. Therein, constraints based on “cluster membership” are imposed for the following optimization steps. For example, the parameter variation in the context of simple gradient optimization scheme can be synchronized with the cluster to speed up convergence and to improve the stability of the optimization. For Newton-like methods, common (or average) Hessian approximations for the cluster members can be calculated to improve the statistical significance of the Hessian particularly for noisy images.
By incorporating the derived constrains the next optimization step is performed, thus leading to an iteration. If the termination criterion, e.g. a predetermined value for the similarity measure, for the optimization is reached, the resulting parameters are stored, otherwise the optimization is continued.
The proposed method can also be used to implement an alternative to current multi-resolution approaches. Rather than reducing the image resolution or the spacing between control points from level to level, the method starts with the large threshold which results in large clusters. Then the threshold is successively set to smaller and smaller values thus implementing a coarse to fine procedure in parameter space rather than in the spatial domain.
For global elastic registration algorithms the proposed method leads to an improvement of the computational efficiency by reducing the dimensionality of search space. This is achieved by varying transformation parameters of clusters of control points rather than optimizing the parameters of each control point individually. Particular large, high resolution data sets resulting from recent developments in CT will benefit from this speed improvement. The robustness of the optimization process is improved as the influence of local deviations due to registration errors is reduced. In contrast to current biomechanical models, no time-consuming segmentation of the images is required. As further advantage more realistic deformation fields will be obtained minimizing the artifacts current schemes produce at tissue boundaries. Anomalies in the deformation pattern that indicate pathology can be detected, and the attention of the user can be directed to image areas where these deviations occur, e.g. by suitable color coding. During template propagation a dynamic clustering procedure will prevent starting estimates from one anatomical region to be propagated to other regions which show a different motion. This will increase the robustness of the method, e.g. at organ surfaces.
It should be noted that the images 10, 10′ may not only be images of the same object but can also be images of different objects, as for instance in the field of inter-patient registration. Further, one image can be an image of an object while the other image can be an image of a corresponding model of the object, as for instance in the field of model-template registration which is often used for the recognition of objects in an image.
While above registration of two images is illustrated, the invention is not limited to the registration of exactly two images. Nowadays there is an increasing tendency to acquire and analyse series of images to study motion patterns. The clustering is then not necessarily done on transformation parameters between two images, but may conceptually also be applied to transformation patterns over time, e.g. cardiac series for wall motion analysis where there is a repetitive motion pattern with the heart cycle (heart embedded in surrounding tissue with different motion characteristics) or lung series for respiration analysis/lung mechanics where there is a repetitive motion pattern with the breathing cycle.
Further, the invention is not limited to elastic registration by means of corresponding point-sets obtained by template registration and an interpolating scheme, which method is illustrated above referring to the figures. Other elastic registration schemes use e.g. a deformable grid rather than corresponding control points. These can also benefit from the present invention: Deformable grids should become more flexible in areas between clusters, otherwise the registration will not converge and finally result in unrealistic transformations. If the critical areas are known, which knowledge can be gained by the present invention, these grids can be made more flexible by either increasing the number of control points on the grid, only in critical areas rather than globally, or by allowing higher bending energies in those areas, i.e. by allowing more flexible deformations.
| Number | Date | Country | Kind |
|---|---|---|---|
| 02079775.9 | Nov 2002 | EP | regional |
| Filing Document | Filing Date | Country | Kind | 371c Date |
|---|---|---|---|---|
| PCT/IB03/04780 | 10/27/2003 | WO | 5/12/2005 |
