Optimal Subpixel Matching of Contour Chains and Segments Bruno Serra  y

Marc Berthod 

INRIA Sophia-Antipolis, Projet Pastis, BP 93, 06902 Sophia-Antipolis Cedex, FRANCE AEROSPATIALE, E/ETRV, 1 Rue Pablo Picasso, 78114 Magny-les-Hameaux, FRANCE

 y

Abstract This paper introduces a new general purpose algorithm that allows the optimal geometric match between contours to be determined, that is the transformation yielding a minimal deformation is obtained. The algorithm relies only on the geometric properties of the contours and does not call for any other constraint, so that it is particularly suitable when no parameterization of the deformation is available or desirable. Contour deformation is explicitly incorporated in the computation, allowing for a thorough use of all geometric information available. Moreover, no discretization is involved in the computation, resulting in two main advantages: rst, the algorithm is robust to dierences in the segmentation of contours and allows the matching of polygonal approximations of contours with very little loss of precision, second, subpixel precision matching can be achieved.

1 Introduction To establish point-to-point correspondences between contours, most methods rely on (1) strong geometric constraints such as epipolar lines or parameterized deformations (2) heuristic criteria such as intensity similarity or local gradient orientation. While the former constraints are not always available (e.g. the motion between the two images might be unknown), the latter prove inapplicable in many cases: changes in luminosity cause the local intensities to vary, perspective deformations aect the gradient orientation. Most existing methods for matching images with weak constraints iteratively distort one image to t the second [2] [3] [7]. The main drawbacks of these approaches are the need of a global preregistration process and their iterative implementation which results in large computational requirements. The approach we propose is aimed at exploiting to its best the geometric information provided by the

contours themselves, using the only reliable information between contours, i.e. shape similarity. The rst step is to dene a deformation measure to compare two contours. This measure should be robust to dierences in the segmentation of contours, especially when dealing with polygonal contours: two contours might be very similar even though not segmented in the same way (dierent number of segments with their endpoints at dierent positions). By dening the correspondence as a continuous function and by using a variational formulation for the deformation measure, we are able to tackle this problem. At the same time, subpixel matching is achieved, since the continuous nature of the correspondence function implies that matches are available for any point on the contours, being an endpoint or not. The second step involves the minimization of the deformation measure over the possible correspondences. This minimization is achieved in two stages. By considering only two linear portions of contours, we rst derive an analytical expression for the correspondence with minimal deformation. This expression is then used as a local cost function in a dynamic programming algorithm which ensures global optimality over the whole contour. The cost function being given as an analytical expression and the correspondence function being continuous, a novel scheme called Continuous Dynamic Programming was devised. This new technique allows an optimal continuous path in a graph to be determined in an ecient fashion.

2 A geometric deformation measure Given two contours dened by their curvilinear coordinates C1(r) and C2(s) (Fig. 1), let = ( rs ) represent a match between point C1 (r) and point C2 (s), and D( ) = C2(s) ? C1 (r) be the displacement vector between the two points. We can associate to

s

C (s) 2

s

nite elements or nite dierences method. However, it may be solved given an analytic expression of D: considering two linear contour portions:

C2

0

D r0

1

r1

C1

C1(r) = a1 r + b1; C2(s) = a2 s + b2; ka1 k = ka2k = 1

C (r) 1

Figure 1: Displacement vector. any match between the two contours a correspondence function (t) = ( rs((tt)) ) dened over an interval [t0 ; t1 ] such that (t0) = ( sr0 ), (t1) = ( sr1 ), and for 1 0 any t, C1(r(t)) is matched with C2 (s(t)). This denition is however not enough to ensure that the function

(t) is unique, since for any function f the function 0 such that 0(t) = (f (t)) will also be consistent with the denition. The following constraint is thus added to ensure uniqueness, by imposing that the correspondence function varies in a uniform fashion over the contours (the matches should not accumulate in some portions of the contours or be sparse in over portions): k _ (t)k = 1 (1) Given two contours and a correspondence function, we can now dene a deformation measure as: Z t1 (2) D = k dD( (t))k2 dt t0

dt

