Hybrid Multiscale Landmark and Deformable Image Registration

An image registration technique is presented for the registration of medical images using a hybrid combination of coarse-scale landmark and B-splines deformable registration techniques. The technique is particularly effective for registration problems in which the images to be registered contain large localized deformations. A brief overview of landmark and deformable registration techniques is presented. The hierarchical multiscale image decomposition of E. algorithm is developed based on combining the multiscale decomposition with landmark and deformable techniques. Successful registration of medical images is achieved by first obtaining a hierarchical multiscale decomposition of the images and then using landmark-based registration to register the resulting coarse scales. Corresponding bony structure landmarks are easily identified in the coarse scales, which contain only the large shapes and main features of the image. This registration is then fine tuned by using the resulting transformation as the starting point to deformably register the original images with each other using an iterated multiscale B-splines deformable registration technique. The accuracy and efficiency of the hybrid technique is demonstrated with several image registration case studies in two and three dimensions. Additionally, the hybrid technique is shown to be very robust with respect to the location of landmarks and presence of noise. 1. Introduction. Image registration is the process of determining the optimal spatial transformation that brings two images into alignment with each other. Image registration is necessary, for example, when images of the same object are taken at different times, from different imaging devices, or from different perspectives. The two images to be registered, called the fixed and moving images, are the input to

our algorithm with respect to the location of landmarks, number of landmarks, and presence of noise.Concluding remarks are given in Section 7.
2. The registration problem.Image registration is the problem of determining the optimal spatial transformation between two given images.The images to be registered are typically (though not necessarily) images of the same object acquired at different times, from different perspectives, or from different modalities.Image registration has many important applications in medicine (e.g., diagnosis, analysis, and treatment planning/evaluation), astrophysics (e.g., alignment of images from different frequencies), military applications (e.g., target recognition), computer vi sion, remote sensing, and many others.See [8], [13] for an overview of applications of image registration.In this paper, we are primarily concerned with registration of medical images.
2.1.The mathematical setting.An n-dimensional gray-scale image f is a map ping which assigns to every point x ∈ Ω ⊂ R n a gray value f (x), called the intensity value of the image f at the point x.Given two images A(x) and B(x), referred to as the fixed and moving images, the goal of the registration problem is to find a spatial transformation φ such that the fixed image A(x) and the transformed moving image B(φ(x)) are similar in some appropriate mathematical way.
A registration algorithm has three components: the transformation model which specifies the way in which the moving image can be transformed to correspond to the fixed image, the similarity metric used to compare the images, and the optimization process that varies the parameters of the transformation model in such a way that the resulting transformation is optimal.Given a distance metric D on L 2 (R 2 ) (for two-dimensional images), and a specified space of transformation models, the general registration problem can be stated as follows: φ = argmin D(A, B(ψ)), where ψ is in the specified space of transformation models.Common metrics used in image registration are mean squares, normalized correlation, and mutual informa tion (for registration of multi-modality images).The space of transformation models used in a particular registration problem depends on the physical and anatomical characteristics of the images to be registered.For example, rigid and affine transfor mations are used when the moving image is assumed to differ from the fixed image by translation, rotation, dilation, and shear.Polynomial and splines translations are used when the transformation between the images is assumed to consist of local ized stretching and non-rigid deformation.See [5], [13] for an overview of the image registration problem.In this paper, we will use a combination of landmark-based and deformable registration techniques.
2.2.Landmark-based registration.Landmark-based registration algorithms use the (manual or automatic) identification of corresponding anatomical structures (or other features) in the images to be registered.For an overview of landmark-based registration, see [20].Given two images A and B to be registered with one another, let F k (A) and F k (B), k = 1, 2, . . .m denote m corresponding features of the im ages.The solution φ of the registration problem is then a map φ : R 2 → R 2 (in the two-dimensional case) such that More generally, if || • || is a norm on the features space, the landmark-based regis tration problem can be stated as the following minimization problem: Eq. ( 2), the transformations ψ are typically chosen to be an element of some finitedimensional space, such as polynomials, splines, or wavelets, and the minimization can be easily solved by expanding D LM (φ) in terms of basis functions of the space and solving the resulting least-squares problem.Landmark-based registration techniques are computationally simple and efficient, but the main drawback of such techniques is that they rely on precise correspondence of features in the images to be registered.Although there has been recent research on automatic identification of landmarks (see [2], [7], [24]), in practice landmarks are typically identified manually.Precise identification of corresponding landmarks is time consuming and tedious, even for a medical expert [20].In addition, there are numerous known examples of cases in which the result of the landmark registration process is a transformation which correctly matches the user-supplied landmarks but is not physically meaningful.See [14, p. 44] for a simple example of such a case.

