- Email: [email protected]
Finding point correspondences using simulated annealing
Pattern Recognition, Vol. 28, No. 2, pp. 231 240, 1995 Elsevier Science Ltd Copyright © 1995 Pattern Recognition Society Printed in Great Britain. All rights reserved 0031 3203/95 $9.50+.00
Pergamon
0031-3203(94) 00087-5
F I N D I N G POINT CORRESPONDENCES USING SIMULATED ANNEALING J. P. PASCUAL STARINKt~ and ERIC BACKERt t Delft University of Technology, Department of Electrical Engineering, P.O. Box 5031, 2628 CD Delft, The Netherlands
(Received 15 September 1993; in revisedform 24 June 1994; receivedfor publication 20 July 1994) Abstract--Identifying corresponding points between two recordings of a point set has always been an important problem in stereo vision applications. We describe this matching problem in terms of cost minimization and present an algorithm to approach the minimal cost mapping using simulated annealing. The algorithm calculates the costs to match all possible point pairs and tries to minimize the sum of the costs of all matched points. Starting from an initial mapping, it uses a random rearrangement scheme to alter the mapping towards the optimal (minimal cost) mapping. Matching Point correspondences Simulated annealing
Search methods
1. INTRODUCTION Stereo vision is concerned with the recording of an object space from slightly different views. Identifying corresponding elements in the views plays a key role. Generally, these relations are not bijective; elements in one set may be absent in the other, and vice versa. This, and the lack of a unifying theory makes the correspondence problem severe and has led to an overwhelming diversity of approaches to solve the problem. ~11To arrive at a structural description of the object space, the problem can be approached in different ways, depending on the representation of the information. From low-level to high-level, these approaches are signal matching, feature matching, and relational matching. In the signal approach, the images are treated as 2D signals consisting of a distribution of grey value intensities. It makes use of the well-explored concepts of signal processing. The elements are usually pixel neighborhoods, which are detected in one image and sought for in the other, for example by grey value correlation or least squares matching. The feature approach was introduced by Marr and coworkers. ~2-4~ They used structureless random dot stereograms to show that the human visual system perceives stereo by detecting elementary tokens rather than matching at a higher, structural level. The elements are local, such as points and edges. They are detected in both images and matched afterwards. In the structural approach, also referred to as relational matching, structures and their relationship to neighboring structures are des-
Stereo vision
cribed by a symbolic representation organized in graphs. 1~'6) Due to the different view points of the recordings, the graphs differ in topology and are therefore matched using inexact graph matchingF I Similarity measures between the detected elements lead to an initial set of possible correspondences. If costs are defined reciprocal to the similarity measures, the minimal cost mapping represents the optimal match. Ullman~8) reformulates his surface model of minimal mapping as a linear programming problem. Also dynamic programming has been used in the minimal cost approach. Otha and Kanade 19~ use dynamic programming for the search both along and across the epipolar lines. In relaxation labeling an element labelling is updated iteratively to achieve a globally consistent interpretation,I~°'t ~) but the final solution depends heavily on the starting point and tends to relax in local optima. "2) Barnard ~t3) bases his (signal) matching approach on simulated annealing to overcome these disadvantages. Point matching belongs basically to the feature approach, but may include flavors of the other two approaches in defining similarity measures.
:~Current address: Max Planck Institute for Biophysical Chemistry, Department of Molecular Biology,P.O. Box 2841, D-37018 Grttingen, Germany. 231
2. MODEL Consider a point set in 3D space. We are given the primary and the secondary projection of this set onto two planes as the points p~, i = 1...Np, and sj, j = 1...N~. Matching is accomplished by selecting a set of correspondences that forms an unambiguous mapping with respect to the desired lay-out. We consider two different lay-out approaches: (i) Single correspondences (L~a ~) (1-to-l, 1-to-0, 0-to-l)
232
J.P.P. STARINK and E. BACKER
Both projection sets probably contain a different number of points, due to occlusion, overprojection, etc. In the single approach, a point Pi is allowed to match either a point sj or the null space, and vice versa. (ii) Multiple correspondences (Z~a,1) (1-to-k, k-to-l, k > 0) In the multiple approach, each point is allowed to match any number of points in the other set. By this, also occluded points can be matched and therefore localized. Reducing the candidate lists is commonly accomplished by setting an upper boundary on the parallax by limiting the allowed distance of the candidates to the epipolar lines. The locations and directions of these lines are determined by the motion parameters. If these parameters are unknown, they can be estimated from a limited number of sure correspondences. For weak perspective and parallel projections, four correspondences of three views are sufficient, la) It is possible to extract the directions and locations of the epipolar lines from four points of two views, ~¢~ but it is impos-
enough so that the system does not get stuck in local minima of some energy function. Originally developed by Kirkpatrick et al., t2°'2 ~ the method has been successfully applied to a variety of hard optimization problems in different fields. Among them, the design of good codes t22~ and attempts to solve the traveling salesman problem. 123)It is also used in general image processing problems,~24) in edge detection,t25~ and stereo matching,t t 3) Here, we are given the finite set ~ of all possible mappings and the energy (cost) function E(J¢). The probability of occurrence of any particular mapping m is proportional to its Boltzman weight, P(m)ocexp(E(m)/T), where T is the temperature. For each state m e J ¢ there is a set N(m) c J//that contains the neighbors of m. In addition, a transition probability matrix R is defined such that R ( m , m ' ) > 0 if and only if m'~N(m). Finally, let T~, T2.... be a sequence of strictly positive numbers such that T~ > T2 > ... and lim T k O. With IX] + = =
k~oo
max (X, 0), the general form of the annealing algorithm is:
Set k = 0 Choose initial mapping X k while Tk # 0 do Choose next state Yk from N(Xk) with probability R(m, m')
f
Set Xk
+~ Set k = k + end.
rk
), Xk
with probability exp F - [ E ( Yk)- E(Xk)] +]
Tk
L
J
otherwise
1
sible to determine the motion parameters. I~s'16~ For a review on perspective projections, see Sabata and Aggarwal.t~ 7~ A cost function c(i, j) is defined by specifying the pair-cost charged when points p~ and sj are matched. This function is of major importance because additional discriminative power between these pair-costs reduces the search time. It should therefore reflect both general and application specific properties, such as grey value neighborhood similarity and relational properties like rigidity.~l sl 2.1. Simulated annealing Simulated annealing is a stochastic optimization algorithm based on the physical analogy of annealing a system of molecules to its ground state. To grow a perfect crystal, one starts at high temperature and then gradually cools the substance, staying as close to equilibrium as is practical. This is directly analogous to what happens in simulated annealing. To bring a synthetic system to equilibrium, the cooling process is simulated using the standard method of Metropolis et al.~191The rate of cooling is controlled by the macroscopic parameter 'temperature', and must be slow
The algorithm produces a discrete-time Markov chain X = {Xklk > 0}. Let J¢* be the set of states in ~ ' at which E attains its minimal value. The idea of the simulated annealing algorithm is to try to achieve lim P [ X k e J / ¢ * ]
= 1
(1)
k~
by letting Tk tend to zero as k leans to infinity. The design of the temperature schedule and the transition distribution R are based on a theorem by Hajek. t26~Before stating this theorem, some additional definitions are necessary: • A state m' is reachable from state m if there is a sequence of states m = m o, ml,.. •, mv = m' for some p > 1 such that ink+ 1eN(mk) for 0 < k < p. • (J¢, N) is irreducible when for any two states m and m', m' is reachable from m. • State m' is reachable at height H from state m if there is a sequence of states m = m o, ml,... ,rap = m' for some p > 1 such that mk÷ 1~N(mk), 0 < k < p, and E(mk) < H for 0 < k < p. • (d/, R, E) is weakly reversible when for any real number H and any pair of states m and m', m is
Finding point correspondences using simulated annealing reachable at height H from m' if m' is reachable at height H from m. • A state m is a local minimum if no state m' with E(m') < E(m) is reachable from m at height E(m). • A state m is a global minimum if no state m' exists for which E(m') < E(m). • The depth of m is + ~ if m is a global minimum, otherwise it is the smallest n u m b e r D _> 0 such that some state m' with E(m') < E(m) can be reached from m at height E(m) + D.
Theorem. Assume that (./¢, R,E) is irreducible and weakly reversible, then the convergence in probability of equation (1) is reached if and only if exp f - d * ' ]
k=1
\
= +oo
(2)
Tk/
where d* is the maximum of depths of all states that are local but not global minima. This result is consistent with the result of other researchers who have established the asymptotic convergence of the simulated annealing algorithm. 127-29) If the temperature schedule assumes the parametric form e
Tk -- - log (k + 1)
(3)
then equation (2), and hence also equation (1), is true if and only ifc _> d*. To limit the n u m b e r of iterations, c should be kept as small as possible. If R is reversible, then, by definition, there is a probability distribution ~ on ~ , such that ~(x)R(x,y)= ~(y)R(y, x) for all x, y e ~ . Necessarily, ~ is the equilibrium distribution for R. A simple transition distribution for which this is true, is given by:
R(x,y) =
233
Table 1. The possible matching combinations between any two correspondences in the single lay-out (S) and the multiple lay-out (M) 0
0 Po p~ Psm p,~
SO
Sss
Ssm
Sm
--
--
S,M
M
--
-S,M M --
S,M S,M M M
S,M S,M M M
M M M M
M M M --
Breaking a correspondence decreases the mapping energy, while creating a new correspondence increases the mapping energy with the corresponding pair-cost. If the change in mapping energy AE is negative, the rearrangement is accepted, otherwise it is accepted with a chance according to the Boltzman probability distribution. The rearrangements are illustrated in Figs 1-4. The configurations left of the arrow represent the situation before, the configuration right of the arrow the situation after the rearrangement. The null space is drawn as a circle on top of the sets. The candidates are grey, p in the left set and s in the right set, If a candidate is the null space, it is indicated with 0. If a point is matched to the null space, it is given the subscript 0. If it is matched to two or more points, the subscript is m (multiple). In case of a single correspondence, there is a two-character subscript, of which the first is s (single). The second character reflects the n u m b e r of elements that its counterpart is matched to. This n u m b e r is either o n e - - t h e candidate itself--or more, represented by s and m, respectively. Table 1 lists the allowed configurations of any two candidates in the two layouts.
if y e N ( x ) otherwise.
(4)
2.2. The matching algorithm To generate an arbitrary, initial mapping, we have used a simple best-first scheme: repeat selecting the minimal cost correspondence from the unmatched points until at least one set is empty. Remaining unmatched points are matched to the null space. Value c = d* is set to the maximal pair-cost in the initial mapping. Although this value is surely too high and results in longer convergence times, it proves to be an adequate guess. Let a rearrangement be a change in the mapping such that the previously unmatched points p~ and st, which are randomized, become a new correspondence, and the resulting mapping m' is part of N(m). The points that are matched to pi and sj are denoted by ~ and gj, respectively. When a multiply matched point is involved, the correspondence to be possibly changed is randomized. When the number of possible rearrangements exceeds one (IN(m) l > 1), a random one is picked according to equation (4).
2.2.1. Single correspondences-- &P~. When a mapping should contain only one-to-one correspondences and unmatched points, there are six possible configurations: 0 - s s, Ps - 0, Po - So, Po - s~, ps - s o and p~ - s s. All can be rearranged in only one way, as illustrated in Fig. 1. The changes in mapping energy due to these rearrangements are given by: AE = - 2c(gi, sj) + c(gj, 0) + c(0, sj)
(5a)
AE = - 2c(pi, Pi) + c(pi, 0) + c(0, 0i)
(5b)
AE = - c(pi, 0) - c(0, sj) + 2c(pi, s) )
(5c)
A E = -c(pi, O)-c(gj, sj) + c(pi, si) + c(gj, O)
(5d)
AE = - c(p/, Pl) - c(0, st) + c(pi, st) + c(0, Pi)
(5e)
AE = - c(pi, p~) - c(gj, sj) + c(p~,sj) + c(~, g2).
(5f)
2.2.2. Multiple correspondences--ZP". To describe the rearrangements in a multiple correspondence mapping, the possible configurations are divided into three sets. The first set consists of the single correspondences, which are described in Section 2.2.1, augmented by four additional ones. These four correspondences play
234
J.P.P. STARINK and E. BACKER
(I,
0 -
,2 - -
m
__
- -
m
m
Pss - 0
(a)
Sss
.~
(b)
C'
PO
" Sss
(d)
Pss
-
So
-7 PO
(e)
-
-7
SO
O
C,
Pss -
Sss
(c)
C
'D
(f)
Fig. 1. The rearrangements for the six possible configurations of the candidates p and s in a single correspondence lay-out.
A E = - c ( g j , sj) + c(pi, sj) AE
=
c(pi, [h) - c(gj, sj) + c(p i, sj) + c(g~,0)
-
AE = - c(p~, Pi) + c(pi, s~) AE = - c(p~, [h) - c(gj, s~) + 2c(E, %). PO
O
-
Sss
(a)
Pss " So
(b)
O
(c)
Pss
-
Sss
(7g) (7h) (7i)
The third set contains six configurations involving multiple-multiple matched points (subscript m). Each can be rearranged in only one way, see Fig. 4. The energy changes are given by: AE
Pss - Sss
(7f)
c(p~, 0) + c(pi, sj)
(8a)
AE = - c(0, sj) + c( E, sj)
(8b)
AE = - 2c(p i, Pi) + c(E, sj) + c(0, Pi)
(8c)
A E = - 2c(gj, si) + c(pi, s~) + c(gj, 0)
(8d)
AE = - c(pi, Pi) + c(pi, s j)
(8e)
AE = - c(gj, sj) + c(pi, sj).
(Sf)
=
-
(d)
Fig. 2. Four rearrangements for the three possible configurations of singly matched candidates p and s resulting in multiple correspondences.
an essential role in the construction of multiple correspondences. The rearrangements of these four are shown in Fig. 2. The related changes in mapping energy are: AE
=
-
c(pi, 0) + c(pi, sj)
(6a)
AE
=
-
c(0, s j) + c(p i, sj)
(6b)
AE = - 2c(pi, Pi) + c(pi, s j) + c(0, Pi, 0)
(6c)
A E = - 2c(gj, sj) + c(pi, s j) + c(gj, 0).
(6d)
The second set holds nine configurations involving the single-multiple matched candidates (subscript sin). Their rearrangements are illustrated in Fig. 3. The changes in mapping energy are:
A E = - c ( p , , pi) + c(O, p~)
(7a)
AE = -c(p~, p~) + c(pl,0)
(7b)
AE = -- c(pi, 0) -- c(gj, sj) + 2c(p~, s j)
(7c)
AE = - c(p~, Pi) - c(0, sj) 6+ 2c(pi, sj)
(7d)
AE = - c ( p i , Pi) - c(gb,s j) + c(E, s i) + c(0, Pi)
(7e)
3. EXPERIMENTALRESULTSAND DISCUSSION In this section, three series of experiments are described to determine the rate of convergence of the matching procedure. In the first series, the mapping is searched within preset cost matrices, in the second within r a n d o m cost matrices. In the third series, the cost matrices are calculated from stereo views of a generated point set. An example from the field of cell biology is included; the localization of DNA. The cost function in the experiments is based on the epipolar constraint only. To obtain the pair-costs, consider the point correspondence (Pi, sj) for i > 0 and j > 0. If the epipolar line of point p/ is 2~ and the positional error of the points is Crpos, then the positional error of ;t~ is also %os. If we assume that the error is normally distributed, then the positional error of sj relative to J,~is ,~-2Crpos. With a probability of 95%, the true correspondence of Pi is located within distance dmax = 1.96 ~/2~po~ ~ 2.8 apo~ of the epipolar line. The cost function is now defined as the square of the distance
235
Finding point correspondences using simulated annealing
O g, C,/~__ C,~D ,D ,D
0o
mT_-_ _
0 - sm
(a)
__
Psm -
m m,~ mxm _
--
)--
__
)__
0
(b)
O C, C~'~ _
~
=×= [] roman'
O O
(d)
O C~
O~,
Psm - Sss
Pss - Ssm
(g)
Psm
(e)
C,
Sss
-
c,,~
(c)
OO
lli~~ _ _
_
Psrn - SO
PO - Ssm
Pss - Ssm
(f)
,D
~
m
m~m m~mX
(h)
Psm
O
- Ssm
O O
(i)
Fig. 3. The nine allowed rearrangements for the seven single/multiple correspondences of the candidates p and s in the multiple correspondence lay-out.
C,/"
,D ,D 'ID~o
_-7
_,
PO - Sm
©
Q O
C~~ • ;, _711 _
_ (a)
Q C C,/D
Pm-
so
:
,D
(b)
,D 0
C,
Srn
(c)
O Q
Q C,
Pss
-
m/m - 'l' m m --
Pm
-
Sss
(d)
Psm
7"
- Sm
(e)
Pm - Ssm
(f)
Fig. 4. The six rearrangements for the configurations of the candidates p and s involding multiply matched points in the multiple correspondence lay-out.
d(s/, 2i) of candidate sj to 2i when it is smaller than d . . . . and +o0 otherwise. For i = 0 or j = 0 , it is set to 1.2 d2ax, so:
d(sj, ,~.i)2, c(i, j) =
+ oo, 1.2 d .2. . .
i > 0, j > 0, d($j,,~i) < dmax i > 0, j > 0, d($j, ~i) ~" dmax
c of the temperature schedule was set to the difference between the m i n i m u m and the maximum of the paircosts in the initial mapping. 3.1. Preset cost matrix
(9)
i = 0 o r j = 0.
If dmax is unknown, the candidate lists may be reduced by selecting t h e / t candidates closest to the epipolar line and assigning + oo to all other pair-costs. For i = 0 or j = O, c(i, j) is set to the most expensive of the /~ candidates. In all experiments, the initial mapping was constructed by matching all the points in both sets to the null space, the maximal cost mapping. The parameter
A cost matrix of 100 x 100 elements was generated such that the true optimal mapping was generated by the best-first search algorithm. Since the optimal mapping equaled the true mapping, it could be recognized when randomized. The positional error was u n k n o w n and could not be used to calculate the pair-costs, therefore a fixed n u m b e r of candidates for each point was selected: 10, 30, 60 and 100 candidates, respectively. The results, plotted in Fig. 5, show that, although not always in the plotted region, both the single and
236
J.P.P. STARINK and E. BACKER Preset Cost Matrix single approach
Preset Cost Matrix multipleapproach
loo
100
"fl
8 == -~ 80 N
20
~,,{,
E 40
:i. '('.~I^
......... :,
20
.
:y, ,----~
rearrangements(xlO0)
l 150
O
200
i
i
r
50
100
150
200
rearrangements (x1130)
Fig. 5. Mapping cost as a function of the number ofaccepted rearrangements in a 100 x 100matrix holding a known optimal mapping, searching for the single and the multiple lay-out.
the multiple correspondence approaches were always able to find the optimal mapping. Also, limiting the number of candidates reduces the number of rearrangements needed to generate the optimal mapping. Furthermore, for a small number of candidates, the multiple approach seems to converge slightly faster than the single approach, while for many candidates this is the other way around. 3.2. Random cost matrix In a 100 x 100 matrix filled with random numbers uniformly distributed between 0 and 100 the optimal mapping was searched. Again, the positional error was unknown, and a fixed number of candidates was selected. Since we were not able to generate and therefore recognize the true mapping, we used the following stop criterion. At first, the mapping will change rapidly to lower energy mappings. Later, the rate of accepted changes towards higher cost mappings will steadily increase until it approximates the rate of changes to lower cost mappings. As temperature increases or the
minimal mapping is approached, changes to higher energy mappings become more unlikely. A reasonable test for equilibrium is when the ratio of changes to higher and lower cost mappings, as measured over a fixed number of accepted changes, is stable. The results, plotted in Fig. 6, show that again both approaches converge to the optimal mapping. Although the single approach seems to converge somewhat faster than the multiple approach, the latter converges to a lower cost mapping. The reason that the multiple approach generates a lower cost mapping than the single approach is that, in certain cases, it can escape local minima easier. The same situation is found in the single approach between including, and not including the null space in the matching. This is illustrated by the example in Fig. 7. The points A and A' form a true correspondence, but are matched to other points. To become a pair, the correspondences (A, B') and (C, A') have to exchange counterparts (Fig. 7a). When the null space is not ineluded and B and C' are very expensive to match, the chance that this rearrangement takes place is very
Random Cost Matrix multiple approach
Random Cost Matrix single approach
.~
100
lO0
80
80
==
eo
E
E
N
40
2o
00
50
100
rearrangements(xlO0)
150
200
o0
5o
100
1~Jo
rearrangement8(xlO0)
Fig. 6. Mapping cost as a function of the number of accepted rearrangements in a 100 × 100 random matrix, while searching for the optimal single and multiple lay-out.
200
Finding point correspondences using simulated annealing
A'
A'
A'
237
A'
l°
A'
b,
Fig. 7. Illustrating the process of matching the true correspondence (A, A') in the single lay-out (a) and in the multiple lay-out (b).
small. Including the null space allows the point A first to be matched (parked) to the null space and then take over element A' from B (Fig. 7b). This solution needs an additional step, which may account for the longer convergence time. The result, however, is a lower cost mapping. 3.3. Stereo views In this experiment, 100 points were randomly distributed in a 100 x 100 x 50 rectangular space, rotated over 30 ° around the x-axis and projected onto the xy-plane. These points were displaced by adding a normally distributed, zero mean vector. Here, ~rr~s was known, and therefore equation (9) could be used to calculate the pair-costs. Figure 8 shows that smaller positional errors lead to lower cost mappings in both approaches, and that the multiple approach generates lower cost mappings, although there is hardly any difference in convergence time between the single and the multiple approach. However, both approaches need more time to converge than in the two experiments described above. The reason for this is that the pair-costs are all smaller than d.... and thus closer to each other. Therefore, to be able to freeze the mapping, the annealing procedure needs a lower temperature.
3.4. Practical example As a practical example we describe a localization study performed on A431 cells. The cells were labeled against DNA with gold tagged antibodies and cut into sections of about 250nm thickness. In the electron microscope, two views of such a section were recorded, one in untilted position and one at a tilt angle of 30°, both at a magnificiation of 34,600 (Fig. 9a). The gold particles were localized in both views using a detection algorithm for circular objects(3°) (Fig. 9b). The projection sets were registered using eight user indicated, sure correspondences. The angle between the x-axis and the tilt axis was estimated as 100.9° + 0.18, the rotation between the two recoding as 0.8 + 0.05, and the scaling difference between both as 1.001 + 0.00. Next, the pair-costs were calculated based on the squared epipolar distance. The mapping was searched for using the multiple approach with an allowed maximum of two counterparts matched to each point at any time (k = 0, 1, 2). The result of the matching is overlayed on the original recordings in Fig. 9(c), and as two perspective views of the 3D positions of the gold markers in Fig. 9(d). The localization step resulted in 105 markers in the untilted view, and 152 markers in the tilted view. Visual inspection showed, that in the untilted view four
Tilted Point Set
Tilted Point Set
single approach
multiple approach
lOO
8o
lOO
3
-i
'1
'1 , \ ~ ',,
8 -~
eo
\
E
E
20
So
211
i 2o
i 40
60
rearrangements (x100)
i
8o
100
0
0
.....
20
. . . . . . . . . . .
'
40
80
80
rearrangements (x100)
Fig. 8. Mapping cost as a function of the number of accepted rearrangements in matching two views (0° and 30°) of 100 points randomly distributed in a 100 x 100 x 50 block. The points were displaced using 0.15, 0.30, 0.45 and 0.60 for ~r,,o,.
100
238
J . P . P . STARINK and E. BACKER
i!
a. original
b. detected markers
c. corresponding markers
×
d. two views of the marker distribution Fig. 9. Localization of immuno gold labeled DNA. Two recordings, one untilted and one tilted at 30 °, are shown in (a). The detected markers are overlayed in (b). The matching result is shown in (e), two 3D views of the markers in (d).
Finding point correspondences using simulated annealing Table 2. Resultsof the matching procedure in the localization study. The number of multiple matches was restricted to 2, resulting in 5 different correspondences Detected
False
Should be
0-to-1
41
3
36
1-to-0
2
1
2
1-to-1 1-to-2 2-to-I
100 3 11
7 1 2
98 3 13
markers were occluded and two were not present in the tilted view. In the tilted view, 9 were occluded and 36 not present in the untilted view. The simulated annealing loop needed about 9000 rearrangements to reach convergence, as decided by the stop criterion described in section 3.2. The matching results are listed in Table 2 for the five different correspondences. They show that for the single correspondences, the matching procedure missed five and made seven mismatches. It is stressed that these figures depend heavily on both the density of the markers in the projection images and on the defined cost function. Using application specific information may change the performance in a positive way. In this example, the cost function could include grey value cross correlation, the distance along the epipolar lines, and lay-out similarity between neighboring markers.
239
point, the multiple approach performs better when the number of candidates p is smaller than about 35% of the total number of points. In practical applications this is probably always the case. All experiments show that the single approach converges slightly faster than the multiple approach. On a SUN SPARCstation 10, the typical difference in convergence time of the tests described in paragraph 3.2 was less than 2 s. However, the multiple approach generates lower cost, and therefore better mappings. Sometimes it is known that each point in one view can be matched only to a limited number of points in the other view. For example, in two views of a cube, a corner point can occlude only one other point. This restriction is easily included in "the multiple approach, as illustrated in the practical example. The multiple approach reduces to the single approach when the number of allowed counterparts is set to one. Our implementation was written in C and matched the gold particles in the DNA example in less than 7 s. It took about 2 min to localize the markers, and an additional 20 s to compute the cost matrix. So, the method presented in this paper runs very fast in comparison to the time required for the additional operations in the matching procedure. Acknowledgements--We would like to thank Dr Bruno
Humbel for providing the biological materials used in the practical example. This work was partially supported by 'The Netherlands Team for Computer Science Research (SPIN)', project 'Three-Dimensional Image Analysis'.
4. CONCLUSIONS In this paper we have presented a point matching method that iteratively changes an arbitrary initial mapping towards the minimal cost mapping based on a simulated annealing minimization scheme. The rearrangements are subject to lay-out constraints to ensure an unambiguous mapping. We have discussed two practical lay-outs: the single lay-out, in which each point can be matched to one point or can be left unmatched, and the multiple lay-out, in which each point can be matched to any number of points. When a point is occluded by another point in one view, but is visible in the other view, the single approach matches one point and leaves the other unmatched. The multiple approach, however, is able to localize both points. Two ways to reduce the candidate lists were described. In the first, the allowed distance of the candidates to the epipolar line is limited, based on the displacement error of the points. In the second, a fixed number of candidates closest to the epipolar line is chosen. Reducing based on the displacement error is preferred, because the chance of recognizing true null correspondences is higher; when all candidates are further away than the allowed distance, the point is automatically left unmatched. Assigning a fixed number of candidates to a point leaves the point unmatched only when all its candidates match cheaper to other points. The experiments with preset cost matrices show, that when choosing a fixed number of candidates for each
REFERENCES
1. M. J. P. M. Lemmens, A survey of matching techniques, in Comm. V of the ISPRS Congress, (Kyoto, Japan) pp. 1-13 (1988). 2. D. Marr, Vision. Freeman and Company, New York (1979). 3. D. Marr and T. Poggio, A theory of human stereo vision, Proc. R. Soc. Lond. Set. B 204, 301-328 (1979). 4. W. E. L. Grimson, From imayes to surfaces: a computational study of the early human visual system. MIT Press, Cambridge, Massachusetts (1981). 5. L. G. Shapiro and R. M. Haralick, Relational matching, Appl. Optics 26, 1845-1851 (1987). 6. K.L. Boyer and A.C. Kak, Structural stereopsis for 3D vision, IEEE Trans. Pattern Anal. Mach. lntell. 10, 144-166 (1988). 7. L. G. Shapiro and R. M. Haralick, Structural descriptions and inexact graph matching, IEEE Trans. Pattern Anal. Mach. Intell. 3, 504-519 (1981). 8. S. Ullman, The Interpretation of Visual Motion. MIT Press, Cambridge, Massachusetts (1979). 9. Y. Otha and T. Kanade, Stereo by intra- and interscanlinesearch usingdynamicprogramming,IEEE Trans. Pattern Anal. Mach. Intell. 7, 139-154 (1985). 10. S.T. Barnard and W. B. Thompson, Disparity analysis of images, IEEE Trans. Pattern Anal. Mach. Intell. 2, 333-340 (1980). 11. R. A. Hummel and S. W. Zucker, On the foundations of relaxation labeling processes, IEEE Trans. Pattern Anal. Mach. lntell. 5, 267-287 (1983). 12. S. Anilyand A. Federgruen,Simulatedannealingmethods, J. Appl. Probab. 28, 657-666 (1987).
240
J.P.P. STARINK and E. BACKER
13. S.T. Barnard, A stochastic approach to stereo vision, Proc. 5th Int. Conf. Artificial Intell. 1, 676-680 (1986). 14. C. H. Lee and T. S. Huang, Finding point correspondences and determining motion of a rigid object from two weak perspective views, Computer Vision, Graphics, and Image Process. 52, 309-327 (1990). 15. L. Aloimonos and C. M. Brown, Perception of structure from motion, in Proc. IEEE Conf. Comput. Vision and Pattern Recognition pp. 22-26 (1986). 16. T. S. Huang and C. H. Lee, Motion and structure from orthographic projections, IEEE Trans. Pattern Anal. Mach. Intell. 11, 536-540 (1989). 17. B. SabataandJ. K. Aggarwal, Estimation of motion from a pair of range images: a review, CVGIP: Image Understanding 54, 309-324 (1991). 18. H. H. Chen and T. S. Huang, Maximal matching of 3D points for multiple-object motion estimation, Pattern Recoffnition 21, 75-90 (1988). 19. N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087-1092 (1953). 20. S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671-680 (1983). 21. S. Kirkpatrick, Optimization by simulated annealing: quantitative studies, J. Star. Phys. 34, 975 (1984).
22. A. A. El Gamal, L. A. Hemachandra, I. Shperling, and V. K. Wei, Using simulated annealing to design good codes, IEEE Trans. Inform. Theory 33, 116-123, (1987), 23. V. t~ern~, Tbermodynamical approach to the traveling salesman problem: an efficient simulation algorithm, J. Optim. Theory Appl. 45, 41-51 (1985). 24. P. Carnevali, L. Coletti and S. Patarnello, Image processing by simulated annealing, IBM J. Res. Develop. 29, 569-579 (1985). 25. H. L. Tan and S. B. Gelfland, A cost minimization approach to edge detection using simulated annealing, IEEE Trans. Pattern Anal. Machine Intell. 14, 3-18 (1992). 26. B. Hajek, Cooling schedules for optimal annealing, Math. Oper. Res. 13, 311-329 (1988). 27. S. Seman and D. Geman, Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images, IEEE Trans. Pattern Anal. Mach. lntell. 6, 721-741 (1984). 28. B. Gidas, Nonstationary markov chains and convergence of the annealing algorithm, J. Star. Phys. 39, 73-131 (1985). 29. D. Mitra, F. Romeo and A. Sangiovanni-Vincentelli, Convergence and finite time behavior of simulated annealing, Adv. Appl. Probab. 18, 747-771 (1986). 30. J. P. P. Starink and I. T. Young, Localization of circular objects, Pattern Reco#nition Lett. 14, 985-905 (1993).
About the Author--PASCUAL STARINK was born in Delft, The Netherlands in 1964. He received the M.Sc. degree in 1987, and the Ph.D. degree in 1993 from Delft University of Technology. His Ph.D. work was on the analysis of electron microscope images. This was done in close collaboration with the department of molecular cell biology of Utrecht University, The Netherlands. Currently he is a Max Planck fellow with the Max Planck Institute for Biophysical Chemistry in G6ttingen, Germany. His work focuses on image processing and analysis, in particular of scanning probe microscope images.
About the Author--ERIC BACKER was born in Soestdijk, The Netherlands in 1940. He received the M.Sc. degree and the Ph.D. degree from Delft University of Technology, in 1969 and 1978, respectively. Hie was appointed as full professor at Delft University of Technology in 1982, and has been visiting professor at Michigan State University in 1978 and 1982, and at the University of South Florida in 1993. Hie is Managing Editor of Pattern Recognition Letters; he has been the general co-chairman of the llICPR, held in The Hague, The Netherlands. His current research interest includes Pattern Recognition, Image Processing, and Machine Intelligence. He has published over 80 technical papers and is (co)-author of six books.
Pergamon
0031-3203(94) 00087-5
F I N D I N G POINT CORRESPONDENCES USING SIMULATED ANNEALING J. P. PASCUAL STARINKt~ and ERIC BACKERt t Delft University of Technology, Department of Electrical Engineering, P.O. Box 5031, 2628 CD Delft, The Netherlands
(Received 15 September 1993; in revisedform 24 June 1994; receivedfor publication 20 July 1994) Abstract--Identifying corresponding points between two recordings of a point set has always been an important problem in stereo vision applications. We describe this matching problem in terms of cost minimization and present an algorithm to approach the minimal cost mapping using simulated annealing. The algorithm calculates the costs to match all possible point pairs and tries to minimize the sum of the costs of all matched points. Starting from an initial mapping, it uses a random rearrangement scheme to alter the mapping towards the optimal (minimal cost) mapping. Matching Point correspondences Simulated annealing
Search methods
1. INTRODUCTION Stereo vision is concerned with the recording of an object space from slightly different views. Identifying corresponding elements in the views plays a key role. Generally, these relations are not bijective; elements in one set may be absent in the other, and vice versa. This, and the lack of a unifying theory makes the correspondence problem severe and has led to an overwhelming diversity of approaches to solve the problem. ~11To arrive at a structural description of the object space, the problem can be approached in different ways, depending on the representation of the information. From low-level to high-level, these approaches are signal matching, feature matching, and relational matching. In the signal approach, the images are treated as 2D signals consisting of a distribution of grey value intensities. It makes use of the well-explored concepts of signal processing. The elements are usually pixel neighborhoods, which are detected in one image and sought for in the other, for example by grey value correlation or least squares matching. The feature approach was introduced by Marr and coworkers. ~2-4~ They used structureless random dot stereograms to show that the human visual system perceives stereo by detecting elementary tokens rather than matching at a higher, structural level. The elements are local, such as points and edges. They are detected in both images and matched afterwards. In the structural approach, also referred to as relational matching, structures and their relationship to neighboring structures are des-
Stereo vision
cribed by a symbolic representation organized in graphs. 1~'6) Due to the different view points of the recordings, the graphs differ in topology and are therefore matched using inexact graph matchingF I Similarity measures between the detected elements lead to an initial set of possible correspondences. If costs are defined reciprocal to the similarity measures, the minimal cost mapping represents the optimal match. Ullman~8) reformulates his surface model of minimal mapping as a linear programming problem. Also dynamic programming has been used in the minimal cost approach. Otha and Kanade 19~ use dynamic programming for the search both along and across the epipolar lines. In relaxation labeling an element labelling is updated iteratively to achieve a globally consistent interpretation,I~°'t ~) but the final solution depends heavily on the starting point and tends to relax in local optima. "2) Barnard ~t3) bases his (signal) matching approach on simulated annealing to overcome these disadvantages. Point matching belongs basically to the feature approach, but may include flavors of the other two approaches in defining similarity measures.
:~Current address: Max Planck Institute for Biophysical Chemistry, Department of Molecular Biology,P.O. Box 2841, D-37018 Grttingen, Germany. 231
2. MODEL Consider a point set in 3D space. We are given the primary and the secondary projection of this set onto two planes as the points p~, i = 1...Np, and sj, j = 1...N~. Matching is accomplished by selecting a set of correspondences that forms an unambiguous mapping with respect to the desired lay-out. We consider two different lay-out approaches: (i) Single correspondences (L~a ~) (1-to-l, 1-to-0, 0-to-l)
232
J.P.P. STARINK and E. BACKER
Both projection sets probably contain a different number of points, due to occlusion, overprojection, etc. In the single approach, a point Pi is allowed to match either a point sj or the null space, and vice versa. (ii) Multiple correspondences (Z~a,1) (1-to-k, k-to-l, k > 0) In the multiple approach, each point is allowed to match any number of points in the other set. By this, also occluded points can be matched and therefore localized. Reducing the candidate lists is commonly accomplished by setting an upper boundary on the parallax by limiting the allowed distance of the candidates to the epipolar lines. The locations and directions of these lines are determined by the motion parameters. If these parameters are unknown, they can be estimated from a limited number of sure correspondences. For weak perspective and parallel projections, four correspondences of three views are sufficient, la) It is possible to extract the directions and locations of the epipolar lines from four points of two views, ~¢~ but it is impos-
enough so that the system does not get stuck in local minima of some energy function. Originally developed by Kirkpatrick et al., t2°'2 ~ the method has been successfully applied to a variety of hard optimization problems in different fields. Among them, the design of good codes t22~ and attempts to solve the traveling salesman problem. 123)It is also used in general image processing problems,~24) in edge detection,t25~ and stereo matching,t t 3) Here, we are given the finite set ~ of all possible mappings and the energy (cost) function E(J¢). The probability of occurrence of any particular mapping m is proportional to its Boltzman weight, P(m)ocexp(E(m)/T), where T is the temperature. For each state m e J ¢ there is a set N(m) c J//that contains the neighbors of m. In addition, a transition probability matrix R is defined such that R ( m , m ' ) > 0 if and only if m'~N(m). Finally, let T~, T2.... be a sequence of strictly positive numbers such that T~ > T2 > ... and lim T k O. With IX] + = =
k~oo
max (X, 0), the general form of the annealing algorithm is:
Set k = 0 Choose initial mapping X k while Tk # 0 do Choose next state Yk from N(Xk) with probability R(m, m')
f
Set Xk
+~ Set k = k + end.
rk
), Xk
with probability exp F - [ E ( Yk)- E(Xk)] +]
Tk
L
J
otherwise
1
sible to determine the motion parameters. I~s'16~ For a review on perspective projections, see Sabata and Aggarwal.t~ 7~ A cost function c(i, j) is defined by specifying the pair-cost charged when points p~ and sj are matched. This function is of major importance because additional discriminative power between these pair-costs reduces the search time. It should therefore reflect both general and application specific properties, such as grey value neighborhood similarity and relational properties like rigidity.~l sl 2.1. Simulated annealing Simulated annealing is a stochastic optimization algorithm based on the physical analogy of annealing a system of molecules to its ground state. To grow a perfect crystal, one starts at high temperature and then gradually cools the substance, staying as close to equilibrium as is practical. This is directly analogous to what happens in simulated annealing. To bring a synthetic system to equilibrium, the cooling process is simulated using the standard method of Metropolis et al.~191The rate of cooling is controlled by the macroscopic parameter 'temperature', and must be slow
The algorithm produces a discrete-time Markov chain X = {Xklk > 0}. Let J¢* be the set of states in ~ ' at which E attains its minimal value. The idea of the simulated annealing algorithm is to try to achieve lim P [ X k e J / ¢ * ]
= 1
(1)
k~
by letting Tk tend to zero as k leans to infinity. The design of the temperature schedule and the transition distribution R are based on a theorem by Hajek. t26~Before stating this theorem, some additional definitions are necessary: • A state m' is reachable from state m if there is a sequence of states m = m o, ml,.. •, mv = m' for some p > 1 such that ink+ 1eN(mk) for 0 < k < p. • (J¢, N) is irreducible when for any two states m and m', m' is reachable from m. • State m' is reachable at height H from state m if there is a sequence of states m = m o, ml,... ,rap = m' for some p > 1 such that mk÷ 1~N(mk), 0 < k < p, and E(mk) < H for 0 < k < p. • (d/, R, E) is weakly reversible when for any real number H and any pair of states m and m', m is
Finding point correspondences using simulated annealing reachable at height H from m' if m' is reachable at height H from m. • A state m is a local minimum if no state m' with E(m') < E(m) is reachable from m at height E(m). • A state m is a global minimum if no state m' exists for which E(m') < E(m). • The depth of m is + ~ if m is a global minimum, otherwise it is the smallest n u m b e r D _> 0 such that some state m' with E(m') < E(m) can be reached from m at height E(m) + D.
Theorem. Assume that (./¢, R,E) is irreducible and weakly reversible, then the convergence in probability of equation (1) is reached if and only if exp f - d * ' ]
k=1
\
= +oo
(2)
Tk/
where d* is the maximum of depths of all states that are local but not global minima. This result is consistent with the result of other researchers who have established the asymptotic convergence of the simulated annealing algorithm. 127-29) If the temperature schedule assumes the parametric form e
Tk -- - log (k + 1)
(3)
then equation (2), and hence also equation (1), is true if and only ifc _> d*. To limit the n u m b e r of iterations, c should be kept as small as possible. If R is reversible, then, by definition, there is a probability distribution ~ on ~ , such that ~(x)R(x,y)= ~(y)R(y, x) for all x, y e ~ . Necessarily, ~ is the equilibrium distribution for R. A simple transition distribution for which this is true, is given by:
R(x,y) =
233
Table 1. The possible matching combinations between any two correspondences in the single lay-out (S) and the multiple lay-out (M) 0
0 Po p~ Psm p,~
SO
Sss
Ssm
Sm
--
--
S,M
M
--
-S,M M --
S,M S,M M M
S,M S,M M M
M M M M
M M M --
Breaking a correspondence decreases the mapping energy, while creating a new correspondence increases the mapping energy with the corresponding pair-cost. If the change in mapping energy AE is negative, the rearrangement is accepted, otherwise it is accepted with a chance according to the Boltzman probability distribution. The rearrangements are illustrated in Figs 1-4. The configurations left of the arrow represent the situation before, the configuration right of the arrow the situation after the rearrangement. The null space is drawn as a circle on top of the sets. The candidates are grey, p in the left set and s in the right set, If a candidate is the null space, it is indicated with 0. If a point is matched to the null space, it is given the subscript 0. If it is matched to two or more points, the subscript is m (multiple). In case of a single correspondence, there is a two-character subscript, of which the first is s (single). The second character reflects the n u m b e r of elements that its counterpart is matched to. This n u m b e r is either o n e - - t h e candidate itself--or more, represented by s and m, respectively. Table 1 lists the allowed configurations of any two candidates in the two layouts.
if y e N ( x ) otherwise.
(4)
2.2. The matching algorithm To generate an arbitrary, initial mapping, we have used a simple best-first scheme: repeat selecting the minimal cost correspondence from the unmatched points until at least one set is empty. Remaining unmatched points are matched to the null space. Value c = d* is set to the maximal pair-cost in the initial mapping. Although this value is surely too high and results in longer convergence times, it proves to be an adequate guess. Let a rearrangement be a change in the mapping such that the previously unmatched points p~ and st, which are randomized, become a new correspondence, and the resulting mapping m' is part of N(m). The points that are matched to pi and sj are denoted by ~ and gj, respectively. When a multiply matched point is involved, the correspondence to be possibly changed is randomized. When the number of possible rearrangements exceeds one (IN(m) l > 1), a random one is picked according to equation (4).
2.2.1. Single correspondences-- &P~. When a mapping should contain only one-to-one correspondences and unmatched points, there are six possible configurations: 0 - s s, Ps - 0, Po - So, Po - s~, ps - s o and p~ - s s. All can be rearranged in only one way, as illustrated in Fig. 1. The changes in mapping energy due to these rearrangements are given by: AE = - 2c(gi, sj) + c(gj, 0) + c(0, sj)
(5a)
AE = - 2c(pi, Pi) + c(pi, 0) + c(0, 0i)
(5b)
AE = - c(pi, 0) - c(0, sj) + 2c(pi, s) )
(5c)
A E = -c(pi, O)-c(gj, sj) + c(pi, si) + c(gj, O)
(5d)
AE = - c(p/, Pl) - c(0, st) + c(pi, st) + c(0, Pi)
(5e)
AE = - c(pi, p~) - c(gj, sj) + c(p~,sj) + c(~, g2).
(5f)
2.2.2. Multiple correspondences--ZP". To describe the rearrangements in a multiple correspondence mapping, the possible configurations are divided into three sets. The first set consists of the single correspondences, which are described in Section 2.2.1, augmented by four additional ones. These four correspondences play
234
J.P.P. STARINK and E. BACKER
(I,
0 -
,2 - -
m
__
- -
m
m
Pss - 0
(a)
Sss
.~
(b)
C'
PO
" Sss
(d)
Pss
-
So
-7 PO
(e)
-
-7
SO
O
C,
Pss -
Sss
(c)
C
'D
(f)
Fig. 1. The rearrangements for the six possible configurations of the candidates p and s in a single correspondence lay-out.
A E = - c ( g j , sj) + c(pi, sj) AE
=
c(pi, [h) - c(gj, sj) + c(p i, sj) + c(g~,0)
-
AE = - c(p~, Pi) + c(pi, s~) AE = - c(p~, [h) - c(gj, s~) + 2c(E, %). PO
O
-
Sss
(a)
Pss " So
(b)
O
(c)
Pss
-
Sss
(7g) (7h) (7i)
The third set contains six configurations involving multiple-multiple matched points (subscript m). Each can be rearranged in only one way, see Fig. 4. The energy changes are given by: AE
Pss - Sss
(7f)
c(p~, 0) + c(pi, sj)
(8a)
AE = - c(0, sj) + c( E, sj)
(8b)
AE = - 2c(p i, Pi) + c(E, sj) + c(0, Pi)
(8c)
A E = - 2c(gj, si) + c(pi, s~) + c(gj, 0)
(8d)
AE = - c(pi, Pi) + c(pi, s j)
(8e)
AE = - c(gj, sj) + c(pi, sj).
(Sf)
=
-
(d)
Fig. 2. Four rearrangements for the three possible configurations of singly matched candidates p and s resulting in multiple correspondences.
an essential role in the construction of multiple correspondences. The rearrangements of these four are shown in Fig. 2. The related changes in mapping energy are: AE
=
-
c(pi, 0) + c(pi, sj)
(6a)
AE
=
-
c(0, s j) + c(p i, sj)
(6b)
AE = - 2c(pi, Pi) + c(pi, s j) + c(0, Pi, 0)
(6c)
A E = - 2c(gj, sj) + c(pi, s j) + c(gj, 0).
(6d)
The second set holds nine configurations involving the single-multiple matched candidates (subscript sin). Their rearrangements are illustrated in Fig. 3. The changes in mapping energy are:
A E = - c ( p , , pi) + c(O, p~)
(7a)
AE = -c(p~, p~) + c(pl,0)
(7b)
AE = -- c(pi, 0) -- c(gj, sj) + 2c(p~, s j)
(7c)
AE = - c(p~, Pi) - c(0, sj) 6+ 2c(pi, sj)
(7d)
AE = - c ( p i , Pi) - c(gb,s j) + c(E, s i) + c(0, Pi)
(7e)
3. EXPERIMENTALRESULTSAND DISCUSSION In this section, three series of experiments are described to determine the rate of convergence of the matching procedure. In the first series, the mapping is searched within preset cost matrices, in the second within r a n d o m cost matrices. In the third series, the cost matrices are calculated from stereo views of a generated point set. An example from the field of cell biology is included; the localization of DNA. The cost function in the experiments is based on the epipolar constraint only. To obtain the pair-costs, consider the point correspondence (Pi, sj) for i > 0 and j > 0. If the epipolar line of point p/ is 2~ and the positional error of the points is Crpos, then the positional error of ;t~ is also %os. If we assume that the error is normally distributed, then the positional error of sj relative to J,~is ,~-2Crpos. With a probability of 95%, the true correspondence of Pi is located within distance dmax = 1.96 ~/2~po~ ~ 2.8 apo~ of the epipolar line. The cost function is now defined as the square of the distance
235
Finding point correspondences using simulated annealing
O g, C,/~__ C,~D ,D ,D
0o
mT_-_ _
0 - sm
(a)
__
Psm -
m m,~ mxm _
--
)--
__
)__
0
(b)
O C, C~'~ _
~
=×= [] roman'
O O
(d)
O C~
O~,
Psm - Sss
Pss - Ssm
(g)
Psm
(e)
C,
Sss
-
c,,~
(c)
OO
lli~~ _ _
_
Psrn - SO
PO - Ssm
Pss - Ssm
(f)
,D
~
m
m~m m~mX
(h)
Psm
O
- Ssm
O O
(i)
Fig. 3. The nine allowed rearrangements for the seven single/multiple correspondences of the candidates p and s in the multiple correspondence lay-out.
C,/"
,D ,D 'ID~o
_-7
_,
PO - Sm
©
Q O
C~~ • ;, _711 _
_ (a)
Q C C,/D
Pm-
so
:
,D
(b)
,D 0
C,
Srn
(c)
O Q
Q C,
Pss
-
m/m - 'l' m m --
Pm
-
Sss
(d)
Psm
7"
- Sm
(e)
Pm - Ssm
(f)
Fig. 4. The six rearrangements for the configurations of the candidates p and s involding multiply matched points in the multiple correspondence lay-out.
d(s/, 2i) of candidate sj to 2i when it is smaller than d . . . . and +o0 otherwise. For i = 0 or j = 0 , it is set to 1.2 d2ax, so:
d(sj, ,~.i)2, c(i, j) =
+ oo, 1.2 d .2. . .
i > 0, j > 0, d($j,,~i) < dmax i > 0, j > 0, d($j, ~i) ~" dmax
c of the temperature schedule was set to the difference between the m i n i m u m and the maximum of the paircosts in the initial mapping. 3.1. Preset cost matrix
(9)
i = 0 o r j = 0.
If dmax is unknown, the candidate lists may be reduced by selecting t h e / t candidates closest to the epipolar line and assigning + oo to all other pair-costs. For i = 0 or j = O, c(i, j) is set to the most expensive of the /~ candidates. In all experiments, the initial mapping was constructed by matching all the points in both sets to the null space, the maximal cost mapping. The parameter
A cost matrix of 100 x 100 elements was generated such that the true optimal mapping was generated by the best-first search algorithm. Since the optimal mapping equaled the true mapping, it could be recognized when randomized. The positional error was u n k n o w n and could not be used to calculate the pair-costs, therefore a fixed n u m b e r of candidates for each point was selected: 10, 30, 60 and 100 candidates, respectively. The results, plotted in Fig. 5, show that, although not always in the plotted region, both the single and
236
J.P.P. STARINK and E. BACKER Preset Cost Matrix single approach
Preset Cost Matrix multipleapproach
loo
100
"fl
8 == -~ 80 N
20
~,,{,
E 40
:i. '('.~I^
......... :,
20
.
:y, ,----~
rearrangements(xlO0)
l 150
O
200
i
i
r
50
100
150
200
rearrangements (x1130)
Fig. 5. Mapping cost as a function of the number ofaccepted rearrangements in a 100 x 100matrix holding a known optimal mapping, searching for the single and the multiple lay-out.
the multiple correspondence approaches were always able to find the optimal mapping. Also, limiting the number of candidates reduces the number of rearrangements needed to generate the optimal mapping. Furthermore, for a small number of candidates, the multiple approach seems to converge slightly faster than the single approach, while for many candidates this is the other way around. 3.2. Random cost matrix In a 100 x 100 matrix filled with random numbers uniformly distributed between 0 and 100 the optimal mapping was searched. Again, the positional error was unknown, and a fixed number of candidates was selected. Since we were not able to generate and therefore recognize the true mapping, we used the following stop criterion. At first, the mapping will change rapidly to lower energy mappings. Later, the rate of accepted changes towards higher cost mappings will steadily increase until it approximates the rate of changes to lower cost mappings. As temperature increases or the
minimal mapping is approached, changes to higher energy mappings become more unlikely. A reasonable test for equilibrium is when the ratio of changes to higher and lower cost mappings, as measured over a fixed number of accepted changes, is stable. The results, plotted in Fig. 6, show that again both approaches converge to the optimal mapping. Although the single approach seems to converge somewhat faster than the multiple approach, the latter converges to a lower cost mapping. The reason that the multiple approach generates a lower cost mapping than the single approach is that, in certain cases, it can escape local minima easier. The same situation is found in the single approach between including, and not including the null space in the matching. This is illustrated by the example in Fig. 7. The points A and A' form a true correspondence, but are matched to other points. To become a pair, the correspondences (A, B') and (C, A') have to exchange counterparts (Fig. 7a). When the null space is not ineluded and B and C' are very expensive to match, the chance that this rearrangement takes place is very
Random Cost Matrix multiple approach
Random Cost Matrix single approach
.~
100
lO0
80
80
==
eo
E
E
N
40
2o
00
50
100
rearrangements(xlO0)
150
200
o0
5o
100
1~Jo
rearrangement8(xlO0)
Fig. 6. Mapping cost as a function of the number of accepted rearrangements in a 100 × 100 random matrix, while searching for the optimal single and multiple lay-out.
200
Finding point correspondences using simulated annealing
A'
A'
A'
237
A'
l°
A'
b,
Fig. 7. Illustrating the process of matching the true correspondence (A, A') in the single lay-out (a) and in the multiple lay-out (b).
small. Including the null space allows the point A first to be matched (parked) to the null space and then take over element A' from B (Fig. 7b). This solution needs an additional step, which may account for the longer convergence time. The result, however, is a lower cost mapping. 3.3. Stereo views In this experiment, 100 points were randomly distributed in a 100 x 100 x 50 rectangular space, rotated over 30 ° around the x-axis and projected onto the xy-plane. These points were displaced by adding a normally distributed, zero mean vector. Here, ~rr~s was known, and therefore equation (9) could be used to calculate the pair-costs. Figure 8 shows that smaller positional errors lead to lower cost mappings in both approaches, and that the multiple approach generates lower cost mappings, although there is hardly any difference in convergence time between the single and the multiple approach. However, both approaches need more time to converge than in the two experiments described above. The reason for this is that the pair-costs are all smaller than d.... and thus closer to each other. Therefore, to be able to freeze the mapping, the annealing procedure needs a lower temperature.
3.4. Practical example As a practical example we describe a localization study performed on A431 cells. The cells were labeled against DNA with gold tagged antibodies and cut into sections of about 250nm thickness. In the electron microscope, two views of such a section were recorded, one in untilted position and one at a tilt angle of 30°, both at a magnificiation of 34,600 (Fig. 9a). The gold particles were localized in both views using a detection algorithm for circular objects(3°) (Fig. 9b). The projection sets were registered using eight user indicated, sure correspondences. The angle between the x-axis and the tilt axis was estimated as 100.9° + 0.18, the rotation between the two recoding as 0.8 + 0.05, and the scaling difference between both as 1.001 + 0.00. Next, the pair-costs were calculated based on the squared epipolar distance. The mapping was searched for using the multiple approach with an allowed maximum of two counterparts matched to each point at any time (k = 0, 1, 2). The result of the matching is overlayed on the original recordings in Fig. 9(c), and as two perspective views of the 3D positions of the gold markers in Fig. 9(d). The localization step resulted in 105 markers in the untilted view, and 152 markers in the tilted view. Visual inspection showed, that in the untilted view four
Tilted Point Set
Tilted Point Set
single approach
multiple approach
lOO
8o
lOO
3
-i
'1
'1 , \ ~ ',,
8 -~
eo
\
E
E
20
So
211
i 2o
i 40
60
rearrangements (x100)
i
8o
100
0
0
.....
20
. . . . . . . . . . .
'
40
80
80
rearrangements (x100)
Fig. 8. Mapping cost as a function of the number of accepted rearrangements in matching two views (0° and 30°) of 100 points randomly distributed in a 100 x 100 x 50 block. The points were displaced using 0.15, 0.30, 0.45 and 0.60 for ~r,,o,.
100
238
J . P . P . STARINK and E. BACKER
i!
a. original
b. detected markers
c. corresponding markers
×
d. two views of the marker distribution Fig. 9. Localization of immuno gold labeled DNA. Two recordings, one untilted and one tilted at 30 °, are shown in (a). The detected markers are overlayed in (b). The matching result is shown in (e), two 3D views of the markers in (d).
Finding point correspondences using simulated annealing Table 2. Resultsof the matching procedure in the localization study. The number of multiple matches was restricted to 2, resulting in 5 different correspondences Detected
False
Should be
0-to-1
41
3
36
1-to-0
2
1
2
1-to-1 1-to-2 2-to-I
100 3 11
7 1 2
98 3 13
markers were occluded and two were not present in the tilted view. In the tilted view, 9 were occluded and 36 not present in the untilted view. The simulated annealing loop needed about 9000 rearrangements to reach convergence, as decided by the stop criterion described in section 3.2. The matching results are listed in Table 2 for the five different correspondences. They show that for the single correspondences, the matching procedure missed five and made seven mismatches. It is stressed that these figures depend heavily on both the density of the markers in the projection images and on the defined cost function. Using application specific information may change the performance in a positive way. In this example, the cost function could include grey value cross correlation, the distance along the epipolar lines, and lay-out similarity between neighboring markers.
239
point, the multiple approach performs better when the number of candidates p is smaller than about 35% of the total number of points. In practical applications this is probably always the case. All experiments show that the single approach converges slightly faster than the multiple approach. On a SUN SPARCstation 10, the typical difference in convergence time of the tests described in paragraph 3.2 was less than 2 s. However, the multiple approach generates lower cost, and therefore better mappings. Sometimes it is known that each point in one view can be matched only to a limited number of points in the other view. For example, in two views of a cube, a corner point can occlude only one other point. This restriction is easily included in "the multiple approach, as illustrated in the practical example. The multiple approach reduces to the single approach when the number of allowed counterparts is set to one. Our implementation was written in C and matched the gold particles in the DNA example in less than 7 s. It took about 2 min to localize the markers, and an additional 20 s to compute the cost matrix. So, the method presented in this paper runs very fast in comparison to the time required for the additional operations in the matching procedure. Acknowledgements--We would like to thank Dr Bruno
Humbel for providing the biological materials used in the practical example. This work was partially supported by 'The Netherlands Team for Computer Science Research (SPIN)', project 'Three-Dimensional Image Analysis'.
4. CONCLUSIONS In this paper we have presented a point matching method that iteratively changes an arbitrary initial mapping towards the minimal cost mapping based on a simulated annealing minimization scheme. The rearrangements are subject to lay-out constraints to ensure an unambiguous mapping. We have discussed two practical lay-outs: the single lay-out, in which each point can be matched to one point or can be left unmatched, and the multiple lay-out, in which each point can be matched to any number of points. When a point is occluded by another point in one view, but is visible in the other view, the single approach matches one point and leaves the other unmatched. The multiple approach, however, is able to localize both points. Two ways to reduce the candidate lists were described. In the first, the allowed distance of the candidates to the epipolar line is limited, based on the displacement error of the points. In the second, a fixed number of candidates closest to the epipolar line is chosen. Reducing based on the displacement error is preferred, because the chance of recognizing true null correspondences is higher; when all candidates are further away than the allowed distance, the point is automatically left unmatched. Assigning a fixed number of candidates to a point leaves the point unmatched only when all its candidates match cheaper to other points. The experiments with preset cost matrices show, that when choosing a fixed number of candidates for each
REFERENCES
1. M. J. P. M. Lemmens, A survey of matching techniques, in Comm. V of the ISPRS Congress, (Kyoto, Japan) pp. 1-13 (1988). 2. D. Marr, Vision. Freeman and Company, New York (1979). 3. D. Marr and T. Poggio, A theory of human stereo vision, Proc. R. Soc. Lond. Set. B 204, 301-328 (1979). 4. W. E. L. Grimson, From imayes to surfaces: a computational study of the early human visual system. MIT Press, Cambridge, Massachusetts (1981). 5. L. G. Shapiro and R. M. Haralick, Relational matching, Appl. Optics 26, 1845-1851 (1987). 6. K.L. Boyer and A.C. Kak, Structural stereopsis for 3D vision, IEEE Trans. Pattern Anal. Mach. lntell. 10, 144-166 (1988). 7. L. G. Shapiro and R. M. Haralick, Structural descriptions and inexact graph matching, IEEE Trans. Pattern Anal. Mach. Intell. 3, 504-519 (1981). 8. S. Ullman, The Interpretation of Visual Motion. MIT Press, Cambridge, Massachusetts (1979). 9. Y. Otha and T. Kanade, Stereo by intra- and interscanlinesearch usingdynamicprogramming,IEEE Trans. Pattern Anal. Mach. Intell. 7, 139-154 (1985). 10. S.T. Barnard and W. B. Thompson, Disparity analysis of images, IEEE Trans. Pattern Anal. Mach. Intell. 2, 333-340 (1980). 11. R. A. Hummel and S. W. Zucker, On the foundations of relaxation labeling processes, IEEE Trans. Pattern Anal. Mach. lntell. 5, 267-287 (1983). 12. S. Anilyand A. Federgruen,Simulatedannealingmethods, J. Appl. Probab. 28, 657-666 (1987).
240
J.P.P. STARINK and E. BACKER
13. S.T. Barnard, A stochastic approach to stereo vision, Proc. 5th Int. Conf. Artificial Intell. 1, 676-680 (1986). 14. C. H. Lee and T. S. Huang, Finding point correspondences and determining motion of a rigid object from two weak perspective views, Computer Vision, Graphics, and Image Process. 52, 309-327 (1990). 15. L. Aloimonos and C. M. Brown, Perception of structure from motion, in Proc. IEEE Conf. Comput. Vision and Pattern Recognition pp. 22-26 (1986). 16. T. S. Huang and C. H. Lee, Motion and structure from orthographic projections, IEEE Trans. Pattern Anal. Mach. Intell. 11, 536-540 (1989). 17. B. SabataandJ. K. Aggarwal, Estimation of motion from a pair of range images: a review, CVGIP: Image Understanding 54, 309-324 (1991). 18. H. H. Chen and T. S. Huang, Maximal matching of 3D points for multiple-object motion estimation, Pattern Recoffnition 21, 75-90 (1988). 19. N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087-1092 (1953). 20. S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671-680 (1983). 21. S. Kirkpatrick, Optimization by simulated annealing: quantitative studies, J. Star. Phys. 34, 975 (1984).
22. A. A. El Gamal, L. A. Hemachandra, I. Shperling, and V. K. Wei, Using simulated annealing to design good codes, IEEE Trans. Inform. Theory 33, 116-123, (1987), 23. V. t~ern~, Tbermodynamical approach to the traveling salesman problem: an efficient simulation algorithm, J. Optim. Theory Appl. 45, 41-51 (1985). 24. P. Carnevali, L. Coletti and S. Patarnello, Image processing by simulated annealing, IBM J. Res. Develop. 29, 569-579 (1985). 25. H. L. Tan and S. B. Gelfland, A cost minimization approach to edge detection using simulated annealing, IEEE Trans. Pattern Anal. Machine Intell. 14, 3-18 (1992). 26. B. Hajek, Cooling schedules for optimal annealing, Math. Oper. Res. 13, 311-329 (1988). 27. S. Seman and D. Geman, Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images, IEEE Trans. Pattern Anal. Mach. lntell. 6, 721-741 (1984). 28. B. Gidas, Nonstationary markov chains and convergence of the annealing algorithm, J. Star. Phys. 39, 73-131 (1985). 29. D. Mitra, F. Romeo and A. Sangiovanni-Vincentelli, Convergence and finite time behavior of simulated annealing, Adv. Appl. Probab. 18, 747-771 (1986). 30. J. P. P. Starink and I. T. Young, Localization of circular objects, Pattern Reco#nition Lett. 14, 985-905 (1993).
About the Author--PASCUAL STARINK was born in Delft, The Netherlands in 1964. He received the M.Sc. degree in 1987, and the Ph.D. degree in 1993 from Delft University of Technology. His Ph.D. work was on the analysis of electron microscope images. This was done in close collaboration with the department of molecular cell biology of Utrecht University, The Netherlands. Currently he is a Max Planck fellow with the Max Planck Institute for Biophysical Chemistry in G6ttingen, Germany. His work focuses on image processing and analysis, in particular of scanning probe microscope images.
About the Author--ERIC BACKER was born in Soestdijk, The Netherlands in 1940. He received the M.Sc. degree and the Ph.D. degree from Delft University of Technology, in 1969 and 1978, respectively. Hie was appointed as full professor at Delft University of Technology in 1982, and has been visiting professor at Michigan State University in 1978 and 1982, and at the University of South Florida in 1993. Hie is Managing Editor of Pattern Recognition Letters; he has been the general co-chairman of the llICPR, held in The Hague, The Netherlands. His current research interest includes Pattern Recognition, Image Processing, and Machine Intelligence. He has published over 80 technical papers and is (co)-author of six books.