This expression measures the variation of the displacement vector along the contours, and thus the deformation from one contour to the other. It is similar to the one used by Hildreth [5] to compute the velocity eld between images, and can be related to the gural continuity and disparity gradient constraints introduced by Mayhew and Frisby [8]. Finding the optimal deformation measure can now be stated as a variational problem where the following functional has to be minimized: Z t1 D( ) = F ( ; t)dt where : F ( ; t) = k dD( (t)) k2 dt

t0

A classical approach involves the resolution of the Euler-Lagrange equation: d F ? F _ = 0 dt

Expanding this equation yields: [D ( (t)) (t) + D

( (t)) _ (t)2]T D ( (t)) = 0 (3) This equation can of course not be solved in the general case, so that it is usually solved by an iterative

eq. 3 reduces to:

(t) = 0 so that the optimal match is a linear function. We can now determine the associated deformation measure (Eq. 2). Since is linear, we can write: r(t) = t + 0 ;

s(t) = t + 0

 and can be determined by using the boundary conditions (t0) and (t1) : =

r1 ? r0 ; t1 ? t0

=

s1 ? s 0 t1 ? t0

and the constraint on (Eq. 1) becomes: (t1 ? t0 )2 = (r1 ? r0 )2 + (s1 ? s0 )2 Substituting the previous equations in the denition of the deformation measure (Eq. 2) nally yields the analytical expression of the optimal deformation: 2 (4) D = k[(ar2 (s?1 ?r s)20)+?(sa1 (?r1s?)2r]01)=k2 1 0 1 0

3 Optimal matching of contours using continuous dynamic programming Given the analytical expression for the optimal matching of linear contour portions, we will now take advantage of the fact that contours can be viewed as chains of linear portions (either contour segments or "micro-segments" between contour points) to derive a global optimum for the matching of whole contours. Considering two contours, the optimal deformation can be expressed as:

D = min

N X i=0

Di

where is piecewise linear and Di is the optimal deformation between two contour segments (eq. 4). This expression can be viewed as the search for an optimal path in a graph constructed by placing on each axis the arc length on each contour (Fig. 2). Any square portion of the graph corresponds to two segments of contour, so that the optimal path is linear in each square of the graph, and we can associate with each path portion the cost function Di given by eq. 4.

contour 2

3 2 point 1 point 1

2

3

contour 1

Figure 2: Matching viewed as a path in a graph. The optimal match is the path with minimal total cost in the graph. It can be entirely dened by its intersections with the lines and columns of the graph, so that only these points need to be determined. Dynamic programming [4] [1] is a method of choice to nd an optimal path in a graph:  it is guaranteed to nd the global optimum;  it is a non-iterative method;  it is computationally ecient. However, classical dynamic programming algorithms are restricted to nding a graph going through the nodes of a graph, whereas in our case the path can cross the lines and columns of the graph anywhere. For this purpose, we have conceived a novel dynamic programming scheme, called Continuous Dynamic Programming, which allows an optimal continuous path to be determined. Discrete dynamic programming proceeds as follows:  A null cost is assigned to all points in the rst column of the graph, allowing the rst contour point to match any point on the second contour.  The graph is explored from left to right in order to nd the best path ending in each node.  For each node in the graph, a simple reasoning yields the best path: for each possible node in the previous column, add the cost of the path ending in that node and the local cost from that node to the current node, and retain the one yielding a minimum cost.  After the last column in the graph has been explored, nd the node in that column with minimal cost, then backtrack in the graph to get the complete path. Continuous dynamic programming proceeds in a similar way, except that instead of associating a numerical cost value to each node, the cost of the optimum path ending in any point between two nodes is given by a parametric expression. The previous point on the path is similarly remembered as a parametric expression. An initial description of this scheme has been described in [9], however with the initial hypothesis that the local cost function is quadratic. In the present case, the local cost function (Eq. 4) is not quadratic, but we will recursively exhibit a linear expression for the total cost function.

Considering a single square of the graph (Fig. 3), we will suppose that the cost of the optimal paths ending on the lower side of the square are linear, and we will show that under that condition the cost of the optimal paths ending on the right side of the square are also linear, so that by recursion the same result holds for any square of the graph (the proof being the same for a path going through the other sides of the square). s (contour 2)

