Composite match index with application of interior deformation field measurement from magnetic resonance volumetric images of human tissues.

Whereas a variety of different feature-point matching approaches have been reported in computer vision, few feature-point matching approaches employed in images from nonrigid, nonuniform human tissues have been reported. The present work is concerned with interior deformation field measurement of complex human tissues from three-dimensional magnetic resonance (MR) volumetric images. To improve the reliability of matching results, this paper proposes composite match index (CMI) as the foundation of multimethod fusion methods to increase the reliability of these various methods. Thereinto, we discuss the definition, components, and weight determination of CMI. To test the validity of the proposed approach, it is applied to actual MR volumetric images obtained from a volunteer's calf. The main result is consistent with the actual condition.

1. Introduction

The physical property is the base of the biological simulation, computer-assisted medical applications, such as clinical diagnosis, and surgical simulation, surgical planning. And estimation of internal deformation field or deformation motion for the biological tissues plays a very significant role in physical parameters estimation. Thus, measuring the internal deformation field of biological tissues is becoming the focus research. Magnetic resonance (MR) imaging (MRI) provides superb anatomic images with excellent spatial resolution and contrasts among soft tissues; thus, it is widely used in computer-assisted medical applications, such as clinical diagnosis, surgery simulation, operation planning, and evaluation of physical characteristics of biological tissues. Increasing number of researchers in medical simulation and medical virtual reality focus on the interior deformation field or motion measurement of biological tissues from MR volumetric images, and it has become one of the significant branches of medical image analysis. Generally, approaches for estimating the deformation of MR volumetric images can be classified into two typical types: elastic deformation model-based and feature matching-based methods.

The elastic deformation model-based method can be classified into either parametric or geometric active models [1]. To obtain the deformation information of an object, the parametric active contours, also called snakes, try to minimize a defined cost function so that the function deforms a given initial contour toward the boundary of the object. This method was first introduced by Kass et al. in 1987 [2] and subsequently developed and used by Lang et al. [3], Cho and Benkeser [4], and Matuszewski et al. [5] to estimate deformation motion of nonrigid objects. In the geometric active model [1, 6–8], the curve and the surface of an object are first detected. Then, the deformation propagation of the curve and the surface is used to track the motion. However, irrespective of what elastic deformation models are employed, disadvantages exist in the deformation estimation; for example, the parametric active model cannot handle changes in the topology of the evolving contours when deformation is performed directly, and often, heuristic topology handling procedures are used [8]. In the geometric active model, when contrast is poor and boundaries are not clear or are continuous in the images, the contours tend to leak through the boundary [9]. The tagged images must have a regular grid pattern in the imaging plane because if the number of tagged points is low, the measurement accuracy would be poor. More important than the former two aspects, regardless of what elastic deformation models are used, they can only handle the deformation at the boundary of nonrigid objects and not the interior deformation.

In recent years, researchers have been increasingly concerned on approaches for matching of nonrigid feature points. Typically, thin-plate spline-robust point matching (RPM) is a famous algorithm for matching non-rigid feature points, which can estimate the joint correspondence and non-rigid transformations between two differently sized point sets. However, optimal processing of the energy function utilized in Chui's method may be trapped in bad local minima [10]. Zheng proposed the RPM-local neighborhood structure (LNS) method of matching non-rigid feature points, based on the supposition that relative distances and orientations among feature points in a neighborhood would be preserved [11]. Lee improved the LNS and presented the topology preserving relaxation labeling (TPRL) algorithm. In the TPRL method, log distance and polar angle bins are utilized to capture the coarse location information of the feature points in a neighborhood. Using shape contexts, Belongie proposed a non-rigid point matching method. In this method, every feature point is represented by a histogram descriptor of the distance and orientation between this feature point and its neighbor feature points [12]. In addition, some other useful methods were also proposed for feature-point matching, such as the coherent point drift matching method of non-rigid points [13, 14] and the preservation of local geometrical characteristics [15]. In these methods, a novel objective function is defined to preserve local image-to-image affine transformations across correspondence. In general, some unsolved problems are involved in the aforementioned matching methods of non-rigid points; for example, the optimal processing of the energy function could be trapped in bad local minima, the topology of the neighboring feature points is not always preserved well, and so on. Most importantly, in these methods, useful information of the feature point is considered singly and lacks a comprehensive approach, which can mix up with the useful and significant information in the point matching of deformation measurement.

