New York University School of Medicine, New York, USA
Received date: July 29, 2014; Accepted date: October 06, 2014; Published date: October 13, 2014
Citation: Huang X (2014) Nonrigid Image Registration Problem Using Fluid Dynamics and Mutual Information. J Biomet Biostat S12:004. doi: 10.4172/2155-6180.S12-004
Copyright: © 2014 Huang X. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are are credited.
Visit for more related articles at Journal of Biometrics & Biostatistics
Image registration is one of the fundamental tasks within image processing. It has wide applications in the fields of medical imaging, computer vision, statistical modeling etc. It is required when one wants to combine valuable statistical information from multiple images, possibly acquired using different modalities, at different time points or from different subjects, or to compare or integrate the data obtained from same or different measures. The subject of this article is nonrigid image registration. In particular, the main focus is on the application of fluid dynamics and mutual information (MI), and the comparison of two different similarity measures, i.e. mutual information and sum of squared differences (SSD). Numerical experiments show that fluid registration using SSD is ideal for mono-modal image registration, while fluid registration using MI does a better job in multi-modal image registration.
Nonrigid image registration; Fluid dynamics; Mutual information; Sum of squared differences; Mathematical and statistical modeling
Image registration is an often encountered problem in many application areas such as, medical imaging, satellite imaging, computer vision, or statistical modeling [1,2]. The task of image registration is to find an optimal geometric transformation that maps one image onto the other image, such that the transformed image is similar to the other image and corresponding points match [3,4].
Image registration adds value to images, e.g. by allowing multiple images acquired using different modalities to be viewed and analyzed in the same coordinate system, and facilitates new uses of images, e.g. to monitor and quantify disease progression over time in the individual or to build statistical models of structural variation in a population.
Though the image registration problem is usually easy to state, it is hard to solve. The main reason is that the problem is ill-posed. Small changes of the template images can lead to completely different registration results. Moreover, the solution may not be unique.
In practice, the concrete types of the geometric transformation as well as the notions of optimal and corresponding depend on specific applications. Each application has its own demand with respect to the meaning of similar and optimal. This is another reason for making image registration such a challenging task. For example, for the registration of X-ray images of bones, one may consider only rigid transformations, whereas for the registration of brain MRI scans one has to consider nonrigid transformations. The way one measures the similarity between images also depends on specific modalities of the images.
To obtain accurate results some applications allow for timeconsuming computations while most applications demand real-time implementation. In this age of Big Data, one always intends to speed up the existent algorithms to improve the performance of the applications. Thus, one also has to compromise between complex modeling and accuracy on one hand and computing time and storage requirements on the other hand. A fast and efficient algorithm is what researchers go in quest of.
A very important application of nonrigid registration is to segment a study image using an electronic atlas. Precise registration allows syntactic and semantic information from the atlas to be transferred to individual images. Another application of registration to an atlas is accumulation of data from images of many patients.
There are basically two factors that influence the classification of image registration methods: the motion model that determines what transformations are allowed and the driving force that drives the motion of the images.
Motion models are usually classified into two categories: rigid and nonrigid transformations. Rigid registration methods allow images of the same object to be registered, which are relatively easy to model and are not covered in this article. Nonrigid registration methods can be used to combine images from either time study of the same object or different objects. It is important to understand that the aim of the transformation is to map one image completely onto the other image in such a way that information from the first image can be applied to the other as well. Practically this means that the transformation must accommodate both high complexity and large deformations.
The possible driving forces are usually related to some similarity measures, so they can be classified into intensity-based, correlationbased and mutual information-based similarity measures.
Bajcsy and Kovacic  were the first to demonstrate volumetric nonrigid registration of medical images. They modeled the image as a linear elastic solid and deformed it using forces derived from an approximation of the local gradient of a correlation-based similarity measure. Miller, Christensen et al. [6,7] also used a global elastic model, but derived the driving force from the derivative of a Gaussian sensor model. These approaches have all suffered from the use of small deformation assumption.
Christensen et al. [8,9] extended their work and proposed the use of a viscous fluid model to control the deformation. It has the advantage of allowing for large curved deformations over the elastic models. Unfortunately, the fluid model requires the repeated solution of Navier-Stokes equations. Christensen solved the equation using the Successive Over-Relaxation (SOR) algorithm, which was computationally expensive.
Bro-Nielsen and Gramkow [10,11] proposed the use of convolution to approximate the solution of Navier-Stokes equations, and presented Christensen's fluid registration algorithm in a multiresolution framework, which sped up the fluid registration at least an order of magnitude. However, the driving force they used was essentially intensity-based, which restricted the use of their model to mono-modality registration. For images that are acquired from several different devices, their intensities cannot be taken directly to measure the image similarity.
In this article, we propose a new approach to do fluid registration. We extend the approach of Christensen's viscous fluid model by applying the convolution of the elastic filter to solve the Navier-Stokes equations and the maximization of mutual information to measure similarity. In order to compute mutual information one has to estimate the intensity distribution of the given images. This is a tricky problem for discrete digital images. We use a Parzen-window estimator by convolving the histogram with a truncated Gaussian function. We also compare the performance of fluid registration with different similarity measures.
The article is structured as follows.
In Section 4, a generic nonrigid image registration problem is described. Stated as a minimization problem, it can be formulated under a unified variational framework and deal with different applications by specific similarity measures and regularizers.
Section 5 introduces the viscous fluid model and discusses its application for nonrigid registration.
Two commonly used similarity measures, sum of squared differences and mutual information, are presented in Section 6. The corresponding driving forces are also formulated there.
Section 7 includes numerical results of fluid registration using sum of squared differences in both two-dimensional and three-dimensional situations, and also compares that with results of fluid registration using mutual information in two-dimensional case. Section 8 concludes with some discussion.
An image can be defined as a two-dimensional or three-dimensional function, and the amplitude of the function at any pair or triplet of spatial coordinates is called the intensity of the image at that point. The term gray level is used often to refer to the intensity of monochrome images. Color images are formed by a combination of individual 2D or 3D images.
In a computer, an image appears as a set of scalar values ordered in a two- or three-dimensional array. We view an image I as a function defined over a two- or three-dimensional bounded domain Ω with smooth boundary The range of an image will be considered to be in the interval [0, 255] throughout the article, in other words we only deal with 8-bit gray-scale (black-and-white) images here.
Image registration aims at finding a suitable spatial transformation such that a transformed image becomes similar to the target image. The task of image registration is to find an optimal geometric transformation between corresponding image data. The registration problem typically occurs if the images of the same object are taken from different perspectives, times or modalities.
Usually, one of the images is viewed as a reference image R and the other one as a deformable template image T. Given R and we are looking for a smooth and invertible transformation such that the deformed image is similar to R.
Obviously similarity between images needs to be measured in an appropriate way. If the images are acquired using the same modality, it makes sense to measure their similarity by their corresponding intensities. However, if the images are taken from different modalities, there might not be a correspondence between the intensities of R(x) and Tφ(x) for an optimal φ. Therefore we may consider other similarity measures which are not related to intensities but based on the intensity pattern of the images.
Classification of image registration
Image registration algorithms can be classified according to the transformation models one uses to relate the template image space to the reference image space.
The first broad category of transformation models is called rigid transformation. It is composed of a rotation followed by a translation. The object in the image therefore retains its shape and form under a rigid transformation. Extensions to the basic rigid transformation include other linear transformations, such as scaling and other more general affine transformations. These transformations are global in nature, thus, they cannot model local geometric differences between images. In most applications a rigid deformation model does not suffice and complex deformations must be estimated. Also rigid registration problems can often be reduced to parametric image registration problems, which are relatively easier to formulate, so we do not discuss these models in this article.
Our main focus is on the second category of transformation models, i.e. nonrigid transformations. These transformations are capable of locally warping the template image to align with the reference image. They are basically nonlinear transformations including elastic models, viscous fluid models etc.. Nonrigid registration aims at recovering a dense field of displacement vectors that maps each pixel individually in one image onto its corresponding pixel in the other, allowing the registration to adapt to local distortions instead of being restricted to global alignment of images. Usually we deal with them as non-parametric registration problems. The idea behind this type of registration is to come up with an appropriate measure for the similarity of images as well as for the likelihood of a non-parametric transformation.
Therefore we can set up a general framework for nonrigid image registration. This framework is based on a variational formulation of the registration problem, and the numerical schemes to be considered are based on the Euler-Lagrange equations which characterize a minimizer.
A general framework of nonrigid image registration
A general approach to image registration is based on a similarity measure and a regularizer. The similarity measure can be viewed as the driving force of the registration while the regularizer controls the transformation.
Given two images, a reference image R and a template image (we may require the boundary ∂Ω to fulfill some regularity constraints e.g. that of being of class C2), we are interested in looking for a diffeomorphic (i.e. differentiable and isomorphic) transformation such that the deformed image
For the following discussion it is also convenient to split the transformation φ into the trivial identity part and the displacement field u, i.e. where we exploit an Eulerian reference frame. Therefore, we are now looking for such a diffeomorphic displacement field u instead of φ. The displacement vector field u is searched for in a set of admissible vector fields such that is a linear subspace of a Hilbert space H, which is endowed with a scalar product
The most intuitive way of approaching the registration problem is to minimize a joint-functional
with a similarity measure D and a regularizer S. The similarity measure D rates the similarity between R and the deformed template The regularizer S penalizes unwanted deformations. The parameter α weights the similarity of the images versus the smoothness of the displacement. In summary, the nonrigid registration problem is defined as the solution of the following minimization problem,
To compute a solution of the minimization problem (4.2) we make use of the fact that the first variation of the joint-functional J has to vanish for a minimizer, i.e. u is a minimizer Since we can express the first variation of each term as follows,
where the force field f(x, u(x)) depends on the particular similarity measure D and is used to drive the flow, and A is a partial differential operator, the minimization problem (4.2) boils down to a problem of solving the following partial differential equations (PDE) system,
Rather than solving the system directly, a gradient descent strategy is used to search for a minimizer of the joint-functional J. Given an initial estimate u0 ∈ H, one can introduce time and a differentiable vector field, also noted as u from the interval [0, T] into H and solve the following initial value problem:
It has been proved (Faugeras et al. ) that if the linear operator A is strong elliptic, invertible and the nonlinear f is bounded and Lipschitz continuous, there exists a unique classical solution to the above initial value problem.
Since image registration is an ill-posed problem, regularization is essential and inevitable. The regularizer is used to pick out the most likely solution. Actually, it is the regularizer that distinguishes the image registration methods. Since elastic models constrain the possible deformation by a compromise between internal and external forces, elastic displacements do not reach the desired deformation. In a viscous fluid model, internal forces disappear over time and the desired deformation can be fully achieved. So the fluid registration method satisfies the general requirements of both complex and large deformations and we therefore regard the fluid registration as the more advanced registration method available.
In terms of the similarity measure, we would like to compare two different methods, sum of squared differences (SSD) and mutual information (MI). The sum of squared differences method is easy to implement and widely applied in a lot of mono-modal image registration, while the mutual information method can successfully deal with multi-modal images.
The intention of fluid registration is to mimic a flow of fluid in a certain sense . Fluid registration is very useful for some medical applications, e.g. the registration of images obtained from different brains . Moreover, this approach allows for smooth and large deformations while maintaining continuity. Although the viscous fluid models seldom have been actual models of the physical phenomena, they have enjoyed great success.
The underlying partial differential equations (PDE) for fluid registration is
Where Δ = ∇2 is the Laplace operator, ∇⋅ is the divergence operator, f is the force field working on a volume element in the interior of the body, v is the velocity field, and μ and λ are the viscosity constants. The term Δv is called the viscous term because it constrains and smooths the velocity field spatially. The term ∇(∇⋅v) limits the gradient of the divergence of the field and allows for contraction or expansion of the fluid. Thus, the viscosity coefficients μ and λ represent the smoothness and the mass injection or compressibility of the fluid. The viscosity constant λ is adjusted to control the rate of growth or shrinkage of local regions within the deforming template, and the viscosity constant μ is adjusted to control shearing between adjacent regions of the template.
In particular, the relationship between the displacement and the velocity in the Eulerian reference frame is given as follows (the material derivative)
Christensen crafted the PDE (5.1) using principles obtained directly from continuum mechanics for deformable bodies. He modified the classical conservation equations of mechanics to account for non-mass-conserving deformations, which allowed for the growth and shrinkage of regions within the continuum.
Now we formulate a variational model for fluid registration which is different from Christensen's approach.
Define the regularizer as
Then the minimization of the energy joint-functional
is equivalent to solving
Here for simplicity we set the coefficient α as 1. For a constant force f this is the PDE for linear elasticity working on the instantaneous velocity field v. Consequently, the equation works by elastically smoothing the instantaneous velocity field of the fluid in the Eulerian reference frame.
Solution of the viscous fluid registration problem requires solving a system of equations determined by the characterizing viscous fluid PDE, the material derivative, and the force field equations. We will discuss the details of the force field equations in the next section.
This system of equations includes nonlinearities in both the force and the material derivative. To solve it, Euler integration is applied over time, using a forward finite difference scheme of the time derivative in (5.5),
The solution can be found by iteratively solving (5.4) for the instantaneous velocity, and integrating over time using (5.7).
From the Euler integration step (5.7), we note that it requires a non-zero Jacobian of the transformation φ
which is the gradient of the transformation φ with respect to the Eulerian reference frame.
Since the entries of this Jacobian matrix are available during the computation, it is easy to check whether this matrix is regular or not. If the matrix is not regular, the transformation φ is not bijective and the step is not admissible. On the other hand, the transformation is required to be homeomorphic. In practice, we require the transformation to be bijective and orientation-preserving or to have a positive Jacobian, i.e.
The transformation becomes singular for large curved transformations. To circumvent this undesirable situation, Christensen  suggested a regridding method, which turned out to be the key remedy to get the wanted deformation.
Solving the linear PDE
The core problem of the fluid registration is solving the linear PDE,
for constant force and time, since it is the most time-consuming part of the fluid registration.
For constant force term f the PDE is linear, and the linear operator is just the Linear elastic operator acting on the velocity field v. linear elastic problems are normally solved using implicit finite element or finite difference methods. However, in the case of images, the size of the problem is huge and in practice unsolvable with these techniques.
We extend the approach proposed by Bro-Nielsen et al. [10,18] to solve the linear PDE with multi-resolution convolution. Using the linearity of the PDE and the superposition principle, a linear elasticity filter can be created as the impulse response of and then applied the filter to f.
Eigenfunction basis of the linear elastic operator:
First we focus on the linear elastic operator
where u denotes the displacement field.
Since mapping the boundary of Ω onto the boundary of Ω is important in our work, the sliding boundary conditions are considered. They are defined for the domain using the Dirichlet boundary conditions,
and Neumann boundary conditions,
The sliding boundary conditions map the unit square onto the unit square, in such a way that boundary points are allowed to slide along the boundary.
With these boundary conditions the eigenfunctions and the eigenvalues have been derived as follows. See  for the complete derivation.
The orthogonal eigenfunctions are given as
The eigenvalues corresponding to the eigenfunctions are
For i=j=0, α1=α2=0. Otherwise they are chosen to ensure that the eigenfunctions are normalized as follows,
3D Eigenfunction basis:
For completeness we also give the eigenfunctions and eigenvalues for the 3D version of the linear elastic operator (5.10)
with the sliding boundary conditions on the domain [0, 1]3.
The orthogonal eigenfunctions are
The eigenvalues corresponding to the eigenfunctions are
For i=j=k=0, α1=α2=α3=0. Otherwise they are given by
Decomposition of forces onto the eigenfunction basis:
In two-dimensional case, to solve (5.9), we now can take advantage of the work Christensen et al. have done for solving (5.10).
We consider the finite truncation of the decomposition of the velocity field under the derived orthonormal basis as follows
Where aijr are the eigenfunction coefficients for the velocity field. Now we apply the linear operator to the field
and take the inner product defined as of the equation with the orthogonal basis
Therefore the eigenfunction coefficients are simply the projection by the inner product of the force vector field onto the eigenfunctions, scaled by the inverse of the eigenvalues.
Convolution filter for linear elasticity:
Now we can derive the convolution filter to solve the linear PDE (5.9) as in . The impulse response of the linear operator is first computed under the linear elasticity eigenfunctions. To get a discrete filter for digital images the impulse response is discretized.
Suppose Ω = [0,1]2 , the impulse response in the direction of x1 is determined as the velocity field corresponding to an impulse force fc = [δ (x-c), 0]T applied at c=[0:5, 0:5]T in the middle of the domain, where δ (x-c) is Dirac's delta function.
Using the previous result, we can get the decomposition coefficients of the impulse response by applying the decomposition projection to the impulse force fc
Where is the x1-coordinate of
With the coefficients known, we obtain the impulse response
This is the response of the linear operator for an impulse force in the x1 direction applied at c. The impulse responses for the other directions are determined by simple rotation of the response for the x1 direction.
Note that the decomposition of the impulse response based on the eigenfunction basis is actually a frequency-based decomposition. Large i and j correspond to high frequencies and small to low frequencies. Hence an ideal low-pass filtering of the impulse response can be implemented by truncating the sequence at number N instead of summing up to infinity.
The sampled filter is defined with dimensions D×D, D odd, in the domain [0, 1]2. The sampling interval is consequently which Shannon's sampling theorem relates to the cut-off frequency f by From equation (5.25) the frequencies corresponding to the summation variables are
and the common truncation point becomes i=j=N=D-1.
Consider a 2D filter of sizeD×D, D odd, and let y=[y1,y2]T, where
The filter implementing the 2D linear elastic operator for the x1 direction is then
A filter for 3D linear elasticity can be derived in a similar fashion . Consider a 3D filter of size D×D×D, D odd, and let y=[y1,y2,y3]T, where The filter implementing the 3D linear elastic operator for the x1direction is then
A few remarks concerning the implementation of fluid registration scheme are listed.
1. Initialization: let i=0 and t0=0;
2. Set u(x, ti)=0;
3. Calculate the force vector field f(x, u(x, ti));
4. If f(x, u(x, ti)) is below the given threshold for all x, STOP;
5. Solve the linear PDE (5.9) for v(x, ti)) by implementing the linear elasticity filter;
6. If the Jacobian J is below the other threshold, regrid the template andgo to 2;
7. Choose the time step so that where is the maximal ow distance allowed in one iteration;
8. Update the displacement field u(x, ti) (Euler integration step);
9. Let i = i + 1, go to step 3.
Every time the Jacobian drops below the given threshold, a new template is generated by applying the current deformation to update the template. In addition, the displacement field is set to zero, while the current velocities remain unchanged. It helps us to evade singular transformation and keep getting the mapping we want. The total deformation becomes the concatenation of the displacement fields, associated with the sequence of propagated templates.
Computation of the forces:
Different similarity measures may provide different force vector fields. The sum of squared differences is suitable for mono-modal image registration because it compares the intensities of the reference and template images to measure their similarity. On the other hand, mutual information is a better similarity measure for multi-modal image registration. Mutual information essentially measures the entropy of the joint density; it is maximal if the images are maximally related. Thus we can find the desired transformation via the maximization of mutual information.
Because of the limited span of the discrete filter, one can implement the viscous fluid registration using the filter in multi-resolution. The fluid registration is first performed on a rough scale. The result of this scale is then propagated to a finer scale and the fluid registration restarts there. This process is continued down to the finest scale of the scale-space, yielding the final registration result.
An image similarity measure quantifies the degree of similarity between intensity patterns in two images. It provides the force vector field in the fluid registration. The force field is the link between the physical model of the viscous fluid and the image data. It is crucial for the success of the registration algorithm.
An image similarity measure can take on many forms, sometimes it depends on the modality of the images to be registered. The sum of squared intensity differences (SSD) is a similarity measure commonly used for registration of images in the same modality, because ideally the intensity of the corresponding pixels in those images should be identical.
A natural extension of SSD is to use the maximum likelihood estimator (MLE) in statistics as the similarity measure, which appears more stable with respect to the parameter α. It is worth noting that the ordinary least squares (OLS) estimator (which is essentially the minimization of SSD) is identical to the MLE under the normality assumption for the error terms. However, there is usually no such direct correspondence between pixel intensities in multi-modal images. Even for images in the same modality, due to different experiment conditions perfect intensity match can rarely be achieved. Alternatively we can measure the similarity between images by intensity patterns instead of intensity values. Mutual information (MI) is the most popular image similarity measure for registration of multi-modal images.
Multi-modal image registration is of greater importance in image processing today, since it allows us to effectively synthesize information from images of different modalities of the same object to fully understand the structure and function of the object. It is also useful to take full advantage of the complementary information coming from multi-modal images. For example, as we know since computerized tomography (CT) and magnetic resonance imaging (MRI) are sensitive to different tissue properties, the appearance of the images obtained with the two techniques differ markedly. By registration we can see clearly both soft tissues and bone structure in the same object such as human bodies.
Sum of squared differences
A straightforward intensity-based approach to similarity measure is based on the minimization of the sum of squared differences.
The measure is defined as L2 norm of the difference of two images as follows,
The force field is determined as
f(x,u(x)) = (R(x) −Tu (x))∇Tu (x). (6.2)
This force field has two factors. The first factor ∇Tu (x) is the gradient of the template image. It reflects all intensity changes on the template image, such as distinct edges and ridges in the image. The second factor R(x) −Tu (x) scales the first factor, such that in positions with large intensity differences, deformation on the template is encouraged and in similar regions deformation is discouraged. The force field is therefore zero in regions that are locally matched.
A general way of comparing the intensity pattern of two images is to use some statistical or information-theoretic similarity measures. Among numerous criteria, the mutual information provides us a good way to measure similarity between images based on intensity distributions. The concept of mutual information is borrowed from information theory , and was introduced in the context of multimodal registration by Viola and Wells III .
We define the mutual informationas the reduction in the entropy of Y given X,
I(X,Y) = H(Y ) - H(Y|X). (6.3)
As Y becomes more dependent on X, H(Y|X) gets smaller and I(X|Y) gets bigger.
The mutual information is positive and symmetric, and measures the amount of information (entropy) that one random variable contains about another random variable. It is the reduction in the uncertainty of one random variable due to the knowledge of the other.
Random variables are considered independent when
H(Y|X) = H(Y );
I(X,Y ) = 0 = H(X) + H(Y ) - H(X,Y ).
In other words, random variables are considered most dependent to each other if their mutual information attains its maximum. This is the theoretical basis on which choosing mutual information as a similarity measure relies.
Mutual information as similarity measure:
In an image registration problem, we now regard the two images as random variables with the joint intensity density and marginal densities
We consider the definition of mutual information based on differential entropy. Mutual information measures how the intensity distributions of two images fail to be independent. Since mutual information essentially measures the entropy of related joint density, it is maximal if the images are maximally related.
By adopting the approach of Hermosillo et al. , we derive the corresponding force field,
The force field also has two factors. It has the gradient of the template image as one factor as well. Therefore it is of interest to interpret the behavior of the other factor.
The first term in FR,Tu (r, t) , namely tries to make the intensity t move closer to a local maximum of It thus tends to cluster On the other hand, the second term, tries to prevent the marginal distribution pTu from becoming too clustered.
Parzen-window density estimation:
In order to compute mutual information we have to estimate the intensity distribution of the given images. The accuracy of this estimation plays a very important role for the outcome of the registration as well as the needed computational effort.
To evaluate the force (6.4) we have to estimate the densities pR,Tu and pTu. We use a Parzen-window estimator as follows,
Where X1,X2, …,XM denote realizations of the random variable X, K is the Parzen-window function and Kh(x)=(1/h)K(x/h). The second equality holds because we always take kernel width parameter h as 1 in image processing.
As Parzen-window function we use a truncated Gaussian
where the standard deviation σ determines the width of the window.
In practice the Parzen-window estimator is much more flexible than other parametric density estimators, because it requires only the smoothness of the density without any assumption about the functional form of the density.
Intuitively the Parzen-window estimator computes a local, or windowed, average of the sample. The window function defines a region centered on x in which sample points contribute to the density estimate. Points that fall outside of the window do not contribute. Getting a reliable estimate depends on having a reasonable number of points fall into the window around the query point x.
2D Fluid registration using SSD
First we did a two-dimensional fluid registration using the sum of squared differences as our similarity measure. The reference image was a square of intensity 255 and the template image was a disk of the same intensity.
We may take it as a synthetic example to show the performance of fluid registration using SSD, since for this case, we may directly determine the similarity between two images by the difference of their intensities.
The numerical result was perfect. We successfully registered the disk to the square (Figure 1).
Then we moved on to a real cell image. The reference image was a cell and the template image was generated by a horizontal translation by 40 pixels and a clockwise rotation by 30 degrees to the reference image. Note that here actually we performed a rigid transformation to the reference image.
Again we may take these two images as acquired from the same modality due to their correspondence of intensity. Therefore it was appropriate to use the sum of squared differences as our similarity measure.
Applying fluid registration which is essentially nonrigid image registration, we could register the template image satisfactorily (Figure 2). It was very impressive for that the internal texture of the cell had been recovered quite accurately.
From Figure 3 we could see that the difference between the registered template image and the reference image was almost negligible. The result was satisfactory as we expected.
3D Fluid registration using SSD
We implemented fluid registration method in three-dimensional situation as well. The biggest problem for three-dimensional registration is the memory for storage. Thus we could only restrict ourselves to a synthetic experiment to show the performance of our algorithm.
We used a cube of intensity 255 as our reference object and a sphere of the same intensity as our template object.
By fluid registration with SSD we successfully deformed the sphere to a decent cube very close to the reference object (Figure 4).
Fluid registration using MI
We tried to do two-dimensional image registration between a disk and a square of different intensities . The disk is of intensity 255 and the square is of 128.
We could take it as a synthetic example of multi-modal image registration. Our main concern now focused on the intensity pattern only. In other words, by deforming a disk we intended to form a square, of whatever the intensity values it would be, as the one in the reference image.
First we applied SSD as the similarity measure in this image registration. We checked the initial step for generating force vector field and then computing velocity vector field. Since the initial displacement field is set as zero vector fields, the velocity field obtained is actually the first deformation to the image.
The initial force and velocity fields were shown in Figure 5. The velocity fields had both positive and negative values on the same layer, which means at the first step some particles would move out and some would move in. Hence there was no hope for the disk to fully expand into a square by SSD.
Figure 5: 2D fluid registration using SSD with μ=10, λ=0: (a) is the reference image of a square of intensity 128; (b) is the template image of a disk of intensity 255; (c),(d) are the horizontal and the vertical components of the force field in the initial iteration; (e),(f) are the corresponding displacement fields. (white is positive value and black is negative).
Then we turned to apply mutual information as the similarity measure. We showed the corresponding initial force and velocity fields in Figure 6. This time by the homogeneous initial velocity field (essentially the first movement), obviously the disk tended to expand outwards, to a square hopefully.
Figure 6: 2D fluid registration using MI with μ=10, λ=0: (a) is the reference image of a square of intensity 128; (b) is the template image of a disk of intensity 255; (c), (d) are the horizontal and the vertical components of the force field in the initial iteration; (e), (f) are the corresponding displacement fields. (white is positive value and black is negative).
The result was shown in Figure 7. Using SSD the disk could at most grow to a bigger disk with a lot of holes inside it. But if we resorted to MI as our similarity measure, a clear shade of a square of the same size as the one in the reference image has been fully recovered. The square was distinguishable but not that obvious because of the sharp contrast between different intensity values in the image.
The advantage of mutual information over sum of squared differences is more clear when we compared them in the same image registration problem (Figure 8). Using SSD the best we could do is to reconstruct a disk connected with four long sharp spikes pointing outwards. Since those spikes were always connected to the disk, we could not remove them as some noises. The pattern appearing in the result would at best be a disk with four corners. However we almost perfectly reconstructed a disk of the same size as in the reference image by MI without any corners.
The research on doing fluid registration with MI on 2D real images and even 3D images is still on-going.
We implemented fluid registration using both sum of squared differences and mutual information. By the synthetic examples in twodimensional and three-dimensional situations, we can conclude that SSD performs very well as the similarity measure in the mono-modal image registration. The fluid registration model also played a key role in the experiment. In the cell image registration, the obscure internal texture had been successfully reconstructed by the fluid dynamics model. Since the force field for SSD is easy to compute, it is an ideal efficient similarity measure for the mono-modal image registration.
But it has no hope to prevail in multi-modal image registration. In the experiment of transforming a disk of 255 to a square of 128, the first movement by the computed displacement was inhomogeneous, which indicated the failure of the process. On the contrary, using MI gave the expected homogeneous movement at the beginning which showed great prospect of success in the experiment.
The comparison was more self-evident in the other experiment of deforming a square of 128 into a disk of 255. With SSD as the similarity measure, the best we can do is to obtain a disk with four long sharp spikes. While using MI we could almost perfectly reconstruct the disk as long as time permitted.
The only drawback for using mutual information as the similarity measure is the computation time. To compute a force field, it involves a lot of computation in estimating the probability density functions of two images, their joint probability density function and the derivatives of these functions. With images of big size or high resolution, the computation would take quite a long time. A fast algorithm is necessary for us to perform the two-dimensional and even three-dimensional fluid registration.