si−1

si

contour 2

ri−1 ri

contour 1

si r i−1 r (contour 1)

Figure 3: A linear path portion. Suppose that the cost D0;i?1 of the optimal path ending in ri?1 is linear: D0;i?1(ri?1) = ri?1 + The optimal path ending in si will be given by: D0;i(si) = min (D0;i?1(ri?1) + Di?1) r ?1

i

Substituting the value of Di?1 (Eq. 4) yields:   2 dr + ds2 ? 2a1 :a2 drds 0 ? dr + D0;i (si) = min dr [dr2 + ds2 ]1=2 with dr = ri ? ri?1; ds = si ? si?1; 0 = + ri This minimization can be simplied as follows: 



0 with : z = dr D0;i(si) = min f ( z ) ds + z ds z2 ? 2a1 :a2 z + 1 p and : f (z) = ? z

1 + z2 Which shows that the optimal path ending in si is linear. We still need to nd the minimum of function f (z) (note that z is simply the inverse of the slope of the path in the graph). Unfortunately, f does not have a single minimum in the general case. However, this diculty can be overcome if we impose the following simple constraint: a1 :a2 and z should have the same sign. This condition is not very restrictive, since it is equivalent to imposing that the angle between the two segments is less than =2. Under that condition, the second derivative of f can easily be proved positive, which implies that there is a single minimum. Furthermore, the minimum is obtained for z = 0 if a1a2  ?2, and for z = 1 if   1. It can be calculated numerically in the other cases. Given the value zmin which minimizes function f , the cost of the optimal path ending in si is simply D0;i (si) = f (zmin)ds + 0 , and the previous point on the path is given by dr = zminds.

However, the following problem remains: considering the case of gure 4, C1 and C3 are linear cost portions corresponding to optimal paths coming from the lower and left sides of the square. In the case of portion C2, the optimal path comes from the lower left corner of the graph, so that dr is xed and the cost function is no longer linear: 2 2 D0;i (si) = dr +[drds2 +?ds2a21]1:a=22 drds ? dr + 0 We deal with this non-linearity by approximating the cost function by successive linear portions, with a xed maximal approximation error.

}C }C

3

2

} C1

Figure 4: Cost portions.

4 Experimental analysis In order to analyze experimentally the precision and robustness to segmentation of the algorithm, a number of experiments have been carried out on two synthetic contour images of a sphere with subpixel edge localization.

4.1 Qualitative and quantitative analysis In order to appreciate qualitatively the matching precision, a 3D reconstruction was done using the output of the algorithm. The result is presented in gure 5, and shows the precision as well as the smoothness of the matching process. For comparison, the result obtained with a discrete algorithm is presented. A quantitative evaluation was then performed: since the exact position of the matches is known on these synthetic images, the mean error between each match and the actual position was calculated, yielding a value of 0.10 pixel. Under the same conditions, the discrete algorithm used for comparison produced a mean error of 0.9 pixel.

Figure 5: Subpixel 3D reconstruction (mean error: 0.10 pixel). vs. discrete algorithm (error: 0.88 pixel). For each approximation, the mean segmentation error was computed, as well as the mean matching error after running the algorithm. Figure 6 shows the graph of the matching error versus the segmentation error. This graph shows the smooth behaviour of the algorithm with respect to segmentation: the increase of the matching error is nearly linear with the segmentation error, which shows that the matching process is little aected by the segmentation of contours and that the algorithm is insensitive to the uneven distribution of points along contours. The polygonal approximation process yields segments whose lengths vary greatly and which are not necessarily similarly distributed in both images, yet the correct correspondences are recovered.

4.3 Computational complexity Fig. 7 shows a plot of the computation time for dierent values of the segmentation tolerance against the number of contour segments in one image. The computational complexity is nearly linear. One must keep in mind that in this experiment, since the size of the contours was xed while the number of segments varied, the number of candidates per segment

4.2 Robustness to segmentation The robustness of the algorithm to segmentation was evaluated by using a number of polygonal approximations of the contours with dierent tolerances.

