This paper describes a method for nonrigid registration of monomodal MRI based on physical laws. The proposed method assumes that the properties of image deformations are like those of viscoelastic matter, which exhibits the properties of both an elastic solid and a viscous fluid. Therefore, the deformation fields of the deformed image are constrained by both sets of properties. After global registration, the local shape variations are assumed to have the properties of the Maxwell model of linear viscoelasticity, and the deformation fields are constrained by the corresponding partial differential equations. To speed up the registration, an adaptive force is introduced according to the maximum displacement of each iteration. Both synthetic datasets and real datasets are used to evaluate the proposed method. We compare the results of the linear viscoelastic model with those of the fluid model on the basis of both the standard and adaptive forces. The results demonstrate that the adaptive force increases in both models and that the linear viscoelastic model improves the registration accuracy.

1. Introduction

A current major topic in computational neuroanatomy is the development of nonrigid image registration algorithms. Nonrigid image registration has many potential applications. It can be applied to building atlases, segmentation, quantifying local morphological differences, comparing the variance of different population, and detecting pathological changes [14]. Various methods have been developed to deal with nonrigid image registration. The methods are usually classified into two categories: feature-based and intensity-based [5]. The former first needs to build a geometric model and identify a number of anatomic characters in the model. These characters include point landmarks, curves, and surfaces [68]. The anatomical characters are then parameterized. The aim of the registration is to find the optimal combinations of the model parameters. These methods are critically dependent on feature extraction quality. However, the anatomical structures are complex, making it difficult to extract them accurately. Usually, human interaction is required during registration, thereby making the process inconvenient and time consuming. The intensity-based method is used to match regional intensity patterns based on mathematical or statistical criteria [913]. This method reduces the need for direct feature extraction or segmentation, can be automatic, and can obtain satisfactory results and is thus widely used. Fluid registration uses physics models and assumes that the deformation obeys fluid mechanics laws. These methods allow flexible deformation with large freedom and are used in many applications.

In the early 1980s, the elastic model was proposed as a means to match images [14]. Brain images are modeled as an elastic solid and the deformations are calculated from elastic mechanic equations [15, 16]. However, this model is only suitable for small deformations. To address this problem, the properties of brain images were assumed to be like those of viscous fluid, and the viscous fluid model was proposed [17], where the deformations are driven by forces that are equal to the gradient of the sum of squared intensity difference (SSD) metrics. The orientation and magnitude of deformation fields are computed using the fluid-dynamical Navier-Stokes equation. This method allows large deformations and serious localized distortions, but with increased likelihood of misregistration [18]. The Navier-Stokes equation is solved by means of the relaxation method and requires much time. To reduce the computation cost, Bro-Nielson proposed a convolution filter method to solve the equation quickly [19]. However, Wollny et al. [20] obtained unsatisfactory results when using a small filter width in the convolution filter methods. However, if the filter width is large, the computational costs are not more advantageous than when using iterative methods. Thus the relaxation method is currently the best method. The brain images are also modeled as diffusions [2123] and have been shown to be similar to the method proposed by Bro-Nielsen et al. [19].

In this study, we try to use the properties of both elastic solids and viscous fluids to register images. Linear viscoelastic matter has these properties. The deformation properties of brain images are assumed to be similar to those of viscoelastic matter and obey viscoelastic laws [24, 25]. As the Maxwell model has the abilities to describe linear viscoelastic deformation [26], this study hence utilizes the Maxwell model to represent and capture large deformation of the brain images. When a force acts on the Maxwell model, the motion of the fluid component relaxes over time, allowing large displacements. The deformation fields are constrained by both elastic and fluid components. To speed up the algorithm, an adaptive force is introduced. Given our aim of monomodal anatomic image registration, the SSD is used as a similarity metric in the registration. Both synthetic and real images are used to demonstrate the performance of the proposed method. The performances of both models (Maxwell and fluid) with both forces (standard and adaptive) are compared with each other. The fluid model with adaptive forces (FMAF) has the fastest registration speed, the Maxwell model with adaptive force (MMAF) is the second, the Maxwell model with standard force (MMSF) is the third, and the fluid model with standard force (FMSF) is the slowest. The ranking of registration accuracy from high to low is as follows: MMAF, MMSF, FMAF, and FMSF.