Deformable registration.
There are many existing approaches to deformable registration: [1], [4], [12], [19], [23].Splines-based transformation models are among the most common and important transformation models used in non-rigid registra tion problems [6].Splines-based registration algorithms use control points in the fixed and moving images and a splines function to define transformations away from these points.The two main spline models used in image registration problems are thin-plate splines and B-splines.Thin-plate splines have the property that each control point has a global influence on the transformation.That is, if the position of one control point is perturbed, then all other points in the image are perturbed as well.This can be a disadvantage because it limits the ability of the transfor mation model to model localized deformations.In addition, the computation time required for a thin-plate spline-based registration algorithm increases significantly as the number of control points increases.See [3] for an overview of thin-plate splines.
In contrast, B-splines are defined locally in the neighborhood of each control point.Thus perturbing the position of one control point affects the transforma tion only in a neighborhood of that point.As a result, B-spline-based registration techniques are more computationally efficient than thin-plate splines, especially for a large number of control points.See [9] and [10] for a detailed description of Bsplines transformation models.The displacement of a node α i,j is specified by a vector x i,j , and the displacement vectors of a collection of nodes characterize the tissue deformation.The displacement at a particular location on the image is de duced by fitting a polynomial, expressed using the basis splines (B-splines), to the grid nodes.The number of control points used determines the number of degrees of freedom of the transformation model, and hence, the computational complexity of the B-splines algorithm.
To define the splines-based deformation model in two dimensions, we let Ω a n x × n y mesh of control points α i,j with uniform spacing δ.Then the B-splines deformation model can be written as the 2-D tensor product of 1-D cubic B-splines: where i = �x/n x � − 1 , j = �y/n y � − 1, and B l represents the l-th basis of the B-splines: Changing the control point α i,j affects the transformation only in a local neighbor hood of α i,j .The control points α act as parameters of the B-splines deformation model, and the degree of non-rigid deformation which can be modeled depends on the resolution of the mesh of control points α.A large spacing of control points allows modeling of global non-rigid deformation, while a small spacing of control points allows modeling of local non-rigid deformations.Additionally, the number of control points determines the number of degrees of freedom of the transformation model, and hence, the computational complexity.
The splines deformation model is illustrated schematically in Figure 1.The deformation is computed explicitly at control points by computing the vector dis placement between the fixed and moving images, and the splines function is used to interpolate the deformation away from the control points.The implementation of the B-splines deformation registration algorithm works in the following way: at each iteration, the distance D between the fixed and moving images is computed; the parameters of the B-splines deformation are computed based on the direction of steepest gradient descent of the metric D, and the resulting transformation φ is applied to the moving image; the distance between the fixed and deformed moving images D(A, B(φ)) is then recomputed, and this process continues until the distance D between the images is optimized.
3. The multiscale decomposition.The multiscale registration techniques to be discussed in this paper are based on the multiscale image representation using the hierarchical (BV, L 2 ) decompositions of [25].This multiscale decomposition will provide a hierarchical expansion of an image that separates the essential features of the image (such as large shapes and edges) from the fine scales of the image (such as details and noise).The decomposition is hierarchical in the sense that it will produce a series of expansions of the image that resolve increasingly finer scales, and hence include increasing levels of detail.The mathematical spaces L 2 , the space of square-integrable functions, and BV , the space of functions of bounded variation, will be used in the decomposition: Generally, images can be thought of as being elements of the space L 2 (R 2 ), while the main features of an image (such as edges) are in the subspace BV (R 2 ).The multiscale image decomposition of [25] interpolates between these spaces, providing a decomposition in which the coarsest scales are elements of BV and the finest scales are elements of L 2 .More precisely, the decomposition is given by the following.Define the J-functional J(f, λ) as follows: where λ > 0 is a scaling parameter that separates the L 2 and BV terms.This functional J(f, λ) was introduced in the context of image processing by Rudin, Osher, and Fatemi [21].Let [u λ , v λ ] denote the minimizer of J(f, λ).The BV component, u λ , captures the coarse features of the image f , while the L 2 component, v λ , captures the finer features of f such as noise.This model is effective in denoising images while preserving edges, though it requires prior knowledge on the noise scaling λ.
Tadmor, et al. proposed in [25] an alternative point of view in which the mini mization of J(f, λ) is interpreted as a decomposition f = u λ +v λ , where u λ extracts the edges of f and v λ extracts the textures of f .This interpretation depends on the scale λ, since texture at scale λ consists of edges when viewed under a refined scale.We refer to v λ = f − u λ as the residual of the decomposition.Upon decomposing f = u λ + v λ , we proceed to decompose v λ as follows: Thus we obtain a two-scale representation of f given by f ∼ = u λ + u 2λ , where now v 2λ = f − (u λ + u 2λ ) is the residual.Repeating this process results in the following hierarchical multiscale decomposition of f .Starting with an initial scale λ = λ 0 , we obtain an initial decomposition of the image f : We then refine this decomposition to obtain After k steps of this process, we have: (5) which is a multiscale image decomposition f ∼ u 0 + u 1 + . . .+ u k , with a residual v k .As k increases, the u k components resolve edges with increasing scales λ k = λ 0 2 k .4. Hybrid Landmark and deformable registration using the multiscale decomposition.In our previous work [15], [16], we used the multiscale decom position of [25] reviewed in Section 3 to develop multiscale rigid and deformable registration algorithms for registration of images that contain high levels of noise.We demonstrated that our multiscale registration techniques enable successful reg istration of images that contain noise levels significantly greater than the levels at which ordinary rigid and deformable registration techniques fail.In this paper, we extend our previous work on multiscale registration to a new registration tech nique based on combining landmark and deformable registration methods with the multiscale decomposition.We shall refer to this technique as a hybrid multiscale landmark and deformable registration algorithm.The main idea behind our new technique is that by incorporating landmark registration, we can use known infor mation about correspondence of bony structures or other anatomical structures to improve the accuracy and efficiency of the registration procedure.
For the general setup, consider two images A (the fixed image) and B (the moving image).We first apply the multiscale hierarchical decomposition to both images.Let m denote the number of hierarchical steps used for the multiscale decompo sitions.For the applications presented in this paper, we use m = 8 hierarchical scales.Upon decomposing both of the images to be registered we use the coarse scales of each image to identify several corresponding anatomical landmarks in each image.Since the coarse scales contain only the main shapes and general features of each image, the identification of landmarks in the coarse scales is a relatively easy task and can potentially be implemented automatically (for example, using the methods of [24]).The determination of which coarse scale to use (e.g., the first or second coarse scale) is image dependent; upon decomposing each image, we choose the coarse scale that contains enough general shapes to identify correspond ing landmarks but contains few details and noise.For the medical image registration applications presented in this paper it is sufficient to use the first coarse scale of each image.
To identify landmarks we must choose anatomical features of the two images that we know must correspond to one another.For example, bony structures (such as the spine or ribs in images of the torso) in one image must correspond to bony structures in another image.We then use the identified landmarks to perform a landmarkbased registration of the coarse scales of each image.The output of the landmark registration algorithm is the transformation φ landmark that optimally brings the first coarse scales C 1 (A) and C 1 (B) into spatial alignment with one another, and is thus a reasonable approximation for the optimal transformation between the original images A and B. Therefore, we use φ landmark as the starting point to deformably register the next coarse scales C 2 (A) and C 2 (B) with one another.We then use the resulting transformation φ 2 as the starting point to register the next scales C 3 (A) and C 3 (B), and iterate this procedure until the final scales are reached.To summarize, the implementation of the hybrid multiscale landmark and deformable registration algorithm is as follows: 1. Decompose the fixed and moving images to be registered using the hierarchical multiscale image decomposition of [25] reviewed in Section 3.