Figure 6: Analysis of the robustness to segmentation: matching error vs. approximation error. Bottom: example of reconstruction of approximated contours.

Figure 7: Complexity analysis : computation time vs. number of contour segments to be matched.

Figure 9: Reconstruction after polygonal approximation of contours (left). Reconstruction with a discrete algorithm (right).

increased with the number of segments (the neighbourhood size for candidate search being constant), so that a quadratic complexity could have been expected.

to the one obtained with the subpixel algorithm, even when the matching is performed on a polygonal approximation of contours.

5 Real images

5.1 Simple scene

The algorithm has been applied to a simple stereo pair of a Rubik's cube. Figure 8 shows a full resolution 3D reconstruction, whereas Figure 9 shows a reconstruction obtained by matching a polygonal approximation of the contours. As can be observed, the quality of the results are comparable, while the computation time was reduced from 207 seconds in the rst case to 16 seconds in the second case, a gain of more than an order of magnitude. Figure 9 also shows the result obtained with a classical discrete algorithm [6]: the 3D reconstruction is strongly aected by discretization noise and the resolution is inferior

Figure 8: 3D reconstruction of a Rubik's cube

5.2 Complex scene This next example demonstrates the high resolution 3D reconstruction that can be achieved using the subpixel matching algorithm on complex real images. Figure 10 shows a view of an indoor scene, the reconstruction obtained with the subpixel matching algorithm, and the reconstruction obtained with a discrete algorithm [6]. As can be observed, the jagged contours produced by the discrete algorithm are much smoother with the subpixel algorithm.

6 Conclusion In this paper, we have described a novel algorithm for the optimal matching of contours. Based on a contour deformation measure, this scheme exploits to a full extent the geometric information contained in the contour. No additional constraints are needed, although any available constraint can be used to reduce the search space and hence the computational complexity. This characteristic yields a general-purpose algorithm, suitable for any task where elastic matching is needed. Furthermore, due to its continuous nature, the algorithm has the ability to produce subpixel results, and can run unchanged on contour segments or contour chains. By feeding the algorithm with a polygonal approximation of contours, the computational complexity can be reduced by an order of magnitude while preserving the quality of the results. An extensive experimental analysis has been carried out, demonstrating the precision, the robustness

to segmentation, and the linear computational complexity of the algorithm. An application to the subpixel 3D reconstruction of stereo images has been presented. Excellent results on complex real images are obtained, and comparisons with the results obtained on the same images with a discrete algorithm show the gains in resolution that can be achieved. Extension of the algorithm to the matching of 3D contours is now underway. The theoretical framework presented in this paper can be applied to the 3D case in a straightforward manner, and preliminary results relative to the elastic matching of 3D medical images have been very encouraging.

References

Figure 10: Oce scene: left image (top), subpixel 3D reconstruction looking from above (middle), discrete 3D reconstruction (bottom)

[1] Amir A. Amini, Terry E. Weymouth, and Ramesh C. Jain. Using dynamic programming for solving variational problems in vision. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(9):855867, sept 1990. [2] Ruzena Bajcsy and Stane Kovacic. Multiresolution elastic matching. Computer Vision, Graphics and Image, 46:121, 1989. [3] D.J. Burr. A dynamic model for image registration. Computer Graphics and Image, 15:102112, 1981. [4] S. Dreyfus. Dynamic Programming and the Calculus of Variations. Academic Press, New-York, 1965. [5] E.C. Hildreth. The Measurement of Visual Motion. The MIT Press, 1983. [6] R. Ma and M. Thonnat. A robust and ecient contour-based stereo matching algorithm. Research Report 1860, INRIA, 1993. [7] Mehran Moshfeghi. Elastic matching of multimodality medical images. CVGIP: Graphical Models and Image Processing, 53(3):271282, may 1991. [8] S. B. Pollard, J. E. W. Mayhew, and J. P. Frisby. PMF - a stereo correspondence algorithm using a disparity gradient limit. Perception, 14(4):449 470, 1985. [9] Bruno Serra and Marc Berthod. Subpixel contour matching using continuous dynamic programming. In Conference on Computer Vision and Pattern Recognition, pages 202207, Seattle, 1994.