2. Materials and Methods

2.1. Maxwell Model

The Maxwell model [26] is made of a spring and a dashpot in series (Figure 1), which is perfectly elastic and viscous. Since the deformation process is assumed to be quasistatic, inertia can be neglected and the force or stress is the same in both parts. The total deformation is the sum deformations of both parts. If the displacement of spring or dashpot is or , the total displacement is If is Young’s modulus of the spring and is the viscosity of the dashpot, and are the stresses of the spring and the dashpot. The stresses are Given that the force on the spring and the depot is equal at any given time, the two parts can be processed independently.

2.2. Reference Frame

Two kinds of reference frames are used to describe deformations in a floating image that is deformed to a target image. One is the Lagrangian reference frame, which describes the deformations by observing changes in the positions and velocities of definite particles. The other is the Eulerian reference frame, which describes the deformations by observing velocity changes at fixed points. The Eulerian reference frame is suitable for large deformations because it does not trace the motion of the particles [27]. Therefore, the Eulerian reference frame is used to track the deformations in our method. Voxel grids are used as the fixed points. A particle at grid position in floating image at time is originated at the position , where is the displacement field. The corresponding velocity field is expressed as where . It comes from the derivative of the displacement field about time. The second term in (3) represents the nonlinearities of the displacement field.

2.3. The Viscoelastic Fluid Algorithm

We extend the Maxwell model to three dimensions. The spring becomes an elastic solid, and the dashpot becomes a viscous fluid. Therefore, the total deformation is similar to that in (1), where and are displacements of the elastic solid part and the viscous fluid part, respectively. The force of the two parts is equal and is expressed as where and are the forces acting on the elastic solid part and the viscous fluid part, respectively.

We used the continuum mechanics method to compute the displacements. The elastic solid displacements are described by the following partial differential equations: where and are Lame’s elastic coefficients and . The velocity of the viscous fluid part is determined using the following equation: where and and are the viscosity constants.

The velocity field of viscous fluid in an Eulerian reference frame can be determined by the following equation: The displacement fields are updated iteratively over time step and are determined as follows: Time step is chosen according to the perturbation of the deformation field; we have The boundary conditions and , and the total displacements on the boundary are set to zero. The elastic equation (5) and fluid equation (6) are solved simultaneously to obtain the total deformation.

2.4. The Adaptive Force

The motivation of the adaptive force is to speed up the registration. In the proposed method, the key parts of the partial differential equations (PDEs) (5) and (6) are the forces that drive the floating image to deform to the target image . The gradient of the SSD metrics is used as these forces. The standard force is defined as where is a constant.

The is the gradient of the floating image at . is the difference in intensity between the deformed floating image and the target image and weighs the . The force is minimized at the location where the floating image and the target image are aligned.

As the registration progresses, the forces become smaller and the corresponding velocities also become smaller which lead to very small deformation in the iteration. Therefore, more iterations are needed to reach the final deformations. To speed up the registration, the forces should increase in the next iteration. Hence, an adaptive force is introduced to solve the problem in the proposed method. The maximum of the displacements should not stay below a specific threshold. When the maximum of the displacements is below the threshold at the current iteration, the forces are adjusted automatically to increase the maximum of the displacements in the next iteration. In our method, an empirical formula is used to define the adaptive force. The adaptive force in the Eulerian reference frame is expressed as where is the next iteration and is the function with respect to the maximum displacements of the current iteration. It is described as where is the maximum of the displacements of the current iteration.

If the maximum displacement is below the threshold , should be larger than one, thereby making larger than . The parameter of the next iteration increases automatically. Therefore, the corresponding forces also increase to prevent the displacement from becoming too small.

2.5. Implementation

When the floating image is deformed by the corresponding deformation field in the registration, the topology of the floating image should be preserved. Keeping all Jacobian of the deformation fields positive can preserve the topology. In the implementation, when the minimum of the Jacobian is below 0.5, the transformation is applied to the floating image to produce a new image and the displacement field is set to zero. The new image is then used as the floating image in the subsequent registration. The process continues as long as the SSD decreases. The pseudocode of the algorithm is as follows.(1)Let and .(2)Calculate the force using (10), (11), and (12).(3)If the SSD stops decreasing or the maximum number of iterations is reached, then stop.(4)Solve PDEs (5) and (6) for displacements and instantaneous velocity , respectively.(5)Choose time step according to (8) and calculate .(6)Calculate the total displacement .(7)If the Jacobian of the transformation is less than 0.5, a new floating image is constructed, and then go to Step 1. Otherwise, update the displacement field according to (4), then set , and go to Step 2.