Coarsest Scale of C1(A)
Image A 2. Register the coarse scales C 1 (A) and C 1 (B) with each other using a standard landmark-based registration algorithm.This step allows the practitioner to incorporate known anatomical information about the images to be registered (such as correspondence of bony structures) into the registration process.Let φ denote the resulting transformation.
computed by the previous scale registration algorithm as the starting point for the current registration.

5.
Results and discussion.In this section, we demonstrate the accuracy, effi ciency, and robustness of the image registration algorithm presented in Section 4 with image registration experiments using both clinical and synthetic deformations.The images used in this section were acquired with a GE Discover-ST Scanner (GE Medical Systems, Miluakee, WI) at the Stanford University Medical Center.We obtained eight sets (phase bins) of images, each consisting of 80 two-dimensional computed tomography (CT) images (slices) of the lungs.Each phase corresponds to a different breathing phase of the respiratory cycle.Each 2D image slice con tains 512 × 512 voxels, and the slice thickness for each phase is 2.5-mm, and the eight breathing phases recorded contain approximately 400 MB of data in DICOM image format.Figure 2 illustrates the CT image slices corresponding to the first and eighth breathing phases of the respiratory cycle.5.1.2D registration of respiratory phases.In this section, we demonstrate the results obtained upon registering two corresponding sample slices (slice 16) from the first and eighth phases of the CT data set using both ordinary deformable registra tion and the multiscale hybrid registration algorithm.The slices are illustrated in Figure 3.For ease of notation, we let P1S16 denote the fixed image (on the left) and P8S16 denote the moving image (on the right).To implement the multiscale hybrid algorithm, we first decompose each of the images into m = 8 hierarchical scales, and identify four corresponding pairs of landmarks in the first coarse scales of each image, as shown in Figure 4.
Figure 5 illustrates the difference between slices before and after registration that we performed using the multiscale hybrid registration technique.For comparison, we also illustrate the difference between the slices after ordinary deformable regis tration.The visual results presented in Figure 5 clearly demonstrate the accuracy of the multiscale hybrid registration technique.The hybrid algorithm successfully recovers the deformation between the two images, and is more accurate than ordi nary deformable registration.Moreover, the hybrid algorithm is computationally more efficient than ordinary deformable registration techniques.Working on a Dell To further quantitatively evaluate the accuracy of the hybrid registration tech nique, we use the multiscale hybrid algorithm to register all 80 slices of the first breathing phase with the corresponding slices of the eighth breathing phase, and for each registration we compute the correlation coefficient ρ between the slices before and after registration.The correlation coefficient ρ(A, B) between two images A and B is given by:

� � (A
¯where A and B are m × n two-dimensional images and A and B represent the mean value of the elements of A and B, respectively.A correlation coefficient of zero indicates a low degree of matching between the images, and a correlation coeffi cient of 1 indicates exact similarity between the images.Correlation coefficients are commonly used representations of similarity between images for the evaluation of deformable registration techniques [17].In Figure 6, we plot the correlation coef ficients between the slices before registration, after hybrid multiscale registration, after multiscale deformable registration (iteratively registering the scales of the im ages with each other without using the initial landmark-based registration) [16], after ordinary deformable registration, and after landmark registration.
The correlation coefficients plot in Figure 6 quantitatively confirms the accuracy of the multiscale hybrid registration technique.For all 80 CT lung slices consid ered, the hybrid technique significantly improves the correlation coefficient between corresponding images, indicating that the algorithm successfully recovers the defor mation between the breathing phases.Moreover, the hybrid multiscale technique is more accurate (based on comparison of correlation coefficients) than all other registration methods considered here.

Large deformation registration.
In this section, we demonstrate the ac curacy of the hybrid technique for registration images that contain large localized deformations.We apply a large splines deformation to the fixed image P1S16 (slice 16, first breathing phase).The corresponding deformation field and the fixed and deformed images are illustrated in Figure 7.
As in Section 5.1, we first decompose the images to be registered into m = 8 hierarchical scales.Using the first coarse scale of each image, we identify four corre sponding pairs of bony structure landmarks, and use landmark-based registration to register the coarse scales of the images with one another.We then apply the result ing transformation to the original deformed image, and use a B-splines deformable registration algorithm to complete the registration process.Figure 8 compares the difference between the images before registration, after ordinary deformable regis tration (for comparison), and after hybrid multiscale registration.
The results presented in Figure 8 demonstrate that the hybrid registration tech nique successfully registers images that contain large localized deformations.For The correlation coefficients between all 80 image slices of the first and eighth breathing phase before registration, after hybrid multiscale registration, after multiscale deformable regis tration, after ordinary deformable registration, and after landmark registration.such images, the multiscale hybrid registration technique is significantly more ac curate than ordinary deformable registration.The improvement of the multiscale hybrid technique over ordinary deformable registration is much more pronounced in this case, indicating that the multiscale hybrid registration algorithm is particularly well-suited for registration problems involving large localized deformations.To vi sualize the deformation between the images we illustrate the deformation computed by the multiscale hybrid registration algorithm applied to a grid image in Figure 9. Table 5.2 compares the correlation coefficients between the fixed and deformed images before registration, after hybrid multiscale registration, after multiscale de formable registration, and after landmark registration.For reference, the correlation coefficient before registration is 0.70.

3D Registration of respiratory phases.
In this section, we demonstrate that our multiscale hybrid registration technique accurately registers three-dimensional images.We combine the 80 two-dimensional images from respiratory phases 1 and 8 (shown in Figure 2) to obtain two three-dimensional CT images, as illustrated in Figure 10.
To register the 3D CT images with one another, we first extend the hierarchical multiscale decomposition of [25] to 3D images.Although the multiscale decompo sition presented in [25] was done in two dimensions only, the hierarchical multiscale   1.The correlation coefficients between the fixed and de formed images after hybrid multiscale registration, after multiscale deformable registration, and after landmark registration for the large deformation test case shown in Figure 7. expansion in equation ( 5) is independent of the image dimensionality.Upon decom posing each 3D image we identify corresponding landmarks in each image.For the purposes of this paper we identify the landmarks slice-by-slice (as in the 2D case) and refer the interested reader to [24] for an algorithm to automatically identify homologous regions in each slice.We then register the 3D coarse scales with one another using the identified landmarks and apply the resulting transformation to the original moving image.This coarse-scale landmark-based registration is imple mented as a three-dimensional registration.Finally, we use an iterated multiscale three-dimensional deformable B-splines registration algorithm to complete the reg istration process.The advantage of using three-dimensional landmark-based and deformable registration procedures to register the volumes is that three-dimensional algorithms allow for transformation in three dimensions (instead of only allowing deformations restricted to a two-dimensional plane).Figure 11 compares the intensity difference between two sample slices before and after multiscale hybrid registration.A comparison of the intensity differences before and after registration demonstrates that the multiscale hybrid registration method indeed recovers the difference between the two images.Similar results were obtained with all other slices.The correlation between the images is 0.8316 before registration and 0.9745 after registration.Working on a Dell Dimension 8400 Intel Pentium 4 CPU (3.40 GHz, 2.00 GB of RAM), the total required computation time for both the 3D multiscale decomposition and multiscale hybrid registration algorithm is on the order of approximately 20 to 40 minutes, depending on the data set; this particular example required 29 minutes.
6. Robustness of the multiscale hybrid registration algorithm.In this sec tion, we demonstrate the robustness of the multiscale hybrid registration algorithm with several image registration experiments.All image registration experiments in this section are performed using the fixed and deformed images shown in Figure 7.The fixed image in each experiment is the image P1S16 (first respiratory phase, slice 16), and the moving image is the image obtained upon applying a large B-splines deformation to the fixed image.We shall consider robustness with respect to the location of the landmarks, number of landmarks, and several types of noise.
6.1.Location of landmarks.In this section, we demonstrate the robustness of the multiscale hybrid registration technique with respect to the location of the bony structure landmarks identified using the coarse-scale representations of the images to be registered.We perform 20 image registration experiments in which we randomly perturb the position of the coarse-scale landmarks in the fixed image only by 10 to 20 mm in each trial (the image sizes are all 512 × 512).In each trial, we place the landmarks in the moving image in exactly the same positions as those illustrated in Figure 4, but we vary the location of each corresponding landmark in the fixed image randomly by 10 to 20 mm, as illustrated in Figure 12.The moving image landmarks are randomly placed in the blue circles shown in Figure 12.
These experiments are designed to determine whether or not the accuracy of the hybrid algorithm is dependent on precise matching between the landmarks.Table 2 presents the mean square differences (MSDs) and correlation coefficients between the images after registration.The MSD between two images A and B is defined as follows: where N is the total number of pixels in each image, A i is the i th pixel of image A, and B i is the i th pixel of image B. Note that the optimum value of the MSD is 0, indicating exact matching between the images, while poor matches between the images result in large MSD values.For reference, MSD before registration is 0.1210 and the correlation coefficient before registration is 0.7924.For reference, we also include the MSD and correlation coefficient after hybrid multiscale registration with exact correspondence of the landmarks (i.e., no intentional perturbation).
The results presented in Table 2 demonstrate that the multiscale hybrid reg istration technique is computationally robust with respect to the location of the landmarks used in the coarse-scale landmark registration phase of the algorithm.Although the precise location of the landmarks differs in each of the 20 trials, the final registration results are essentially the same.This is because the coarse-scale landmark registration is only the first step in the hybrid registration algorithm; the multiscale deformable registration step fine tunes the result of the landmark regis tration, thus correcting any anomalies introduced in the landmark selection process.This robustness indicates that the hybrid algorithm is particularly well suited for This robustness to the location of the landmarks is one of the most notable features of our hybrid multiscale landmark and deformable registration algorithm.Most landmark-based registration algorithms require precise identification of cor responding landmarks and are thus tedious, time-consuming, and difficult to im plement, even for a medical expert.However, since the accuracy of our hybrid algorithm is not dependent on precise correspondence of the landmarks in the fixed and moving images, the identification of the landmarks can be completed quickly or implemented automatically.

Number of landmarks.
In this section, we consider the effect of the number of pairs of landmarks used in the landmark-based registration step on the accuracy and efficiency of the multiscale hybrid registration algorithm.Table 3 presents the correlation coefficients between the images after registration and total registration time (including the multiscale decomposition, landmark-based registration of the coarse scales, and deformable registration) for image registration experiments using a varying number of landmarks.The locations of the landmarks are illustrated in Figure 13.The results presented in Table 3 indicate that accurate registration can be obtained using four pairs of corresponding landmarks.6.3.Noise.In this section, we demonstrate the robustness of the multiscale hybrid registration algorithm with respect to the presence of noise in the images to be registered.Our multiscale registration algorithms [15], [16] were initially developed in the context of registration of noisy medical images.We consider three types of noise models: Gaussian, multiplicative, and impulse noise.2. The MSDs and correlation coefficients between the fixed and deformed images after multiscale hybrid registration.Each trial number represents a different perturbation of the landmarks in the fixed image.3. The correlation coefficients between the fixed and deformed images after hybrid registration and total computation time required for hybrid registration for image registration experiments using a varying number of landmarks.
where n(x) is uniformly distributed random noise with mean µ and variance v. Figure 14 illustrates the fixed and moving images with added Gaussian noise of mean 0 and variance 0.1 Table 4 presents the correlation coefficients between the noisy images shown in Figure 14 before registration, after ordinary deformable registration, and after hybrid registration.6.3.2.Multiplicative noise.Multiplicative noise, or speckle noise, is commonly ob served in medical imaging.It is defined by the following model, where s(x) denotes the actual image and f (x) denotes the observed (noisy) image: where η(0, δ) is uniformly distributed random noise of mean 0 and variance δ.
where n(x) is uniformly distributed random noise with mean µ and variance v. Figure 15 illustrates the fixed and moving images with added multiplicative noise of mean 0 and variance 0.  5.The correlation coefficients between the noisy fixed and moving images shown in Figure 15 before registration, after ordi nary deformable registration, and after multiscale hybrid registra tion.
Table 5 we presents the correlation coefficients between the noisy images shown in Figure 15 before registration, after ordinary deformable registration, and after hybrid registration.6.3.3.Impulse noise.Impulse noise, or salt and pepper noise, is noise that resembles salt and pepper granules randomly distributed over the image.Impulse noise is typically defined by the following model, where s(x) denotes the actual image and f (x) denotes the observed (noisy) image: Table 6.The correlation coefficients between the noisy fixed and moving images shown in Figure 16 before registration, after ordi nary deformable registration, and after multiscale hybrid registra tion.
� s(x), with probability 1 − δ, f (x) = (8) η(x), with probability δ, where η(x) is an identically distributed, independent random process which sets corrupted pixels alternatively to zero (black) or one (white); unaffected pixels re main unchanged.An arbitrary pixel is affected by noise with probability δ, and not affected with probability 1 − δ.We refer to δ as the impulse noise density, since adding impulse noise of density δ to an image f (x) affects approximately δ • size(f ) pixels. Figure 16 illustrates the fixed and moving images with added impulse noise of density 0.2 Table 6 presents the correlation coefficients between the noisy images shown in Figure 16 before registration, after ordinary deformable registration, after multiscale deformable registration, and after hybrid registration, as well as the total computation time required for each registration.6.4.Discussion.The results presented in this section demonstrate the accuracy, efficiency, and robustness of the hybrid multiscale landmark and deformable reg istration technique for registration of medical images in both two and three di mensions.In particular, the hybrid technique is shown to be more accurate than ordinary landmark registration, ordinary deformable registration, and multiscale deformable registration, especially in the case of registering images that contain large localized deformations and/or noise.Because the accuracy of the hybrid technique is not sensitive to precise identification of the location of correspond ing landmarks, the hybrid registration method is a significant improvement over ordinary landmark-based registration methods, which are known to be dependent on exact spatial correspondence of landmarks.Additionally, the hybrid technique accurately registers deformed images that contain noise levels large enough that ordinary registration techniques fail to produce meaningful results.
The development of our hybrid multiscale technique was motivated by the ob servation that the correspondence between some of the regions in the moving image (such as bony structures) can be easily identified visually, while the correspondence between other regions (such as those that contain tissue deformation, breathing movement, or lack of distinct image features) is less obvious.Thus, rather than approach the mapping of the two types of regions equally, as with ordinary image registration algorithms, mapping should be approached separately for each region.Our proposed technique combines landmark registration and deformable registra tion so that prior knowledge about the correspondence between the images, such as visual identification of corresponding landmarks, can be incorporated into the deformable registration process.The use of the coarse scales (which contain only the main shapes and general features of the images) obtained via the hierarchical multiscale image decomposition of [25] enables quick and easy identification of cor responding landmarks, and the deformable registration component of the algorithm fine-tunes the registration result produced by the coarse-scale landmark registration.7. Conclusions.In this paper, we proposed a novel hybrid two-step image reg istration technique combining the hierarchical multiscale image decomposition of [25] with landmark and deformable registration methods.The hybrid approach has several major advantages over ordinary registration algorithms.It allows the practi tioner to incorporate a priori knowledge of corresponding bony or other anatomical structures into the registration process using a coarse-scale landmark registration.Because the coarse scales of an image contain only its main shapes, the coarse-scale landmark registration phase of the hybrid algorithm is easy to implement and is computationally efficient.The transformation produced by the coarse-scale land mark registration is used as the starting point for a deformable registration, which significantly decreases the computation time required for convergence of the de formable registration algorithm.The hybrid method was applied to both two-and three-dimensional deformable image registration problems, and our results demon strated that the hybrid registration is more accurate that ordinary landmark and deformable registration methods, and that the technique significantly improves the accuracy of registration of images that contain large localized deformations.The hybrid algorithm is very robust with respect to the location of the landmarks used in the coarse-scale landmark registration phase of the algorithm, indicating that the accuracy of the technique is not dependent on precise spatial correspondence of landmarks.Compared to ordinary landmark registration, this method significantly reduces the time required for the difficult and tedious task of landmark selection.The hybrid method is also robust with respect to the presence of several different types of noise, and greatly improves the accuracy of registration of images that contain noise or other artifacts.
Appendix A. Implementation of the multiscale decomposition.As de scribed in [25], the initial scale λ 0 should capture the smallest oscillatory scale in f , given by ) 2λ 0 λ 0 where W −1,∞ is the Sobolev space with norm given by: � However, in practice we may not be able to determine the size of ||f || W −1,∞ , so we determine the initial choice of λ 0 experimentally.Fol lowing [25], for the applications presented in this paper, we will use λ 0 = 0.01 and λ j = λ 0 2 j .
We follow the numerical algorithm of [25] for the construction of our hierarchical decomposition.In each step we use finite-difference discretization of the Euler-Lagrange equations associated with the J(v j , λ j+1 ) to obtain the next term, u j+1 , in the decomposition of the image f .Due to the singularity when |�u λ | = 0, we replace J(f, λ) by the regularized functional u+v=f Ω and at each step, we find the minimizer u λ of J � .The Euler-Lagrange equation for with the Neumann boundary conditions: where ∂Ω is the boundary of the domain Ω and n is the unit outward normal.We k thus obtain an expansion f ∼ � u j , where the u j are constructed as approximate j=0 solutions of the recursive relation given by the following elliptic PDE: To numerically implement the method, we cover the domain Ω with a grid (x i := ih, y j := jh), and discretize the elliptic PDE of Eq. ( 12) as follows:.