Therefore, to improve further the proposed feature-matching-based approach and improve the robustness of the matching result, this paper proposes a composite match index (CMI). In Section 2, we introduce some previous work, in Section 3, we describe the concept and definition of CMI, and Section 4 introduces the CMI application on feature matching of image pairs from non-rigid objects. In Section 5, examples and preliminary experimental results are given, and discussion and conclusions are presented in the final section.

2. Previous Work

Feature matching plays a significant role in human visual perception, recognition, and computer vision. In medical imaging, most existing feature matching-based research has focused on non-rigid registration and internal deformation field measurement. The general idea of these works is first to extract enough feature points or markers from medical images acquired from non-rigid objects on natural and deformed states, respectively. Next, the feature matching algorithm is applied on extracted feature-point sets to establish robust corresponding pairs. Finally, corresponding pairs are used as control points in non-rigid registration and are used to calculate sparse deformation fields in internal deformation field measurements. Therefore, finding robust corresponding pairs is a vital problem in the present work. We surveyed existing works on feature-point matching in computer vision. Relaxation is a valid technique to disambiguate matches and improve the robustness of matches. Finding a globally optimal or reasonably good suboptimal solution in relaxation is a difficult task, and such matching techniques in non-rigid medical image processing have been rarely addressed. However, a potential advantage is that harder matching problems can be solved using global optimization techniques.

Papademetris et al. [16] presented a method for the integration of feature and intensity information for non-rigid registration. In this case, a distance-based robust point matching framework was proposed to estimate feature-point correspondences. A disadvantage of the algorithm is that it estimates transformation using weighted least squares, which affects the strength of matching.

Zhang et al. [17] introduced a feature matching-based algorithm and considered the problem of 2D deformation field measurement as an example. Matching strengths are measured using correlation and relative distance between two feature points. Relaxation by the optimization algorithm is deductive of the function of matching strength. In later research [18], after slight revision, the algorithm was extended to a 3D situation because the intensity in a magnetic resonance (MR) image is the information of tissue mapped on an image. Thus, the correlation intensity of regions between two points in matching and relaxation can effectively use the properties of tissue.

The work [19] proposed a local geometric preserving algorithm to find corresponding feature pairs from given feature points set in MR volumes acquired from an object on natural and deformed states, respectively. The main contribution of the algorithm to feature matching is that for a non-rigid tissue, when an outside force is applied on it, the deformation magnitude and orientation are different in different regions. However, for a local region on the object, the difference is actually very slight and can sometimes be ignored.

Problems in image feature-point matching remain as great challenges for medical image processing. Thus, the accuracy of feature matching needs to be further improved. Typically, single factors, such as intensity and distance, are effective in matching algorithm for specific areas. However, total accuracy cannot be improved. The integration of multifactors to form a composite approach can make use of the advantages of each factor to improve total accuracy. The present work proposes a composite framework that can pose multicomponents in a single cost function with associated weights to find corresponding feature pairs.

3. CMI

CMI-based feature-point matching approach was proposed to address the fusion of different operator types and to improve the reliability of results from single operators. Here, CMI is a scalar quantity that describes the matching possibility of point pairs. Let be the vector of component value and its corresponding weight, respectively. Then, according to the linear weighting method, CMI is defined as subject to where ξ represents the CMI of the i-th pair, c represents the value of the kth component consisting of the CMI, w represents the weight of the kth component, and K is the number of components in the CMI. Here, component is a factor that can be used to evaluate feature-point pair similarities. Weight w is used to measure the significance of a component for CMI. Various weighting methods have been reported for different research fields. In this case, to consider the independence of each component, the correlation weighted method is used to determine the weight of each component. Let r = [r 1,1, r 1,2,…, r 1,] be the correlation vector consisting of correlation score of the component 0 and k. Then, the weight of the t + 1 time w is defined by with where c is the value of the i-th feature-point pair and N is the total number of match pairs in the potential matching set obtained at time t.

Since feature-point pairs within the potential matching set obtained at time t are used as samples to compute the weight w of the kth component in t + 1 times iteration, the pairs in potential matching set are different at each time. Thus, w values are also different at different times, keeping iterations in the matching process.

CMI is an effective way to fuse multifeature matching algorithm. CMI takes full advantage of all the considered factors to generate a more robust feature matching approach and obtain more accurate matching results. Thus, the feature-point matching algorithm, which decides the strength of matching via a similarity judge function, can theoretically be integrated as a CMI component. In this case, the local geometric persistence (LGP), local intensity similarity (LIS), and local correlation score (LCS) between regions around participants are selected as the components to compute the CMI of a match pair (p , p ) and demonstrate the validity of CMI. The following section will discuss how to compute LCS, LGP, and LIS.