As it has been proved that relaxation is currently the best method [20], we solve PDEs (6) and (7) by means of successive overrelaxation [28].

2.6. Evaluation

The performance of the proposed method is evaluated on the basis of three analyses. The first analysis uses the golden deformation field . We then compare the recovered deformation field by the root mean square (RMS) error over all voxels: where is the number of total voxels in the image.

The second analysis uses the mean of the SSD. They are defined as where is the number of total voxels in the image, and are floating and target images, respectively. While is the deformed floating image.

The last analysis uses tissue overlaps, which are defined as where is the volume of the tissues.

If the floating image completely matches the target image, the value would be one, and the RMS or SSD would be minimized. If there is no overlap between the two images, the value would be zero, and the RMS or SSD would be maximized.

3. Experiments and Results

Four experiments are conducted to demonstrate the proposed method. The first two experiments are about 2D data and the rest 3D volumes. The method is implemented in C and complies with VC++ [29]. The whole image is modeled using a single set of material parameters for simplification purposes. The parameters , , are all set to one and is set to zero. The parameters and are both set to 1, and is set to 0.8 voxels. The maximum iteration is set to 200. These parameters are used in all the experiments.

3.1. 2D Synthetic Datasets

The experiment shows that the proposed method can deal with large deformation well. The image sizes are pixels, as shown in Figure 2. Figure 2(a) is the floating image, with a rectangular image, and Figure 2(b) is the target image, a C-shape image. The results of FMSF, MMSF, FMAF, and MMAF are all successful to deform the rectangular image to C-shape image. Figure 2(c) shows the results of FMAF. The computing costs are listed in the second row of Table 1. The computing times of FMSF, FMAF, MMSF, and MMAF are 36, 12, 26, and 18 seconds, respectively. The ranking of speed from the fastest to the slowest is as follows: FMAF, MMAF, MMSF, and FMSF.

3.2. 2D Brain MRI Datasets

The second experiment shows the effectiveness of the proposed method when it is applied to brain MRI. The floating image size is 256 256, as shown in Figure 3(a). This image is registered to a selected image by the finite element method [30] and obtains the deformation fields , which is used as the golden standard. The known deformation fields are applied to the floating image to obtain a target image, as shown in Figure 3(b). The known field is shown in Figure 3(c) using the following equation: where , are the known deformation fields in the and directions, respectively.

The fluid model and the Maxwell model with the standard forces and the adaptive forces are applied to the images. The computing time is listed in the third row of Table 1. The FMAD only costs 23 seconds, which is the fastest. By contrast, the MMAF and the MMSF cost 32 and 41 seconds, respectively, and they are the second and the third in terms of speed. The FMSF is the slowest and costs 60 seconds.

All of these methods can successfully deform the floating image to the target image. However, the matching accuracy is different. Table 2 lists the RMS and SSD acquired by the various methods. The RMS of FMSF, FMAF, MMSF, and MMFA are 0.3783, 0.3066, 0.2742, and 0.2412. According to RMS, MMAF has obtained the best result, followed by the MMSF, then FMAF, and FMSF last.

Figure 4 shows the difference among the known deformation fields obtained using various methods using the following equation: where is the difference of the known deformation fields. , are the deformation fields of and directions, respectively. We find that regardless of which forces act, the differences of the known deformation fields with that of the fluid model (shown in Figures 4(a) and 4(b)) are much larger than that with the Maxwell model (shown in Figures 4(c) and 4(d)). The Maxwell model has obtained better results. Among them, the Maxwell model with the adaptive forces has obtained the best result, whereas the fluid model with the standard forces has the worst results. The SSD of FMSF, FMAF, MMSF, and MMFA are 0.0381, 0.0293, 0.0367, and 0.0293. Based on the SSD, the results of the MMAF and FMAF have the same rank, whereas the results of the MMSF and the FMSF are ranked as third and fourth. It indicates that the adaptive force is superior to the standard force and the Maxwell model is more robust.