�
where D + , D − , and D 0 denote the forward, backward, and centered divided dif ferences, respectively.To solve the discrete regularized Euler-Lagrange equations (13), we use the Gauss-Siedel iterative method to obtain: To satisfy the Neumann boundary conditions (11), we first reflect f outside Ω by 0 adding grid lines on all sides of Ω.As the initial condition, we set u We i,j = f i,j .n ∞ iterate this numerical scheme for n = 0, 1, . . .N until ||u − u n ∞ −1 || is less than n ∞ some preassigned value so that u i,j is an accurate approximation of the fixed point steady solution u λ .
n ∞ Finally, we denote the final solution u λ := {u } i,j .To obtain the hierarchical i,j multiscale decomposition, we reiterate this process, each time updating f and λ in the following way: That is, at each step, we apply the J(f current − u λ , 2λ) minimization to the residual f current − u λ of the previous step.Taking λ j = λ 0 2 j , we obtain after k steps a hierarchical multiscale decomposition f = u λ 0 + u λ 1 + . . .+ u λ k + v λ k , where we write u λj = u j .We call the u j , j = 1, 2, . . ., k the components of f , and the v k the residuals.For ease of notation, given an image f , we let C k (f ) denote the k th scale of the image f , k = 1, . . ., m: Thus C k (A) will denote the k th scale of the image A, and C k (B) will denote the k th scale of image B.
A.1.Three-dimensional implementation.To implement the iterated multiscale decomposition in three dimensions, we cover the image domain Ω with a grid (x i := ih, y j := jh, z k := kh), and let D + , D − , and D 0 denote the forward, backward, and centered divided differences, respectively.Then the 3D extension of the iterated multiscale decomposition given by equation ( 14) in Section A is:

Figure 1 .
Figure 1.Schematic visualization of the splines deformation model.

3 .
Use φ landmark as the starting point to deformably register the next scales C 2 (A) and C 2 (B) with each other.Let φ 2 denote the resulting transformation.4. Use φ 2 as a starting point to deformably register C 3 (A) and C 3 (B) with each other.Let φ 3 denote the transformation obtained upon registering C 3 (A) with C 3 (B).Iterate this method, at each stage using the transformation

Figure 2 .
Figure 2. The CT image slices (80 per phase) corresponding to the first and eighth breathing phases of the respiratory cycle.

Figure 3 .
Figure 3. Two corresponding sample slices from two breathing phases of the same patient.The image on the left (denoted P1S16) is from the first breathing phase, and the image on the right (de noted P8S16) is the corresponding slice from the eighth breathing phase.

Figure 4 .
Figure 4.The coarse-scale landmarks used to implement the landmark-based registration portion of the hybrid algorithm.

Figure 5 .
Figure 5.The difference between the CT respiratory phase slices P1S16 and P8S16 before registration, after ordinary deformable registration, and after hybrid multiscale registration.

Figure 7 .Figure 8 .
Figure 7.The deformation field (left) corresponding to the de formation transformation between the fixed (center) and deformed (right) images to be registered with one another.

