Frontiers of Information Technology & Electronic Engineering www.zju.edu.cn/jzus; engineering.cae.cn; www.springerlink.com ISSN 2095-9184 (print); ISSN 2095-9230 (online)
E-mail: [email protected]
Feature matching using quasi-conformal maps
∗Chun-xue WANG, Li-gang LIU†‡
(School of Mathematical Sciences, University of Science and Technology of China, Hefei 230026, China )
†E-mail: [email protected]
Received Nov. 26, 2015; Revision accepted Mar. 24, 2016; Crosschecked Apr. 22, 2017
Abstract: We present a fully automatic method for finding geometrically consistent correspondences while discarding outliers from the candidate point matches in two images. Given a set of candidate matches provided by scale-invariant feature transform (SIFT) descriptors, which may contain many outliers, our goal is to select a subset of these matches retaining much more geometric information constructed by a mapping searched in the space of all diffeomorphisms. This problem can be formulated as a constrained optimization involving both the Beltrami coefficient (BC) term and quasi-conformal map, and solved by an efficient iterative algorithm based on the variable splitting method. In each iteration, we solve two subproblems, namely a linear system and linearly constrained convex quadratic programming. Our algorithm is simple and robust to outliers. We show that our algorithm enables producing more correct correspondences experimentally compared with state-of-the-art approaches.
Key words: Feature correspondence; Quasi-conformal map; Splitting method http://dx.doi.org/10.1631/FITEE.1500411 CLC number: TP391
1 Introduction
Finding geometrically consistent correspon- dences between pairs of images is one of the most fundamental operations in graphics and vision and it forms the basis of some problems such as feature tracking, image classification or retrieval, object de- tection, and shape matching. Given two images of the objects in different poses or even different ob- jects of the same class, the task is to select a set of corresponding points, one from each image, such that the two points in a pair correspond to the same location. Identifying corresponding points is a chal- lenging problem, as the shape of objects and their pose can significantly change across images.
Different matching strategies have been used among different image cues. Over the past decade, correspondence techniques have evolved significantly
‡Corresponding author
*Project supported by the National Natural Science Foundation of China (Nos. 61672482 and 11626253) and the One Hundred Talent Project of the Chinese Academy of Sciences
ORCID: Chun-xue WANG, http://orcid.org/0000-0002-2118- 3016
cZhejiang University and Springer-Verlag Berlin Heidelberg 2017
in the graphics and vision literature. We refer the interested reader to Heider et al. (2011) and Tuyte- laars and Mikolajczyk (2008) for a survey of the state-of-the-art methods in feature matching. These methods, such as the scale-invariant feature trans- form (SIFT) algorithm, create a collection of candi- date correspondences by matching local signatures.
However, as they consider only local intensity, many methods are globally inconsistent.
To filter the geometrically inconsistent corre- spondences, a low-dimensional parameterized model of deformations was proposed, such as in Chui and Rangarajan (2003), which performs feature- based nonrigid registration with the thin-plate spline (TPS). However, the resulting mapping may contain many incorrect pairs and only a few correct ones, as there is no guarantee that it is bijective. Focusing on the quality of the mapping, Lipmanet al. (2014) selected a large subset of correspondences that are aligned by a global deformation of bounded distor- tion with an alternative approach, but the subset would filter more correct pairs while selecting more incorrect ones as the bounded distortion constraints
are set to be same everywhere. In this study, by bal- ancing between these approaches, we focus mainly on the space of all diffeomorphisms by optimizing the deformation distance and keeping bijectivity.
Given several suitable landmarks, most existing algorithms can find registration in the space of all diffeomorphisms accurately and efficiently only by fixing the boundary vertex onto squares or circles (Zeng and Gu, 2011; Luiet al., 2012; Lam and Lui, 2014). However, the registration problem becomes challenging when all points are set to be landmarks in the feature-matching problem. In this case, bi- jectivity is hard to keep in the obtained registration as the initial correspondences contain many outliers.
To tackle this problem, we formulate an optimiza- tion problem of certain energy functionals over the space of all diffeomorphisms. Specifically, we mini- mize an energy functional involving both a Beltrami coefficient (BC) term and a deformation term, where the first term measures the distortion of the quasi- conformal map and the second measures the geomet- ric distance between corresponding pairs. Different from Lipman et al. (2014), we find an orientation- preserving diffeomorphism ˜f : {pi}Ni=1 → {qi}Ni=1 from a low-dimensional deformation space, where {pi}Ni=1 and {qi}Ni=1 are two input correspondences to be matched.
We propose an efficient splitting algorithm to perform this minimization in which we optimize vari- ables alternatively. Using the proposed algorithm, diffeomorphisms (1-1 and onto) between two input feature data can be effectively obtained, even with the large number of incorrect matches. Numerical results show that our algorithm performs well in a number of experiments of synthetic data and real images over state-of-the-art approaches. The contri- butions of our work can be summarized as follows:
1. We are the first to introduce quasi-conformal maps into the feature-matching problem, which is formulated as a constrained optimization involving both the BCs and quasi-conformal map and solved by an efficient iterative algorithm based on vari- able splitting. Experiments show that our method can select more geometrically consistent correspon- dences compared with state-of-the-art approaches, even when the input contains many outliers.
2. Our method is not sensitive to the choice of the parameters and is robust to outliers. We demon- strate it on a series of synthetic data.
2 Related work
Finding meaningful mapping or matching be- tween corresponding data that optimizes certain kinds of energy functionals has been extensively studied (Montagnatet al., 2001; Nealen et al., 2006;
van Kaick et al., 2011; Zhao et al., 2012). In this section, we review several works on matching point features of nonrigid deformations. We also discuss some previous works closely related to registration via quasi-conformal maps.
1. Low-dimensional deformation spaces. With a set of candidate pairs of correspondences as input, which may include many outliers, several works aim to find a subset belonging to some low-dimensional deformation spaceD. RANSAC (Fischler and Bolles, 1981) finds a large subset of pairs aligned by a global deformation in D up to an ε-deviation denoted as Dε, which is formulated precisely as
minf∈D
N i=1
f(pi)− qi0.
Lipmanet al. (2014) performed a similar formulation to find the maximal subset of bounded distortion denoted asFK, and the problem is stated as follows:
f∈FminK
N i=1
f(pi)− qi0,
where FK is the collection of all deformations de- cided by a distortion value less than or equal to K, and both ε and K decide the deformation space to be low- or large-dimensional. Similarly, our method produces a low-dimensional deformation space rep- resented by an optimal BC, denoted asDb.
2. Space-deformation based methods. Space de- formation has been widely used in feature match- ing, and maximizes the quality of the matching and minimizes the complexity of the deformation simul- taneously. Chui and Rangarajan (2003) minimized the deformation energy based on thin-plate splines, while Belongie et al. (2002) used shape context lo- cal descriptors to optimize a similar energy. Hinton et al. (1991) minimized deformation energy includ- ing a spline-based deformation cost and a generative model of appearance with an elastic matching algo- rithm. Jian et al. (2005) performed a nonrigid reg- istration between these local frequency maps using the Riesz transform. However, these methods may be sensitive to initialization and outliers.
3. Quasi-conformal maps. We review mainly the work related to ours, including conformal maps, quasi-conformal maps, and Beltrami flows. Various conformal geometric methods have been presented for nonrigid surface matching and registration (Yezzi and Mennucci, 2005; Wanget al., 2007; Zeng et al., 2009). Taimouri and Hua (2014) introduced a new quasi-conformal metric which measures the curva- ture changes at each vertex of each pose during the deformation. Given boundary constraints, quasi- conformal maps are widely used to parameterize meshes onto disk domains and obtained by solving the Beltrami equation. Several methods have been proposed to deal with simple domains in the complex plane (Mastin and Thompson, 1984; Daripa, 1991;
1992; Lui and Ng, 2015). Ho and Lui (2016) proposed an algorithm to compute the quasi-conformal param- eterization for a connected 2D domain or surface. A method called the Beltrami holomorphic flow was used to obtain the optimal BC associated with the registration (Luiet al., 2012). Matching landmarks consistently with quasi-conformal mapping has also been proposed. Given several features on the sur- faces, Zeng and Gu (2011) proposed a method to register 3D surfaces with large deformations using quasi-conformal curvature flow. However, bijectiv- ity of the mapping is difficult to guarantee, especially when all feature pairs are set to be landmarks. Our approach can provide a diffeomorphism by setting each energy term a weight parameter properly.
3 Background
In this work, we apply quasi-conformal (QC) maps to find geometrically consistent correspon- dences between pairs of images. In this section, we review some basic theories on quasi-conformal geom- etry. We refer the interested reader to Gardiner and Lakic (2000) and Lehtoet al. (1973) for details.
Given two surfacesM and N, a map f : M → N is conformal if it preserves the surface metric up to a multiplicative factor called the conformal factor.
A generalization of conformal maps is the quasi- conformal maps, which are orientation-preserving homeomorphisms between Riemann surfaces with bounded conformality distortion, in the sense that their first-order approximations take small circles to small ellipses of bounded eccentricity (Gardiner and Lakic, 2000). Supposef : D → D, where D and D
are two domains in C. Then f is quasi-conformal, provided that it satisfies the Beltrami equation
∂f
∂z = μ(z)∂f
∂z (1)
for some complex-valued Lebesgue measurable μ : C → C, satisfying μ∞ < 1 (Bers, 1977). μ is called the Beltrami coefficient, which measures how far the map at each point deviates from a conformal map. Eq. (1) admits a geometrical interpretation.
EquipD with the metric tensor
ds2= Ω(z)2|dz + μ(z)d¯z|2,
where Ω(z) > 0. Then, f satisfies Eq. (1) pre- cisely when it is a conformal transformation from D equipped with this metric to domain D equipped with the standard Euclidean metric, and f is called the μ-conformal. In particular, map f is conformal around the neighborhood ofz when μ(z) = 0.
Consider the effect of the pullback underf of the usual Euclidean metric. The resulting metric is then given by
∂f
∂z
2|dz + μ(z)d¯z|2,
which is relevant to the background Euclidean metric dzd¯z, and has eigenvalues
(1 +|μ|)2
∂f
∂z
2,
and
(1− |μ|)2
∂f
∂z
2.
The eigenvalues represent the squared length of the major and minor axes of the ellipse obtained by pulling back along f the unit circle in the tangent plane (Fig. 1). Accordingly, the dilatation off at a pointz is defined by
K(z) =1 +|μ(z)|
1− |μ(z)|. The dilatation off is given by
K(f ) = sup
z∈D|K(z)| = 1 +μ∞ 1− μ∞.
A simple calculation leads to the following property:
Theorem 1 If f : C → C satisfies μ(f )∞ <
1, then f is a diffeomorphism (Gardiner and Lakic, 2000).
í_ȝ(p _ _ȝ(p _ arg(ȝ(p))/2
Fig. 1 Illustration of how the Beltrami coefficient controls conformality distortion
Theorem 1 plays the fundamental role of obtain- ing a diffeomorphism. That is, given a BCμ : C → C with μ∞ < 1, we can find a one-to-one quasi- conformal mapping fromC onto itself satisfying the Beltrami equation (Eq. (1)).
4 Problem formulation
Given a set of candidate pairs of correspon- dences in the complex plane (pi, qi) ∈ C × C (i = 1, 2, · · · , N ), our aim is to extract a geometrically consistent subset {(pil, qil)}, l = 1, 2, · · · , n ≤ N matched by ˜f ∈ Db(Dbis the space of all diffeomor- phisms). We can set up our problem as
f = arg min˜
f∈Db
N i=1
f(pi)− qi0,
where the0-norm measures the Euclidean distance between f (pi) and qi; i.e., f(pi)− qi0 = 1 if f (pi) = qi, and f(pi)− qi0 = 0 otherwise. In this way, we select matches with small conformal dis- tortion. Unfortunately, it is impossible to construct space Db by adding several constraints, such as in Lipmanet al. (2014). Going back to Eq. (1), we solve Beltrami’s equation by involving anLp-minimization of the Beltrami energy:
EB(f ) =
∂f
∂z − μ∂f
∂z
p. (2)
Note that, if we setp large enough, the optimization problem gives a good approximation. In particular, when p = 2, the result of the least-square problem (2) is called the least-square quasi-conformal map associated to μ (Gu and Yau, 2008; Weber et al., 2012). In this study, we choosep = 2 and solve it with an alternating minimization algorithm.
Going back to Theorem 1, to obtain a diffeo- morphic registration associated with an optimal BC, which is guaranteed to be bijective, our optimization can generally be written as finding ˜f : S1→ S2 that
satisfies f = arg min˜
f E(f)
=
|∇μf|2dS + α
|μf|2dS
+ γ
|f¯z− μffz|2dS + β
N i=1
f(pi)− qi0 s.t. μf∞< 1,
(3)
where μf is the BC associated withf ∈ Db, and we write
S1 as
for short.
The first term ofE is a regularization term en- suring the smoothness of f . The second term of E is minimized to control the conformality distortion off. The third term of E aims to minimize the dis- cretization of the least-square Betrami energy, and the last term is designed to align as many of pi’s as we can withqi’s.
The optimization problem (3) poses two main challenges:
1. The problem involves the optimization of both quasi-conformal mapsf and μf, which depend highly on each other, leading to a highly nonlinear optimization.
2. The L2,0 functional in problem (3) is not smooth and is nonconvex.
We tackle the problem by combining the split- ting method and iterative reweighted least-square (IRLS) method. More specifically, following Lipman et al. (2014), we replace the L2,0 term by a smooth functional to approximate its minima. Sinceμf and f are highly related and the energy is highly nonlin- ear, we introduce a new variable ν and reformulate problem (3) to be an equality-constrained problem.
Thus, the problem can be further formulated as look- ing for an optimal BC ˜ν : S1 → C, which is the BC of some diffeomorphism ˜f : S1→ S2, minimizing the following energy functionalE(ν, f):
(˜ν, ˜f ) = arg min
f,ν E(ν, f)
=
|∇ν|2dS +α
|ν|2dS +γ
|f¯z− νfz|2dS
+ β
N i=1
f(pi)− qi0 s.t. ν∞< 1, ν = μ(f).
(4) Problems (3) and (4) are equivalent to each other.
Based on problem (4), we apply the splitting method
and the IRLS method to solve the constrained opti- mization problem, by replacing theL2,0 functional term and equation constraints, as follows:
(˜ν, ˜f ) = arg min
f,ν Eδsplit(ν, f )
=
| ∇ν|2dS + α
|ν|2dS
+ σ
|ν − μ(f)|2dS + γ
|f¯z− νfz|2dS
+ β
N i=1
(f(pi)− qi2+ δ)p2 s.t. ν∞< 1,
(5)
where p is a small constant, the sequence δ → 0 is a parameter updated during the iterations, and σ is the weight of the penalty term required to in- crease to infinity theoretically. We iteratively min- imize Eδsplit(ν, f ) subject to the constraints. More precisely, given an initial matchf(0), set ν(0) = 0, δ = δ(0). Suppose we have obtained ν(k) and δ(k) at the kth iteration. Fixing ν(k), we obtain f(k) by minimizing Eδsplit(k)(ν(k), f ) over f . Once f(k) is updated, similarly, by fixing f(k), we ob- tainν(k+1) by minimizingEδsplit(k)(ν, f(k)) over ν. If
Eδsplit(k+1)(ν(k+1), f(k+1))− Eδsplit(k)(ν(k), f(k)) < ε, up- dateδ(k).
5 Numerical implementation
5.1 Discretization
We use a Delaunay triangulation of the pla- nar point set {pi} denoted as T = (V, F ), where V = {pi}Ni=1 is the set of vertices and F = {fj}Tj=1 is the set of oriented faces. We choose f from the space of continuous piecewise linear (CPL) map- pings defined by the values at vertices. The com- plex derivatives f¯z and fz are naturally defined on each triangle. Givenei,ej,ek as edge vectors ofT with opposite vertices (i, j, k), the gradient of f is (fiti+ fjtj+ fktk)/(2AT), where AT is area of the triangle, ti = e⊥i , and fi is the value at vertex pi. Then per-triangle derivatives are given by
fz= 1
4AT(fi¯ti+ fj¯tj+ fkt¯k), fz= 1
4AT(fiti+ fjtj+ fktk),
where ti = txi + i tyi is the complex form of ti = (txi, tyi).
The BC μ is a piecewise constant function de- fined on each triangle ofT . Different from Lam and Lui (2014) and Lui and Ng (2015), who computed μ on each vertex and converted the value to each face by averaging the values on each vertex, we re- strict the value ofμ on triangle τ as μτ and compute it by μτ = f¯z|τ/fz|τ (defined above) for each tri- angle, so μ = (μ0, μ1, · · · , μT −1) ∈ RT. For any μ1, μ2, μ ∈ RT, we define the inner product and norm as
(μ1, μ2)RT =
τ
μ1τμ2τAτ,
μRT =
(μ, μ)RT.
For any μ ∈ RT, we define the jump ofμ over an edgee as
[μ]e=
e≺τμτsgn(e, τ ), e ∂T,
0, e ⊆ ∂T,
where e ≺ τ denotes that e is an edge of triangle τ , and sgn(e, τ ) defines the relative orientation of edge e to triangle τ . More specifically, we set all the triangles with an anticlockwise orientation, and all edges are with fixed orientations chosen randomly.
For an edgee, if the orientation of e is consistent with the orientation of τ , then sgn(e, τ ) = 1; otherwise, sgn(e, τ ) = −1. In the piecewise constant function space, the gradient operator can be defined as
∇ : μ → ∇μ, ∇μe= [μ]e. Thus, the term
|∇μ|2dS can be discretized as
ele[μ]2e, where le denotes the length of edge e. Since μτ∞ < 1 is hard to optimize, we have simplified it intoμxτ∞<√
2/2 and μyτ∞<√ 2/2 for allτ ∈ F .
To summarize, we discretize energy problem (5) with complex variables fi to present the map and per-face complex variablesντ for BC:
(˜ν, ˜f ) = arg min
f,ν Esplit(ν, f )
= 1 ne
e
le[ν]2e+α T
τ
|ντ|2+σ T
τ
Aτ|ντ−μτ(f )|2
+ γ T
τ
Aτf¯z|τ − ντfz|τ2 + β
N
N i=1
f(pi)− qi2+ δ p2
s.t. ντj∞<
√2
2 , τ ∈ F, j = 0, 1,
(6)
where ne denotes the number of edges, and Aτ the area of triangleτ .
5.2 Optimization
Energy problem (6) is highly nonlinear. Stan- dard methods, e.g., Newton methods or ADMM methods (Boydet al., 2011; Li et al., 2015), do not perform well for optimizing this energy. Besides, it is quadratic with respect to both vector variablesf = [f (p1), f (p2), · · · , f (pN)] and ν = [ν1, ν2, · · · , νT].
This suggests the use of an alternating-descent algo- rithm to optimizef and ν in an alternating fashion.
More specifically, if f is fixed, we can obtain ν by solving a linearly constrained convex quadratic prob- lem. Similarly, we can easily obtainf by solving a linear system whenν is fixed. While there is no guar- antee of the convergence to a global minimum, this algorithm is relatively stable with consistent practi- cal behavior (Fig. 7).
1.f-subproblem. We discuss the minimization of Esplit(ν(k), f ) over f by fixing ν(k). Note that for each p > 0 and δ > 0, the standard family of functionals
Ep,δ(f) =N
i=1
f(pi)− qi2+ δ p2
is smooth, and we can take a sequenceδ → 0 while p is a small constant. Let us denoteEδ(k)as the energy at thekth iteration:
Eδ(k)=
N i=1
ω(k−1)i f(pi)− qi2. (7)
Since μ(f ) characterizes the distortion of the map (scale, stretch, etc.), we set a small weight to the region of large distortion, that is,
ωi(k−1)= exp
−μ(f(k−1))vi2 2m
· di, (8)
where
di =
f(k−1)(pi)− qi2+ δ(k)
p2−1 ,
μ(·)vi denotes the value of μ(·) on the ith vertex computed by averaging the values of its neighboring faces, and we setm = 0.005 in our implementation.
Then the energy is given by f(k)= arg min
f
γ T
τ
Aτf¯z|τ − ντ(k)fz|τ2 + β
N
N i=1
ω(k−1)i f(pi)− qi2
,
(9)
which is a quadratic problem inf whose solution can be easily obtained by solving a linear system.
2.ν-subproblem. Once f(k)is obtained, we min- imizeEsplitp,δ
ν, f(k)
overν while fixing f(k). That is, we obtainν(k+1)by solving the following problem:
minν Esplit
ν, f(k)
= 1 ne
e
le[ν]2e+ σ T
τ
Aτντ − μτ(f(k))2 +α
T
τ
Aτ|ντ|2+ γ T
τ
f¯z(k)|τ − ντfz(k)|τ2 s.t. ντj∞<
√2
2 , τ ∈ F, j = 0, 1,
(10) whereμ(f ) is the BC associated with f .
Theoretically, the conventional penalty method requires that σ should increase to infinity. In our experiments, we found that setting σ to be a large constant can also give satisfactory results. Obvi- ously, this is a linearly constrained convex quadratic programming problem, optimized over ν with Mat- lab’s quadratic programming routine. The two-step optimization is summarized in Algorithm 1.
6 Experimental results
We tested our algorithm on both synthetic data and real images, and compared it with several exist- ing algorithms, including BD (Lipmanet al., 2014), Tensor (Duchenne et al., 2011), and RANSAC (Fis- chler and Bolles, 1981) algorithms, both visually and quantitatively. For comparison, we tested the Ten- sor method with a cost function to penalize changes in the angles between triplets of feature points and compute three RANSAC models, affine (denoted as RANSAC-AFF), epipolar (denoted as RANSAC- EPI), and projective (denoted as RANSAC-PRO) (using Matlab’s estimateGeometricTransform and estimateFundamentalMatrix function, respectively), with suitable parameters demonstrated in Lipman et al. (2014). In our experiments we chose constants p = 0.001 and σ = 20. We decreased δ by half in
Algorithm 1 Feature matching using quasi- conformal maps
Require: candidate correspondences (pi, qi) ∈ R2×R2, i = 1, 2, · · · , N
Ensure: subset of pairs (pil, qil), l = 1, 2, · · · , n ≤ N , and the bijective mapf
1: Initialization: Triangulation T = {V, F } from {pi}Ni=1
f(0)(pi) = pi,ν(0)= 0
δ(0)≈ diam{pi}, ω0i from Eq. (8)
// diam{·} gives the minimum radius of a circle that contains all the pointspi
k = 0, δmin= 1 × 10−8,kmax= 50 2: while δ(k)> δminandk < kmax do
3: while Eδsplit(k)(ν(k), f(k)) − Esplitδ(k−1)(ν(k−1), f(k−1))
> ε do 4: k = k + 1 5: δ(k)= δ(k−1)
6: Setωi(k−1)according to Eq. (8)
7: Optimize f subproblem via problem (9) with ν(k)fixed
8: Optimize ν subproblem via problem (10) with f(k)fixed
9: end while 10: δ(k)= δ(k)/2 11: end while
12: Return all pairs(pi, qi) for (pi, qi) < ε
each iteration, and set the lower bound as 1× 10−8. Experiments showed that our algorithm is not sen- sitive to these choices. To measure whether f (pi) reaches the location closely to the targetqi, we use the distance threshold of five pixels.
6.1 Synthetic data
To compare our algorithm with other state-of- the-art methods, we first tested it on synthetic data.
We sampled points of a square randomly and uni- formly on the 2D plane, and produced another set of points by a random radial basis function (RBF) mapping of the first one. We produced six RBFs.
For each RBF, we produced n = 64 inlier pairs of points sampled uniformly from an 8× 8 grid p1, p2, · · · , pn, chose 10 control points randomly de- noted as{ck}10k=1, and then obtained the correspond- ing RBF as
f (p) =
10 i=1
φ(d(p, pi))ci+P1(p),
whereφ(·) is an RBF, and we used φ(r) = 1 − 30r2− 10r3+ 45r4− 6r5− 60r3· log r. P1is a vector valued
linear polynomial. Next, following Lipman et al.
(2014), N − n outlier pairs (pi, qi) were added by first creating pi uniformly in the square, i = n + 1, n+2, · · · , N , and then choosing qiaccording to the distribution of outliers estimated from real images.
To evaluate the performance of our algorithm, we carried out 100 trials randomly with an outlier fraction in the range [0, 0.95] for each RBF and com- puted theF -measure for each trial:
F = 2 · precision· recall precision + recall,
where precision is the fraction of retrieved instances that are relevant, and recall is the fraction of relevant instances that are retrieved (Sasaki, 2007). We plot- ted the F -measure of the six RBFs as a function of the fraction of outliers over all trials in Fig. 2, which shows that our algorithm gives a higher F -measure even for a large fraction of outliers.
Fraction outlier 0 0.2 0.4 0.6 0.8 1.0
Fraction outlier 0 0.2 0.4 0.6 0.8 1.0
Fraction outlier 0 0.2 0.4 0.6 0.8 1.0
Fraction outlier 0 0.2 0.4 0.6 0.8 1.0
Fraction outlier 0 0.2 0.4 0.6 0.8 1.0
Fraction outlier 0 0.2 0.4 0.6 0.8 1.0
F-measure
0 0.5 1.0
F-measure
0 0.5 1.0
F-measure
0 0.5 1.0
F-measure
0 0.5 1.0
F-measure
0 0.5 1.0
F-measure
0 0.5 1.0
BD Ours RANSAC-AFF
RANSAC-EPI Tensor
(a) (b)
(c) (d)
(e) (f)
Fig. 2 Performance curves on six RBFs using our method and other methods: (a) RBF 1; (b) RBF 2;
(c) RBF 3; (d) RBF 4; (e) RBF 5; (f ) RBF 6
6.2 Real images
Results on real images are shown in Figs. 3–
5. We tested our algorithm on different kinds of image pairs including images with large deforma- tions, images from different viewpoints, and images
(a) Input images (b) SIFT (c) QC (our algorithm)
(d) BD (K =3) (e) BD (K =20) (f) Tensor
(g) RANSAC-AFF (h) RANSAC-EPI (i) RANSAC-PRO
Fig. 3 Experiments on real images with large deformations. There are nine image pairs: original images (a), SIFT correspondences as input (b), and several methods for filtering the SIFT correspondences (c)–(i).
Disks of the same color denote correspondences in an image pair, and the color varies linearly according to the horizontal position of points in the left image. On the right image of each pair, red lines show the errors between selected correspondences and ground-truth correspondences. The black disks are used to visualize the pairs selected by our method but missed by other methods (References to color refer to the online version of this figure)
of different animals of the same species. We used SIFT descriptors in the VLFeat software package to obtain the initial correspondences in all experi- ments (Vedaldi and Fulkerson, 2010). We compared RANSAC, Tensor, and BD algorithms with distor- tion bound parameterK = 3 and K = 20 character- izing different deformation spaces.
In Figs. 3–5, disks of the same size and color indicate matching pairs to visualize matching results clearly. In the left image of each image pair, we set the color varying according to the horizontal position of the points. For all real images, we used ground- truth correspondences marked by a human observer,
and showed the correspondence error by connecting a red line between each correspondence chosen and the respective ground-truth correspondence in each right image. Our goal is to select as many SIFT cor- respondences compatible with the ground-truth or with small positional deviations (indicated by short red lines), while filtering the ones with large posi- tional deviations (indicated by long red lines). In addition, we used black disks to visualize the pairs selected by our method but missed by other meth- ods. The results show that our method selects more consistent correspondences than the others.
(a) Input images (b) SIFT (c) QC (our algorithm)
(d) BD (K =3) (e) BD (K =20) (f) Tensor
(g) RANSAC-AFF (h) RANSAC-EPI (i) RANSAC-PRO
(j) Input images (k) SIFT (l) QC (our algorithm)
(m) BD (K =3) (n) BD (K =20) (o) Tensor
(p) RANSAC-AFF (q) RANSAC-EPI (r) RANSAC-PRO
Fig. 4 Experiments on real images from different viewpoints: red lines show the errors between selected correspondences and ground-truth correspondences and the black disks in the results are used to visualize the pairs selected by our method but missed by other methods (see caption in Fig. 3 for details). The images were taken with permission from Lazebniket al. (2004; 2005) (References to color refer to the online version of this figure)
(a) Input images (b) SIFT (c) QC (our algorithm)
(d) BD (K =3) (e) BD (K =20) (f) Tensor
(g) RANSAC-AFF (h) RANSAC-EPI (i) RANSAC-PRO
(j) Input images (k) SIFT (l) QC (our algorithm)
(m) BD (K =3) (n) BD (K =20) (o) Tensor
(p) RANSAC-AFF (q) RANSAC-EPI (r) RANSAC-PRO
Fig. 5 Experiments on real images of different animals of the same species: red lines show the errors between selected correspondences and ground-truth correspondences and the black disks in the results are used to visualize the pairs selected by our method but missed by other methods (see caption in Fig. 3 for details). The images were taken with permission from Lazebniket al. (2004; 2005) (References to color refer to the online version of this figure)
6.2.1 Images with large deformations
For images with large deformations, the SIFT descriptors are less discriminative and contain more incorrect correspondences, and it is difficult to filter correct correspondences. Fig. 3 shows a scene where two objects are close to each other in the left im- age while away from each other in the right one. In this case, the distances are preserved between fea- ture points within each object while changing sig- nificantly across the two objects, and our aim is to
capture the large deformation while discarding the outliers. Fig. 3 shows that previous methods cannot capture correct pairs near the eyes since the input SIFT matches contain many incorrect ones. In con- trast, our method selects more good correspondences even though the SIFT matches are of poor quality.
6.2.2 Images from different viewpoints
We tested our method on the images from dif- ferent viewpoints (Fig. 4). It is a challenging test case for images from different viewpoints since some
feature points are captured in the first image, while they are in shade in the second one. In addition, the location of one object with respect to another or background may change significantly. Exam- ples in Fig. 4 demonstrate that previous works ei- ther miss more good pairs (indicated by black disks) and introduce more false ones (indicated by long red lines) or fail to recognize symmetry and repetitions (RANSAC-AFF). Our method can capture repeated structures or breaks correctly and choose more pairs fitting the input with fewer outliers.
6.2.3 Images of different animals of the same species Fig. 5 shows the performance of our algorithm over animal images from the natural world. We tested the images of different animals of the same species. Since the SIFT descriptor is unstable, many outliers are included, but our algorithm can still per- form well. The first example in Fig. 5a shows that RANSAC-PRO and BD (K = 3) lose more correct correspondences and other algorithms always select more incorrect ones, while we can balance them well especially on the body of the bird. For the other example, both RANSAC-EPI and RANSAC-PRO lose more correct correspondences on the right wing, while RANSAC-AFF cannot filter the outliers effi- ciently. However, we can locate the left and right wings correctly with more consistent matches.
7 Discussion
7.1 Quality
To quantify the correspondences retrieved, we present theF -measure in Table 1 on the examples in Figs. 3–5 of the tested algorithms. For each exam- ple, since there might be some discrepancy between
a human observer and the geometrically consistent SIFT correspondences, we define ground-truth as an acceptable distance threshold of five pixels from the human observer. From Table 1, we can find that our algorithm achieves largerF -measures than the other algorithms.
7.2 Bijectivity
Following Section 3, we investigated the bijec- tivity of deformationf by computing the norm of the BC on each face. Fig. 6 shows the histograms of the BC norms of examples from Figs. 3–5. All the BC norms are less than 1.
7.3 Running time
We have implemented the alternating-descent algorithm using MATLAB on a laptop with Pentium IV, 2.16 GHz CPU, and 2 GB RAM. We solved f us- ing a sparse linear system andν using quadratic pro- gramming in Matlab. Table 2 gives the model statis- tic and running time of the examples from Figs. 3–5.
Our method takes more time than RANSAC and BD algorithms, as it solves a quadratic program and a linear system in each iteration.
7.4 Convergence
We solve energy problem (6) by the alternating- descent algorithm discussed in Section 5. Although it is difficult to prove the convergence of the splitting algorithm, the energy decreases consistently through iterations, similar to the coordinate descent algo- rithms (Wright, 2015). We further examine the con- vergence of our algorithm via numerical experiments and plot the energy curves of examples in Figs. 3–5 (Fig. 7). The results show that the energy always decreases during the iterations.
Table 1 Comparison ofF-measure of real images
Method F -measure
Fig. 3 Figs. 4a–4i Figs. 4j–4r Figs. 5a–5i Figs. 5j–5r
BD (K = 3) 0.919 0.943 0.796 0.896 0.847
BD (K = 5) 0.899 0.922 0.773 0.854 0.806
BD (K = 20) 0.835 0.808 0.564 0.714 0.728
RANSAC-AFF 0.910 0.938 0.773 0.862 0.786
RANSAC-EPI 0.851 0.848 0.782 0.705 0.561
RANSAC-PRO 0.924 0.890 0.733 0.902 0.776
Tensor 0.859 0.812 0.664 0.744 0.813
Ours 0.945 0.972 0.858 0.962 0.959
Number of faces Number of faces (a)
(d) (e)
(b) (c)
600
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 500 400 300 200 100 0
1200 1000 800 600 400 200 0
Number of faces
140 120 100 80 60 40 20 0
Number of faces
120 100 80 60 40 20 0
Number of faces
400
300
200
100
0
Norm 0 0.10.20.30.40.50.60.70.80.9
Norm 0 0.10.20.30.40.50.60.70.80.9 Norm
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 Norm
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 Norm
Fig. 6 Histograms of the BC norm of examples in Fig. 3 (a), Figs. 4a–4i (b), Figs. 4j–4r (c), Figs. 5a–5i (d), and Figs. 5j–5r (e)
Table 2 Comparison of running time
Method
Running time (s)
Fig. 3 Figs. 4a–4i Figs. 4j–4r Figs. 5a–5i Figs. 5j–5r
(|v|=1905, |f|=3980) (|v|=2925, |f|=6064) (|v|=495, |f|=908) (|v|=533, |f|=980) (|v|=902, |f|=1690)
BD 215.491 441.196 52.367 60.543 105.553
RANSAC-AFF 0.028 0.151 0.411 0.257 0.092
RANSAC-EPI 10.436 10.543 7.145 7.361 7.537
RANSAC-PRO 0.023 0.790 0.533 0.522 0.572
Tensor 792.466 1827.547 132.693 169.544 321.593
Ours 242.911 410.278 54.736 74.501 122.583
0 10 20 30 40 50
0.4 0.5 0.6 0.7 0.8 0.9 1.0
Number of iterations
Energy
Fig. 3c Fig. 4c Fig. 4l Fig. 5c Fig. 5l
Fig. 7 Energy curves of examples in Figs. 3–5. These curves show that the energy decreases during the it- erations and yields a stable result
7.5 Choice of parameters
We have tested our algorithm with different pa- rameter settings on RBF 1. Specifically, we have ap- plied our algorithm with different values of σ in the interval (10, 30) on data RBF 1 andε in the interval (2, 8) to see the average performance with different outlier fractions of data RBF 1. Fig. 8 shows that our algorithm is not sensitive to all of these choices.
8 Conclusions
In this work, we considered the problem of find- ing a geometrically consistent set of correspondences between two input images. Given a set of candidate matches provided by SIFT descriptors, which may include many outliers, we selected a subset of these correspondences that can be aligned perfectly from
2 3 4 5 6 7 8 0.80
0.85 0.90 0.95 1.00
Distance threshold
)íPHDVXUH Precision Recall
0 0.2 0.4 0.6 0.8 1.0
0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
)UDFWLRQRXWOLHU
FíPHDVXUH
ı=10 ı=15 ı=20 ı=25 ı=30 (a)
(b)
Fig. 8 The performance of our method with different choices of parameterσ (a) and the effect of varying the distance thresholdε (b). The performance is stable for a wide range of parameter values
the space of all diffeomorphisms. We formulated this problem as a constrained optimization problem in- volving both the BC term and quasi-conformal map term, and solved it using a splitting method. In each iteration of our algorithm, we solved thef sub- problem by solving a linear system and ν subprob- lem by linearly constrained convex quadratic pro- gramming. We further examined the convergence of our algorithm via numerical experiments. Exper- iments showed that our algorithm is not sensitive to the parameter choice and is robust to outliers, and produces excellent results on synthetic data and real images in comparison to the state-of-the-art ap- proaches.
We dealt with the one-to-one correspondence of geometrically consistent features in this work, but there are still many issues to be solved. Given two sets of features, which may contain different numbers of features, we hope that we can generalize this algo- rithm to find the overall consistent correspondence mapping (one-to-one or one-to-many). In addition,
since SIFT may not be stable when the deforma- tions and lighting changes are large, how to build a descriptor capturing more geometric information would be an interesting topic for future work.
References
Belongie, S., Malik, J., Puzicha, J., 2002. Shape matching and object recognition using shape contexts. IEEE Trans. Patt. Anal. Mach. Intell., 24(4):509-522.
http://dx.doi.org/10.1109/34.993558
Bers, L., 1977. Quasiconformal mappings, with applications to differential equations, function theory and topology.
Bull. Am. Math. Soc., 83(6):1083-1100.
http://dx.doi.org/10.1090/S0002-9904-1977-14390-5 Boyd, S., Parikh, N., Chu, E., et al., 2011. Distributed
optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach.
Learn., 3(1):1-122.
http://dx.doi.org/10.1561/2200000016
Chui, H., Rangarajan, A., 2003. A new point matching algorithm for non-rigid registration. Comput. Vis.
Image Understand., 89(2-3):114-141.
http://dx.doi.org/10.1016/S1077-3142(03)00009-2 Daripa, P., 1991. On a numerical method for quasi-conformal
grid generation. J. Comput. Phys., 96(1):229-236.
http://dx.doi.org/10.1016/0021-9991(91)90274-O Daripa, P., 1992. A fast algorithm to solve nonhomoge-
neous Cauchy-Reimann equations in the complex plane.
SIAM J. Sci. Stat. Comput., 13(6):1418-1432.
http://dx.doi.org/10.1137/0913080
Duchenne, O., Bach, F., Kweon, I.S., et al., 2011. A tensor- based algorithm for high-order graph matching. IEEE Trans. Patt. Anal. Mach. Intell., 33(12):2383-2395.
http://dx.doi.org/10.1109/TPAMI.2011.110
Fischler, M.A., Bolles, R.C., 1981. Random sample consen- sus: a paradigm for model fitting with applications to image analysis and automated cartography. Commun.
ACM, 24(6):381-395.
http://dx.doi.org/10.1145/358669.358692
Gardiner, F.P., Lakic, N., 2000. Quasiconformal Teichmüller Theory. American Mathematical Society, Providence, USA. http://dx.doi.org/10.1090/surv/076
Gu, X.D., Yau, S.T., 2008. Computational Conformal Ge- ometry. International Press, Somerville, MA, USA.
Heider, P., Pierre-Pierre, A., Li, R., et al., 2011. Local shape descriptors, a survey and evaluation. Eurographics Workshop on 3D Object Retrieval, p.1-8.
http://dx.doi.org/10.2312/3DOR/3DOR11/049-056 Hinton, G.E., Williams, C.K.I., Revow, M.D., 1991. Adap-
tive elastic models for hand-printed character recogni- tion. 4th Int. Conf. on Neural Information Processing Systems, p.512-519.
Ho, K.T., Lui, L.M., 2016. QCMC: quasi-conformal pa- rameterizations for multiply-connected domains. Adv.
Comput. Math., 42(2):279-312.
http://dx.doi.org/10.1007/s10444-015-9424-1
Jian, B., Vemuri, B.C., Marroquin, J.L., 2005. Robust nonrigid multimodal image registration using local fre- quency maps. Biennial Int. Conf. on Information Processing in Medical Imaging, p.504-515.
http://dx.doi.org/10.1007/11505730_42
Lam, K.C., Lui, L.M., 2014. Landmark and intensity-based registration with large deformations via quasi-conformal maps. SIAM J. Imag. Sci., 7(4):2364-2392.
http://dx.doi.org/10.1137/130943406
Lazebnik, S., Schmid, C., Ponce, J., 2004. Semi-local affine parts for object recognition. British Machine Vision Conf., p.779-788. http://dx.doi.org/10.5244/C.18.98 Lazebnik, S., Schmid, C., Ponce, J., 2005. A maximum
entropy framework for part-based texture and object recognition. ICCV, p.832-838.
http://dx.doi.org/10.1109/ICCV.2005.10
Lehto, O., Virtanen, K.I., Lucas, K.W., 1973. Quasiconfor- mal Mappings in the Plane. Springer New York.
Li, Y., Xie, X., Yang, Z., 2015. Alternating direction method of multipliers for solving dictionary learning. Commun.
Math. Stat., 3:37-55.
http://dx.doi.org/10.1007/s40304-015-0050-5
Lipman, Y., Yagev, S., Poranne, R., et al., 2014. Fea- ture matching with bounded distortion. ACM Trans.
Graph., 33(3):26.
http://dx.doi.org/10.1145/2602142
Lui, L.M., Ng, T.C., 2015. A splitting method for diffeomor- phism optimization problem using Beltrami coefficients.
J. Sci. Comput., 63(2):573-611.
http://dx.doi.org/10.1007/s10915-014-9903-4
Lui, L.M., Wong, T.W., Zeng, W., et al., 2012. Optimization of surface registrations using Beltrami holomorphic flow.
J. Sci. Comput., 50(3):557-585.
http://dx.doi.org/10.1007/s10915-011-9506-2
Mastin, C.W., Thompson, J.F., 1984. Quasiconformal map- pings and grid generation. SIAM J. Sci. Stat. Comput., 5(2):305-310. http://dx.doi.org/10.1137/0905022 Montagnat, J., Delingette, H., Ayache, N., 2001. A review
of deformable surfaces: topology, geometry and defor- mation. Image Vis. Comput., 19(14):1023-1040.
http://dx.doi.org/10.1016/S0262-8856(01)00064-6 Nealen, A., Müller, M., Keiser, R., et al., 2006. Physically
based deformable models in computer graphics. Com- put. Graph. For., 25(4):809-836.
http://dx.doi.org/10.1111/j.1467-8659.2006.01000.x Sasaki, Y., 2007. The Truth of the F -measure. School of
Computer Science, University of Manchester.
Taimouri, V., Hua, J., 2014. Deformation similarity measure- ment in quasi-conformal shape space. Graph. Models, 76(2):57-69.
http://dx.doi.org/10.1016/j.gmod.2013.12.001
Tuytelaars, T., Mikolajczyk, K., 2008. Local invariant fea- ture detectors: a survey. Found. Trends Comput.
Graph. Vis., 3(3):177-280.
http://dx.doi.org/10.1561/0600000017
van Kaick, O., Zhang, H., Hamarneh, G., et al., 2011. A survey on shape correspondence. Comput. Graph.
For., 30(6):1681-1707.
http://dx.doi.org/10.1111/j.1467-8659.2011.01884.x Vedaldi, A., Fulkerson, B., 2010. Vlfeat: an open and
portable library of computer vision algorithms. Proc.
18th ACM Int. Conf. on Multimedia, p.1469-1472.
http://dx.doi.org/10.1145/1873951.1874249
Wang, S., Wang, Y., Jin, M., et al., 2007. Conformal geometry and its applications on 3D shape matching, recognition, and stitching. IEEE Trans. Patt. Anal.
Mach. Intell., 29(7):1209-1220.
http://dx.doi.org/10.1109/TPAMI.2007.1050
Weber, O., Myles, A., Zorin, D., 2012. Computing ex- tremal quasiconformal maps. Comput. Graph. For., 31(5):1679-1689.
http://dx.doi.org/10.1111/j.1467-8659.2012.03173.x Wright, S.J., 2015. Coordinate descent algorithms. Math.
Program., 151(1):3-34.
http://dx.doi.org/10.1007/s10107-015-0892-3
Yezzi, A., Mennucci, A., 2005. Conformal metrics and true
“gradient flows” for curves. ICCV, p.913-919.
http://dx.doi.org/10.1109/ICCV.2005.60
Zeng, W., Gu, X.D., 2011. Registration for 3D surfaces with large deformations using quasi-conformal curvature flow. CVPR, p.2457-2464.
http://dx.doi.org/10.1109/CVPR.2011.5995410 Zeng, W., Hua, J., Gu, X., 2009. Symmetric conformal
mapping for surface matching and registration. Int. J.
CAD/CAM, 9(1):103-109.
Zhao, Z., Feng, X., Teng, S., et al., 2012. Multiscale point correspondence using feature distribution and frequency domain alignment. Math. Probl. Eng., 2012:382369.
http://dx.doi.org/10.1155/2012/382369