3.3. The Internet Brain Segmentation Repository (IBSR) Database

High-resolution 3D MR images are used to evaluate the proposed method. The MRI data are downloaded from the IBSR [31] and include 20 normal MR brain datasets and the skulls are all stripped. The column and row are 256, and the slice is from 58 to 64. The voxel size is . We have randomly selected 30 couples from the data to test the methods. The mean time is listed in the fourth row of Table 1. The computing times of FMSF, FMAF, MMSF, and MMAF are 2891, 1143, 2025, and 1793. The computation time of the fluid model with the adaptive forces is the fastest, the Maxwell model with the adaptive forces is the second, the model with the standard forces is the third, and the fluid model with the standard force is the slowest.

The mean of SSD is listed in the second row of Table 3. The SSD of the FMSF, FMAF, MMSF, and MMAF are 0.0553, 0.0490, 0.0504, and 0.0484, respectively. According to the SSD value, the results obtained by the MMAF and the FMAF are ranked as the first and the second, respectively, and those acquired by the MMSF are the third. The result obtained using the FMSF is the worst. However, finding the difference using visual inspections in the results is difficult. An example is shown in Figure 5.

The average tissue overlap values are listed in the second row of Table 4. The overlap values of FMSF, FMAF, MMSF, and MMAF are 0.8813, 0.8879, 0.8872, and 0.8917. The overlap values of the Maxwell model are larger than that of the fluid model, and the model with the adaptive forces performs better than that with the standard forces.

3.4. Real Datasets

The real datasets are acquired from the local hospital. The scans are acquired using a SIEMENS TRIO 3 Tesla scanner installed at the Institute of Biophysics of the Chinese Academy of Sciences. These scans are T1 sagittal images (TR = 1730 ms, TE = 3.93 ms, thickness = 1.0 mm, no gap, in-plane resolution = , slice = 192, and flip angle = 15). The scans are resampled as and the voxel size is . Thirty couples are randomly selected from the datasets. The mean time is listed in the fifth row of Table 1, and the SSD is listed in the third row of Table 3. The computing times of FMSF, FMAF, MMSF, and MMAF are 1256, 532, 1077, and 625, respectively. The SSD of FMSF, FMAF, MMSF, and MMAF are 0.0328, 0.0293, 0.0301, and 0.0279. The computation costs and matching accuracy of the real datasets are similar to those of the IBSR datasets. The overlap values are listed in the third row of Table 4. The overlap values of FMSF, FMAF, MMSF, and MMAF are 0.8823, 0.8912, 0.8892, and 0.8920. The results are similar to those obtained in Section 3.3. The proposed method is also compared with a method using fluid model and mutual information (FMMI) [32]. As the second column of Table 4 shows, the overlaps of the IBSR dataset and the real dataset from FMMI are the smallest, respectively. This indicates that the SSD is better than mutual information in monomodal images.

4. Conclusions

The proposed method is driven by the fluid and elastic models [1517], which uses the Maxwell model, a linear viscoelastic model that combines the properties of elastic and fluid models, to represent the image deformation. The proposed method introduces an adaptive force to speed up the registration.

The performances of the elastic and fluid models are compared in [17]. Therefore, we only compare the proposed method with the fluid method in this paper. The successive over relaxation method is used to solve the corresponding PDE, which is not the fastest but the most accurate among the evaluated methods. The computational cost can be reduced if PDEs are solved quickly, such as when using filter convolution [19] with a small filter width and parallel computing. However, the relative computation costs should be the same as those obtained in this paper. Actually, the fluid model is a special case of the linear viscoelastic model. When and , the linear viscoelastic model becomes a fluid model.

Our experimental results show that the linear viscoelastic model has several potential applications and that adaptive force can greatly reduce the registration time. The proposed method includes many parameters, and these should be analyzed further. We would also like to analyze the characterization of the transformation and how to obtain the optimal parameters for the corresponding transformation in the future.

Conflict of Interests

The authors declare that they have no conflict of interests.


This research was supported by the National Basic Research Program of China (2010CB732505 and 2013CB328806), the Key Projects in the National Science & Technology Pillar Program (2013BAI01B01), the National Hi-Tech Research and Development Program (2013AA013703), and the National Natural Science Foundation of China (61272360).