Figure 9 .
Figure 9.The original and deformed grid image.The deformed grid on the right is obtained by applying the deformation computed by the multiscale hybrid registration algorithm to the grid image on the left.

Figure 10 .
Figure 10.The 3D CT volumes obtained upon combining the 80 2D CT phase-binned images corresponding to the first (left) and eighth (right) respiratory phases.

Figure 11 .
Figure 11.The difference between two corresponding slices of the 3D CT images before and after multiscale hybrid registration.

Figure 12 .
Figure 12.The perturbed locations of the coarse-scale landmarks in the moving images.In each trial, the landmarks in the moving image are randomly placed in the blue circles.

6. 3 . 1 .
Gaussian noise.Gaussian noise is an independent additive noise model in which the observed (noisy) image f (x) is the sum of the true image s(x) and the noise n(x):

Figure 13 .
Figure 13.The locations of the landmarks used in the experi ments presented in Table3.

Figure 14 .
Figure 14.The fixed and deformed images with Gaussian noise of mean 0 and variance 0.1.

Figure 16 .
Figure 16.The fixed and deformed images with impulse noise of density 0.2.

Table 4 .
2The correlation coefficients between the noisy fixed and moving images shown in Figure14before registration, after ordinary deformable registration, and after multiscale hybrid registration.Figure 15.The fixed and deformed images with multiplicative noise of mean 0 and variance 0.2.