- Email: [email protected]
Optimal matching of convex polygons
Pattern Recognition Letters 9 (1989) 327 334
June 1989
North-Holland
Optimal matching of convex polygons Pedro COX, Henri M A I T R E Laboraloire Image, Ecole Natianale Sup~;rieure des T~;h;communicalions, 46 rue Barrault, Paris 75634 Cedex 13, France
Michel M I N O U X Lahoratoire MASI, University; Paris VI. Place Jussieu, Paris 75005, France
Celso RIBEIRO* Department qf Eh,ctrical Engineering, Catholic University ~[' Rio de Janeiro, G~h'ea, Caixa Postal 38063, Rio de Janeiro 22452, Brazil
Received 20 April 1988 Revised 24 November 1988
Abstract: We describe in this paper the problem of matching two series of observations of a 2-D contour, in which each sampling of the contour has different hidden parts. This problem is formulated as the minimization of the distance between the convex hulls of the two samplings and a solution method, based on nonlinear programming techniques, is presented. Computational results are reported to illustrate the practical applicability of the approach. Key words." Contour matching, contour fitting, image matching, pattern recognition, 2-D recognition.
1. Introduction In the domain of robot vision, the use of laser range-finders to acquire a description of the environment is becoming more and more popular; see e.g. (Faugeras and Hebert, 1986). Similar systems are also of increasing use for quality control, 3-D data basis formation and CAD/CAM; see e.g. (Schmitt et al., 1985). For these applications, the object to be measured is usually placed in a structured light beam and the coordinates of the lit points of the surface are measured by triangulation. One of such data acquisition systems is described in details by Lelandais (1984) and Brunet (1987). When only one laser source is used, the surface measurement is often incomplete because of the existence of hidden parts. For this reason, the experi-
ment has to be repeated with a second laser beam in a different position. As a matter of fact, matching these two sets of measurements arises as an important stage in the process of providing a single, more extensive description of the object. This matching problem is addressed in this paper in a simplified form: we limit ourselves to the two-dimensional case, where the observation of the object surface is reduced to one planar closed contour line. Nevertheless, our approach is currently being extended to the three-dimensional case and will be dealt with in a forthcoming paper by Bloch-Boulanger et al. (1988). Besides structural methods (Tsai and Yu, 1985) and contour coding techniques (Pavlidis, 1978), which can not be used here because they give no information about the relative position of the two
* The work of the last author was sponsored by FINEP (research contract grant number 4.3.86.0689-00), CNPq (research contract grant number 302281-85/EE), IBM Brasil and the French Ministry of Foreign Relations. 0167-8655,'89/$3.50 (~ 1989, Elsevier Science Publishers B.V. (North-Holland)
327
Volume 9, Number 5
PATTERN RECOGNITION LETTERS
shapes, most of the existing work rely on an estimation of the Euclidean distance between the shapes, called 'cross correlation', 'template matching' (Pratt, 1978) or 'Hough transform' (Sklansky, 1978). These methods require the computation of a distance measure of the two shapes for each relative position. To reduce the computational cost, 'chamfer' techniques, which precompute the distance of each pixel of the field to one of the shapes, were proposed; see e.g. (Barrow et al., 1977). Later, Borgefors (1986) introduced efficient approximations of the Euclidean distance to compute the chamfer. Burr (1981) introduced a local measure of the distance between the two shapes, based on a 'weighted' Euclidean distance, providing a perfect matching after a rubber-like distortion of the plane. Widrow (1973) and Kass et al. (1987) have proposed mechanical models to express the distances between contours. Unfortunately, for all the above mentioned techniques, no efficient strategy guaranteed to converge to a minimum distance solution is available. Several techniques were proposed to match shapes when the correspondence between the points of the two shapes is known beforehand: an exact technique in 2-D and an interactive one in 3-D (Lin et al., 1986), an exact one in 3-D using quaternions (Faugeras and Hebert, 1983) and another exact one in 3-D based on the singular value decomposition of a 3 x 3 matrix (Arun et al., 1987). But, in many cases, finding the appropriate correspondence is by far the most difficult part of the problem. The matching technique described in the next sections does not rely on such a restrictive assumption.
2. Matching two observations of a 2-D contour: problem statement We consider a 2-D contour (not necessarily convex) from which we extract discretized samples by the laser range-finder system mentioned in the introduction. Depending on various parameters of the data acquisition system (such as the angle of observation and the position of the rotation axis), the hidden parts are not the same for each sampling of the contour: it may be the case that some point oh328
June 1989
tained in a first sample corresponds to a hidden point for a second sample. We study from now on the problem of how to exploit two distinct samples of the same contour in order to get a more extensive representation of the whole contour (i.e., with fewer non-observed points). Our approach proceeds as follows: (i) For each sample, we construct the convex hull of the points in the sample and we select the subset of points which are the extreme points of the convex hull (the reader is referred to (Preparata and Shamos, 1985) for a survey of algorithms for convex hull calculations), (ii) Once the convex hulls of the two samples have been determined, we solve the problem of matching in an optimal way (according to a criterion to be defined below) the two convex polygons whose extreme points are those selected in (i). In the remainder, we concentrate on the design of an efficient way of solving part (ii).
3. Optimal matching of two convex polygons The two convex polygons P l and P2 are given by specifying the following data: (i) N~ (resp. N2), the number of extreme points of polygon P1 (resp. P2); (ii) the coordinates a s, bj~ (resp. a~, b~) of the jth extreme point of polygon P~ (resp. P2), for j = 1..... N 1 (resp. j = 1..... N2); without loss of generality the origin of the axes is supposed to coincide with the center of gravity of polygon P2; and (iii) the numbering of the extreme points is supposed to be such that for any j~{1 ..... N1} (resp. j 6 { 1..... N21 ), if we define
j,
IJ + 1 ,
ifj
l, ,
ifj = N1 (resp. j = N2),
then the segment joining the extreme points numberedj andj' is an edge of polygon P~ (resp. P2)In order to have a measure of how well the two polygons match, we shall need an objective function. This proximity measure should be designed in such a way as to be the smallest, the closest the two boundaries of the polygons are to each other. For any point with coordinates (u,v), we define its distance to polygon PI as
Volume 9, Number 5
PATTERN R E C O G N I T I O N LETTERS
{ ll(u,v)- Projp,(u,v)l I, @t(u, v) =
otherwise,
where Proje~ (u, v) denotes the point of P~ achieving minimum euclidean distance to (u,v) and int(P0 denotes the interior of P~. The distance to polygon P2 is defined in a similar way. With this notation, the proximity measure corresponding to the initial relative position of the two polygons will be taken as NI =
N2
2-2 j=l
4. Minimizing the proximity measure
if (u, v) ¢ int(Pt), 0,
D
June 1989
, 1 b l,
°P2~aj,
J)-{-
2 ,j=l
6z ta2 bZ~
P1 ~" J" J*"
Note that very small values of D are attained only if, for each of the polygons, each extreme point is either in the interior of the other polygon, or very close to it. In order to get such a desirable situation, we see that we have to minimize the proximity measure over all possible relative positions of the two polygons. To generate all possible relative positions, it is enough to apply to one of the polygons (P2 in this case) all possible displacements on the plane, an arbitrary displacement being specified by (i) a translation vector (x,y): and (ii) a rotation angle 0 around the origin of the axes.
Applying a displacement with parameters (x,y, O) to any point of polygon P2 with coordinates (:
Fc°sO
If
ksin 0
-sinO](~) cos
oJ\IU"
In a natural way, we will denote by T(x,y,O,P 2) the polygon resulting from the application of a displacement with parameters (x,y,O) to each of the extreme points of P2. With this notation, looking for values of x, y and 0 resulting in a minimum proximity measure can be reduced to the (unconstrained) minimization of the objectivc function (in three continuous variables only): N1
F(x, .v,O) = ~, (~T(x,y,O, -2 P2)\~aj~, b~ jl j= 1 N2
+ ~ a~ r(x,y,O,a~,b~). j=l
The objective function F(x,y,O) has a number of properties which have to be taken into account in the design of an algorithm for the minimization problem. In particular: (i) F is continuous differentiable with respect to all three variables x, y, and 0: (ii) F is convex with respect to x and y, for fixed O: and (iii) F is not convex with respect to the variable O, for fixed x and y. In view of the nonconvexity with respect to the variable O, we propose here to decompose the minimization by considering separately the variables corresponding to translation (x and y) and the variable corresponding to rotation (0).
4.I, Minimizing the ohjective fimction for lixed 0 The objective function F(x, y, O) being convex and differentiable for 0 fixed to some value [), many unconstrained optimization methods can be applied. In our experiments we tried the steepest descent algorithm (SD) and the conjugate gradient algorithm (CG) (see e.g. (Luenberger, 1973)) and observed no significant differences in the behavior of the two: in each case the number of iterations required to achieve sufficient accuracy was consistently very low (typically less than five iterations). This can be explained by the fact that F, as a sum of distance functions, is very well-conditioned. It is well known that both SD and CG algorithms require, at each step, the computation of: (i) the value of the function F a t the current point; (ii) the gradient g = (g.~,gy) of F at the current point: and (iii) the minimum of F along the chosen descent direction d = (d,., d~.).
(a) Computing the current value o f F at (x,y).[br 0 = OJixed Computing the projection of any point (u, v) on polygon PI or deciding that it lies in cony[P1] can be done in time O(N0. Similarly, for projecting on P2, this computation requires O(N2) elementary operations. Therefore, the total complexity of computing the value of F a t any point (x,y,O) would be 329
Volume 9, N u m b e r 5
PATTERN RECOGNITION LETTERS
O(N~ • Nz). However, it should be pointed out that this computation can be done in time complexity O(N1 + Nz), as suggested by Aggarwal (1988). As a result of this computation, we get for every 1 vertex (a j,1 bj) of P1 the coordinates of its projection on conv[P2], which will be denoted by (~/J,/~J), and for every vertex (a],by) of P2 the coordinates (~/~,/~2) of its projection on conv[P1]. Note that ~/), /~), ~/2, a n d / ~ implicitly depend on x, y and 0. With this notation, the value of F(x, y, O) can be rewritten as
r(x,y,O)
N1
=
[(a) -
+ (b) - t;))
j=l N2
+ E [ ( a ] - - t i ] ) 2 + ( b 2 - @ ) 2 3. j=l
(b) Computing the gradient o f F at (x,y) for 0 = fixed The gradient of a sum of differentiable functions being the (vector) sum of their gradients, the gradient of F at (x,y) for 0 = 0 fixed can be easily computed with linear complexity O(N1 + N2) as N1
V F(x,y,O) = 2 • (a) - aj,b~-I 1 _ ~)) j:l N2
+22(4
-
a),b)
-
j=l
(c) Minimizing the objective function along the descent direction (dx, dy) There are many possible methods for minimizing a given function of a single variable along a (descent) direction. One may use, for instance, methods based on interval reduction (see e.g. (Minoux, 1986)). However, when dealing with a convex differentiable function such that computing the gradient is computationally cheap once the function value has already been computed, the so-called secant method features good global and asymptotic convergence properties. Each step of this secant method requires only O(N 1 + N2) time.
4.2. Finding the best value./or 0 If, for each value 0 of 0, we apply the algorithm outlined in (a), (b) and (c) above, we get optimal values for x and y which are implicit functions of 0. 330
June 1989
We denote these implicit functions respectively by )~(0) and .9(0). Thus, what we want is to find 0* minimizing
z(O) = F(£(O), y(O), O) over all possible angles 0. However, due to the nonconvexity of F(x,y,O), we cannot expect z(O) to be convex. Therefore, neither gradient methods nor interval reduction methods are guaranteed to produce a global minimum (because of non-unimodality). In view of this, the following technique has been used. We assume that a good approximation 0o to the global minimum 0* can be computed (subsection 4.3 below describes a way of getting such an approximate solution). We then define a discretization step eo (typically, eo = I degree) and look for a local minimum in the vicinity of 0 o by the following simple method:
Procedure Approx (0o, ~ ) Step 0: Set 0 c = 0 o, 01eft : 0 0 - - g O , and 0right : 00 + to. Compute Z(0o), z(01m), and Z(0right ). Step 1: Set 0m = 0~eft, if Z(01m)~< Z(0right); 0m = 0right, otherwise. If Z(0m)~> z(Oc) go to Step 4. If 0 m : 01eft g o to Step 2, if 0m = 0right go to Step 3. Step 2: Set 0right = 0c, 0 c : 01eft , and 01m = 0¢ eo. Compute z(01m) and go back to Step 1. Step 3: Set 01eft = 0c, 0 c : 0right , and 0right : 0 c + eo. Compute Z(0right) and go back to Step I. Step 4: Output 0~ as the result of Procedure Approx (0o, eo). Note that in the above procedure each time a value z(O) for some 0 needs to be computed, a complete optimization process on x and y (subsection 4.1) has to be carried out. Computational experience (see Section 5 below) shows that, if e o is not chosen too small, this procedure requires on the average the computation of z(O) for very few values of 0, in order to get an eo-approximation to the local minimum closest in some sense to the starting point 0o. If high accuracy were required, it would always be possible to iteratively apply Procedure Approx with smaller and smaller values of eo. A possible strategy could be the following: (a) Set 01 = 0 0 and e I = ~0" (b) Compute 02 -- Approx (01,el). If el is smaller
Volume 9. Number 5
PATTERN RECOGNITION LETTERS
than some prescribed accuracy, then stop with 02 being the final solution obtained. Otherwise take ~2 = q ' P (where p is some constant reduction factor, e.g. p = 0.1). (c) Replace el by g2~ 01 by Oz and go back to (b).
June 1989
4.3. Starting solution The procedure described above may be made to converge faster if good starting values Xo, Yo and 0o for x, y, and 0 are available. If exact coincidence of
~
P.01
a_
j -
.02
a_
b
b
C
C
Figure I, Illustration for 5 test problems: (a) original relative position of P I and P> (b) initial solution obtained by matching the principal axes of inertia of Pj and P2, and (c) optimal solution. 331
Volume 9, Number 5
PATTERN RECOGNITION LETTERS
P .05
June 1989
P .06
a
a
b
C Figure 1 (continued).
the two contours could be achieved, such a coincidence would be obtained by shifting P2 to the center of gravity of P1 and by rotating P2 so as to get the same inertia axes for both polygons. The coordinates (cx, cy) of the center of gravity of 332
P1 can be computed in time complexity O(N1). Determining the principal axes of inertia of PI (resp. P2) requires O ( N 0 time (resp. O(N2) time). Let be the angle between the principal axes of inertia of P1 and P2. Then, we take as starting solution:
Volume 9. Numbcr 5
PATTERN RECOGNITION LETTERS
P .10
June 1989
ment to P2 such that its major axis of inertia coincides with that of P1.
5. Computational results
a
b
We present in Table l the results obtained by the application of our approach to some test problems. For each of them we present the number of nodes N~ and N 2 of polygons Pt and P2, respectively: the proximity measure Do for the starting solution Xo, Yo, 0o; the final value D* obtained for the proximity measure after application of the optimization process: the relative value D*/Do of the optimal solution over the initial solution: the number n t of computations of z(0): the average number n2 of iterations of the steepest descent method for each fixed value 0 = {~: and the CPU time (in seconds) needed to solve the problem on a DEC Microvax computer. The computational results presented in Table 1 confirm the practical applicability of our approach. In particular, the number of iterations is seen to be always quite low, and this is mainly due to the quality of the starting solution. Figure 1 illustrates the results obtained for 5 test problems, and allows one to check visually the quality of the initial solution (obtained by the method proposed in subsection 4.3) and that of the final optimal solution. Apart from other applications in image processing and pattern analysis, the ideas presented in this paper are currently being extended to the three-dimensional case (see (Bloch-Boulanger et al.. 1988)) and will be incorporated into the 3-D data acquisition system mentioned in the Introduction (see (Lelandais, 1984), (Brunet, 1987)).
Acknowledgements
C Figure I (continued).
(Xo,Yo) = (cx, cy), 0, if r ( x o, Yo, [9) <~ F(xo, Y0, b + rt),
0° =
0 + n, otherwise,
which corresponds to the application of a displace-
We would like to acknowledge lsabelle BlochBoulanger for her active participation in the computational experiments.
References Aggarwal, A. (1988). Private communication. Arun. K.S., T.S. Huang and S.D. Blostein(1987). Least-squares fitting of two 3-D point sets. IEEE Trans. Pattern Anal. Machine Intell. 9(5), 698 700. 333
June 1989
PATTERN RECOGNITION LETTERS
Volume 9, Number 5 Table 1 Computational results
Problem
NI
N2
DO
D*
D*/D o
nI
n2
CPU(sec)
P.01 P.02 P.03 P.04 P.05 P.06 P.07 P.08 P.09 P. 10
10 15 20 20 30 20 50 60 11 18
9 15 15 20 25 20 40 50 14 17
11.22 49.83 60.87 215.51 79.17 44.39 284.91 401.20 1.17 5.70
0.78 26.71 45.01 150.59 32.95 19.19 121.41 361.16 0.57 3.70
7% 54% 74% 70% 42% 43% 43% 90% 49% 65%
5 5 5 11 6 4 5 5 18 II
3 3 3 4 3 3 4 6 2 3
1.39 2.41 4.91 13.38 7.11 3.18 16.71 37.91 3.36 0.74
Barrow, H.J., J.M. Tenenbaum, R.C. Bolles and H.C. Wolf (1977). Parametric correspondence and the chamfer matching: two new techniques for image matching. Proc. 5th International Joint Conference on Artificial Intelligence, Cambridge, 659-663. Bloch-Boulanger, I., H. Maitre, M. Minoux and C. Ribeiro (1988). Recalage de polygones et polyh6dres convexes. Rapport de recherches ENST-88D004, Ecole Nationale Sup6rieure des T616communications, Paris. Borgefors, G. (1986). Distance transformations in digital images. Computer Vision. Graphics, and Image Processing 34, 344-371. Brunet, M. (1987). Laser-based 3-D digitizing. Computer Graphics World, November 1987, 103-104. Burr, D.J. ( 1981 ). Elastic matching of line drawing. IEEE Trans. Pattern Anal. Machine Intell. 3(6), 703-713. Faugeras, O.D. and M. Hebert (1983). A 3-D recognition and positioning algorithm using geometrical matching between primitive surfaces. Proc. 8th International Joint Conference on Artificial Intelligence, Karlsruhe, 996 1002. Faugeras, O.D. and M. Hebert (1986). The representation, recognition and positioning of 3-D shapes from range data. In: A. Rosenfeld, Ed., Techniques./br 3-D Machine Perception. North-Holland, Amsterdam, 13 52, Kass, M., A. Witkin and D. Terzopoulos (1987). Snakes: active contour models. Proc. 1st International Conference on Computer Vision, London, 259 268. Lelandais, S. (1984). R6alisation et premi6re exploration d'un
334
syst6me de numbrisation de formes tridimensionelles. These de Doctorat, Universit6 de Compiegne, France. Lin, Z.-C., T.S. Huang, S.D. Blostein, H. Lee and E.A. Margerum (1986). Motion estimation from 3-D point sets with and without correspondences. Proc, IEEE Conference on Computer Vision and Pattern Recognition, Miami Beach, 194- 201. Luenberger, D.G. (1973). Introduction to Linear and Nonlinear Programming. Addison-Wesley, Reading, MA. Minoux, M. (1986). Mathematical Programming. Wiley, New York. Pavlidis, T. (1978). A review of algorithms for shape analysis. Computer Graphics and Image Processing 7(2), 211 242. Pratt, T.K. (1978). Digital Image Proeessing. Wiley, New York. Preparata, F.P. and M.I. Shamos (1985). Computational Geometry: An Introduction. Springer, New York. Schmitt, F., H. Maitre, A. Clainchard and J. Lopez-Krahe (1985). Acquisition and representation of real object surface data. 2nd International Technical Symposium on Optical and Electro-Optical Applied Science and Engineering, BioStereoMetrics "85 Conference, Cannes, SPIE Proc. 602. Sklansky, J. (1978). On the Hough technique for curve detection. 1EEE Trans. Computers 27 (100), 923-926. Tsai, W.H. and S.S. Yu (1985). Attributed string matching with merging for shape recognition. IEEE Trans. Pattern Anal. Machine lntell. 7(4), 453~460. Widrow, B. (1973). The rubber mask technique I. Pattern Recognition 5, 175 211.
June 1989
North-Holland
Optimal matching of convex polygons Pedro COX, Henri M A I T R E Laboraloire Image, Ecole Natianale Sup~;rieure des T~;h;communicalions, 46 rue Barrault, Paris 75634 Cedex 13, France
Michel M I N O U X Lahoratoire MASI, University; Paris VI. Place Jussieu, Paris 75005, France
Celso RIBEIRO* Department qf Eh,ctrical Engineering, Catholic University ~[' Rio de Janeiro, G~h'ea, Caixa Postal 38063, Rio de Janeiro 22452, Brazil
Received 20 April 1988 Revised 24 November 1988
Abstract: We describe in this paper the problem of matching two series of observations of a 2-D contour, in which each sampling of the contour has different hidden parts. This problem is formulated as the minimization of the distance between the convex hulls of the two samplings and a solution method, based on nonlinear programming techniques, is presented. Computational results are reported to illustrate the practical applicability of the approach. Key words." Contour matching, contour fitting, image matching, pattern recognition, 2-D recognition.
1. Introduction In the domain of robot vision, the use of laser range-finders to acquire a description of the environment is becoming more and more popular; see e.g. (Faugeras and Hebert, 1986). Similar systems are also of increasing use for quality control, 3-D data basis formation and CAD/CAM; see e.g. (Schmitt et al., 1985). For these applications, the object to be measured is usually placed in a structured light beam and the coordinates of the lit points of the surface are measured by triangulation. One of such data acquisition systems is described in details by Lelandais (1984) and Brunet (1987). When only one laser source is used, the surface measurement is often incomplete because of the existence of hidden parts. For this reason, the experi-
ment has to be repeated with a second laser beam in a different position. As a matter of fact, matching these two sets of measurements arises as an important stage in the process of providing a single, more extensive description of the object. This matching problem is addressed in this paper in a simplified form: we limit ourselves to the two-dimensional case, where the observation of the object surface is reduced to one planar closed contour line. Nevertheless, our approach is currently being extended to the three-dimensional case and will be dealt with in a forthcoming paper by Bloch-Boulanger et al. (1988). Besides structural methods (Tsai and Yu, 1985) and contour coding techniques (Pavlidis, 1978), which can not be used here because they give no information about the relative position of the two
* The work of the last author was sponsored by FINEP (research contract grant number 4.3.86.0689-00), CNPq (research contract grant number 302281-85/EE), IBM Brasil and the French Ministry of Foreign Relations. 0167-8655,'89/$3.50 (~ 1989, Elsevier Science Publishers B.V. (North-Holland)
327
Volume 9, Number 5
PATTERN RECOGNITION LETTERS
shapes, most of the existing work rely on an estimation of the Euclidean distance between the shapes, called 'cross correlation', 'template matching' (Pratt, 1978) or 'Hough transform' (Sklansky, 1978). These methods require the computation of a distance measure of the two shapes for each relative position. To reduce the computational cost, 'chamfer' techniques, which precompute the distance of each pixel of the field to one of the shapes, were proposed; see e.g. (Barrow et al., 1977). Later, Borgefors (1986) introduced efficient approximations of the Euclidean distance to compute the chamfer. Burr (1981) introduced a local measure of the distance between the two shapes, based on a 'weighted' Euclidean distance, providing a perfect matching after a rubber-like distortion of the plane. Widrow (1973) and Kass et al. (1987) have proposed mechanical models to express the distances between contours. Unfortunately, for all the above mentioned techniques, no efficient strategy guaranteed to converge to a minimum distance solution is available. Several techniques were proposed to match shapes when the correspondence between the points of the two shapes is known beforehand: an exact technique in 2-D and an interactive one in 3-D (Lin et al., 1986), an exact one in 3-D using quaternions (Faugeras and Hebert, 1983) and another exact one in 3-D based on the singular value decomposition of a 3 x 3 matrix (Arun et al., 1987). But, in many cases, finding the appropriate correspondence is by far the most difficult part of the problem. The matching technique described in the next sections does not rely on such a restrictive assumption.
2. Matching two observations of a 2-D contour: problem statement We consider a 2-D contour (not necessarily convex) from which we extract discretized samples by the laser range-finder system mentioned in the introduction. Depending on various parameters of the data acquisition system (such as the angle of observation and the position of the rotation axis), the hidden parts are not the same for each sampling of the contour: it may be the case that some point oh328
June 1989
tained in a first sample corresponds to a hidden point for a second sample. We study from now on the problem of how to exploit two distinct samples of the same contour in order to get a more extensive representation of the whole contour (i.e., with fewer non-observed points). Our approach proceeds as follows: (i) For each sample, we construct the convex hull of the points in the sample and we select the subset of points which are the extreme points of the convex hull (the reader is referred to (Preparata and Shamos, 1985) for a survey of algorithms for convex hull calculations), (ii) Once the convex hulls of the two samples have been determined, we solve the problem of matching in an optimal way (according to a criterion to be defined below) the two convex polygons whose extreme points are those selected in (i). In the remainder, we concentrate on the design of an efficient way of solving part (ii).
3. Optimal matching of two convex polygons The two convex polygons P l and P2 are given by specifying the following data: (i) N~ (resp. N2), the number of extreme points of polygon P1 (resp. P2); (ii) the coordinates a s, bj~ (resp. a~, b~) of the jth extreme point of polygon P~ (resp. P2), for j = 1..... N 1 (resp. j = 1..... N2); without loss of generality the origin of the axes is supposed to coincide with the center of gravity of polygon P2; and (iii) the numbering of the extreme points is supposed to be such that for any j~{1 ..... N1} (resp. j 6 { 1..... N21 ), if we define
j,
IJ + 1 ,
ifj
l, ,
ifj = N1 (resp. j = N2),
then the segment joining the extreme points numberedj andj' is an edge of polygon P~ (resp. P2)In order to have a measure of how well the two polygons match, we shall need an objective function. This proximity measure should be designed in such a way as to be the smallest, the closest the two boundaries of the polygons are to each other. For any point with coordinates (u,v), we define its distance to polygon PI as
Volume 9, Number 5
PATTERN R E C O G N I T I O N LETTERS
{ ll(u,v)- Projp,(u,v)l I, @t(u, v) =
otherwise,
where Proje~ (u, v) denotes the point of P~ achieving minimum euclidean distance to (u,v) and int(P0 denotes the interior of P~. The distance to polygon P2 is defined in a similar way. With this notation, the proximity measure corresponding to the initial relative position of the two polygons will be taken as NI =
N2
2-2 j=l
4. Minimizing the proximity measure
if (u, v) ¢ int(Pt), 0,
D
June 1989
, 1 b l,
°P2~aj,
J)-{-
2 ,j=l
6z ta2 bZ~
P1 ~" J" J*"
Note that very small values of D are attained only if, for each of the polygons, each extreme point is either in the interior of the other polygon, or very close to it. In order to get such a desirable situation, we see that we have to minimize the proximity measure over all possible relative positions of the two polygons. To generate all possible relative positions, it is enough to apply to one of the polygons (P2 in this case) all possible displacements on the plane, an arbitrary displacement being specified by (i) a translation vector (x,y): and (ii) a rotation angle 0 around the origin of the axes.
Applying a displacement with parameters (x,y, O) to any point of polygon P2 with coordinates (:
Fc°sO
If
ksin 0
-sinO](~) cos
oJ\IU"
In a natural way, we will denote by T(x,y,O,P 2) the polygon resulting from the application of a displacement with parameters (x,y,O) to each of the extreme points of P2. With this notation, looking for values of x, y and 0 resulting in a minimum proximity measure can be reduced to the (unconstrained) minimization of the objectivc function (in three continuous variables only): N1
F(x, .v,O) = ~, (~T(x,y,O, -2 P2)\~aj~, b~ jl j= 1 N2
+ ~ a~ r(x,y,O,a~,b~). j=l
The objective function F(x,y,O) has a number of properties which have to be taken into account in the design of an algorithm for the minimization problem. In particular: (i) F is continuous differentiable with respect to all three variables x, y, and 0: (ii) F is convex with respect to x and y, for fixed O: and (iii) F is not convex with respect to the variable O, for fixed x and y. In view of the nonconvexity with respect to the variable O, we propose here to decompose the minimization by considering separately the variables corresponding to translation (x and y) and the variable corresponding to rotation (0).
4.I, Minimizing the ohjective fimction for lixed 0 The objective function F(x, y, O) being convex and differentiable for 0 fixed to some value [), many unconstrained optimization methods can be applied. In our experiments we tried the steepest descent algorithm (SD) and the conjugate gradient algorithm (CG) (see e.g. (Luenberger, 1973)) and observed no significant differences in the behavior of the two: in each case the number of iterations required to achieve sufficient accuracy was consistently very low (typically less than five iterations). This can be explained by the fact that F, as a sum of distance functions, is very well-conditioned. It is well known that both SD and CG algorithms require, at each step, the computation of: (i) the value of the function F a t the current point; (ii) the gradient g = (g.~,gy) of F at the current point: and (iii) the minimum of F along the chosen descent direction d = (d,., d~.).
(a) Computing the current value o f F at (x,y).[br 0 = OJixed Computing the projection of any point (u, v) on polygon PI or deciding that it lies in cony[P1] can be done in time O(N0. Similarly, for projecting on P2, this computation requires O(N2) elementary operations. Therefore, the total complexity of computing the value of F a t any point (x,y,O) would be 329
Volume 9, N u m b e r 5
PATTERN RECOGNITION LETTERS
O(N~ • Nz). However, it should be pointed out that this computation can be done in time complexity O(N1 + Nz), as suggested by Aggarwal (1988). As a result of this computation, we get for every 1 vertex (a j,1 bj) of P1 the coordinates of its projection on conv[P2], which will be denoted by (~/J,/~J), and for every vertex (a],by) of P2 the coordinates (~/~,/~2) of its projection on conv[P1]. Note that ~/), /~), ~/2, a n d / ~ implicitly depend on x, y and 0. With this notation, the value of F(x, y, O) can be rewritten as
r(x,y,O)
N1
=
[(a) -
+ (b) - t;))
j=l N2
+ E [ ( a ] - - t i ] ) 2 + ( b 2 - @ ) 2 3. j=l
(b) Computing the gradient o f F at (x,y) for 0 = fixed The gradient of a sum of differentiable functions being the (vector) sum of their gradients, the gradient of F at (x,y) for 0 = 0 fixed can be easily computed with linear complexity O(N1 + N2) as N1
V F(x,y,O) = 2 • (a) - aj,b~-I 1 _ ~)) j:l N2
+22(4
-
a),b)
-
j=l
(c) Minimizing the objective function along the descent direction (dx, dy) There are many possible methods for minimizing a given function of a single variable along a (descent) direction. One may use, for instance, methods based on interval reduction (see e.g. (Minoux, 1986)). However, when dealing with a convex differentiable function such that computing the gradient is computationally cheap once the function value has already been computed, the so-called secant method features good global and asymptotic convergence properties. Each step of this secant method requires only O(N 1 + N2) time.
4.2. Finding the best value./or 0 If, for each value 0 of 0, we apply the algorithm outlined in (a), (b) and (c) above, we get optimal values for x and y which are implicit functions of 0. 330
June 1989
We denote these implicit functions respectively by )~(0) and .9(0). Thus, what we want is to find 0* minimizing
z(O) = F(£(O), y(O), O) over all possible angles 0. However, due to the nonconvexity of F(x,y,O), we cannot expect z(O) to be convex. Therefore, neither gradient methods nor interval reduction methods are guaranteed to produce a global minimum (because of non-unimodality). In view of this, the following technique has been used. We assume that a good approximation 0o to the global minimum 0* can be computed (subsection 4.3 below describes a way of getting such an approximate solution). We then define a discretization step eo (typically, eo = I degree) and look for a local minimum in the vicinity of 0 o by the following simple method:
Procedure Approx (0o, ~ ) Step 0: Set 0 c = 0 o, 01eft : 0 0 - - g O , and 0right : 00 + to. Compute Z(0o), z(01m), and Z(0right ). Step 1: Set 0m = 0~eft, if Z(01m)~< Z(0right); 0m = 0right, otherwise. If Z(0m)~> z(Oc) go to Step 4. If 0 m : 01eft g o to Step 2, if 0m = 0right go to Step 3. Step 2: Set 0right = 0c, 0 c : 01eft , and 01m = 0¢ eo. Compute z(01m) and go back to Step 1. Step 3: Set 01eft = 0c, 0 c : 0right , and 0right : 0 c + eo. Compute Z(0right) and go back to Step I. Step 4: Output 0~ as the result of Procedure Approx (0o, eo). Note that in the above procedure each time a value z(O) for some 0 needs to be computed, a complete optimization process on x and y (subsection 4.1) has to be carried out. Computational experience (see Section 5 below) shows that, if e o is not chosen too small, this procedure requires on the average the computation of z(O) for very few values of 0, in order to get an eo-approximation to the local minimum closest in some sense to the starting point 0o. If high accuracy were required, it would always be possible to iteratively apply Procedure Approx with smaller and smaller values of eo. A possible strategy could be the following: (a) Set 01 = 0 0 and e I = ~0" (b) Compute 02 -- Approx (01,el). If el is smaller
Volume 9. Number 5
PATTERN RECOGNITION LETTERS
than some prescribed accuracy, then stop with 02 being the final solution obtained. Otherwise take ~2 = q ' P (where p is some constant reduction factor, e.g. p = 0.1). (c) Replace el by g2~ 01 by Oz and go back to (b).
June 1989
4.3. Starting solution The procedure described above may be made to converge faster if good starting values Xo, Yo and 0o for x, y, and 0 are available. If exact coincidence of
~
P.01
a_
j -
.02
a_
b
b
C
C
Figure I, Illustration for 5 test problems: (a) original relative position of P I and P> (b) initial solution obtained by matching the principal axes of inertia of Pj and P2, and (c) optimal solution. 331
Volume 9, Number 5
PATTERN RECOGNITION LETTERS
P .05
June 1989
P .06
a
a
b
C Figure 1 (continued).
the two contours could be achieved, such a coincidence would be obtained by shifting P2 to the center of gravity of P1 and by rotating P2 so as to get the same inertia axes for both polygons. The coordinates (cx, cy) of the center of gravity of 332
P1 can be computed in time complexity O(N1). Determining the principal axes of inertia of PI (resp. P2) requires O ( N 0 time (resp. O(N2) time). Let be the angle between the principal axes of inertia of P1 and P2. Then, we take as starting solution:
Volume 9. Numbcr 5
PATTERN RECOGNITION LETTERS
P .10
June 1989
ment to P2 such that its major axis of inertia coincides with that of P1.
5. Computational results
a
b
We present in Table l the results obtained by the application of our approach to some test problems. For each of them we present the number of nodes N~ and N 2 of polygons Pt and P2, respectively: the proximity measure Do for the starting solution Xo, Yo, 0o; the final value D* obtained for the proximity measure after application of the optimization process: the relative value D*/Do of the optimal solution over the initial solution: the number n t of computations of z(0): the average number n2 of iterations of the steepest descent method for each fixed value 0 = {~: and the CPU time (in seconds) needed to solve the problem on a DEC Microvax computer. The computational results presented in Table 1 confirm the practical applicability of our approach. In particular, the number of iterations is seen to be always quite low, and this is mainly due to the quality of the starting solution. Figure 1 illustrates the results obtained for 5 test problems, and allows one to check visually the quality of the initial solution (obtained by the method proposed in subsection 4.3) and that of the final optimal solution. Apart from other applications in image processing and pattern analysis, the ideas presented in this paper are currently being extended to the three-dimensional case (see (Bloch-Boulanger et al.. 1988)) and will be incorporated into the 3-D data acquisition system mentioned in the Introduction (see (Lelandais, 1984), (Brunet, 1987)).
Acknowledgements
C Figure I (continued).
(Xo,Yo) = (cx, cy), 0, if r ( x o, Yo, [9) <~ F(xo, Y0, b + rt),
0° =
0 + n, otherwise,
which corresponds to the application of a displace-
We would like to acknowledge lsabelle BlochBoulanger for her active participation in the computational experiments.
References Aggarwal, A. (1988). Private communication. Arun. K.S., T.S. Huang and S.D. Blostein(1987). Least-squares fitting of two 3-D point sets. IEEE Trans. Pattern Anal. Machine Intell. 9(5), 698 700. 333
June 1989
PATTERN RECOGNITION LETTERS
Volume 9, Number 5 Table 1 Computational results
Problem
NI
N2
DO
D*
D*/D o
nI
n2
CPU(sec)
P.01 P.02 P.03 P.04 P.05 P.06 P.07 P.08 P.09 P. 10
10 15 20 20 30 20 50 60 11 18
9 15 15 20 25 20 40 50 14 17
11.22 49.83 60.87 215.51 79.17 44.39 284.91 401.20 1.17 5.70
0.78 26.71 45.01 150.59 32.95 19.19 121.41 361.16 0.57 3.70
7% 54% 74% 70% 42% 43% 43% 90% 49% 65%
5 5 5 11 6 4 5 5 18 II
3 3 3 4 3 3 4 6 2 3
1.39 2.41 4.91 13.38 7.11 3.18 16.71 37.91 3.36 0.74
Barrow, H.J., J.M. Tenenbaum, R.C. Bolles and H.C. Wolf (1977). Parametric correspondence and the chamfer matching: two new techniques for image matching. Proc. 5th International Joint Conference on Artificial Intelligence, Cambridge, 659-663. Bloch-Boulanger, I., H. Maitre, M. Minoux and C. Ribeiro (1988). Recalage de polygones et polyh6dres convexes. Rapport de recherches ENST-88D004, Ecole Nationale Sup6rieure des T616communications, Paris. Borgefors, G. (1986). Distance transformations in digital images. Computer Vision. Graphics, and Image Processing 34, 344-371. Brunet, M. (1987). Laser-based 3-D digitizing. Computer Graphics World, November 1987, 103-104. Burr, D.J. ( 1981 ). Elastic matching of line drawing. IEEE Trans. Pattern Anal. Machine Intell. 3(6), 703-713. Faugeras, O.D. and M. Hebert (1983). A 3-D recognition and positioning algorithm using geometrical matching between primitive surfaces. Proc. 8th International Joint Conference on Artificial Intelligence, Karlsruhe, 996 1002. Faugeras, O.D. and M. Hebert (1986). The representation, recognition and positioning of 3-D shapes from range data. In: A. Rosenfeld, Ed., Techniques./br 3-D Machine Perception. North-Holland, Amsterdam, 13 52, Kass, M., A. Witkin and D. Terzopoulos (1987). Snakes: active contour models. Proc. 1st International Conference on Computer Vision, London, 259 268. Lelandais, S. (1984). R6alisation et premi6re exploration d'un
334
syst6me de numbrisation de formes tridimensionelles. These de Doctorat, Universit6 de Compiegne, France. Lin, Z.-C., T.S. Huang, S.D. Blostein, H. Lee and E.A. Margerum (1986). Motion estimation from 3-D point sets with and without correspondences. Proc, IEEE Conference on Computer Vision and Pattern Recognition, Miami Beach, 194- 201. Luenberger, D.G. (1973). Introduction to Linear and Nonlinear Programming. Addison-Wesley, Reading, MA. Minoux, M. (1986). Mathematical Programming. Wiley, New York. Pavlidis, T. (1978). A review of algorithms for shape analysis. Computer Graphics and Image Processing 7(2), 211 242. Pratt, T.K. (1978). Digital Image Proeessing. Wiley, New York. Preparata, F.P. and M.I. Shamos (1985). Computational Geometry: An Introduction. Springer, New York. Schmitt, F., H. Maitre, A. Clainchard and J. Lopez-Krahe (1985). Acquisition and representation of real object surface data. 2nd International Technical Symposium on Optical and Electro-Optical Applied Science and Engineering, BioStereoMetrics "85 Conference, Cannes, SPIE Proc. 602. Sklansky, J. (1978). On the Hough technique for curve detection. 1EEE Trans. Computers 27 (100), 923-926. Tsai, W.H. and S.S. Yu (1985). Attributed string matching with merging for shape recognition. IEEE Trans. Pattern Anal. Machine lntell. 7(4), 453~460. Widrow, B. (1973). The rubber mask technique I. Pattern Recognition 5, 175 211.