For convenient descriptions, several definitions are first clarified as follows

Initial feature set p , feature-point set extracted from the MR volume acquired from the object at a natural state.

Deformed feature set p , feature-point set extracted from the MR volume acquired from the object at a deformed state.

PMS, a potential feature match set composed of a match pair (p , p ) if and only if the best match of p is p and conversely p is also the best match of p .

p represents the feature point i in the initial feature set, and p represents the feature point i in the deformed feature set.

3.1. LGP

Let c and c be the moment center computed using the initial feature set and its mapping in the deformed feature set, respectively, let p be the ith point in the initial feature set, and let the mapping in the deformed feature set be p . Based on the consistent deformation in a local region, the distance ratio of a potential match pair in a local region far from their moment center is equivalent and thus yields μ : where μ is the distance ratio of the jth potential match pair (p , p ) in the local region around pair (p , p ), d(p , c ) is the Euclidian distance between p and c , d(p , c ) is the Euclidian distance between p and c , and J is the number of potential match pairs in the local region. Ideally, μ should be a constant in the local region.

Moreover, d = [d(p , c ), d(p , c ),…, d(p , c )] and d = [d(p , c ), d(p , c ),…, d(p , c )] are the distance sets of the potential pairs within a local region around pair (p , p ), respectively. Based on the definition of mathematical expectation, we yield Thus, if p in the deformed feature set is the best match of a given feature p in the initial feature set, then, the geometric deformation of potential match pair (p , p ) within a local region around pair (p , p ) is defined as subject to In a small local region, all the g (j = 1,2,…, J) should be approximately identical and go to zero; the smaller the value of g , the better the geometric persistence of a potential match pair (p , p ). This is called geometric persistence in this case. Thus, the impact factor of the j-th feature pair for the LGP within a small local region is The geometric property within a local region is approximately consistent in the initial and deformed states. If a pair is the best match for each other, then the correlation of potential matches within a local region around the pair must be a strong one. The correlated score gc(p , p ) of the geometric persistence of PMS in a small local region around (p , p ) can represent the LGP of feature pair (p , p ), specifically: where J is the number of potential matches within a local region. In (11), if λ is large, the pair (p , p ) may be a strong match pair; thus, its weight must also be large. In addition, the value range of gc(p , p ) should be [−1,1]. Normalizing gc(p , p ) yields normalized LGP as

3.2. LIS

LIS is used to describe the intensity difference between regions around a feature-point pair in the initial and deformed volumes. As mentioned earlier, the tissue within a local region is the same in the initial and deformed states. Thus, based on the MRI principle, the intensity difference is small. The inner product between two regions has the same properties with the invariance of rotation, zoom in, and zoom out. The normalized inner product between regions around (p , p ) is adopted to define the similarity of two regions. Thus, where X is the region in the initial volume centered at feature p and X is the mapping region of X centered at feature p .

3.3. LCS

Let I(p ) and I(p ) be the intensity of the m-th voxel within the region centered at p and p in the initial and deformed MR volumes, respectively. Let 𝒪 be the local cubic region with a size of w × h × l. The local correlation score between local cubic regions around feature p in the initial MR volume and its candidate match feature p in the deformed MR volume is defined as where Here, σ 2(I(p )) and σ 2(I(p )) are the standard derivation of the local region 𝒪 around feature p and p , respectively. They are given by where and are the averaged intensity in the neighborhood of feature p and p , respectively.

4. Application in Feature Matching

This section describes the measurement of internal deformation fields using CMI. First, the cost function is given to obtain optimal feature pairs iteratively. Then, the actual feature matching algorithm is described. Finally, the internal deformation fields are measured using optimal feature pairs.

4.1. Cost Function

CMI is an index that measures the strength between a given feature and its candidate matches in feature matching. In theory, for a given reference feature, its potential match must have the strongest CMI among all the candidates. Thus, for an optimal potential matching set, its whole CMI will also be the strongest. Based on this idea, we yield where S is the cost function in iteration and N represents the total number of match pairs in the PMS obtained at time t.

4.2. Actual Matching Algorithm

The objective of the feature matching algorithm is to obtain an optimal PMS ultimately. The idea of PMS optimization is to maximize the aforementioned cost function S iteratively. In each iterative step, the current PMS strength is evaluated by all candidate matches within PMS using the defined cost function S. The iterative steps will stop until S no longer increases or is subjected to stop conditions. Specifically, the inputs are two feature-point sets obtained from MR volumetric images of an object under natural and deformed states, respectively. The output is an optimal PMS. The specific process of the algorithm is summarized as follows.

