- Email: [email protected]
Matching point features under small nonrigid motion
Pattern Recognition 34 (2001) 2353}2365
Matching point features under small nonrigid motion Senthil Kumar *, Maha Sallam, Dmitry Goldgof Bell Laboratories, 101 Crawfords Corner Road, Room 2K-523, Holmdel, NJ 07733, USA Intelligent Systems M.D., Tampa, FL, USA University of South Florida, Tampa, FL, USA Received 21 October 1999; accepted 19 October 2000
Abstract This paper describes a method for matching point features between images of objects that have undergone small nonrigid motion. Feature points are assumed to be available and, given a properly extracted set of feature points, a robust matching is established under the condition that the local nonrigid motion of each point is restricted to a circle of radius , where is not too large. This is in contrast to other techniques for point matching which assume either rigid motion or nonrigid motion of a known kind. The point matching problem is viewed in terms of weighted bipartite graph matching. In order to account for the possibility that the feature selector can be imprecise, we incorporate a greedy matching strategy with the weighted graph matching algorithm. Our algorithm is robust and insensitive to noise and missing features. The resulting matching can be used with image warping or other techniques for nonrigid motion analysis, image subtraction, etc. We present our experimental results on sequences of mammograms, images of a deformable clay object and satellite cloud images. In the "rst two cases we provide quantitative comparison with known ground truth. 2001 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. Keywords: Point correspondence problem; Nonrigid motion analysis; Graph matching
1. Introduction One of the "rst steps in many image registration and nonrigid motion analysis algorithms is the estimation of point correspondences between two images. The problem is to extract a set of n points from each image and to match each point in the "rst image to its corresponding point in the second. There are two primary issues that have to be addressed in solving the point correspondence problem in nonrigid motion: feature selection and feature matching. The "rst issue, namely, feature selection involves extracting a valid set of features that appear in both images. Because of the large variety of nonrigid motions that occur in the real world, feature extraction is very much a problem-dependent issue. For instance, * Corresponding author. Tel.: #1-732-949-9376; fax: #1732-949-0399. E-mail addresses: [email protected] (S. Kumar), sallam@ ismd.com (M. Sallam), [email protected] (D. Goldgof).
a texture-based feature extractor may be suitable for digital mammograms but it would be totally unsuitable for the images of a human face. Also, in the presence of occlusions and large 3D motions it becomes very di$cult to extract the same set of features from both images. The second issue, namely, feature matching is relatively easy in the case of rigid motion where feature points satisfy certain positional and angular constraints. In the case of nonrigid motion, such constraints are not satis"ed and hence feature matching becomes very di$cult even when reliable features are available. Selecting and matching point correspondences in nonrigid motion is a very di$cult problem and it is unlikely that any single algorithm will provide a solution general enough to address all applications. In this paper, we address the issue of feature matching. We assume that a valid set of feature points is available and we provide a technique for matching those feature points under certain constraints. Extensive research has gone into solving this problem for the case of rigid motion [1}5] and recently some work has been done for the case of nonrigid motion. For
0031-3203/01/$20.00 2001 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 1 - 3 2 0 3 ( 0 0 ) 0 0 1 6 6 - 7
2354
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
instance, in Ref. [6], the authors determine point matches by assuming that the surface under study is undergoing conformal motion and, in Ref. [7], they estimate matches by curvature analysis under the assumption that the local nonrigid motion is very small. In Ref. [8], the authors establish correspondences in nonrigid motion using the `smoothness of motion assumptiona. In this technique, which requires at least three consecutive time frames, it is assumed that the motion is smooth and that there are no drastic changes between two time frames. Such assumptions are not always satis"ed. In digital mammography, for example, one normally compares just two mammograms. Also, since the mammograms are taken on the order of one year apart, the smoothness of motion assumption does not hold. In this paper, we present an algorithm that matches features under the constraint that the local nonrigid motion of each point is restricted to a circle of radius , where need not be very small. We have used this method successfully in our work on digital mammograms and it has other applications as well. The application domain consists of those problems where the issue of feature extraction can be addressed reliably. Feature extraction is very much an application dependent issue and we address it brie#y for the special case of digital mammograms. We designed our algorithm for addressing the matching problem in digital mammography and the results (as measured both visually and in terms of the goodness of a warping that is based on the match) are very encouraging. We feel that, the matching strategy is applicable to other problems as well, provided a good feature extractor is available. Our matching scheme is based on the fact that corresponding points in two images would be close to each other in some N-dimensional (N-D) feature space. Let us denote the set of feature points derived from the "rst image by P"p , p ,2, p and that from the second L image by Q"q , q ,2, q . In general, m need not be K equal to n. To each of these points we can assign a feature value or a feature vector (the feature measures used depend on the application). If p and q are corresponding G G points, then the distance between p and q would be very G G small in our N-D feature space. Loosely, a matching can be viewed as an optimal pairing of points in P to those in Q such that the sum of the squared distances between the paired points is a minimum compared with all possible pairings. We solve the correspondence problem by using results from graph theory. A number of researchers have used graph-theoretic concepts in solving the point correspondence problem for rigid motion [1,2,5]. The basic idea behind all these techniques is subgraph isomorphism. Some of these techniques "rst build two di!erent graphs, one for each image, and look for subgraph isomorphism between the two graphs. Others construct separate graphs for the two images while maintaining
graph isomorphism. While these techniques work well for rigid motion, they are not suitable for nonrigid motion because distances between vertices and angles between edges are not preserved in nonrigid motion. In fact in some cases, such as X-ray images and images of clouds, even the relative positioning of feature points may not be preserved throughout an image sequence. We view our problem in terms of weighted bipartite graph matching in which a single graph is used to represent the feature points derived from both images. Kim and Kak [3] have used unweighted bipartite graph matching along with discrete relaxation to match scenes of rigid objects to those stored in a model library. In their technique, bipartite matching is used to prune the search space. The use of weighted bipartite matching would prune their search space even further. A weighted bipartite graph is a weighted graph whose vertex set can be divided into two disjoint sets ; and < such that all edges run between ; and <. A matching in this graph is a maximal subset M of edges such that each vertex of the graph is incident upon at most one edge of M. The correspondence between bipartite graph matching and the point matching is established in Section 2. A greedy algorithm for solving the correspondence problem is discussed in Section 3. This algorithm tries to minimize an error measure in a greedy fashion but it does not guarantee that an overall minimum error is always reached. A well-known polynomial-time algorithm, called the Hungarian Method, for solving the weighted matching problem is discussed in Section 4. The exact problem formulation in terms of bipartite graph matching is not necessarily appropriate in a realworld situation where feature points could be missing and feature values could be erroneous. We address this issue in Section 5, where we combine the greedy algorithm and the Hungarian method to produce a hybrid method for extracting a reliable set of point correspondences between the two images. We feel that the hybrid method is very useful in practical situations. Problems that can arise due to multiple minima are addressed in Section 6. Selection of feature points in digital mammograms is addressed in Section 7. Section 8 discusses the applications of this technique to the analysis of mammograms, a deformable clay "gure and cloud sequences along with experimental results. Quantitative comparison with known ground truth is presented in the "rst two cases. Section 9 provides some conclusions.
2. Point matching and weighted bipartite graphs A weighted bipartite graph is a weighted graph whose vertex set can be divided into two disjoint sets ; and < such that all edges run between ; and <. Formally, such a graph is denoted by G"(;, <, E, =) where ; and < are disjoint sets of vertices, E-;;< is a set of edges
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
2355
Fig. 1. An example of a weighted bipartite graph and a minimum weight matching. (a) Bipartite graph. (b) A minimum weight matching.
and =:EPR is a function that associates a real number called a weight to each edge. For simplicity, let us assume that ;"<"n and that E";;< (i.e., every vertex in ; is connected to every vertex in <). Given any subset A of edges, the weight of A is de"ned to be the sum of the weights of the edges in A. We de"ne a matching to be a set M of edges such that (1) M"n and (2) each vertex in the graph is incident on exactly one edge in M. If an edge (u, v) belongs to M then vertex u is said to be matched to vertex v and vice versa. A matching with a minimum weight is said to be a minimum weight matching or simply a minimum matching. Fig. 1 shows a weighted bipartite graph with n"3 and a minimum weight matching of weight 11. There exists an elegant algorithm called the Hungarian method for "nding a minimum weight matching in a weighted bipartite graph. This algorithm has a time complexity of O(n). We give an overview of this algorithm in Section 4. We now map the point matching problem to the minimum weight graph matching problem. Suppose we are given a set of feature points from each image. Each feature point has an associated feature value or a vector of feature values (intensity, texture density, edge strength/direction, etc.) such that corresponding points in the two images would be close to each other in the feature space. Point matching could then be viewed as the pairing of points from the two images such that the sum of the distances between respective pairs of feature vectors is a minimum. This is, of course, under the assumption that non-corresponding points would have dissimilar feature values. When the feature points are distributed over the entire image, accidental similarities can occur between noncorresponding points. To reduce the error that might arise due to such accidental similarities, we impose a locality constraint: given a point p in one image, the search for its corresponding point q is
conducted only within a neighborhood around the coordinates of p. (We have assumed that global translation and rotation, if present, have been approximately accounted for.) The size of the neighborhood is determined based on our a priori knowledge about the amount of nonrigid motion and the errors in our estimates of the global motion. Let us denote the set of feature points derived from the "rst image by P"p , p ,2, p and that from the sec J ond image by Q"q , q ,2, q . We now construct K a bipartite graph G"(;, <, E, =) where the nodes in ; correspond to points in P and the nodes in < correspond to points in Q (here after, we will treat nodes and points synonymously). In order to simplify analysis, we set ;"< by adding a suitable number of dummy nodes to the smaller set. There is an edge between every node in ; and every node in <. The weight of an edge (u, v) between nondummy nodes u3; and v3< is the squared distance between the corresponding feature vectors. If v does not lie in the -neighborhood of u, then the weight of the edge (u, v) is set to a very large number, BIG. If either v or u is a dummy node then the weight is set to BIG. With the above construction, the point matching problem can be viewed as a problem of determining a minimum weight matching in G: given two sets of feature points, we construct a bipartite graph as mentioned above and use the Hungarian method to determine a minimum weight matching M. If an edge (u, v) belongs to M and if the weight of the edge is not BIG then u and v are corresponding points. The above formulation gives a matching of maximum cardinality. Such a formulation is not always the desired formulation due to the fact that the feature selector might be inaccurate. In any practical situation, the second image will be missing some features detected in the "rst and it might have new features that were not detected in the
2356
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
"rst. Moreover, the feature values themselves will have an uncertainty associated with them. We will discuss this in detail in Section 5, where we combine a greedy algorithm with the Hungarian method to determine a more robust matching. The greedy algorithm and the Hungarian method are discussed in the following two sections.
a with d. This match has a weight of 16 as opposed the other possible match whose weight is 3. Despite the fact that the greedy algorithm does not always provide a minimum weight matching, it is useful in combination with the Hungarian method. The Hungarian method is discussed in the next section.
3. A greedy algorithm for point matching
4. The Hungarian method * an overview
In this section, we present a greedy algorithm for solving the matching problem in a weighted bipartite graph G(;, <, E, =). For each node in ; this algorithm picks up a match in a greedy fashion. The "nal matching provided by this algorithm is not guaranteed to be of minimum weight. However, in conjunction with the Hungarian method, it is quite useful for extracting a robust set of matches. The greedy method by itself was "rst proposed and used by us to match feature points in mammograms [9]. It works as follows. Given the graph G, scan the set of edges and "nd an edge (u, v) with the smallest weight. Match u with v and remove the edges incident on u and those incident on v. Repeat the above step until every node in ; is matched to some node in <. The algorithm is given below.
The Hungarian method is a well-known algorithm for "nding a minimum weight matching in a bipartite graph. In this section, we present an overview of the algorithm. For a detailed discussion, the reader is referred to Refs. [10}12]. Let us denote the nodes in ; by u , u ,2, u and those L in < by v , v ,2, v . De"ne the variables x i, j3 L GH 1, 2,2, n such that x "1 if u is matched with v GH G H and 0 otherwise. Let w denote the weight of the edge GH (u , v ). Then the problem of "nding a minimum weight G H matching can be posed as the following linear programming problem:
Greedy Algorithm Input: A bipartite graph G(;, <, E, =) Output: A matching M 1. MQ 2. Repeat ; times (a) Find an edge (u, v)3E with the smallest weight. (b) MQM (u, v) (c) Remove from E the edges (u, k) and (k, v) k"1, 2,2, ; if they have not been removed already. In the example of Fig. 1(a), the greedy algorithm will "rst select the edge (b, d), followed by the edge (c, f ) and "nally (a, e). This is, in fact, a minimum matching. The greedy algorithm does not always produce a matching of minimum weight. For example, for the graph shown in Fig. 2 the greedy algorithm would match b with c and
Fig. 2. The greedy algorithm would match b with c and a with d which has a weight of 16 rather than the other possible match which has a weight of 3.
L L Minimize w x GH GH G H subject to the constraints
(1)
L x "1, i"1, 2,2, n, (2) GH H L x "1, j"1, 2,2, n, (3) GH G x "0 or 1, i, j"1, 2,2, n. (4) GH Eq. (2) imposes the constraint that each node u 3; is G matched to exactly one node in < and Eq. (3) imposes the constraint that each node v 3<, is matched to exactly G one node in ;. The above linear program can be solved e$ciently by the Hungarian method. The algorithm is given below. Hungarian Method [12] 1. Form an n;n matrix W whose (i, j)th entry is the weight w . GH 2. For each row in W, subtract the smallest entry in the row from all entries in the row. 3. For each column in W, subtract the smallest entry in the column from all entries in the column. (In steps 2 and 3, the smallest entry can be a zero in which case no subtraction is necessary.) 4. (a) Identify a row or column in W with the smallest number of zeros (ignoring rows and columns that were already identi"ed). Ties are settled arbitrarily. Pick a zero in the row(column) selected. If there are two or more zeros, pick one arbitrarily. If what we identi"ed is a row, draw a vertical line through this
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
zero. If what we identi"ed is a column, draw a horizontal line through this zero. (b) Repeat the step 4(a) above, ignoring zeros in rows and columns that were identi"ed earlier and zeros with lines through them, until either each zero has exactly one line through it or n lines have been drawn. (c) If the number of lines drawn is n, then we are done. The zeros identi"ed for the drawing of lines constitute a solution. In other words, if a line was drawn through a zero at row i and column j, then match u with v . If the number of lines drawn is less than G H n, go to step 5. 5. (a) Find the entries with no lines through them lying in unidenti"ed rows/columns and "nd their minimum. Subtract this minimum from all these entries. (b) Add the number subtracted in step 5(a) above to each entry in W which has two lines (both horizontal and vertical) through it. (c) Remove all lines and go to step 4. The above algorithm has a time complexity of O(n) and is guaranteed to "nd a minimum weight matching. There can be more than one minimum weight matching and exactly which one of these is found depends on the rule used for breaking ties in step 4(a) of the algorithm. In the next section, we take a careful look at our formulation of the point matching problem and propose a combination of the greedy and Hungarian methods that provides a reliable solution to the problem.
5. A hybrid method We formulated the point matching problem in terms of minimum weight matching in a bipartite graph. In a real world situation, this is not always appropriate. To appreciate this statement, consider the hypothetical case shown in Fig. 3. In this "gure, the positions of feature
Fig. 3. Feature points from the two images: a & b are from the "rst image; c & d are from the second image. Feature values are given in parentheses.
2357
points extracted from two corresponding images are shown. The two points from the "rst image are denoted by a and b while those from the second image are denoted by c and d. The numbers shown in parentheses are the feature values associated with these points. The associated bipartite graph and the minimum matching are shown in Fig. 4. The feature values are shown in parentheses next to the nodes. The weight of the edge joining two nodes is the square of di!erence between the corresponding feature values. The minimum matching has matched a with c and b with d. This matching has a weight of 8 compared with a weight of 16 for the other possible matching. In this case, it is not clear if the minimum matching is, in fact, the correct matching. One might argue that since b and c both have the same feature value of 100, we should match them together and completely ignore a and d. The reasoning would be that it is possible that a's actual corresponding point was not extracted from the second image and d's correct corresponding point was not extracted from the "rst image. The problem here is that, by de"nition, a minimum matching matches all the nodes of the bipartite graph and this is not always necessary or correct. Due to the imprecise nature of the feature selector and image noise, some feature points in one image may not have a corresponding point in the next image. In the above example, the feature values form such a tight cluster in the feature space that it is too close to call one way or the other and it is best to ignore all 4 points and instead concentrate on other feature points that can be derived from the images. Here we are following the philosophy that it is better to have a smaller number of matches than to have wrong matches. So, what we need is a matching between a subset ; of ; and a subset < of < such that the matching is of minimum weight in relation to other possible matchings between ; and < and the matching ignores nodes like the ones discussed above. It is di$cult to formulate the above observation mathematically. At "rst thought, it looks as though we should ask for a matching of any cardinality such that the average weight per edge is a minimum. This is not the correct formulation as can be seen from Fig. 5. In this "gure, we have shown a bipartite graph and two matches. The minimum weight matching of cardinality 1 matches a with c with an average weight of 0. The minimum weight matching of cardinality 2 matches a with c and b with d with an average weight of 0.5. Even though the second matching has a higher average weight, it is, in fact, a good matching. Due to the above reasons, we keep the problem formulation rather intuitive. We need a matching between subsets of the feature points such that the matching is of minimum weight as far as these subsets are concerned. We would like the subsets to be as large as possible while ignoring tight clusters of points in the feature space that
2358
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
Fig. 4. A weighted bipartite graph for matching points given in Fig. 3 and a minimum weight matching for these points. (a) Weighted bipartite graph. (b) A minimum weight matching.
Fig. 5. Smaller average weight may not necessarily indicate a better matching. A matching of cardinality 1 of point a with point c has an average weight of 0 while a matching of cardinality 2 of a with c, and b with d has an average weight of 0.5. (a) A bipartite graph. (b) A matching.
could lead to confusion. Such a matching strategy would be highly robust and insensitive to noise and errors in feature measurements. Where there is a potential con#ict, we ignore feature points completely and concentrate on other feature points, thereby providing a reliable match. Towards this end, we combine the Hungarian method with a greedy algorithm. Referring back to Fig. 4 again, we note that the greedy algorithm would give the matching (b, c),(a, d) where the matching of a with d is incorrect. Recall that the Hungarian method provided the mapping (a, c),(b, d). The results of the two methods are in con#ict because of the tight cluster formed by these points in the feature space. The greedy method always matches a point u3; to that unmatched point v3< that is closest to u in the feature space. On the other hand, the Hungarian method chooses the matches more carefully so as to minimize the total weight. Their results will be in con#ict at tight clusters in the feature space. Feature points forming tight clusters in the feature space can be avoided by a simple combination of the greedy method and the Hungarian method. In this combined algorithm, we "rst apply the Hungarian method to get a matching M and then apply the greedy algorithm to get a matching M and form their intersection
M"M M . This matching M combines the best of both worlds: (1) it is a minimum weight matching between the matched points because it is part of the Hungarian matching (proved below), (2) each point u3; is matched to an unmatched point v3< that is closest to u in the feature space (because of the greedy matching) and (3) feature points forming tight clusters in the feature space are ignored whenever the greedy algorithm con#icts with the Hungarian method. Note that we are not guaranteeing that all such clusters will be removed. We now prove that the matching M does, in fact, provide a minimum weight matching between the matched points. De5nition. Given a bipartite graph G(;, <, E, =), its restriction to a node set S-; < is the bipartite graph G(;, <, E, =) such that ;"; S, <"< S, E" ;;< and =: EPR such that =(u, v)"=(u, v). Given a bipartite graph G and a subset of nodes S, to construct the restriction of G to S, remove from G any node that is not in S along with any edge that is incident on it. Lemma. Let M be a minimum weight matching for a given bipartite graph G(;, <, E, =) and M-M. Let S be the set of nodes that have an edge in M incident on them. Then M is a minimum weight matching for G(;, <, E, =) where G is the restriction of G to S. Proof. Let x "1 if (u , v )3M and 0 otherwise. Then the GH G H total weight of M is L L ¹" x w . (5) GH GH G H By the construction of the Hungarian method, ¹ is a minimum. Let I"1, 2,2, n. Denote the nodes in ; by u and G the nodes in < by v , i3I. Also, let A"i3I u 3S, G G B"j3I v 3S, AM "IA and BM "IB. Then ¹ can be H
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
written as
and
¹" x w # x w # x w GH GH GH GH GH GH GZ HZ GZ M HZ GZ HZ M
¹ (¹.
# x w . (6) GH GH GZ M HZ M Now consider the restriction G of G to S and the corresponding matching M. If the numbering of nodes in G is the same as in G, then the weight of the match given by M is ¹" x w . (7) GH GH GZ HZ Suppose there exists another matching M whose weight is less than M. De"ne y by y "1 if (u , v )3M and GH GH G H 0 otherwise. Then the weight of M is ¹ " y w GH GH GZ HZ
(8)
2359
(9)
De"ne z by GH y if i3A and j3B, (10,11) z " GH GH x otherwise. GH Then, the weight ¹ " L L z w is less than G H GH GH ¹ (this follows from Eqs. (6), (8) and (9)) and the z 's GH satisfy the constraints given by Eqs. (2), (3) and (4). This implies the existence of a matching M for G whose weight is less than M, a contradiction. Hence, M is of minimum weight in G. 䊐
An example of using the combined algorithm is given in Fig. 6. Fig. 6(a) shows a bipartite graph. The numbers in parentheses denote the feature values. The weights of an edge between two nodes is the square of the di!erence between the corresponding feature values. Edges with
Fig. 6. Using a hybrid algorithm which combines the Hungarian and greedy methods produces a smaller set of more reliable matches. (a) A bipartite graph. (b) Matching produced by hungarian method. (c) Matching produced by greedy algorithm. (d) Combined matching.
2360
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
a weight of BIG are not shown. Fig. 6(b) shows the matching M produced by the Hungarian method. Fig. 6(c) shows the matching M produced by the greedy algorithm. The combined matching M"M M is shown in Fig. 6(d). The hybrid technique has produced matches that are reliable and ignored others. The hybrid algorithm is very useful in practical situations. Instead of trying to "nd a match which matches all the feature points, the hybrid algorithm gives a match between subsets of points. This match is more reliable than those given by the greedy algorithm or the Hungarian method. We tested this method both quantitatively (by comparing with known ground truth) and qualitatively on di!erent types of images. We will describe these results later.
6. Multiple minima As mentioned before, there can be more than one minimum weight matching in a weighted bipartite graph. Given a bipartite graph, theoretically, we can assign weights to the edges in an in"nite number of ways so as to create more than one minimum matching. For example, if we assign equal weights to all the edges, then there are n! matches with the same minimum weight. In our case, however, weights are based on feature values at the nodes. Multiple minima can still occur but, under rather restrictive conditions. Given any vector >"[y , y ,2, y ], a permutation L of > is another vector >N which is obtained from > by an arbitrary rearrangement of its elements. A matching between node sets ; and < can be seen as a permutation denote the vector of feature values of <. Suppose weights are assigned such that the weight of each edge is equal to the squared diwerence between the feature values at nodes incident on the edge. Let > be a permutation of > such that the corresponding nodes of the corresponding elements of X and > are matched by a minimum matching M . Then there exists another minimum weight matching M if and only if there exists another permutation > of > such that X2> "X2> . Proof. The weight of M is w " L (x !y ). G G G By our hypothesis w is a minimum. The weight of M corresponding to the permutation > is w " L (x !y ). M would also be a minimum G G G matching if w "w . i.e., L L (x !y )" (x !y ) (12) G G G G G G
L L N x y " x y G G G G G G NX2> "X2> .
(13) 䊐
(14)
The above lemma states the conditions under which multiple minima can occur and it is comforting to know that such minima can occur only under restrictive conditions. They can still occur and at this time we do not have any suggestions for avoiding them other than to recommend a careful selection of feature values or feature vectors. The one and only constraint that we imposed on the motion, namely, the restriction of local deformation to a circle of radius , was imposed to avoid multiple minima that can occur due to similarities between features that lie in di!erent regions of the image. A careful selection of feature vectors along with a good value for will reduce the possibility of multiple minima. For example, if we are tracking a moving line segment of constant intensity, a bad choice for the feature value would be intensity while a better choice would be distance from one end of the segment. In some applications, it is possible to use relative positioning of points to validate our matches and remove subsets of matches that fail our validation criteria. But this amounts to assuming more restrictions on the motion. Our purpose here is to present a general strategy onto which application dependent ideas can be added. (For one such application, please refer to our work on mammograms [9].)
7. Feature extraction The selection of feature values or feature vectors depends on the type of application at hand. The key is to select points which are distinct from their local surroundings and which remain distinct throughout the image sequence. If the application deals with images which have clearly de"ned edges or curves, then corners and points of maximum curvature along the curves might be useful. Other applications may require matching points in textured images, in which case the use of interest operators (such as those that extract points of maximum local variance) may be appropriate. In this section we describe the feature extractor that we used on digital mammograms. Our interest points are based on the Moravec interest operator [13,14]. This operator selects points which have a maximum directional local variance. To select a point, the operator "rst calculates a measure of local variance by calculating the average squared di!erence between adjacent pixels in each one of four di!erent directions (horizontal, vertical, and two diagonal directions) in a window of size s centered around the point. Each point is then assigned a value I which is the minimum of these
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
averages. A point is selected as an interest point if (1) it has the maximum I in a window of size d and (2) I is above a certain threshold ¹ . In our implementaKMP?TCA tion, s is "xed to a value of 7 while d and ¹ are KMP?TCA varied to control the number of points produced to be matched. Each interest point is then assigned a four dimensional feature vector. The feature vector is derived using four of Laws' texture energy measures (see Ref. [15, pp. 467}469]). The masks we use are ¸52E5, ¸52S5, R52S5, and R52R5, using Laws' notation. Essentially, each feature vector captures texture information over a small region surrounding the corresponding feature point. Note that this feature extractor is certainly not the best possible operator because it uses only the 2D information available from the mammogram whereas the mammogram itself is a 2D projection of a 3D object. Despite this limitation, we have found this feature extractor to be very e!ective in providing good matches as judged by the quality of image registration that is based on these matches. More details are presented in the experimental section. Regardless of what type of images we consider, moving towards satisfying the following ideal conditions will improve the "nal matching: 1. Every extracted feature point in one image should have its corresponding point extracted from the second. This ideal condition would not be satis"ed in most applications. However, if we extract a large number of points from both images, we should increase the number of corresponding pairs that are extracted. 2. The points extracted are well distributed throughout the feature space so that no two points have the same features. This will decrease the amount of `confusiona caused by close clusters of points and hence will improve the performance of the matching algorithm. The point distribution in feature space is dependent on both the points we select and the types of features we associate with them. 3. The points and their features are selected such that the values of the features do not change as a result of changes in the imaging process and local/global motion of objects from one image to the next. Perfect satisfaction of the above conditions is not possible in most cases. Our hybrid method is meant to overcome much of the di$culty introduced by the gap between the ideal and the practical.
8. Experimental results In this section we discuss the results of applying the hybrid algorithm to mammogram images, images of a deformable clay "gure and clouds in a hurricane formation.
2361
Although each application deserves a careful selection of feature extractor, we used the same feature extractor (described in the previous section) for all these applications. A variable parameter in our implementation is the size of the window within which we allow each point in one image to be matched in the second image. The size of this window is determined by the amount of local nonrigid motion possible and the amount of error in the global positioning. There is a trade-o! in selecting a smaller size window to increase the robustness of the matching process and selecting a larger size window to ensure that the actual match does indeed fall within the window. 8.1. Mammogram pair In the "rst experiment we use a pair of mammograms, one of which was obtained from the other by using a synthetic warping technique. Since the warping is synthetic, the exact point correspondences are known and it is possible to quantitatively evaluate our technique. This experiment provides an excellent quantitative test to measure the accuracy of our algorithm. We compare the matching produced by the Hungarian method and the greedy algorithm, each acting alone, with those produced by the hybrid method. We use the mammogram pair to validate our proposition that using the hybrid method produces a better subset of matches than either the Hungarian or the greedy methods. The warping used simulates a complex transformation of the image which allows for local deformations. The transformation is performed by selecting a set of evenly distributed points along the breast region boundary, displacing the points by varying amounts and then using a thin-plate spline warping algorithm given in Ref. [16] to warp the entire image. Since the warping function is known, we know the exact point correspondences between the two images. We selected interest points from the "rst image and used the known warping to compute their locations in the next image. We then calculated the feature values at these points in both images. These feature points and values were input to the three methods and the resulting matches were compared with the actual ones. Fig. 7 shows the set of interest points which was selected in one image and then mapped onto the other (top). It also shows the matching produced by the hybrid method (bottom). Corresponding points are indicated by the same label in both images. Here, each image is of size 388;665 and "141 pixels (approximately a third of the image width). Table 1 shows the total points matched, total incorrect matches, and percentage of incorrect matches for the greedy, Hungarian and hybrid methods. It is clear from this table that while the hybrid method matched a smaller set of points, the percentage of correct matches is signi"cantly higher than that given by the other two methods.
2362
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
Fig. 7. Interest points used (top) and the matching produced by the hybrid method on a mammogram image and a warped version of the same image.
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365 Table 1 Results for the mammogram pair: although the Hybrid method matched a smaller number of points, it produced the least number of incorrect matches and the smallest percentage of incorrect matches Method
Greedy Hungarian Hybrid
Total mactched
Incorrect matches
%Incorrect
139 138 127
11 14 4
7.91 10.14 3.14
Table 2 Results for the mammogram pair with missing feature points Method
Total matched
Incorrect matched
%Incorrect
Greedy Hungarian Hybrid
115 109 105
11 8 4
9.56 7.33 3.81
To test the e!ects of missing points, we repeated the matching process after removing every "fth point in the "rst image while retaining all the feature points in the second image. The results are shown in Table 2. As can be seen, the percentage of correct matches remains almost the same.
2363
Table 3 The hybrid method produces a smaller percentage of mismatched points in the dinosaur images Method
Total matched
Incorrect matches
%Incorrect
Greedy Hungarian Hybrid
227 225 175
59 70 29
25 31 16
obtained gray-scale images of the model under di!erent deformations. To evaluate our matching algorithm, interest points were selected automatically in a pair of images and a human operator was allowed to select the correct match for every interest point in the "rst image by clicking a mouse close to one of several marked potential matching interest points in the second image. The operator was allowed to leave a point in the "rst image unmatched if the corresponding point was missed by the interest point selector in the second image. This established the correct correspondence to which the results of the Hungarian, greedy, and hybrid methods were compared. Fig. 8 shows the image sequence of the model. To avoid cluttering we have labeled only a small subset of the matched points. Table 3 provides a comparison of the three matching methods. Just as in the case of the mammogram images, the hybrid method produces a smaller percentage of incorrect matches than either of the other two methods. Here, the images are of size 320;240 and "81 pixels.
8.2. Deformable clay model 8.3. Cloud images In the second experiment, we constructed a deformable clay model using textured modeling clay on which interest points were marked using small black beads. We then
In this experiment we applied our algorithm to the images of clouds in a hurricane formation. The motion
Fig. 8. The matching produced by the hybrid method on two deformation instances of a clay dinosaur. The second object was obtained from the "rst by slightly bending the tail, the neck and the body. To avoid cluttering, only a small subset of the matched points is shown.
2364
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
Fig. 9. Interest points extracted (top) and the matching produced by the hybrid method on a sequence of hurricane images.
Fig. 10. Corresponding regions from two hurricane images.
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
involved is that of wind-trapped clouds and is composed of complex 3D movements. Since we are working with 2D images and our feature extractor is also two dimensional, we cannot hope to match the actual cloud particles. This experiment is intended to study the e!ectiveness of our algorithm on a time varying 2D textured scene where, two points are considered to be corresponding to each other if their feature vectors (based on the 2D scene) are identical. We want to stress that this experiment was conducted more out of curiosity * a serious analysis of hurricane images would require much more than our simple 2D feature extractor. Fig. 9 shows the interest points selected in two images of hurricane Hugo that were taken 30 min apart (top). The "gure also shows the matches produced by the hybrid method (bottom). Here, the images are of size 512;512 and "101 pixels. Portions of the original images are shown in Fig. 10. 9. Conclusions In this paper we presented a hybrid algorithm for the estimation of point correspondences under small nonrigid motion. Our algorithm presents a general strategy under the assumption that the local nonrigid deformation is limited to a neighborhood of size . Our design philosophy was to extract a reliable set of correspondences from possible inaccurate measurements. As we mentioned in the introduction, the application domain for our algorithm consists of those applications where the issue of feature extraction can be addressed reliably. Feature extraction is very much an application dependent issue the success of the algorithm depends on a careful selection of feature values. While it is possible to envision applications where such a selection is not possible, we believe that our algorithm is applicable to a wide range of problems. References [1] H.H. Chen, T.S. Huang, Maximal matching of 3D points for multiple-object motion estimation, Pattern Recognition 214 (2) (1988) 75}90.
2365
[2] A. Goshtasby, G.C. Stockman, Point pattern matching using convex hull edges, IEEE Trans. Systems Man Cybernet. 15 (5) (1985) 631}637. [3] W. Kim, A.C. Kak, 3-D object recognition using bipartite matching embedded in discrete relaxation, IEEE Trans. Pattern Recognition Mach. Intell. 13 (3) (1991) 224}251. [4] K.S. Arun, T.S. Huang, S.D. Blostein, Least square "tting of two point sets, IEEE Trans. Pattern Recognition Mach. Intell. 9 (5) (1987) 698}700. [5] S. Umeyama, Parametrized point pattern matching and its applications to recognition of object families, IEEE Trans. Pattern Recognition Mach. Intell. 15 (2) (1993) 136}144. [6] C. Kambhamettu, D. Goldgof, Point correspondence recovery in nonrigid motion, Proceedings of the International Conference on Computer Vision and Pattern Recognition, 1992, pp. 222}227. [7] C. Kambhamettu, D. Goldgof, Curvature-based approach to point correspondence recovery in nonrigid motion, Proceedings of the International Conference on Computer Vision and Pattern Recognition, July 1994, pp. 222}227. [8] R. Jain, I.K. Sethi, Establishing Correspondence of Non-Rigid Objects using Smoothness of Motion, IEEE, New York, 1984, pp. 83}87. [9] M. Sallam, K.W. Bowyer, Registering time sequences of mammograms using a two-dimensional image unwarping technique, Proceedings of the Second International Workshop on Digital Mammography, July 1994, pp. 121}130. [10] C. Papadimitriou, K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity, Prentice-Hall, Englewood Cli!s, NJ, 1982. [11] G. Strang, Introduction to Applied Mathematics, Wellesley Cambridge Press, Wellesley, MA, 1995. [12] L.R. Foulds, Combinatorial Optimization for Undergraduates, Springer, Berlin, 1984. [13] D. Ballard, C. Brown, Computer Vision, Prentice-Hall, Englewood Cli!s, NJ, 1982. [14] H.P. Moravec, Rover visual obstacle avoidance, Proceedings of the International Joint Conference on A.I., 1981, pp. 785}790. [15] R. Haralick, L. Shapiro, Computer and Robot Vision, Addison-Wesley, Reading, MA, 1992. [16] F.L. Bookstein, Principal warps: thin plate splines and the decomposition of deformations, IEEE Trans. Pattern Recognition Mach. Intell. 11 (6) (1989) 567}585.
Matching point features under small nonrigid motion Senthil Kumar *, Maha Sallam, Dmitry Goldgof Bell Laboratories, 101 Crawfords Corner Road, Room 2K-523, Holmdel, NJ 07733, USA Intelligent Systems M.D., Tampa, FL, USA University of South Florida, Tampa, FL, USA Received 21 October 1999; accepted 19 October 2000
Abstract This paper describes a method for matching point features between images of objects that have undergone small nonrigid motion. Feature points are assumed to be available and, given a properly extracted set of feature points, a robust matching is established under the condition that the local nonrigid motion of each point is restricted to a circle of radius , where is not too large. This is in contrast to other techniques for point matching which assume either rigid motion or nonrigid motion of a known kind. The point matching problem is viewed in terms of weighted bipartite graph matching. In order to account for the possibility that the feature selector can be imprecise, we incorporate a greedy matching strategy with the weighted graph matching algorithm. Our algorithm is robust and insensitive to noise and missing features. The resulting matching can be used with image warping or other techniques for nonrigid motion analysis, image subtraction, etc. We present our experimental results on sequences of mammograms, images of a deformable clay object and satellite cloud images. In the "rst two cases we provide quantitative comparison with known ground truth. 2001 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. Keywords: Point correspondence problem; Nonrigid motion analysis; Graph matching
1. Introduction One of the "rst steps in many image registration and nonrigid motion analysis algorithms is the estimation of point correspondences between two images. The problem is to extract a set of n points from each image and to match each point in the "rst image to its corresponding point in the second. There are two primary issues that have to be addressed in solving the point correspondence problem in nonrigid motion: feature selection and feature matching. The "rst issue, namely, feature selection involves extracting a valid set of features that appear in both images. Because of the large variety of nonrigid motions that occur in the real world, feature extraction is very much a problem-dependent issue. For instance, * Corresponding author. Tel.: #1-732-949-9376; fax: #1732-949-0399. E-mail addresses: [email protected] (S. Kumar), sallam@ ismd.com (M. Sallam), [email protected] (D. Goldgof).
a texture-based feature extractor may be suitable for digital mammograms but it would be totally unsuitable for the images of a human face. Also, in the presence of occlusions and large 3D motions it becomes very di$cult to extract the same set of features from both images. The second issue, namely, feature matching is relatively easy in the case of rigid motion where feature points satisfy certain positional and angular constraints. In the case of nonrigid motion, such constraints are not satis"ed and hence feature matching becomes very di$cult even when reliable features are available. Selecting and matching point correspondences in nonrigid motion is a very di$cult problem and it is unlikely that any single algorithm will provide a solution general enough to address all applications. In this paper, we address the issue of feature matching. We assume that a valid set of feature points is available and we provide a technique for matching those feature points under certain constraints. Extensive research has gone into solving this problem for the case of rigid motion [1}5] and recently some work has been done for the case of nonrigid motion. For
0031-3203/01/$20.00 2001 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 1 - 3 2 0 3 ( 0 0 ) 0 0 1 6 6 - 7
2354
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
instance, in Ref. [6], the authors determine point matches by assuming that the surface under study is undergoing conformal motion and, in Ref. [7], they estimate matches by curvature analysis under the assumption that the local nonrigid motion is very small. In Ref. [8], the authors establish correspondences in nonrigid motion using the `smoothness of motion assumptiona. In this technique, which requires at least three consecutive time frames, it is assumed that the motion is smooth and that there are no drastic changes between two time frames. Such assumptions are not always satis"ed. In digital mammography, for example, one normally compares just two mammograms. Also, since the mammograms are taken on the order of one year apart, the smoothness of motion assumption does not hold. In this paper, we present an algorithm that matches features under the constraint that the local nonrigid motion of each point is restricted to a circle of radius , where need not be very small. We have used this method successfully in our work on digital mammograms and it has other applications as well. The application domain consists of those problems where the issue of feature extraction can be addressed reliably. Feature extraction is very much an application dependent issue and we address it brie#y for the special case of digital mammograms. We designed our algorithm for addressing the matching problem in digital mammography and the results (as measured both visually and in terms of the goodness of a warping that is based on the match) are very encouraging. We feel that, the matching strategy is applicable to other problems as well, provided a good feature extractor is available. Our matching scheme is based on the fact that corresponding points in two images would be close to each other in some N-dimensional (N-D) feature space. Let us denote the set of feature points derived from the "rst image by P"p , p ,2, p and that from the second L image by Q"q , q ,2, q . In general, m need not be K equal to n. To each of these points we can assign a feature value or a feature vector (the feature measures used depend on the application). If p and q are corresponding G G points, then the distance between p and q would be very G G small in our N-D feature space. Loosely, a matching can be viewed as an optimal pairing of points in P to those in Q such that the sum of the squared distances between the paired points is a minimum compared with all possible pairings. We solve the correspondence problem by using results from graph theory. A number of researchers have used graph-theoretic concepts in solving the point correspondence problem for rigid motion [1,2,5]. The basic idea behind all these techniques is subgraph isomorphism. Some of these techniques "rst build two di!erent graphs, one for each image, and look for subgraph isomorphism between the two graphs. Others construct separate graphs for the two images while maintaining
graph isomorphism. While these techniques work well for rigid motion, they are not suitable for nonrigid motion because distances between vertices and angles between edges are not preserved in nonrigid motion. In fact in some cases, such as X-ray images and images of clouds, even the relative positioning of feature points may not be preserved throughout an image sequence. We view our problem in terms of weighted bipartite graph matching in which a single graph is used to represent the feature points derived from both images. Kim and Kak [3] have used unweighted bipartite graph matching along with discrete relaxation to match scenes of rigid objects to those stored in a model library. In their technique, bipartite matching is used to prune the search space. The use of weighted bipartite matching would prune their search space even further. A weighted bipartite graph is a weighted graph whose vertex set can be divided into two disjoint sets ; and < such that all edges run between ; and <. A matching in this graph is a maximal subset M of edges such that each vertex of the graph is incident upon at most one edge of M. The correspondence between bipartite graph matching and the point matching is established in Section 2. A greedy algorithm for solving the correspondence problem is discussed in Section 3. This algorithm tries to minimize an error measure in a greedy fashion but it does not guarantee that an overall minimum error is always reached. A well-known polynomial-time algorithm, called the Hungarian Method, for solving the weighted matching problem is discussed in Section 4. The exact problem formulation in terms of bipartite graph matching is not necessarily appropriate in a realworld situation where feature points could be missing and feature values could be erroneous. We address this issue in Section 5, where we combine the greedy algorithm and the Hungarian method to produce a hybrid method for extracting a reliable set of point correspondences between the two images. We feel that the hybrid method is very useful in practical situations. Problems that can arise due to multiple minima are addressed in Section 6. Selection of feature points in digital mammograms is addressed in Section 7. Section 8 discusses the applications of this technique to the analysis of mammograms, a deformable clay "gure and cloud sequences along with experimental results. Quantitative comparison with known ground truth is presented in the "rst two cases. Section 9 provides some conclusions.
2. Point matching and weighted bipartite graphs A weighted bipartite graph is a weighted graph whose vertex set can be divided into two disjoint sets ; and < such that all edges run between ; and <. Formally, such a graph is denoted by G"(;, <, E, =) where ; and < are disjoint sets of vertices, E-;;< is a set of edges
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
2355
Fig. 1. An example of a weighted bipartite graph and a minimum weight matching. (a) Bipartite graph. (b) A minimum weight matching.
and =:EPR is a function that associates a real number called a weight to each edge. For simplicity, let us assume that ;"<"n and that E";;< (i.e., every vertex in ; is connected to every vertex in <). Given any subset A of edges, the weight of A is de"ned to be the sum of the weights of the edges in A. We de"ne a matching to be a set M of edges such that (1) M"n and (2) each vertex in the graph is incident on exactly one edge in M. If an edge (u, v) belongs to M then vertex u is said to be matched to vertex v and vice versa. A matching with a minimum weight is said to be a minimum weight matching or simply a minimum matching. Fig. 1 shows a weighted bipartite graph with n"3 and a minimum weight matching of weight 11. There exists an elegant algorithm called the Hungarian method for "nding a minimum weight matching in a weighted bipartite graph. This algorithm has a time complexity of O(n). We give an overview of this algorithm in Section 4. We now map the point matching problem to the minimum weight graph matching problem. Suppose we are given a set of feature points from each image. Each feature point has an associated feature value or a vector of feature values (intensity, texture density, edge strength/direction, etc.) such that corresponding points in the two images would be close to each other in the feature space. Point matching could then be viewed as the pairing of points from the two images such that the sum of the distances between respective pairs of feature vectors is a minimum. This is, of course, under the assumption that non-corresponding points would have dissimilar feature values. When the feature points are distributed over the entire image, accidental similarities can occur between noncorresponding points. To reduce the error that might arise due to such accidental similarities, we impose a locality constraint: given a point p in one image, the search for its corresponding point q is
conducted only within a neighborhood around the coordinates of p. (We have assumed that global translation and rotation, if present, have been approximately accounted for.) The size of the neighborhood is determined based on our a priori knowledge about the amount of nonrigid motion and the errors in our estimates of the global motion. Let us denote the set of feature points derived from the "rst image by P"p , p ,2, p and that from the sec J ond image by Q"q , q ,2, q . We now construct K a bipartite graph G"(;, <, E, =) where the nodes in ; correspond to points in P and the nodes in < correspond to points in Q (here after, we will treat nodes and points synonymously). In order to simplify analysis, we set ;"< by adding a suitable number of dummy nodes to the smaller set. There is an edge between every node in ; and every node in <. The weight of an edge (u, v) between nondummy nodes u3; and v3< is the squared distance between the corresponding feature vectors. If v does not lie in the -neighborhood of u, then the weight of the edge (u, v) is set to a very large number, BIG. If either v or u is a dummy node then the weight is set to BIG. With the above construction, the point matching problem can be viewed as a problem of determining a minimum weight matching in G: given two sets of feature points, we construct a bipartite graph as mentioned above and use the Hungarian method to determine a minimum weight matching M. If an edge (u, v) belongs to M and if the weight of the edge is not BIG then u and v are corresponding points. The above formulation gives a matching of maximum cardinality. Such a formulation is not always the desired formulation due to the fact that the feature selector might be inaccurate. In any practical situation, the second image will be missing some features detected in the "rst and it might have new features that were not detected in the
2356
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
"rst. Moreover, the feature values themselves will have an uncertainty associated with them. We will discuss this in detail in Section 5, where we combine a greedy algorithm with the Hungarian method to determine a more robust matching. The greedy algorithm and the Hungarian method are discussed in the following two sections.
a with d. This match has a weight of 16 as opposed the other possible match whose weight is 3. Despite the fact that the greedy algorithm does not always provide a minimum weight matching, it is useful in combination with the Hungarian method. The Hungarian method is discussed in the next section.
3. A greedy algorithm for point matching
4. The Hungarian method * an overview
In this section, we present a greedy algorithm for solving the matching problem in a weighted bipartite graph G(;, <, E, =). For each node in ; this algorithm picks up a match in a greedy fashion. The "nal matching provided by this algorithm is not guaranteed to be of minimum weight. However, in conjunction with the Hungarian method, it is quite useful for extracting a robust set of matches. The greedy method by itself was "rst proposed and used by us to match feature points in mammograms [9]. It works as follows. Given the graph G, scan the set of edges and "nd an edge (u, v) with the smallest weight. Match u with v and remove the edges incident on u and those incident on v. Repeat the above step until every node in ; is matched to some node in <. The algorithm is given below.
The Hungarian method is a well-known algorithm for "nding a minimum weight matching in a bipartite graph. In this section, we present an overview of the algorithm. For a detailed discussion, the reader is referred to Refs. [10}12]. Let us denote the nodes in ; by u , u ,2, u and those L in < by v , v ,2, v . De"ne the variables x i, j3 L GH 1, 2,2, n such that x "1 if u is matched with v GH G H and 0 otherwise. Let w denote the weight of the edge GH (u , v ). Then the problem of "nding a minimum weight G H matching can be posed as the following linear programming problem:
Greedy Algorithm Input: A bipartite graph G(;, <, E, =) Output: A matching M 1. MQ 2. Repeat ; times (a) Find an edge (u, v)3E with the smallest weight. (b) MQM (u, v) (c) Remove from E the edges (u, k) and (k, v) k"1, 2,2, ; if they have not been removed already. In the example of Fig. 1(a), the greedy algorithm will "rst select the edge (b, d), followed by the edge (c, f ) and "nally (a, e). This is, in fact, a minimum matching. The greedy algorithm does not always produce a matching of minimum weight. For example, for the graph shown in Fig. 2 the greedy algorithm would match b with c and
Fig. 2. The greedy algorithm would match b with c and a with d which has a weight of 16 rather than the other possible match which has a weight of 3.
L L Minimize w x GH GH G H subject to the constraints
(1)
L x "1, i"1, 2,2, n, (2) GH H L x "1, j"1, 2,2, n, (3) GH G x "0 or 1, i, j"1, 2,2, n. (4) GH Eq. (2) imposes the constraint that each node u 3; is G matched to exactly one node in < and Eq. (3) imposes the constraint that each node v 3<, is matched to exactly G one node in ;. The above linear program can be solved e$ciently by the Hungarian method. The algorithm is given below. Hungarian Method [12] 1. Form an n;n matrix W whose (i, j)th entry is the weight w . GH 2. For each row in W, subtract the smallest entry in the row from all entries in the row. 3. For each column in W, subtract the smallest entry in the column from all entries in the column. (In steps 2 and 3, the smallest entry can be a zero in which case no subtraction is necessary.) 4. (a) Identify a row or column in W with the smallest number of zeros (ignoring rows and columns that were already identi"ed). Ties are settled arbitrarily. Pick a zero in the row(column) selected. If there are two or more zeros, pick one arbitrarily. If what we identi"ed is a row, draw a vertical line through this
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
zero. If what we identi"ed is a column, draw a horizontal line through this zero. (b) Repeat the step 4(a) above, ignoring zeros in rows and columns that were identi"ed earlier and zeros with lines through them, until either each zero has exactly one line through it or n lines have been drawn. (c) If the number of lines drawn is n, then we are done. The zeros identi"ed for the drawing of lines constitute a solution. In other words, if a line was drawn through a zero at row i and column j, then match u with v . If the number of lines drawn is less than G H n, go to step 5. 5. (a) Find the entries with no lines through them lying in unidenti"ed rows/columns and "nd their minimum. Subtract this minimum from all these entries. (b) Add the number subtracted in step 5(a) above to each entry in W which has two lines (both horizontal and vertical) through it. (c) Remove all lines and go to step 4. The above algorithm has a time complexity of O(n) and is guaranteed to "nd a minimum weight matching. There can be more than one minimum weight matching and exactly which one of these is found depends on the rule used for breaking ties in step 4(a) of the algorithm. In the next section, we take a careful look at our formulation of the point matching problem and propose a combination of the greedy and Hungarian methods that provides a reliable solution to the problem.
5. A hybrid method We formulated the point matching problem in terms of minimum weight matching in a bipartite graph. In a real world situation, this is not always appropriate. To appreciate this statement, consider the hypothetical case shown in Fig. 3. In this "gure, the positions of feature
Fig. 3. Feature points from the two images: a & b are from the "rst image; c & d are from the second image. Feature values are given in parentheses.
2357
points extracted from two corresponding images are shown. The two points from the "rst image are denoted by a and b while those from the second image are denoted by c and d. The numbers shown in parentheses are the feature values associated with these points. The associated bipartite graph and the minimum matching are shown in Fig. 4. The feature values are shown in parentheses next to the nodes. The weight of the edge joining two nodes is the square of di!erence between the corresponding feature values. The minimum matching has matched a with c and b with d. This matching has a weight of 8 compared with a weight of 16 for the other possible matching. In this case, it is not clear if the minimum matching is, in fact, the correct matching. One might argue that since b and c both have the same feature value of 100, we should match them together and completely ignore a and d. The reasoning would be that it is possible that a's actual corresponding point was not extracted from the second image and d's correct corresponding point was not extracted from the "rst image. The problem here is that, by de"nition, a minimum matching matches all the nodes of the bipartite graph and this is not always necessary or correct. Due to the imprecise nature of the feature selector and image noise, some feature points in one image may not have a corresponding point in the next image. In the above example, the feature values form such a tight cluster in the feature space that it is too close to call one way or the other and it is best to ignore all 4 points and instead concentrate on other feature points that can be derived from the images. Here we are following the philosophy that it is better to have a smaller number of matches than to have wrong matches. So, what we need is a matching between a subset ; of ; and a subset < of < such that the matching is of minimum weight in relation to other possible matchings between ; and < and the matching ignores nodes like the ones discussed above. It is di$cult to formulate the above observation mathematically. At "rst thought, it looks as though we should ask for a matching of any cardinality such that the average weight per edge is a minimum. This is not the correct formulation as can be seen from Fig. 5. In this "gure, we have shown a bipartite graph and two matches. The minimum weight matching of cardinality 1 matches a with c with an average weight of 0. The minimum weight matching of cardinality 2 matches a with c and b with d with an average weight of 0.5. Even though the second matching has a higher average weight, it is, in fact, a good matching. Due to the above reasons, we keep the problem formulation rather intuitive. We need a matching between subsets of the feature points such that the matching is of minimum weight as far as these subsets are concerned. We would like the subsets to be as large as possible while ignoring tight clusters of points in the feature space that
2358
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
Fig. 4. A weighted bipartite graph for matching points given in Fig. 3 and a minimum weight matching for these points. (a) Weighted bipartite graph. (b) A minimum weight matching.
Fig. 5. Smaller average weight may not necessarily indicate a better matching. A matching of cardinality 1 of point a with point c has an average weight of 0 while a matching of cardinality 2 of a with c, and b with d has an average weight of 0.5. (a) A bipartite graph. (b) A matching.
could lead to confusion. Such a matching strategy would be highly robust and insensitive to noise and errors in feature measurements. Where there is a potential con#ict, we ignore feature points completely and concentrate on other feature points, thereby providing a reliable match. Towards this end, we combine the Hungarian method with a greedy algorithm. Referring back to Fig. 4 again, we note that the greedy algorithm would give the matching (b, c),(a, d) where the matching of a with d is incorrect. Recall that the Hungarian method provided the mapping (a, c),(b, d). The results of the two methods are in con#ict because of the tight cluster formed by these points in the feature space. The greedy method always matches a point u3; to that unmatched point v3< that is closest to u in the feature space. On the other hand, the Hungarian method chooses the matches more carefully so as to minimize the total weight. Their results will be in con#ict at tight clusters in the feature space. Feature points forming tight clusters in the feature space can be avoided by a simple combination of the greedy method and the Hungarian method. In this combined algorithm, we "rst apply the Hungarian method to get a matching M and then apply the greedy algorithm to get a matching M and form their intersection
M"M M . This matching M combines the best of both worlds: (1) it is a minimum weight matching between the matched points because it is part of the Hungarian matching (proved below), (2) each point u3; is matched to an unmatched point v3< that is closest to u in the feature space (because of the greedy matching) and (3) feature points forming tight clusters in the feature space are ignored whenever the greedy algorithm con#icts with the Hungarian method. Note that we are not guaranteeing that all such clusters will be removed. We now prove that the matching M does, in fact, provide a minimum weight matching between the matched points. De5nition. Given a bipartite graph G(;, <, E, =), its restriction to a node set S-; < is the bipartite graph G(;, <, E, =) such that ;"; S, <"< S, E" ;;< and =: EPR such that =(u, v)"=(u, v). Given a bipartite graph G and a subset of nodes S, to construct the restriction of G to S, remove from G any node that is not in S along with any edge that is incident on it. Lemma. Let M be a minimum weight matching for a given bipartite graph G(;, <, E, =) and M-M. Let S be the set of nodes that have an edge in M incident on them. Then M is a minimum weight matching for G(;, <, E, =) where G is the restriction of G to S. Proof. Let x "1 if (u , v )3M and 0 otherwise. Then the GH G H total weight of M is L L ¹" x w . (5) GH GH G H By the construction of the Hungarian method, ¹ is a minimum. Let I"1, 2,2, n. Denote the nodes in ; by u and G the nodes in < by v , i3I. Also, let A"i3I u 3S, G G B"j3I v 3S, AM "IA and BM "IB. Then ¹ can be H
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
written as
and
¹" x w # x w # x w GH GH GH GH GH GH GZ HZ GZ M HZ GZ HZ M
¹ (¹.
# x w . (6) GH GH GZ M HZ M Now consider the restriction G of G to S and the corresponding matching M. If the numbering of nodes in G is the same as in G, then the weight of the match given by M is ¹" x w . (7) GH GH GZ HZ Suppose there exists another matching M whose weight is less than M. De"ne y by y "1 if (u , v )3M and GH GH G H 0 otherwise. Then the weight of M is ¹ " y w GH GH GZ HZ
(8)
2359
(9)
De"ne z by GH y if i3A and j3B, (10,11) z " GH GH x otherwise. GH Then, the weight ¹ " L L z w is less than G H GH GH ¹ (this follows from Eqs. (6), (8) and (9)) and the z 's GH satisfy the constraints given by Eqs. (2), (3) and (4). This implies the existence of a matching M for G whose weight is less than M, a contradiction. Hence, M is of minimum weight in G. 䊐
An example of using the combined algorithm is given in Fig. 6. Fig. 6(a) shows a bipartite graph. The numbers in parentheses denote the feature values. The weights of an edge between two nodes is the square of the di!erence between the corresponding feature values. Edges with
Fig. 6. Using a hybrid algorithm which combines the Hungarian and greedy methods produces a smaller set of more reliable matches. (a) A bipartite graph. (b) Matching produced by hungarian method. (c) Matching produced by greedy algorithm. (d) Combined matching.
2360
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
a weight of BIG are not shown. Fig. 6(b) shows the matching M produced by the Hungarian method. Fig. 6(c) shows the matching M produced by the greedy algorithm. The combined matching M"M M is shown in Fig. 6(d). The hybrid technique has produced matches that are reliable and ignored others. The hybrid algorithm is very useful in practical situations. Instead of trying to "nd a match which matches all the feature points, the hybrid algorithm gives a match between subsets of points. This match is more reliable than those given by the greedy algorithm or the Hungarian method. We tested this method both quantitatively (by comparing with known ground truth) and qualitatively on di!erent types of images. We will describe these results later.
6. Multiple minima As mentioned before, there can be more than one minimum weight matching in a weighted bipartite graph. Given a bipartite graph, theoretically, we can assign weights to the edges in an in"nite number of ways so as to create more than one minimum matching. For example, if we assign equal weights to all the edges, then there are n! matches with the same minimum weight. In our case, however, weights are based on feature values at the nodes. Multiple minima can still occur but, under rather restrictive conditions. Given any vector >"[y , y ,2, y ], a permutation L of > is another vector >N which is obtained from > by an arbitrary rearrangement of its elements. A matching between node sets ; and < can be seen as a permutation
L L N x y " x y G G G G G G NX2> "X2> .
(13) 䊐
(14)
The above lemma states the conditions under which multiple minima can occur and it is comforting to know that such minima can occur only under restrictive conditions. They can still occur and at this time we do not have any suggestions for avoiding them other than to recommend a careful selection of feature values or feature vectors. The one and only constraint that we imposed on the motion, namely, the restriction of local deformation to a circle of radius , was imposed to avoid multiple minima that can occur due to similarities between features that lie in di!erent regions of the image. A careful selection of feature vectors along with a good value for will reduce the possibility of multiple minima. For example, if we are tracking a moving line segment of constant intensity, a bad choice for the feature value would be intensity while a better choice would be distance from one end of the segment. In some applications, it is possible to use relative positioning of points to validate our matches and remove subsets of matches that fail our validation criteria. But this amounts to assuming more restrictions on the motion. Our purpose here is to present a general strategy onto which application dependent ideas can be added. (For one such application, please refer to our work on mammograms [9].)
7. Feature extraction The selection of feature values or feature vectors depends on the type of application at hand. The key is to select points which are distinct from their local surroundings and which remain distinct throughout the image sequence. If the application deals with images which have clearly de"ned edges or curves, then corners and points of maximum curvature along the curves might be useful. Other applications may require matching points in textured images, in which case the use of interest operators (such as those that extract points of maximum local variance) may be appropriate. In this section we describe the feature extractor that we used on digital mammograms. Our interest points are based on the Moravec interest operator [13,14]. This operator selects points which have a maximum directional local variance. To select a point, the operator "rst calculates a measure of local variance by calculating the average squared di!erence between adjacent pixels in each one of four di!erent directions (horizontal, vertical, and two diagonal directions) in a window of size s centered around the point. Each point is then assigned a value I which is the minimum of these
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
averages. A point is selected as an interest point if (1) it has the maximum I in a window of size d and (2) I is above a certain threshold ¹ . In our implementaKMP?TCA tion, s is "xed to a value of 7 while d and ¹ are KMP?TCA varied to control the number of points produced to be matched. Each interest point is then assigned a four dimensional feature vector. The feature vector is derived using four of Laws' texture energy measures (see Ref. [15, pp. 467}469]). The masks we use are ¸52E5, ¸52S5, R52S5, and R52R5, using Laws' notation. Essentially, each feature vector captures texture information over a small region surrounding the corresponding feature point. Note that this feature extractor is certainly not the best possible operator because it uses only the 2D information available from the mammogram whereas the mammogram itself is a 2D projection of a 3D object. Despite this limitation, we have found this feature extractor to be very e!ective in providing good matches as judged by the quality of image registration that is based on these matches. More details are presented in the experimental section. Regardless of what type of images we consider, moving towards satisfying the following ideal conditions will improve the "nal matching: 1. Every extracted feature point in one image should have its corresponding point extracted from the second. This ideal condition would not be satis"ed in most applications. However, if we extract a large number of points from both images, we should increase the number of corresponding pairs that are extracted. 2. The points extracted are well distributed throughout the feature space so that no two points have the same features. This will decrease the amount of `confusiona caused by close clusters of points and hence will improve the performance of the matching algorithm. The point distribution in feature space is dependent on both the points we select and the types of features we associate with them. 3. The points and their features are selected such that the values of the features do not change as a result of changes in the imaging process and local/global motion of objects from one image to the next. Perfect satisfaction of the above conditions is not possible in most cases. Our hybrid method is meant to overcome much of the di$culty introduced by the gap between the ideal and the practical.
8. Experimental results In this section we discuss the results of applying the hybrid algorithm to mammogram images, images of a deformable clay "gure and clouds in a hurricane formation.
2361
Although each application deserves a careful selection of feature extractor, we used the same feature extractor (described in the previous section) for all these applications. A variable parameter in our implementation is the size of the window within which we allow each point in one image to be matched in the second image. The size of this window is determined by the amount of local nonrigid motion possible and the amount of error in the global positioning. There is a trade-o! in selecting a smaller size window to increase the robustness of the matching process and selecting a larger size window to ensure that the actual match does indeed fall within the window. 8.1. Mammogram pair In the "rst experiment we use a pair of mammograms, one of which was obtained from the other by using a synthetic warping technique. Since the warping is synthetic, the exact point correspondences are known and it is possible to quantitatively evaluate our technique. This experiment provides an excellent quantitative test to measure the accuracy of our algorithm. We compare the matching produced by the Hungarian method and the greedy algorithm, each acting alone, with those produced by the hybrid method. We use the mammogram pair to validate our proposition that using the hybrid method produces a better subset of matches than either the Hungarian or the greedy methods. The warping used simulates a complex transformation of the image which allows for local deformations. The transformation is performed by selecting a set of evenly distributed points along the breast region boundary, displacing the points by varying amounts and then using a thin-plate spline warping algorithm given in Ref. [16] to warp the entire image. Since the warping function is known, we know the exact point correspondences between the two images. We selected interest points from the "rst image and used the known warping to compute their locations in the next image. We then calculated the feature values at these points in both images. These feature points and values were input to the three methods and the resulting matches were compared with the actual ones. Fig. 7 shows the set of interest points which was selected in one image and then mapped onto the other (top). It also shows the matching produced by the hybrid method (bottom). Corresponding points are indicated by the same label in both images. Here, each image is of size 388;665 and "141 pixels (approximately a third of the image width). Table 1 shows the total points matched, total incorrect matches, and percentage of incorrect matches for the greedy, Hungarian and hybrid methods. It is clear from this table that while the hybrid method matched a smaller set of points, the percentage of correct matches is signi"cantly higher than that given by the other two methods.
2362
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
Fig. 7. Interest points used (top) and the matching produced by the hybrid method on a mammogram image and a warped version of the same image.
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365 Table 1 Results for the mammogram pair: although the Hybrid method matched a smaller number of points, it produced the least number of incorrect matches and the smallest percentage of incorrect matches Method
Greedy Hungarian Hybrid
Total mactched
Incorrect matches
%Incorrect
139 138 127
11 14 4
7.91 10.14 3.14
Table 2 Results for the mammogram pair with missing feature points Method
Total matched
Incorrect matched
%Incorrect
Greedy Hungarian Hybrid
115 109 105
11 8 4
9.56 7.33 3.81
To test the e!ects of missing points, we repeated the matching process after removing every "fth point in the "rst image while retaining all the feature points in the second image. The results are shown in Table 2. As can be seen, the percentage of correct matches remains almost the same.
2363
Table 3 The hybrid method produces a smaller percentage of mismatched points in the dinosaur images Method
Total matched
Incorrect matches
%Incorrect
Greedy Hungarian Hybrid
227 225 175
59 70 29
25 31 16
obtained gray-scale images of the model under di!erent deformations. To evaluate our matching algorithm, interest points were selected automatically in a pair of images and a human operator was allowed to select the correct match for every interest point in the "rst image by clicking a mouse close to one of several marked potential matching interest points in the second image. The operator was allowed to leave a point in the "rst image unmatched if the corresponding point was missed by the interest point selector in the second image. This established the correct correspondence to which the results of the Hungarian, greedy, and hybrid methods were compared. Fig. 8 shows the image sequence of the model. To avoid cluttering we have labeled only a small subset of the matched points. Table 3 provides a comparison of the three matching methods. Just as in the case of the mammogram images, the hybrid method produces a smaller percentage of incorrect matches than either of the other two methods. Here, the images are of size 320;240 and "81 pixels.
8.2. Deformable clay model 8.3. Cloud images In the second experiment, we constructed a deformable clay model using textured modeling clay on which interest points were marked using small black beads. We then
In this experiment we applied our algorithm to the images of clouds in a hurricane formation. The motion
Fig. 8. The matching produced by the hybrid method on two deformation instances of a clay dinosaur. The second object was obtained from the "rst by slightly bending the tail, the neck and the body. To avoid cluttering, only a small subset of the matched points is shown.
2364
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
Fig. 9. Interest points extracted (top) and the matching produced by the hybrid method on a sequence of hurricane images.
Fig. 10. Corresponding regions from two hurricane images.
S. Kumar et al. / Pattern Recognition 34 (2001) 2353}2365
involved is that of wind-trapped clouds and is composed of complex 3D movements. Since we are working with 2D images and our feature extractor is also two dimensional, we cannot hope to match the actual cloud particles. This experiment is intended to study the e!ectiveness of our algorithm on a time varying 2D textured scene where, two points are considered to be corresponding to each other if their feature vectors (based on the 2D scene) are identical. We want to stress that this experiment was conducted more out of curiosity * a serious analysis of hurricane images would require much more than our simple 2D feature extractor. Fig. 9 shows the interest points selected in two images of hurricane Hugo that were taken 30 min apart (top). The "gure also shows the matches produced by the hybrid method (bottom). Here, the images are of size 512;512 and "101 pixels. Portions of the original images are shown in Fig. 10. 9. Conclusions In this paper we presented a hybrid algorithm for the estimation of point correspondences under small nonrigid motion. Our algorithm presents a general strategy under the assumption that the local nonrigid deformation is limited to a neighborhood of size . Our design philosophy was to extract a reliable set of correspondences from possible inaccurate measurements. As we mentioned in the introduction, the application domain for our algorithm consists of those applications where the issue of feature extraction can be addressed reliably. Feature extraction is very much an application dependent issue the success of the algorithm depends on a careful selection of feature values. While it is possible to envision applications where such a selection is not possible, we believe that our algorithm is applicable to a wide range of problems. References [1] H.H. Chen, T.S. Huang, Maximal matching of 3D points for multiple-object motion estimation, Pattern Recognition 214 (2) (1988) 75}90.
2365
[2] A. Goshtasby, G.C. Stockman, Point pattern matching using convex hull edges, IEEE Trans. Systems Man Cybernet. 15 (5) (1985) 631}637. [3] W. Kim, A.C. Kak, 3-D object recognition using bipartite matching embedded in discrete relaxation, IEEE Trans. Pattern Recognition Mach. Intell. 13 (3) (1991) 224}251. [4] K.S. Arun, T.S. Huang, S.D. Blostein, Least square "tting of two point sets, IEEE Trans. Pattern Recognition Mach. Intell. 9 (5) (1987) 698}700. [5] S. Umeyama, Parametrized point pattern matching and its applications to recognition of object families, IEEE Trans. Pattern Recognition Mach. Intell. 15 (2) (1993) 136}144. [6] C. Kambhamettu, D. Goldgof, Point correspondence recovery in nonrigid motion, Proceedings of the International Conference on Computer Vision and Pattern Recognition, 1992, pp. 222}227. [7] C. Kambhamettu, D. Goldgof, Curvature-based approach to point correspondence recovery in nonrigid motion, Proceedings of the International Conference on Computer Vision and Pattern Recognition, July 1994, pp. 222}227. [8] R. Jain, I.K. Sethi, Establishing Correspondence of Non-Rigid Objects using Smoothness of Motion, IEEE, New York, 1984, pp. 83}87. [9] M. Sallam, K.W. Bowyer, Registering time sequences of mammograms using a two-dimensional image unwarping technique, Proceedings of the Second International Workshop on Digital Mammography, July 1994, pp. 121}130. [10] C. Papadimitriou, K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity, Prentice-Hall, Englewood Cli!s, NJ, 1982. [11] G. Strang, Introduction to Applied Mathematics, Wellesley Cambridge Press, Wellesley, MA, 1995. [12] L.R. Foulds, Combinatorial Optimization for Undergraduates, Springer, Berlin, 1984. [13] D. Ballard, C. Brown, Computer Vision, Prentice-Hall, Englewood Cli!s, NJ, 1982. [14] H.P. Moravec, Rover visual obstacle avoidance, Proceedings of the International Joint Conference on A.I., 1981, pp. 785}790. [15] R. Haralick, L. Shapiro, Computer and Robot Vision, Addison-Wesley, Reading, MA, 1992. [16] F.L. Bookstein, Principal warps: thin plate splines and the decomposition of deformations, IEEE Trans. Pattern Recognition Mach. Intell. 11 (6) (1989) 567}585.