Abstract
Nonrigid image registration is widely used to estimate tissue deformations in highly deformable anatomies. Among the existing methods, nonparametric registration algorithms such as optical flow, or Demons, usually have the advantage of being fast and easy to use. Recently, a diffeomorphic version of the Demons algorithm was proposed. This provides the advantage of producing invertible displacement fields, which is a necessary condition for these to be physical. However, such methods are based on the matching of intensities and are not suitable for registering images with different contrast enhancement. In such cases, a registration method based on the local phase like the Morphons has to be used. In this paper, a diffeomorphic version of the Morphons registration method is proposed and compared to conventional Morphons, Demons, and diffeomorphic Demons. The method is validated in the context of radiotherapy for lung cancer patients on several 4D respiratory-correlated CT scans of the thorax with and without variable contrast enhancement.
1. Introduction
In the context of image-based medical diagnostics and treatment, highly deformable anatomies are a problem for multiple-time imaging analysis along the course of treatment. Indeed, a precise tracking of organs is made difficult because of shape and position variations. Nonrigid registration may be used to compute a displacement vector for each voxel of an image [1], enabling the estimation of the spatial variations of the anatomy. The displacement vectors are computed as pointing to the best corresponding location of the voxels in another image according to a metric which is a measure of the image matching and under some constraints on global properties of the resulting deformation, such as invertibility and smoothness.
Several registration methods have been used in the past years to estimate deformations in highly deformable anatomies [2–5]. Many efforts have been made to improve the quality of displacement estimates and also to reduce the amount of required preprocessing or modeling and improve registration speed [6–8]. Besides, the choice of a registration method for medical application depends on the characteristics (e.g., modality) of the images to be registered [1]. The existing methods [9] can be divided into parametric, or model-based, methods (B-splines [10], thin-plate splines [11], radial basis functions [12], linear elastic FEM [13], etc.) and nonparametric methods (viscous fluid [14], optical flow [15], etc.). In this second category, the algorithm called Demons [16, 17] is fast, efficient, and easy to use, as it requires no particular preprocessing nor patient-specific modeling. This method aims at calculating a regular displacement field which produces a good matching of the intensities in both images by minimizing a metric, such as the sum of squared differences (SSDs) [18] or the mutual information (MI) [8] between images along with a measure of the field regularity.
In a growing number of applications, the displacement fields resulting from registration are used to deform images from other modalities or other spatial distribution maps (e.g., the dose map associated to CT scans in radiotherapy [19, 20]). Therefore, the matching of structures in images based on their intensities is not a sufficient constraint for producing realistic anatomical deformation estimations [21]. This is the reason why a priori information on the physical characteristics of anatomical deformations has to be included in the registration process. Diffeomorphism is a necessary condition for displacement fields to be physical [22]. Indeed, organs can be compressed and deformed, but cannot undergo noninvertible spatial transformations, for example, showing mirror effects. A method has been proposed in [23] to limit the displacement fields computed by the Demons to a set of diffeomorphic transformations, using diffeomorphic flows and Lie algebra.
In several medical protocols, contrast agents are used in order to facilitate interpretation. This makes the registration problem incompatible with the hypothesis of intensity conservation. Furthermore, an histogram equalization is often not able to correct for contrast agent variability, as different regions will be enhanced in different ways inside the image. Therefore, simple metrics, such as SSD or cross-correlation, are not suitable for matching those images, and methods that are suitable for registering variable contrast images have to be investigated [24, 25].
A method similar to Demons but using a phase-based approach was first proposed in [26] and was called Morphons. The principle of the method is to match transitions (between dark and bright zones) rather than intensities, by looking locally at the spatial oscillations in intensities. This method uses Gaussian smoothing as regularization of the displacement field and additive accumulation during the iterative process. This is nevertheless not sufficient to ensure the invertibility of the deformation [22, 27].
In this paper, a Morphons registration using a diffeomorphic accumulation step is proposed and its accuracy is assessed in the case of thorax image registration, also in presence of different contrast enhancements, and compared to the Demons. The paper is organized as follows. In Section 2, the main mathematical concepts and definitions are presented. Then, in Section 3 a generic nonparametric registration process is presented, and its particularization to Morphons and to diffeomorphisms is proposed in Section 4. In Section 5, different registrations are applied on images of the thorax, without contrast enhancement in the first experiment and with contrast enhancement in the second. The results of these experiments are eventually discussed in Section 6.
2. Mathematical Framework
For the sake of clarity, let us introduce some key mathematical concepts used throughout this paper.
2.1. Images and Deformation Fields
In this paper, we always denote 3D images by lower case letters. For instance, in the process of estimating a displacement field, the fixed and the moving images are written and , respectively. We consider them as real valued functions on the volume of points , that is, . Most of the time, these functions, and also the continuous operations performed on them, such as convolutions or integrals, must be understood as approximated on the discrete voxel grid , omitting the treatment of volume boundaries. In this study, image convolutions were performed using zero-padding outside the boundaries.
A displacement field on is a vectorial field . It is associated to the “deformation” operation , that is, , with the identity deformation: . The operation , and by extension its vector field , is said to be diffeomorphic if it is invertible, differentiable, and its inverse is differentiable. For the transformation to be invertible, its Jacobian must not vanish in any point , that is, if , with . Moreover, it has to be positive (). Indeed, a transformation with negative Jacobians does not correspond to physical deformations (as the mirror operation).
Mathematically, given the images and , we will see that our global objective of our study is to estimate such that the warping of by is “close” to , that is, with the common function composition. We will use sometimes the notation to insist on the warping action of on . By extension, this warping symbol can also be used on vector fields themselves, for example, for two displacement fields and , .
In practice, the warping is applied on discrete images. The transformation might therefore need to be truncated (on the volume boundaries) to the closest point inside the volume in order to avoid extrapolation of the images to be warped.
2.2. Compositive Accumulation
In this paper, we promote a particular way to combine, or accumulate, properly two displacement fields and . Adding them to form (as performed by many nonparametric registration methods; see Section 3) is of course computationally efficient, but it breaks the consistency with the composition of the corresponding spatial transformations, as illustrated in Figure 1.
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
| (f) |
| (g) moving image |
| (h) |
| (i) |
Indeed, one can clearly see in Figure 1(h) that the warping of the image in Figure 1(g) by the sum of two diffeomorphic fields, and does not correspond to the successive warping of this image by and then by , which is represented in Figure 1(i).
However, the compositive operation, denoted by , solves this issue. It is simply defined as By construction, the deformation operation linked to the deformation field is therefore . If both displacement fields are diffeomorphic, their composition is also diffeomorphic [28].
The operation has some interesting and useful properties. First, the neutral accumulation is of course obtained with the null displacement field, that is, . Second, it is easy to prove the associative relations for three displacement fields , , and . And finally, and are linked through the simple relation meaning that the displacement field is equivalent to summing the field with the field warped by . This is illustrated in Figure 1: the vectors in Figure 1(i), corresponding to the successive warping by and then , are the sum of the vectors in Figures 1(a) and 1(c), as shown in Figure 1(f).
2.3. Diffeomorphic Flow and Exponentiation
An important notion used in Section 4.2 is the concept of (continuous) diffeomorphic flow [27, 29, 30]. Given a point and a smooth vector field , the flow is the dynamic solution of the following (autonomous) ordinary differential equation: At a given “time” , the position is simply a point on the trajectory following tangentially from the initialization on (see Figure 2). Following [27], the exponential of a vector field , that is, , is the nonlinear deformation operation obtained by the flow of at time , that is, . Interestingly, this exponential map acts as the common scalar-valued exponential, that is, for , and it is invertible by simply considering the inverted vector field, that is, . In addition, for differentiable , is also a diffeomorphism on . In other words, modifies the 3D coordinates with no intersection between the motions of points. Indeed, such a possibility would induce a point with two different motion vectors, a situation that is forbidden by (4) since is uniquely defined.
2.4. Scaling and Squaring
A numerical scheme exists to compute approximately but efficiently when belongs to a regular grid of voxels . Indeed, when the field is close enough to zero (i.e., ), the exponential of the field can be approximated using the first-order Taylor expansion , that is, by the transformation itself. On the other hand, the solution of the flow equation (4) in can be approximated by “discretizing” between 0 and 1. Indeed, as (where the exponent expresses the number of times the deformation operation is combined with itself), one can use the scaling and squaring strategy for computing the exponential [31]. If one chooses such that the field is close enough to zero, the first-order approximation can be used to estimate (based on the Padé approximant near the origin). Then, the solution of the flow equation is computed by performing recursive compositions of the field by itself, given that such compositions are computationally affordable. Notice that taking is equivalent to the simple first-order approximation. The scaling and squaring steps for field exponentiation [22] are depicted hereafter.(i)Scaling. Divide by a factor such that is small enough, for example, when = voxels.(ii)Exponentiation. Compute first-order explicit integration of the flow: .(iii)Squaring. Perform recursive squarings (using field composition) of the flow at time in order to obtain the flow at time 1, which is the field exponential. In other words, starting with , do times the computation , in order to get .
We see that using this method, only compositions (and therefore interpolations) are needed for estimating the exponential. Compared to standard estimation of the flow over a regular discretization of the time interval in steps, the scaling and squaring method limits the numerical errors due to composition of vector fields, but it does not decrease the amplification of the error due to the field estimation at time .
3. Generic Registration Pipeline
Nonrigid registration methods can be divided into parametric and nonparametric methods. Parametric (or model-based) methods aim at calculating the parameters of a deformation model in a high-dimensional space in order to optimize a global objective function that takes into account image similarity and transformation regularity [10]. In this case, the a priori information is included in the modelization and regularity criteria of the nonrigid transformation. For example, the harmonic energy of transformation can be explicitely included in the objective function [32].
On the other hand, nonparametric methods make it possible to decouple similarity optimization from regularization by directly acting on the displacement field. The a priori information has then to be included in the optimization process by using proper regularization techniques. Decoupled optimization makes the registration computationally efficient [8], mainly because the computation of each displacement vector is independent from others, but it prevents us from easily including more complex regularization constraints in the process, for example, such as in volume preserving registrations [33, 34].
3.1. Multiscale Nonparametric Registration
Most nonparametric registrations are based on an iterative process which is composed of 3 steps: (i) field computation, (ii) field accumulation, and (iii) field regularization. The idea is to progressively build a proper displacement field by iteratively improving the matching between the fixed image and the moving image warped by this displacement field, according to a certain metric. Note that, depending on the nature of the displacement one tries to model, the regularization is applied either on the increment field or on the accumulated field. Regularizing the field increment corresponds to a viscous fluid modeling, while regularizing the global transformation corresponds to an elastic solid modeling [14]. Only the second is considered in this study.
In this paper, our general nonparametric registration framework (e.g., valid for Demons and Morphons) adopts a multiscale approach; that is, the displacement field estimation is stabilized by decomposing the fixed and the moving images in several scales, for example, using a simple smoothing and downsampling procedure [35].
The three steps mentioned above are then applied a certain number of times (until the algorithm reaches a certain stopping criterion) to each scale separately from coarse to fine scales (Figure 3). The general explanation of these three basic blocks is given hereafter. The way they are iteratively applied at each scale is described in Section 3.2.
3.1.1. Field Computation
At each iteration of the registration process, an update displacement field () is first computed as a function () of the fixed image () and the moving image () warped by the displacement field resulting from previous iterations (): where and denote the deformation operations linked to and , respectively.
Depending on the nature of the images to be registered, this local displacement estimation can be based on different local image metrics, such as SSD [17], mutual information computed on blocks of voxels [8, 36], and local phase [26].
3.1.2. Field Accumulation
After the field computation, the total displacement must be increased by the update field
This accumulation operation is sometimes implemented as a simple addition of accumulated and update fields (as in [18, 37, 38]). However, as explained in Section 2.2, this accumulation is perhaps computationally efficient but is not consistent with the composition of the corresponding spatial transformations. The solution is, therefore, to replace it by the compositive accumulation introduced earlier. The accumulation of the displacement fields with is then compatible with the way is estimated. Indeed, since is computed from , the accumulation of and must modify by , a process intrinsically integrated by the operation . Moreover, the associativity of validates the compositive accumulation of displacement fields over several iterations, as illustrated in Figure 3.
3.1.3. Field Regularization
Eventually, the field is regularized in order to get a smoother transformation and reduce the impact of image noise on the registration output: This operation is achieved by applying a low-pass filter on each component of the displacement field. We assume it to be a Gaussian smoothing with a size of a few voxels, which tends to reduce the harmonic energy of the transformation [32].
It is always possible to produce invertible fields by performing a very strong Gaussian smoothing. This, however, may reduce significantly the accuracy of the estimated displacement by limiting the solution to excessively smooth displacement fields. On the other hand, by preventing the displacement field from being noninvertible, the diffeomorphic accumulation acts in some way as a regularization, allowing the estimation of invertible fields while performing only moderate smoothing.
3.2. Registration Algorithm
Let us explain now the whole multiscale nonparametric registration algorithm relying on the three specific procedures defined in Section 3.1.
The algorithm takes as inputs the fixed and the moving images and , some parameters described below, and outputs the estimated transformation such that . The whole procedure described in Table 1 and depicted in Figure 3 involves computations on different scales , from coarse () to fine (). Each scale is associated to a subsampled grid of voxels , where is the subsampling factor (e.g., in this study) between scale and scale . The functions and , defined on the initial grid , are downsampled (after antialiasing smoothing) at any scale by the operation Downj. An upsampling operator Up, implemented as a simple linear interpolation, is used to transfer any displacement field defined on a grid to the finer grid using as upsampling factor. For each scale , the accumulated displacement field is iteratively updated until one reaches a particular stopping criterion (e.g., based on the convergence of or on the SSD, as precised in Section 5).
4. Diffeomorphic Morphons
Our paper adapts the global registration method explained in the previous section to Morphons [26, 39] by taking care of the invertibility of the accumulated displacement field, that is, by introducing diffeomorphic field accumulations.
As already mentioned above, the particularity of Morphons, compared to other nonparametric methods, is that the field computation (function in (5)) is based on the local phase rather than intensity difference. In other words, knowing the phase difference between periodic signals of the same frequency allows the estimation of the spatial shift between them. Therefore, under the assumption that images can locally be considered as a sum of periodic signals, the computation of the local phase difference is equivalent to the estimation of the local displacement between images. This procedure is stabilized by the multiscale approach described in Section 3.2. Besides, Morphons combines the estimation of displacement vectors with a measure of the confidence we have in these estimations, resulting in a certainty map. Therefore, for Morphons, given two images and , the displacement field estimation is actually split into two quantities that is, respectively, an update of the displacement field along with an update of the certainty map. A similar split is also performed on subsequent operations and .
Here are the details about the three steps of the pipeline of Section 3 for this specific registration, including our contribution to the field accumulation step.
4.1. Displacement Field Calculation
In Morphons, a displacement field is estimated thanks to the dephasing between the local phases of the fixed and the moving images. This local phase can be probed at a certain frequency and in a particular direction using quadrature filters [40]. More precisely, Morphons method uses a quadrature filter of direction (also called loglets [40]) defined in frequency by the polar separable function where is the frequency vector, if and 0 else, , is the unit vector supporting , and is a radial function centered on and defined as for .
Since their support corresponds to the half volume and since (with the angle separating and ), loglets can be seen as the analytic counterparts of the steerable filters introduced by Freeman and Adelson [41]. As a matter of fact, only a limited number of orientations are necessary to cover the whole frequency plane. Typically, in 2D, these directions are taken as with for , and in 3D, is taken as the 6 normal vectors to the faces of a hemi-icosahedron [42, 43]. Notice also that each filter in the spatial domain is centered around the origin with a typical width given by .
Morphons take advantage of the following behavior. Given an image , defining the filtering with the common convolution operation and the shorthand , we can write since . Therefore, by processing the warped image similarly, the local phase difference can be computed as with the complex conjugation and the local dephasing between and in direction .
An important observation is that the nonnegative value represents also the correlation between and the translated filter ; that is, the filter translated on . If the image was perfectly represented by the latter, that is, if we had locally for any and some constant , a displacement of by a displacement field approximately constant over the support of would induce a dephasing since the frequency vector of is . An important implicit assumption is nevertheless that since the dephasing is known up to modulo . Moreover, only and not can be determined, as another manifestation of the blank wall problem [44].
In practice, for most of , is not perfectly represented by one filter but by a linear combination of them where the amplitude measures the adequacy of the fit between and . Consequently, the local update displacement linking and in each is estimated by solving the weighted least square optimization where the are the certainty map of the filter . As explained above, reflects for each voxel how reliable the field estimation is; that is, how contrasted the bandpass-filtered images are.
Numerically, the optimization in (13) is a standard weighted least square minimization; that is, it corresponds the minimization of the energy , using the diagonal matrix , the matrix , and the vector . An easy computation shows that the solution of (13) is then given by the Moore-Penrose pseudoinverse , that is, with arbitrary set to 0 when is not invertible.
Jointly to the estimation (13), a global certainty map associated to the quality of the estimation of is defined as [43] that is, the sum of all certainty measures for each quadrature filter. This update of the certainty map must then be combined with an accumulated certainty computed from previous iterations (see Section 4.2).
In the multiscale approach described in Section 3.2, using the same quadrature filters at decreasing scales is equivalent to estimating the phase of the bandpass-filtered image around increasing cutoff frequencies, that is, with each time . This sustains the coarse-to-fine displacement estimation, that is, the computation of and on different scale bands and of and .
Convolutions with quadrature filters can be implemented efficiently in the Fourier domain thanks to the FFT and the convolution theorem. However, since the spatial extent of those filters is small, it is also possible to use efficient spatial convolutions with truncated kernels, as done in this study. As the local phase is invariant to local intensity scaling, the Morphons procedure is suitable for registering images with various contrast enhancements. Besides, some studies indicate that the phase extraction allows a fast convergence and a subvoxel precision in displacement estimation (e.g., see [39]).
4.2. Field Accumulation
In the original Morphons method, the accumulated field is computed as a weighted sum of the update field and the previous accumulated field, as used in damped optimization schemes. The weights are given by the certainty on the update field (, as computed from ) and the accumulated certainty map (). As the certainty map must also be accumulated in order to reflect the confidence in all previous displacement computations, the accumulation step must be divided into two operations (field accumulation) and (certainty accumulation): where in the last formula, similar to the field accumulation, the certainty map is updated by its own certainty [43].
However, as it was explained before, the addition of displacement fields is not really appropriate for accumulating spatial transformations, in contrast to composition. The compositive accumulation may also be damped using the certainty as a weighting factor
The (SSD-based) Demons registration is a nonparametric algorithm which performs the optimization of the SSD between images. In [27], a diffeomorphic field accumulation is proposed as improvement of the Demons method. The idea is to use an adaptation of the optimization method to Lie groups [45] in order to limit the possible solutions to diffeomorphic transformations. In practice, this is done by replacing the accumulation step of the Demons by an accumulation using the diffeomorphic flow introduced in Section 2. This accumulation reads then where the field exponential can be efficiently estimated using a small number of recursive compositions of the field by itself. Consequently, the displacement field is linked to the deformation operation .
In the case of the Morphons, the accumulation step can be achieved in the same way. This will produce smoother fields than the traditional addition or composition. However, the accumulation step in the Morphons method involves a damping based on the certainty. Therefore, we propose the following accumulation step for diffeomorphic Morphons: Since for any vector field , the accumulation fades away when . The accumulation of the certainty map remains as explained previously in (17).
Notice that, because the field is discretized on a grid of voxels, interpolation is needed for computing the composition of two diffeomorphisms. Therefore, errors due to successive interpolations could potentially lead to noninvertible transformations. However, such problems were not observed in practical experiments using reasonable smoothing of the field.
4.3. Field Regularization
During the displacement estimation step, the relevance of local phase computation is estimated and used as weight for the accumulation. This certainty map may also be used for a smart regularization of the displacement field. Regularization is performed using a normalized convolution [46] of the field by a Gaussian kernel, taking into account the certainty map in order to put greater importance to high certainty locations. The certainty is also regularized in the same way as displacement field components in order to preserve the correspondence between the displacement vectors and their corresponding certainty.
Mathematically, given a positive function and a filter (typically a Gaussian kernel of variance ), the normalized convolution of a (scalar) function by as involved by the normalization is This operation does not increase the maximum amplitude of the filtered function. Indeed, for a nonnegative kernel , we show easily that , with . The accumulated displacement field and subsequently the certainty map are therefore regularized thanks to this operation using for normalization the certainty map , that is, Notice that, for computing , the normalized convolution is performed separately on all components of the vector field.
This operation tends to propagate the displacement field from high certainty areas to areas which show less significant filter responses. Besides, by setting to zero the certainty outside the volume boundaries, normalized convolution cancels the influence of the padding strategy. This step produces a smooth version of the accumulated field that may reduce the accuracy of image matching resulting from the displacement estimation step, as it limits the possible solutions to smooth displacement fields.
However, if the iterative algorithm is to converge, the solution will be regular and invertible (except for large numerical errors), thanks to accumulation and regularization constraints, but it will also be (at least locally) optimal in terms of local phase difference. Indeed, as the phase is monotonic and smooth, a mismatch between local structures will automatically lead to nonzero field update with a high certainty value, which will tend to improve the displacement estimate and fit the structures together.
The Jacobian of the displacement field may be used as a criterion for validating the physical behavior of the deformation. Indeed, the Jacobian gives for each voxel the change in volume this voxel encounters during deformation. Jacobian indicates expansion when it is greater than 1, and compression when it is smaller than 1. A negative Jacobian means that the voxel is “inverted” (getting a negative volume), which is incompatible with the mass-preservation principle.
In the following, the diffeomorphic version of Demons and Morphons are denoted respectively D-Demons and D-Morphons.
5. Experiments and Results
The methods were first compared for several simple 2D virtual situations in order to demonstrate the interest in chosing the accurate registration method with respect to the images to be registered.
For the clinical validation, Morphons and D-Morphons registrations were first validated on a 10-phase point-validated pixel-based breathing thorax model (POPI-model) from the Léon Bérard Cancer Center, Lyon, France [4], in order to compare the D-Morphons to Morphons, Demons, and D-Demons in the case of intensity conservation between images. Then, it was applied to lung images with different contrast enhancements, in order to illustrate the benefit of a phase-based approach compared to traditional SSD-based registration methods in the case where intensities are not conserved between the images to be registered.
All simulations were performed using Linux, on a single processor Intel Core 2 (2.4 GHz). Our MATLAB implementation used for the prototyping of the methods was also used for simulation. Notice that no efforts were made for achieving good performances in terms of computational cost and memory requirements in the implementations used in this study. The local phase estimation was performed using convolutions with quadrature filters. Less than 1 GB of RAM was required for registering two volumes of voxels using all registrations. The time required for registering such images, using the parameters presented hereafter, was around 6 minutes for Demons, 42 for Morphons, 7 for D-Demons and 43 for D-Morphons. However, preliminary results based on a C++ implementation of the Morphons, which uses operations in the Fourier domain instead of convolutions (as done in our matlab implementation) and using 4 threads on a quad-core CPU, allowed a division of the computation time by 50, leading to Morphons registrations taking about one minute for such a typical image size.
5.1. Illustrative Virtual Experiments
Two 2D virtual experiments were performed. The first experiment, illustrated in Figure 4, is based on a virtual disk image after blurring. Two images of the same disk were created, the only difference being the scale of intensities (multiplication by 0.75). This experiment shows the interest in using a phase-based method (conventional Morphons in this example) while registering identical shapes with different contrasts, compared to an intensity-based method (conventional Demons).
The second virtual experiment is based on two images of a disk (see Figure 5). In the fixed image, a disk of radius was created, and a hole (disk of radius ) was added in its center. In the moving image, a disk of radius was created with the same intensity scaling as in the fixed image. This example illustrates the case where a structure is missing in one image compared to the other, as it may occur in practice (e.g., the problem of bowel gas in CT images of the abdomen). This experiment illustrates how the diffeomorphic version of the Morphons algorithm can prevent from producing negative volumes after registration, without increasing the smoothing by using a larger Gaussian regularization kernel.
5.2. Accuracy Assessment on a Breathing Thorax Model
The POPI model [4] is composed of 10 volumes reconstructed from a 4D respiration-correlated CT scan (RCCT) of the thorax, each volume corresponding to a particular phase of an average breathing cycle. 41 landmarks were identified by medical experts in each of the 10 images for registration validation.
Conventional Morphons, D-Demons, and D-Morphons were applied between a reference phase and the 9 others. For all methods, the number of scales was set to , with final resolution of . The variance of the Gaussian kernel used for regularization was empirically set to twice the voxel size ( voxels). For this experiment, a minimum of 10 and a maximum of 20 iterations was used at each scale. In between, the iterative process was stopped if the changes, measured in terms of SSD, were inferior to 0.01%. Such a convergence criterion was usually reached before the 20th iteration, supporting the fact that both Demons and Morphons behave like optimization methods.
The results were then compared with each other and with the results from a conventional Demons algorithm as used in [4]. The comparisons were achieved in terms of error in landmark position, SSD between images, harmonic energy, and minimum Jacobian.(i)The landmark position error evaluates the ability of the registration in finding the physical motion of organs. (ii)The SSD between fixed and deformed images is a measure of the image matching according to the assumption of intensity conservation. It is computed as .(iii)The harmonic energy [29, 32] of the displacement field indicates how regular the field is and is computed as .(iv)The Jacobian of the field indicates the volume change of each voxel. Recall that negative values of the Jacobian correspond to inverted volumes, which is not acceptable in a physical point of view. The Jacobian is computed as , with , where is Kronecker's delta ( if , 0 else) and is the component of the displacement field. In practice, the partial derivatives can be computed using centered finite difference approximations.
The comparisons of landmark position errors (expressed in ) resulting from the different registrations can be seen in Table 2 with, from left to right, the error in landmark position (norm of the difference) before registration, using Demons (values from the POPI website), Morphons, D-Demons, and D-Morphons. Position errors are noted as follows: mean/std (max). On average, for Morphons, D-Demons, and D-Morphons, the error in landmark position was equal or inferior to 1 mm, which is half the size of the voxels at the finest scale of the registration process.
Results showed that all registrations greatly improved the matching of intensities. The SSD between fixed and deformed image was similar for Morphons, D-Demons, and D-Morphons (see Figure 6). The harmonic energy of the fields resulting from these registrations was also comparable (see Figure 6).
(a)
(b)
(c)
The matching and the harmonic energy obtained by Demons (as presented by the authors of [4] on the POPI website) was slightly less good than for the 3 other methods. However, this is most likely due to the parameters used for registration (e.g., the number of scales, the variance for smoothing, etc.). In particular, for very similar images (first 2 phases of the RCCT), the algorithm was not able to find a smooth displacement field that reduced the SSD.
The minimum Jacobian of the displacement fields resulting from conventional methods gets down to for both Demons and Morphons (see Figure 6), as, respectively, 67 and 460 voxels were inverted for the corresponding phase when applying the field on the moving image (which is composed of almost 6 mega voxels). However, when using diffeomorphic accumulation, the minimum Jacobian was raised to 0.2 for the Demons and 0.1 for the Morphons, showing that the diffeomorphic accumulation step prevented the field from inverting voxels.
5.3. Application to Images of the Thorax with and without Lodine Contrast Agent
The breathing-correlated motion of tumor is a typical feature of lung cancer that has to be dealt with in radiotherapy planning. RCCT images provide information about the tumor motion throughout the breathing cycle. From the different respiratory phases, an adequate margin around the tumor (the ITV, i.e., the Internal Target Volume) can be estimated, integrating thus all tumor positions through the respiratory cycle [20].
However, the lack of contrast enhancement, as well as the high noise level and the presence of artifacts that characterize 4D RCCT, may significantly impair the accurate delineation of the target volumes on these images. More particularly, the iodine contrast agent is of prime importance to help at differentiating tumor extents from vascular structures in the centrally located lung tumors. In this context, the acquisition of a conventional contrast-enhanced CT (CE-CT) acquired during free breathing should be considered for the delineation task, while the 4D RCCT is used to estimate the motion range of the tumor during breathing. To automatize this process, the delineated tumor volume at the CE-CT can be deformed on the various respiratory phase images from the 4D RCCT using nonrigid registration to finally get the ITV, as illustrated on Figure 7.
The purpose of this experiment is to compare Demons and Morphons algorithms (conventional and diffeomorphic versions) for the registration between images with and without contrast enhancement, while keeping the same setting as for the POPI experiment.
A CE-CT scan of 3 lung cancer patients was acquired as well as a 4D RCCT scan at another time point. The first CT scan was taken in free breathing using an iodine contrast agent. The 4D RCCT scan was acquired without any contrast agent and was reconstructed into 10 phases. Histogram equalization was not able to correct for localized contrast differences between the CE-CT and RCCT phase images. For all 3 patients, Demons, Morphons, D-Demons, and D-Morphons were applied between each of the 10 RCCT images and the CE-CT, with the same registration parameters as for the POPI simulation.
The displacement fields resulting from these registrations were compared in terms of harmonic energy and minimum Jacobian (see Figure 8). The resulting images were compared in terms of SSD and mutual information.
(a)
(b)
The harmonic energy of displacement fields resulting from Demons and D-Demons was quite higher than that with the Morphons and D-Morphons, and the minimum Jacobian of the displacement fields was positive only for registrations using the diffeomorphic accumulation. In the worst case, 7455 and 1114 voxels were inverted using respectively Demons and Morphons without diffeomorphic accumulation (on an image of 5 mega voxels). An example of area leading to bad transformations (with negative Jacobians) using conventional methods is depicted in Figure 9. D-Morphons lead to the smoothest transformation, with minimum Jacobian values around 0.2. These quite low values, however, were very sporadic within the image volume.
We noticed that, unlike the results obtained with the POPI simulation, the SSD resulting from the Morphons and D-Morphons was a bit higher than the SSD resulting from Demons and D-Demons. However, as illustrated in the example of Figure 4, the SSD does not reflect the matching in variable contrast areas. On the other hand, no significant differences in terms of mutual information were observed between images resulting from the different registrations. This is likely due to the very low contrasts in the noncontrasted images within the regions corresponding to contrast-enhanced tissues in the other image whereas the main differences in terms of displacement field were located in these regions, as illustrated in Figure 10.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
In order to illustrate the effect of the registration on contrast-enhanced tissues, one phase of the RCCT scan of one of the 3 patients was chosen as example. For this patient, the tumor was located close to contrasted tissues. The tumor and the blood vessels were delineated by a physician, on the contrast-enhanced scan and on one phase of the RCCT scan. The delineations on the phase image were deformed according to the fields resulting from the different registrations. The results are illustrated in Figure 10.
The change in volume due to warping was computed, as well as the harmonic energy inside the delineated stuctures and the difference between the center of mass of the tumor with and without registration.
The change in volume was very small when using a phase-based field computation for both the vessels (around 0%) and for the tumor delineations (around 1%), while it rose up to 23% for the vessels and to 6% for the tumor while using the Demons. In the same way, the harmonic energy and the error on the center of mass of the tumor were much smaller for the phase-based registration methods. These results are summarized in Table 3. One can notice that the diffeomorphic accumulation of the field in the Morphons did not change the results in terms of harmonic energy and volume changes compared to conventional Morphons. This is due to the fact that the displacement of the considered organs is small and smooth.
6. Discussion
The first medical experiment showed that D-Morphons and D-Demons lead to similar matching of both image intensities and anatomical landmarks. This shows that for monomodal registration of lung CT scans, the phase difference has an efficiency comparable to the efficiency of the SSD metric. Furthermore, the D-Morphons produced displacement fields as smooth as those obtained with D-Demons. In opposition to conventional Demons and Morphons, both diffeomporphic methods produced invertible displacement fields which are physically meaningful.
The second medical experiment illustrates the limitations in registering images with various levels of contrast enhancement with the Demons method. Indeed, the intensity matching resulting from Demons was better than that from Morphons, but the field was obviously wrong, as the Demons results in a global shrinking of the contrasted tissues (arteries) that does not reflect a proper anatomical behavior, but that is due to the fact that the Demons registration is based on the minimization of the SSD, which produces an improper displacement estimation when the intensities of identical tissues are different in the fixed and moving images. This mismatch between registered anatomical structures is clearly visible on Figure 10. As illustrated in the example of Figure 4, the field produced by Demons tries to match structures of same intensity, which do not correspond to identical anatomical structures because of the difference in contrast agent concentration. Therefore, the field resulting from Demons (see the field on the left part of Figure 10) is far less smooth than it should be and can lead to wrong deformation estimations as it illustrated in the example (see Table 3). In this case, the difference in intensity between the images with and without contrast enhancement lead to important volume changes for vessels and tumor by using Demons or D-Demons, while almost no changes in volume were observed for these tissues when using a phase-based approach. Besides, the harmonic energy inside these tissues shows that the field is much more smooth using the phase-based registration. It is important to notice that these effects are mostly limited by the regularization of the displacement field during the Demons and D-Demons registrations, and that they will still be worse if less regularization is used (smaller variance of the Gaussian kernel used for smoothing the displacement field). This is not the case for the fields produced by the Morphons and diffeomorphic Morphons, which are much smoother and preserve the anatomical topology even with contrast variations between images (see Figure 10). Notice that the reduction of the smallest segmentation that can be observed in the Morphons results is mostly due to interslices motion, as confirmed by the Jacobian close to 1 in this area that shows that there is no important volume changes within this segmented region. Finally, one can see that the invertibility of the displacement field is observed with both diffeomorphic registrations.
These results can be summarized by classifying the different registration strategies according to the smoothness (harmonic energy) and the invertibility (minimum Jacobian) of the resulting displacement fields (see Table 4) for the variable contrast experiment.
One can notice that the D-Morphons algorithm combines both advantages: the field is invertible and smooth, which suggests that it is likely a better estimation of the real transformation which is known to be smooth in this area.
7. Conclusion
The D-Morphons is a multiresolution registration algorithm which computes a diffeomorphic displacement field based on the minimization of the local intensity phase. The method managed to estimate the deformations in a breathing thorax, with an accuracy comparable to the accuracy of the D-Demons, and leads to the same requisite property of invertibility of the field. Moreover, the D-Morphons managed to accurately estimate the deformations between images with variable contrast, while the conventional SSD-based methods led to misalignment of anatomical structures affected by the contrast variation.
Acknowledgments
The authors thank the Medical Informatics Laboratory of Linköping University for sharing their implementation of the Morphons. G. Janssens and J. O. de Xivry are funded by an FRIA grant. L. Jacques is a Postdoctoral Researcher funded by the Belgian FRS-FNRS.