Compute LCS and LIS. For each given pair (p , p ) consisting of features in initial and deformed volumes, we use a local region (size of 9 × 9 × 3 in this case) centered at features to compute LCS and LIS according to (14) and (13), respectively.

Form initial PMS. The LCS is used as the initial CMI of each match pair in the step of initial PMS formation. In other words, LCS is the only criterion of this step.

Compute LGP. For each given pair (p , p ), we first search for neighbor potential matches within a small window (size of 17 × 17 × 3 in this case) centered at p . The potential matches contained within the window are participants in the LGP computation using the approach in Section 3.1.

Compute w. Compute the weight for each CMI component using potential matches in current PMS as samples. The specific computing method can be seen in (4).

Update the CMI of each pair. For each given pair (p , p ), its corresponding CMI is updated through the weighting sum of the components LCS, LIS, and LCP, which are computed in (0) and (2).

Form PMS and compute the cost function S. The updated CMI of each pair forms new PMS. The cost function in (17) is then computed using potential matches in the current PMS.

Repeat (2) to (5) until S no longer increases.

Return the current PMS.

Although candidate sets LCS and LIS of each pair are constant, PMS is dynamic because of the varying LGP and w of the component at t + 1 times iteration. Thus, the match strength index of CMI is varied. Dynamic cost function will move potential matches into or out of the PMS. The best candidate of a feature-point may also change.

4.3. Measuring Density Deformation Fields

After obtaining the optimal PMS, the internal density deformation fields of non-rigid objects are then obtained. In this study, the method proposed in our previous work [20] is used to obtain the internal density deformation fields. In summary, the internal density deformation fields are interpolated by sparse deformation fields using a finite element model. In detail, the magnitude of the sparse deformation field is first computed by its corresponding pair in PMS using Euclidian distance. The start and end points of a field direction are defined by the points of the corresponding pair. Next, a non-rigid object is reconstructed using tetrahedra, whose nodes are points in the PMS. The density deformation fields can then be interpolated using the finite element method.

Let P be an arbitrary volume voxel at x = (x, y, z) within a tetrahedron ◊P P P P consisting of nodal points P , P , P , and P . Its displacement may be approximated by weighting the finite element node displacements u (x) using their shape function [20]: where u (x) is the displacement of nodal i, and the shape function N (x) on tetrahedron ◊P P P P is given by where ♦ P P P P and ♦ P P P P are the volume of tetrahedron ◊P P P P and ◊P P P P , respectively.

5. Experiments and Results

Our approach consists of four steps: feature extraction, affine transformation, feature matching, and deformation field measurement. Extracting sufficient features from the initial and deformed volumes is necessary to find enough homologous feature pairs. In this study, high-curvature 3D points were preextracted as features from MR volumetric images. In this case, the two-dimensional Harris operator [21] was extended to a 3D operator by extracting features from the MR volumetric images [22].

Some experiments were designed to demonstrate the performance of the proposed approach. All experiments were performed using our own tool developed with Visual C++, which runs on Microsoft Windows XP. All described experimental results were obtained on a Lenovo Portable PC with a 2.20 GHz Intel(R) Core(TM) 2 Duo CPU T6600 and 4 GB of RAM.

In the experiment, the MR images were acquired from a volunteer's calf (Figure 1(a)) using an MRI scanner at natural state and deformed states (initial and under forcing), respectively. In both cases, the FOV was 20 × 20 cm, and the slice gap was 2 mm. Some slices (Figures 1(b) and 1(c)) placed at the middle section of the calf were selected to form the MR volumes. As a result, initial and deformed volumes with size of 512 × 512 × 57 voxels were generated for the experiment.

Figure 1

Acquired MR volumes. (a) Place of acquired volume; (b) MR slices in volume obtained at natural state; (c) MR slices in volume obtained at deformed state.

First, 500 and 800 features were extracted from the volume acquired on the natural and deformed states, respectively. Next, the proposed CMI-based feature match approach was applied on the two feature-point sets to obtain the optimal PMS. As the result, a PMS with 245 potential match pairs was obtained. The sparse and density deformation fields were computed using the method mentioned in Section 4.3. Figure 2 shows 50000 internal density deformation fields, with large deformation at the bottom of the calf. This result is consistent with the actual situation.

Figure 2

Density deformation fields. Deformation fields generated using PMS obtained using the (a) CMI-based feature match algorithm and the (b) RPFM algorithm).

To prove the validity of the proposed CMI-based feature match algorithm, we compared it with a robust point feature matching (RPFM) algorithm proposed by Chen [23]. In the present study, we applied the RPFM algorithm to the same feature-point sets, which resulted in a PMS with 316 potential match pairs.

We selected 12 landmarks in the slice (z = 40) of deformed MR volume to test the accuracy of the measured internal deformation fields, as shown in the middle picture of Figures 3 and 4. Then, the landmarks were subjected to reverse moving using the internal deformation fields measured through the CMI-based algorithm and RPFM algorithm. The results on the MR volume acquired at natural state were projected to check the accuracy of the deformation fields. Figures 3 and 4 show the reverse moving results of the landmarks.

Figure 3

Reverse moving result of the landmarks using deformation fields measured through the CMI-based approach.

Figure 4

The reverse moving result of the landmarks using deformation fields measured through the RPFM approach.

In Figures 3 and 4, the center of each red rectangle in the middle picture (z = 40) gives the landmark position. Slices that lie on the left and right sides (the middle layer) give the reverse moving result of the landmarks and the z value of different slices, respectively. The outer layer is the zoom in for the reverse moving result of each landmark. In the middle and outer layers, the red rectangles represent the reverse moving position of the landmarks, the green rectangles are actual position of landmarks, and the yellow rectangles represent the reverse moving positions and actual position consistency. From Figures 3 and 4, we note the accuracy of the reverse moving position of landmarks using deformation fields calculated by PMS obtained using CMI-based approach obviously is higher than that of RPFM, that is, the reverse moving position of landmarks 0, 1, 2, 3, and 5. Table 1 shows the quantitative accuracy of the reverse moving results of the landmarks using internal deformation fields obtained by PMS via CMI and RPFM.

Table 1

Comparison of the accuracy in different directions.

Approach	Error	Number of landmarks												RMSE
		0#	1#	2#	3#	4#	5#	6#	7#	8#	9#	10#	11#
	x-error	− 12	− 8	− 7	− 13	0	− 4	− 6	0	0	− 7	0	0	6.626965
CMI	y-error	8	− 6	7	− 7	0	1	6	0	0	5	0	0	4.654747
	z-error	0	− 2	2	2	− 1	0	− 1	3	− 1	2	3	2	1.848423
	x-error	− 21	− 23	− 9	− 12	0	− 11	0	0	− 1	6	0	0	10.61838
RPFM	y-error	20	− 8	26	− 10	0	− 13	0	0	− 7	4	0	0	11.08302
	z-error	2	2	2	− 2	1	3	0	0	3	4	3	1	2.254625

RMSE: root mean square error; x-error: error in the x-direction; y-error: error in the y-direction; z-error: z-error in the z-direction.

As shown in Table 1, regardless of the direction (i.e., x-, y-, and z-directions), the accuracy of the deformation fields measured through PMS obtained using the CMI-based approach is better than that using the RPFM algorithm.

The number of potential matches in optimal PMS obtained using the CMI-based feature matching algorithm is fewer than that of RPFM because the CMI-based approach is combined with the multifeatures in feature matching, whereas RPFM is a single-feature approach. In other words, the match requirements of CMI are stricter compared with those of RPFM. The reliability of optimal PMS obtained using the CMI-based algorithm is higher because it has more accurate deformation fields than the RPFM algorithm. This conclusion is supported by the reverse moving results of the landmarks.

6. Conclusions

In this work, a new method called CMI is presented for the integration of feature-based internal deformation field measurements. In general, feature match algorithms using a single property are highly accurate in specific aspects. However, the overall accuracy is limited because the full advantages of different properties in feature-point matching are not fully used. Fusion multialgorithms offer the use of advantages in algorithms to improve accuracy. Such a fusion is necessary for feature matching in non-rigid objects, where the improvement will be more obvious. In addition, the most advantage of the proposed approach is to provide a feasible option to integrate various feature matching algorithms. Each feature matching algorithm can act as the component of the CMI, and if the appropriate weight can be assigned to the component, then, one can obtain more reliable potential matches. Obviously, the effect of the component weight should also be considered. Thus, (1) investigating an approach to determine the appropriate weights should be the focus of future research; (2) the imaging mechanism of MRI should be further considered in component of the CMI to remove the aberrance of machine to improve the accuracy of feature-point matching as possible.