- Email: [email protected]
Comparison of smooth polynomial functional surfaces for use in alignment
ELSEVIER
Image and Vision Computing
Comparison of smooth polynomial D.S. Drpurtmmt
15 (1997) 529-534
functional surfaces for use in alignment
Meek,
D.J.
Walton
of Computer Science, University of Munitoha, Received 26 April 1996; accepted
Winnipeg, Munitoba,
Canclda
27 November 1996
Abstract Two functional surfaces can be compared through an integral that gives a measure of matching between the two surfaces. If each surface is represented as a polynomial function of two variables, the integral can be done mathematically. This leads to a dramatlc decrease in the computing time required for the comparison of two surfaces. The alignment of two surfaces is the process of transforming one of the surfaces until the measure of matching between the transformed surface and the other surface is optimized. The reduction in the time for comparison leads to a reduction in the time required for alignment. 0 1997 Elsevier Science B.V. Keywords:
Comparison
of smooth surfaces; Alignment;
Registration
1. Introduction Two functional integral
that gives
surfaces can be compared through an a measure of matching between the two
surfaces. If each surface is represented as a polynomial function of two variables, the integral can be done mathematically. This leads to a dramatic decrease in the computing time, compared to using numerical methods, required for the comparison of two surfaces (see Table 2). The alignment of two surfaces is the process of applying rigid transformations to one of the surfaces until it matches the other surface as well as possible; that is, until the measure of matching between the transformed surface and the other surface is optimized. The alignment can be done automatically using a program that minimizes a function of several variables. The comparison of two surfaces is generally the most time-consuming part of the alignment process, so a faster comparison method is very important. The use of exact formulas results in a smoothly varying measure-of-matching function. This is important because the n-variable minimization process works better on a smoothly varying function. A good survey of methods for aligning or registering images, including a discussion on comparing images, has appeared recently [l]. Many of the methods for comparing images apply directly to the comparison of surfaces. The methods used here differ from two commonly used methods for surface comparison. One method for comparing surfaces depends on the identification of intrinsic features such as sharp edges or points of maximum curvature [2-81. Another 0262.8856/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved PII SO262-8X56(96)00002-9
method for comparing surfaces depends on the fact that identifiable marks appear on both surfaces and those marks give a set of corresponding points [9]. The method
proposed here makes a comparison without searching for features and without using identifiable marks. This is similar to ideas proposed recently, but differs in that here continuous surfaces are aligned rather than surfaces represented by dense range data [ 10-121, and here translations and rotations are considered rather than just translations 1131. The following four sections detail the way in which alignment is done in this paper. Section 2 describes the type of surfaces that are to be aligned. Section 3 describes the transformations that will be applied to one of the surfaces. Section 4 describes the method of comparing two surfaces. Section 5 describes the alignment process and gives an example. Conclusions are presented in Section 6. An expanded version of this paper is given in a Technical Report [ 141.
2. Surface representation The two surfaces to be aligned are of the form z = A(x,y) and z = B&y), where both A(x,y) and B&y) are uniform polynomial B-splines. Assume that the boundaries of the patches of each surface, projected onto the xy-plane, form an integral grid of n X n squares. The theory and examples presented here use biquadratic patches, although any degree can be used. The integral square grid of patch boundary lines simplifies the problem of determining which patches in one surface lie above or below a given patch in the other
530
D.S. Meek, D.J. Walton/image
and Vision Computing 15 (1997) 529-534
surface. It also simplifies the organization and many of the computations; for example, the index of a patch of surface A is the same as the (x,y) coordinates of the lower left corner of its square boundary. Without loss of generality, the patch boundary grid of surface A is aligned with the x- and y-axes; the boundary lines are parallel to, and an integral distance from, the axes, and the lower left corner of the grid is at the origin. The boundary grid of surface B is not necessarily aligned with the x- and y-axes. Both B-spline surfaces are expressed in terms of their B-spline control points. Let the (n + 2)* control points for surface A be (n + 0.5, - 0.5, (- 0.5, - OS,C,,e), (0.5, - OS,C,,e) )...) C n+,,O), (0.5,0.5,CO,,), . . . . (n + 0.5,0.5,C,,l,1),..., ). The (ij) patch boundary square (n + 0.572 + o.5,c,+l,n+l of surface A, 0 9 i, j 5 n - 1, is the square in the xy-plane whose corners have (x, y) coordinates (ij), (i + 1j), (ij + l), (i + lj + 1). The (ij) patch of surface A is represented parametrically as
with a different polynomial over each boundary square. Because the surfaces are biquadratic B-splines, neighbouring polynomials will join with C’ continuity. The surfaces are represented in terms of a common coordinate system before comparison. For the (ij) patch of surfaceA,letx=i+sandy=j+t,orreplacesbyx-iand replace t by y -j in the above formula (1) for Ai,Js,t). In terms of (x,y), the (ij) patch of.surface A is the polynomial Aij(x - i, y -j). For the (k, I) patch of surface B, let x = TX + (k + s)cos 6 - (1+ t)sin 0 and let y = TY + (k + S) sin 0 + (1+ t)cos 0. Solving for s and t in terms of x and y,
The (k,l) patch of surface B can be written as a polynomial in terms of (x,y) by substituting eqn (3) into eqn (2) for Bk,J(s,t).
3. Surface transformation
i+s Surface A is not transformed; surface B is transformed as follows: it is rotated about its local origin by 8, then its local origin is translated to (T,, T,) (see Fig. 1). The rotation 0 is restricted to the range (- r/2,&2). This simplifies the matching by making it easier to find which boundary squares of B overlap a given boundary square of A.
,OSs,t
j+t
Ai,j(s, t) =
( Ai,j(s, t) i where Ai,j(s, t) = i(( 1 - s)~, 1 + 2~ - 2~~7s2) ci,j+
ci,j X
1
(1 - t)2
ci,ji-2
4. Surface comparison
ci+l.j
ci+l,j+I
Ci+l,j+2
1+ 2t - 2t2
ci+2,j
ci-t2,j+1
Ci+2,j+2
t2
and the Ci,j are the z-values of the control points of surface A. Assume the patch boundary grid of surface B has been rotated by an angle of 8 about its local origin and its local origin has been translated to (T,, T,). The (k,Z) patch of the surface B is represented parametrically as
To compare two surfaces, it is useful to define a meaning of corresponding points between the two surfaces. Here points are said to correspond if they are vertically aligned; that is, point (x,y,A(x,y)) in surface A corresponds to point (x,y,B(x,y)) in surface B. This is a fairly good way to
T, + (k + s)cos 13- (1+ t)sin 0
Bk,L-Qt) =
TY + (k + s)sin 0 + (1+ t)cos 6
&, lb, t)
, OSs,t<
1,
/
where B,,,(s, t) = ;((l
- S)2, 1 + 2.9- 2S2, S2) C k,l+l
ck,i
X
ck, ck+2,1
1.1
ck,l,2
C k+l,[+l
ck,
ck+2,1+1
ck,2,1,2
1,1+2
(1 - t)2 1 + 2t - 2t2 t2
(2) and the Ckl are the z-values of the control points of surface B. Notice that each surface consists of rr* patches generally
(T,,TJ Fig.
1. Overlapping boundary squares.
D.S. Meek, D.J. Walton/Image
and Vision Compufing
associate points between two surfaces whose normals are pointing roughly in the z-direction. The comparison of two surfaces will be summarized by a single number called the measure of matching. The measure of matching used here is the integral over the overlapping area of the square of the differences in the z-values of the two surfaces. This measure is large when the surfaces are far apart and small when they are close, so the goal will be to minimize the measure of matching. The measure of matching is normalized by dividing by the areaofoverlap so that aclose fit with a tiny overlap does not appear to be a good overall match. Generally, the patch boundary grids of the two surfaces will not be aligned (see Fig. 1). This complicates the calculation of the measure of matching because the regions over which the difference of the two piecewise surfaces is given by non-piecewise, mathematical formulae have to be identified. These regions are the polygons that are the intersections of a boundary square of surface A and a boundary square of surface B. For example, the polygonal region that is the intersection of boundary square (1,l) of A and boundary square (2,l) of B is shaded in Fig. 1. The measure of matching over this polygonal region is II
15 (1997) 529-534
4.1.1. Case I: Assuming (5)
531
0 2 0,
x cos e + (7) x sin e gives (i - T,)cos 8 + Cj - T,)sin r3< k + S,
(6) X cos 0 + (8) X sin 0 gives k + s < (i t- 1 - lrr)cos e + lj + 1 - T,)sin 8, (7) x cos t9+ (6) x (- sin 0) gives (‘j- T,,)cos e - (i + 1 - T~r)sin 0 < 1+ t, and (8)Xcos13+(5)X(-sin@givesI+t<(j+l--T,)cose - (i - T,)sin 8.
4.1.2. Case 2: assuming 0 < 0, (5) x cos 8 + (8) x sin 0 gives (i - T,)COS e + fj + 1 - T,)sin 0 < k + S,
(A I I (~7 Y) - &I (~9 y:02 dx dy,
where the double integral is over the shaded polygonal region, and A ,,(x,y) and B2,(x,y) are the appropriate nonpiecewise polynomials from surfaces A and B. The integral can be found by exact mathematical formulae as shown below and in the Appendix.
+ (j - T,)sin 8, (7) x cos e + (5)
x ( - sin e) gives 0’ - T,)COS e
- (i - T,)sin 0 < 1+ t, and (8) X cos 0 + (6) x ( - sin 0) gives 1+ t < (j + 1 -
4.1. Finding which patch boundary squares overlap
-
Listing all the polygonal regions that are the intersection of a boundary square ofA and a boundary square of B can be difficult. One organized procedure, given a boundary square of A, is to look through all the boundary squares of B to see which ones overlap it. A more efficient way to find which boundary squares of B overlap boundary square (ij) of A follows. The (k,Z) boundary square of surface B is the set of points (;)
+(k+s)(;oss”)
-_(I+))(
+;;),
OSs,t<
1.
The integral values of (k,l) such that i
(3
T,+(k+s)cos0-(1+r)sinf3
(6)
j
(7)
and TY+ (k + s)sin 8 + (1+ t)cos e < j + 1
T,)COS
8
(i + 1 - TJsin 8.
Considering the first two inequalities in case 1, the bounding values on k + s can be thought of as real numbers consisting of integer and fraction parts. If k + s equals one of the bounds, k will be the integer part and s will be the fraction part. The values that k may have are floor[(i - TJcos 0 + 0’ - Ty)sin f3] 5 k 5 floor[(i + 1 - T,) x cos 8 + tj + 1 - T,)sin e]. The range of 1 values for case 1 and the range of (k,l) values for case 2 can be found in a similar way. 4.2. Integrating squares
over the intersection
of two patch boundary
Suppose the (ij) boundary square of surface A is being considered and it was determined from the above bounds that the (k,l) boundary square of surface B overlaps it. There are two main cases: the boundaries of the grid square of B are parallel to the axes, and the boundaries are not parallel to the axes.
(8)
are required. Assume that -a/2 < i9 < 1rl2, so that cos 0 > 0. There are two cases: 0 2 0 (so that sin B 2 0), and e < 0 (so that sin 0 < 0).
4.2.1. Case I: parallel, 8 = 0 Let the coordinates of the (k,l) boundary square of B be (L,, BY), (R,,B,), (L,,T,), (R,,T,), where L stands for left, R
D.S. Meek, D.J. Walton/Image and Vision Computing 15 (1997) 55?-534
532
T = (T,,T,)
(LvTy)
(R,vTY)
1
I L
(L,vB,)
1
I
(R,,By)
I
1
B = (B,.By)I
Fig. 2. Overlap when the grid of surface B is parallel and not parallel to the axes.
stands for right, T stands for top, and B stands for bottom (see Fig. 2). The measure of matching is the integral under the top line (L,,T,) to (R,,T,) minus the integral under the bottom line (L,,B,) to (K,,B$, where the (x,y) values are limited to the (ij) boundary square of surface A. 4.2.2. Case 2: not parallel, 8 # 0 For the (k,Z) boundary square of surface B, identify the left vertex L = (L,,LY), right vertex R = (Rx,R,), bottom vertex B = (B,,BJ, and the top vertex T = (T,,T,) (see Fig. 2). The measure of matching is the integral under line LT plus the integral under line TR minus the integral under line LB minus the integral under line BR, where the (x,y) values are limited to the (ij) boundary square of surface A. 4.3. Integrating
under a line
Care must be taken in integrating to keep the (x,y) points in the (ij) boundary square of surface A. This gives rise to many cases. Let y = L(x) be the line that joins a starting point S = (S,,S,) to an ending point E = (E,,E,). Define xsi i
intersects the line y = j + 1. If x, I I,, the measure of matching is s: si/’ ’ P(x, y) dy dx. Otherwise, Z, < x, and the measure of matching is Ji; SJ’” P(x, y) dy dx + c J’(x) P(xt Y) dy dx ( see Fig. 3 for the above cases). 4.3.2. Case 2: j < L(x,) < j + I If the slope of the line y = L(x) is 0, the measure of matching is s: s/L@,)P(x, y) dy dx. If the slope of the line y = L(x) is > 0, find the point (Z,j + 1) where y= L(x) intersects the line y = j + 1. If x, % I,, the measure of matching is J;Te&‘@)P(x,y) dy dx. Otherwise, Z, < x, and the measure of matching is sk l”” P(x, y) dy dx + c J’ ’ P(x, y) dy dx. If the slope of the line y = L(x) is < 0, find the point (Z,j) where y = L(x) intersects the line Z,, the measure of matching is y = j. If x, I s’ &‘(‘)P(x, y) dy dx. Otherwise, Z, < x, and the measure of matching is si; l”“’ P(x, y) dy Q. 4.3.3. Case 3: L(x,) 5 j If the slope of the line y = L(x) is 5 0, the measure of matching is zero. If the slope of the line y = L(x) is > 0, find the point (ZJ) where y = L(x) intersects the line y =j. If x, 5 I,, the measure of matching is zero. Otherwise, I, < x, and the measure of matching is t 6’” P(x, y) dy dx.
i+lsx
then x, = bound(S,) and x, = bound(E,) are the limits of integration along the x-axis. For convenience, let P(x, y) = (Ajj(x, y) - Bk&x,y))2, the square of the differences in z-values of the two surfaces. The position of L(x,) gives three cases. 4.3.1. Case 1: j + I 5 L(x,) If the slope of the line y = L(x) B 0, the measure of matching is j’ 6j+’ P(x y) dy dx. If the slope of the line y = L(x) is < 0, find the point (Z,j + 1) where y= L(x)
5. Alignment of surfaces The transformation of surface B can be thought of as a function of three variables: rotation about the z-axis, the translation in the x-direction, and the translation in the y-direction. Consequently, the measure of matching of surface A to the transformed surface B can be thought of as a function of three variables. The goal is to find the transformation that minimizes the measure of matching. Powell’s method of
xs xe Fig. 3. Integral under a line.
D.S. Meek, D. .I. Waltonflmagr
and Vi.sion Computing
15
(I YY7)
533
529-534
was created (the same values as in surface A except that the middle value is increased by 1). This surface was then rotated by 0.2 radians and translated by 0.5 in the X- and y-directions to give surface B. Surface B is represented by light grid and contour lines in Fig. 4. The minimization program automatically chooses transformations in an attempt to minimize the measure of matching between the two surfaces. One iteration of Powell’s method is considered to be the result of one linear search in each of the n search directions. The results of the alignment are given in Table I. Fig. 5 shows the tinal position of the best match. Timings were done to compare the speed of alignment when exact integrals are used and when a numerical method (rectangle rule) is used for the integrals. The results in Table 2 show that the numerical method (assuming 60 function values between grid lines in Fig. 4, which gives 180 X 180 function values) takes about one hundred times as long as the exact integration. Using numerical methods, the integration part of the alignment process is 98.3% of the total time. Fig. 4. Original
positions:
surface
A dark,
surface
B light.
6. Conclusions minimizing a function of n-variables [ 151 was used to vary the transformation of surface B until the best possible match between surface A and the transformed surface B was found. Suppose surface A has control points with the following r-values: 1
2
3
I
1
3
8
7
5
2
2
3
4
3
2
1
3
3
2
1
Polynomial surfaces can be compared through exact integration. The comparison with exact integration can be considerably faster than with numerical methods, even though the computer program with exact integration is more complicated. Alignment can be performed much faster using comparison with exact integration. It is expected that the optimization algorithm also works better with exact integration because the function to optimize is smoother.
Acknowledgements
Surface A is represented by dark grid and contour lines in Fig. 4. A surface which has control points with the following z-values:
The authors acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada for this research. The authors also acknowledge the helpful comments of two anonymous referees.
12311 3
8
7
5
2
3
5
7
4
3
Appendix
2
3
4
3
2
1
3
3
2
1
This appendix gives formulae for the exact double integration of a polynomial over the type of region that is
Table Results Iteration
I
I of the optimization
program Translation
0.000
in x
Translation
0.000
2
-0.256
-0.714
in y
Angle
Measure
of matching
0.000
1.1296
X
0.016
3.9203
X IO -’
0.004
IO”
3
-0.290
-0.739
3.6228
X IO-’
4
-0.433
-0.576
-0.124
2.9244
X IO-’
5
-0.507
-0.569
-0.164
2.3975
X IO-’
6
-0.516
-0.567
-0.169
2.3918
x
IOY’
D.S. Meek, D.J. Waltodlmage and Vision Computing 15 (1997) 529-534 li,j =
-Li+j+2
Xl
x
[x’+‘(ax+by’+l]~~+O.+l)b
{ = &{x;+$xl
x’(ux+b)idx I X0
+b)i+‘-X;+l(aXo+b)i+’
l)bzi,j-l}
+G+
An initial value for this recurrence
642)
is Ii,_1
Xi+1 _,;+I zj, _ 1 =
Fig. 5. Final positions;
surface A dark, transformed
surface B light.
Table 2 Timings for optimization
38.18 4032.92 69.72
required by this paper. Let P&y)
= 9 F* AiixiY’. i=Oj=O
The double integral
where XI Z. .= ‘3J
X’(ax+b)i+‘dr. I X0
A recurrence for Ii,, i B 0, j 2 integration by parts of
i+l
One can define Zi,_z = 0 and calculate Ii,_, from (A2). The values I,r,~,I.r,l,...,Zi,n that are needed for the inner sum of (Al) can be calculated from (A2).
References Time (s)
Time for exact integrals Time for numerical integrals (60 values per grid unit) Time for the program other than integration
l
- 1, is obtained
from
[l] L.G. Brown, A survey of registration techniques, Computing Surveys 29 (1992) 325-376. [2] Y. Chen, G. Medioni, Object modelling by registration of multiple range images, Image and Vision Computing 10 (1992) 145-155. [3] P. Gerlot-Chiron, Y. Bizais, Registration of multimodality medical images using a region overlap criterion, CVGIP: Graphical Models and Image Proc. 54 (1992) 396-406. [4] S.Z. Li, Matching: invariant to translations, rotations, and scale changes, Pattern Recognition 25 (1992) 583-594. [5] G. Medioni, A. Huertas, M. Wilson, Automatic registration of color separation films, Machine Vision and Applications 4 (1991) 3 l-5 1. [6] G. Medioni, R. Nevatia, Matching images using linear features, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. PAMI(1984) pp. 675-685. [7] J. Weng, N. Ahuja, T.S. Huang, Matching two perspective views, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. PAMI(1992) pp. 806-825. [8] Y. Zhang, J.J. Gerbrands, Method for matching general stereo planar curves, Image and Vision Computing 13 (1995) 645-655. [9] R. DeLong, M.C. Bramwell, Matching an initial surface to a worn surface subject to rotation and translation, in: R.R. Martin (Ed.) Mathematics of Surfaces II, Clarendon, Oxford, 1987, pp. 137-149. [IO] M. Herbin, A. Venot, J.Y. Devaux, E. Walter, J.F. Lebruchec, L. Dubertret, J.C. Roucayrol, Automated registration of dissimilar images: application to medical imagery, Computer Vision, Graphics and Image Processing 47 (1989) 77-88. [ 1 l] B. Kamgar-Pa&, J.L. Jones, A. Rosenfeld, Registration of multiple overlapping range images without distinctive features, IEEE Computer Sot. Conf. on Computer Vision and Pattern Recognition, 4-8 June 1989, San Diego, CA, pp. 2X2-290. [12] D.S. Meek, D.J. Walton, Surface alignment or registration of multiple range images, Congressus Numerantium 99 ( 1994) 17 I- 177. [13] K. Deguchi, Registration technique for partially covered image sequence, IEEE Proc. 8th Int. Conf. on Pattern Recognition, Paris, October 27-31, 1986, pp. 1186-1189. [ 141 D.S. Meek, D.J. Walton, Comparison of smooth polynomial functional surfaces for use in alignment, Technical Report 96/I 1, Department of Computer Science, University of Manitoba, Winnipeg, Canada, 1996. [15] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes in C, 2nd Edn., Cambridge University Press, Cambridge, UK, 1992.
Image and Vision Computing
Comparison of smooth polynomial D.S. Drpurtmmt
15 (1997) 529-534
functional surfaces for use in alignment
Meek,
D.J.
Walton
of Computer Science, University of Munitoha, Received 26 April 1996; accepted
Winnipeg, Munitoba,
Canclda
27 November 1996
Abstract Two functional surfaces can be compared through an integral that gives a measure of matching between the two surfaces. If each surface is represented as a polynomial function of two variables, the integral can be done mathematically. This leads to a dramatlc decrease in the computing time required for the comparison of two surfaces. The alignment of two surfaces is the process of transforming one of the surfaces until the measure of matching between the transformed surface and the other surface is optimized. The reduction in the time for comparison leads to a reduction in the time required for alignment. 0 1997 Elsevier Science B.V. Keywords:
Comparison
of smooth surfaces; Alignment;
Registration
1. Introduction Two functional integral
that gives
surfaces can be compared through an a measure of matching between the two
surfaces. If each surface is represented as a polynomial function of two variables, the integral can be done mathematically. This leads to a dramatic decrease in the computing time, compared to using numerical methods, required for the comparison of two surfaces (see Table 2). The alignment of two surfaces is the process of applying rigid transformations to one of the surfaces until it matches the other surface as well as possible; that is, until the measure of matching between the transformed surface and the other surface is optimized. The alignment can be done automatically using a program that minimizes a function of several variables. The comparison of two surfaces is generally the most time-consuming part of the alignment process, so a faster comparison method is very important. The use of exact formulas results in a smoothly varying measure-of-matching function. This is important because the n-variable minimization process works better on a smoothly varying function. A good survey of methods for aligning or registering images, including a discussion on comparing images, has appeared recently [l]. Many of the methods for comparing images apply directly to the comparison of surfaces. The methods used here differ from two commonly used methods for surface comparison. One method for comparing surfaces depends on the identification of intrinsic features such as sharp edges or points of maximum curvature [2-81. Another 0262.8856/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved PII SO262-8X56(96)00002-9
method for comparing surfaces depends on the fact that identifiable marks appear on both surfaces and those marks give a set of corresponding points [9]. The method
proposed here makes a comparison without searching for features and without using identifiable marks. This is similar to ideas proposed recently, but differs in that here continuous surfaces are aligned rather than surfaces represented by dense range data [ 10-121, and here translations and rotations are considered rather than just translations 1131. The following four sections detail the way in which alignment is done in this paper. Section 2 describes the type of surfaces that are to be aligned. Section 3 describes the transformations that will be applied to one of the surfaces. Section 4 describes the method of comparing two surfaces. Section 5 describes the alignment process and gives an example. Conclusions are presented in Section 6. An expanded version of this paper is given in a Technical Report [ 141.
2. Surface representation The two surfaces to be aligned are of the form z = A(x,y) and z = B&y), where both A(x,y) and B&y) are uniform polynomial B-splines. Assume that the boundaries of the patches of each surface, projected onto the xy-plane, form an integral grid of n X n squares. The theory and examples presented here use biquadratic patches, although any degree can be used. The integral square grid of patch boundary lines simplifies the problem of determining which patches in one surface lie above or below a given patch in the other
530
D.S. Meek, D.J. Walton/image
and Vision Computing 15 (1997) 529-534
surface. It also simplifies the organization and many of the computations; for example, the index of a patch of surface A is the same as the (x,y) coordinates of the lower left corner of its square boundary. Without loss of generality, the patch boundary grid of surface A is aligned with the x- and y-axes; the boundary lines are parallel to, and an integral distance from, the axes, and the lower left corner of the grid is at the origin. The boundary grid of surface B is not necessarily aligned with the x- and y-axes. Both B-spline surfaces are expressed in terms of their B-spline control points. Let the (n + 2)* control points for surface A be (n + 0.5, - 0.5, (- 0.5, - OS,C,,e), (0.5, - OS,C,,e) )...) C n+,,O), (0.5,0.5,CO,,), . . . . (n + 0.5,0.5,C,,l,1),..., ). The (ij) patch boundary square (n + 0.572 + o.5,c,+l,n+l of surface A, 0 9 i, j 5 n - 1, is the square in the xy-plane whose corners have (x, y) coordinates (ij), (i + 1j), (ij + l), (i + lj + 1). The (ij) patch of surface A is represented parametrically as
with a different polynomial over each boundary square. Because the surfaces are biquadratic B-splines, neighbouring polynomials will join with C’ continuity. The surfaces are represented in terms of a common coordinate system before comparison. For the (ij) patch of surfaceA,letx=i+sandy=j+t,orreplacesbyx-iand replace t by y -j in the above formula (1) for Ai,Js,t). In terms of (x,y), the (ij) patch of.surface A is the polynomial Aij(x - i, y -j). For the (k, I) patch of surface B, let x = TX + (k + s)cos 6 - (1+ t)sin 0 and let y = TY + (k + S) sin 0 + (1+ t)cos 0. Solving for s and t in terms of x and y,
The (k,l) patch of surface B can be written as a polynomial in terms of (x,y) by substituting eqn (3) into eqn (2) for Bk,J(s,t).
3. Surface transformation
i+s Surface A is not transformed; surface B is transformed as follows: it is rotated about its local origin by 8, then its local origin is translated to (T,, T,) (see Fig. 1). The rotation 0 is restricted to the range (- r/2,&2). This simplifies the matching by making it easier to find which boundary squares of B overlap a given boundary square of A.
,OSs,t
j+t
Ai,j(s, t) =
( Ai,j(s, t) i where Ai,j(s, t) = i(( 1 - s)~, 1 + 2~ - 2~~7s2) ci,j+
ci,j X
1
(1 - t)2
ci,ji-2
4. Surface comparison
ci+l.j
ci+l,j+I
Ci+l,j+2
1+ 2t - 2t2
ci+2,j
ci-t2,j+1
Ci+2,j+2
t2
and the Ci,j are the z-values of the control points of surface A. Assume the patch boundary grid of surface B has been rotated by an angle of 8 about its local origin and its local origin has been translated to (T,, T,). The (k,Z) patch of the surface B is represented parametrically as
To compare two surfaces, it is useful to define a meaning of corresponding points between the two surfaces. Here points are said to correspond if they are vertically aligned; that is, point (x,y,A(x,y)) in surface A corresponds to point (x,y,B(x,y)) in surface B. This is a fairly good way to
T, + (k + s)cos 13- (1+ t)sin 0
Bk,L-Qt) =
TY + (k + s)sin 0 + (1+ t)cos 6
&, lb, t)
, OSs,t<
1,
/
where B,,,(s, t) = ;((l
- S)2, 1 + 2.9- 2S2, S2) C k,l+l
ck,i
X
ck, ck+2,1
1.1
ck,l,2
C k+l,[+l
ck,
ck+2,1+1
ck,2,1,2
1,1+2
(1 - t)2 1 + 2t - 2t2 t2
(2) and the Ckl are the z-values of the control points of surface B. Notice that each surface consists of rr* patches generally
(T,,TJ Fig.
1. Overlapping boundary squares.
D.S. Meek, D.J. Walton/Image
and Vision Compufing
associate points between two surfaces whose normals are pointing roughly in the z-direction. The comparison of two surfaces will be summarized by a single number called the measure of matching. The measure of matching used here is the integral over the overlapping area of the square of the differences in the z-values of the two surfaces. This measure is large when the surfaces are far apart and small when they are close, so the goal will be to minimize the measure of matching. The measure of matching is normalized by dividing by the areaofoverlap so that aclose fit with a tiny overlap does not appear to be a good overall match. Generally, the patch boundary grids of the two surfaces will not be aligned (see Fig. 1). This complicates the calculation of the measure of matching because the regions over which the difference of the two piecewise surfaces is given by non-piecewise, mathematical formulae have to be identified. These regions are the polygons that are the intersections of a boundary square of surface A and a boundary square of surface B. For example, the polygonal region that is the intersection of boundary square (1,l) of A and boundary square (2,l) of B is shaded in Fig. 1. The measure of matching over this polygonal region is II
15 (1997) 529-534
4.1.1. Case I: Assuming (5)
531
0 2 0,
x cos e + (7) x sin e gives (i - T,)cos 8 + Cj - T,)sin r3< k + S,
(6) X cos 0 + (8) X sin 0 gives k + s < (i t- 1 - lrr)cos e + lj + 1 - T,)sin 8, (7) x cos t9+ (6) x (- sin 0) gives (‘j- T,,)cos e - (i + 1 - T~r)sin 0 < 1+ t, and (8)Xcos13+(5)X(-sin@givesI+t<(j+l--T,)cose - (i - T,)sin 8.
4.1.2. Case 2: assuming 0 < 0, (5) x cos 8 + (8) x sin 0 gives (i - T,)COS e + fj + 1 - T,)sin 0 < k + S,
(A I I (~7 Y) - &I (~9 y:02 dx dy,
where the double integral is over the shaded polygonal region, and A ,,(x,y) and B2,(x,y) are the appropriate nonpiecewise polynomials from surfaces A and B. The integral can be found by exact mathematical formulae as shown below and in the Appendix.
+ (j - T,)sin 8, (7) x cos e + (5)
x ( - sin e) gives 0’ - T,)COS e
- (i - T,)sin 0 < 1+ t, and (8) X cos 0 + (6) x ( - sin 0) gives 1+ t < (j + 1 -
4.1. Finding which patch boundary squares overlap
-
Listing all the polygonal regions that are the intersection of a boundary square ofA and a boundary square of B can be difficult. One organized procedure, given a boundary square of A, is to look through all the boundary squares of B to see which ones overlap it. A more efficient way to find which boundary squares of B overlap boundary square (ij) of A follows. The (k,Z) boundary square of surface B is the set of points (;)
+(k+s)(;oss”)
-_(I+))(
+;;),
OSs,t<
1.
The integral values of (k,l) such that i
(3
T,+(k+s)cos0-(1+r)sinf3
(6)
j
(7)
and TY+ (k + s)sin 8 + (1+ t)cos e < j + 1
T,)COS
8
(i + 1 - TJsin 8.
Considering the first two inequalities in case 1, the bounding values on k + s can be thought of as real numbers consisting of integer and fraction parts. If k + s equals one of the bounds, k will be the integer part and s will be the fraction part. The values that k may have are floor[(i - TJcos 0 + 0’ - Ty)sin f3] 5 k 5 floor[(i + 1 - T,) x cos 8 + tj + 1 - T,)sin e]. The range of 1 values for case 1 and the range of (k,l) values for case 2 can be found in a similar way. 4.2. Integrating squares
over the intersection
of two patch boundary
Suppose the (ij) boundary square of surface A is being considered and it was determined from the above bounds that the (k,l) boundary square of surface B overlaps it. There are two main cases: the boundaries of the grid square of B are parallel to the axes, and the boundaries are not parallel to the axes.
(8)
are required. Assume that -a/2 < i9 < 1rl2, so that cos 0 > 0. There are two cases: 0 2 0 (so that sin B 2 0), and e < 0 (so that sin 0 < 0).
4.2.1. Case I: parallel, 8 = 0 Let the coordinates of the (k,l) boundary square of B be (L,, BY), (R,,B,), (L,,T,), (R,,T,), where L stands for left, R
D.S. Meek, D.J. Walton/Image and Vision Computing 15 (1997) 55?-534
532
T = (T,,T,)
(LvTy)
(R,vTY)
1
I L
(L,vB,)
1
I
(R,,By)
I
1
B = (B,.By)I
Fig. 2. Overlap when the grid of surface B is parallel and not parallel to the axes.
stands for right, T stands for top, and B stands for bottom (see Fig. 2). The measure of matching is the integral under the top line (L,,T,) to (R,,T,) minus the integral under the bottom line (L,,B,) to (K,,B$, where the (x,y) values are limited to the (ij) boundary square of surface A. 4.2.2. Case 2: not parallel, 8 # 0 For the (k,Z) boundary square of surface B, identify the left vertex L = (L,,LY), right vertex R = (Rx,R,), bottom vertex B = (B,,BJ, and the top vertex T = (T,,T,) (see Fig. 2). The measure of matching is the integral under line LT plus the integral under line TR minus the integral under line LB minus the integral under line BR, where the (x,y) values are limited to the (ij) boundary square of surface A. 4.3. Integrating
under a line
Care must be taken in integrating to keep the (x,y) points in the (ij) boundary square of surface A. This gives rise to many cases. Let y = L(x) be the line that joins a starting point S = (S,,S,) to an ending point E = (E,,E,). Define xsi i
intersects the line y = j + 1. If x, I I,, the measure of matching is s: si/’ ’ P(x, y) dy dx. Otherwise, Z, < x, and the measure of matching is Ji; SJ’” P(x, y) dy dx + c J’(x) P(xt Y) dy dx ( see Fig. 3 for the above cases). 4.3.2. Case 2: j < L(x,) < j + I If the slope of the line y = L(x) is 0, the measure of matching is s: s/L@,)P(x, y) dy dx. If the slope of the line y = L(x) is > 0, find the point (Z,j + 1) where y= L(x) intersects the line y = j + 1. If x, % I,, the measure of matching is J;Te&‘@)P(x,y) dy dx. Otherwise, Z, < x, and the measure of matching is sk l”” P(x, y) dy dx + c J’ ’ P(x, y) dy dx. If the slope of the line y = L(x) is < 0, find the point (Z,j) where y = L(x) intersects the line Z,, the measure of matching is y = j. If x, I s’ &‘(‘)P(x, y) dy dx. Otherwise, Z, < x, and the measure of matching is si; l”“’ P(x, y) dy Q. 4.3.3. Case 3: L(x,) 5 j If the slope of the line y = L(x) is 5 0, the measure of matching is zero. If the slope of the line y = L(x) is > 0, find the point (ZJ) where y = L(x) intersects the line y =j. If x, 5 I,, the measure of matching is zero. Otherwise, I, < x, and the measure of matching is t 6’” P(x, y) dy dx.
i+lsx
then x, = bound(S,) and x, = bound(E,) are the limits of integration along the x-axis. For convenience, let P(x, y) = (Ajj(x, y) - Bk&x,y))2, the square of the differences in z-values of the two surfaces. The position of L(x,) gives three cases. 4.3.1. Case 1: j + I 5 L(x,) If the slope of the line y = L(x) B 0, the measure of matching is j’ 6j+’ P(x y) dy dx. If the slope of the line y = L(x) is < 0, find the point (Z,j + 1) where y= L(x)
5. Alignment of surfaces The transformation of surface B can be thought of as a function of three variables: rotation about the z-axis, the translation in the x-direction, and the translation in the y-direction. Consequently, the measure of matching of surface A to the transformed surface B can be thought of as a function of three variables. The goal is to find the transformation that minimizes the measure of matching. Powell’s method of
xs xe Fig. 3. Integral under a line.
D.S. Meek, D. .I. Waltonflmagr
and Vi.sion Computing
15
(I YY7)
533
529-534
was created (the same values as in surface A except that the middle value is increased by 1). This surface was then rotated by 0.2 radians and translated by 0.5 in the X- and y-directions to give surface B. Surface B is represented by light grid and contour lines in Fig. 4. The minimization program automatically chooses transformations in an attempt to minimize the measure of matching between the two surfaces. One iteration of Powell’s method is considered to be the result of one linear search in each of the n search directions. The results of the alignment are given in Table I. Fig. 5 shows the tinal position of the best match. Timings were done to compare the speed of alignment when exact integrals are used and when a numerical method (rectangle rule) is used for the integrals. The results in Table 2 show that the numerical method (assuming 60 function values between grid lines in Fig. 4, which gives 180 X 180 function values) takes about one hundred times as long as the exact integration. Using numerical methods, the integration part of the alignment process is 98.3% of the total time. Fig. 4. Original
positions:
surface
A dark,
surface
B light.
6. Conclusions minimizing a function of n-variables [ 151 was used to vary the transformation of surface B until the best possible match between surface A and the transformed surface B was found. Suppose surface A has control points with the following r-values: 1
2
3
I
1
3
8
7
5
2
2
3
4
3
2
1
3
3
2
1
Polynomial surfaces can be compared through exact integration. The comparison with exact integration can be considerably faster than with numerical methods, even though the computer program with exact integration is more complicated. Alignment can be performed much faster using comparison with exact integration. It is expected that the optimization algorithm also works better with exact integration because the function to optimize is smoother.
Acknowledgements
Surface A is represented by dark grid and contour lines in Fig. 4. A surface which has control points with the following z-values:
The authors acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada for this research. The authors also acknowledge the helpful comments of two anonymous referees.
12311 3
8
7
5
2
3
5
7
4
3
Appendix
2
3
4
3
2
1
3
3
2
1
This appendix gives formulae for the exact double integration of a polynomial over the type of region that is
Table Results Iteration
I
I of the optimization
program Translation
0.000
in x
Translation
0.000
2
-0.256
-0.714
in y
Angle
Measure
of matching
0.000
1.1296
X
0.016
3.9203
X IO -’
0.004
IO”
3
-0.290
-0.739
3.6228
X IO-’
4
-0.433
-0.576
-0.124
2.9244
X IO-’
5
-0.507
-0.569
-0.164
2.3975
X IO-’
6
-0.516
-0.567
-0.169
2.3918
x
IOY’
D.S. Meek, D.J. Waltodlmage and Vision Computing 15 (1997) 529-534 li,j =
-Li+j+2
Xl
x
[x’+‘(ax+by’+l]~~+O.+l)b
{ = &{x;+$xl
x’(ux+b)idx I X0
+b)i+‘-X;+l(aXo+b)i+’
l)bzi,j-l}
+G+
An initial value for this recurrence
642)
is Ii,_1
Xi+1 _,;+I zj, _ 1 =
Fig. 5. Final positions;
surface A dark, transformed
surface B light.
Table 2 Timings for optimization
38.18 4032.92 69.72
required by this paper. Let P&y)
= 9 F* AiixiY’. i=Oj=O
The double integral
where XI Z. .= ‘3J
X’(ax+b)i+‘dr. I X0
A recurrence for Ii,, i B 0, j 2 integration by parts of
i+l
One can define Zi,_z = 0 and calculate Ii,_, from (A2). The values I,r,~,I.r,l,...,Zi,n that are needed for the inner sum of (Al) can be calculated from (A2).
References Time (s)
Time for exact integrals Time for numerical integrals (60 values per grid unit) Time for the program other than integration
l
- 1, is obtained
from
[l] L.G. Brown, A survey of registration techniques, Computing Surveys 29 (1992) 325-376. [2] Y. Chen, G. Medioni, Object modelling by registration of multiple range images, Image and Vision Computing 10 (1992) 145-155. [3] P. Gerlot-Chiron, Y. Bizais, Registration of multimodality medical images using a region overlap criterion, CVGIP: Graphical Models and Image Proc. 54 (1992) 396-406. [4] S.Z. Li, Matching: invariant to translations, rotations, and scale changes, Pattern Recognition 25 (1992) 583-594. [5] G. Medioni, A. Huertas, M. Wilson, Automatic registration of color separation films, Machine Vision and Applications 4 (1991) 3 l-5 1. [6] G. Medioni, R. Nevatia, Matching images using linear features, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. PAMI(1984) pp. 675-685. [7] J. Weng, N. Ahuja, T.S. Huang, Matching two perspective views, IEEE Trans. Pattern Analysis and Machine Intelligence, Vol. PAMI(1992) pp. 806-825. [8] Y. Zhang, J.J. Gerbrands, Method for matching general stereo planar curves, Image and Vision Computing 13 (1995) 645-655. [9] R. DeLong, M.C. Bramwell, Matching an initial surface to a worn surface subject to rotation and translation, in: R.R. Martin (Ed.) Mathematics of Surfaces II, Clarendon, Oxford, 1987, pp. 137-149. [IO] M. Herbin, A. Venot, J.Y. Devaux, E. Walter, J.F. Lebruchec, L. Dubertret, J.C. Roucayrol, Automated registration of dissimilar images: application to medical imagery, Computer Vision, Graphics and Image Processing 47 (1989) 77-88. [ 1 l] B. Kamgar-Pa&, J.L. Jones, A. Rosenfeld, Registration of multiple overlapping range images without distinctive features, IEEE Computer Sot. Conf. on Computer Vision and Pattern Recognition, 4-8 June 1989, San Diego, CA, pp. 2X2-290. [12] D.S. Meek, D.J. Walton, Surface alignment or registration of multiple range images, Congressus Numerantium 99 ( 1994) 17 I- 177. [13] K. Deguchi, Registration technique for partially covered image sequence, IEEE Proc. 8th Int. Conf. on Pattern Recognition, Paris, October 27-31, 1986, pp. 1186-1189. [ 141 D.S. Meek, D.J. Walton, Comparison of smooth polynomial functional surfaces for use in alignment, Technical Report 96/I 1, Department of Computer Science, University of Manitoba, Winnipeg, Canada, 1996. [15] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes in C, 2nd Edn., Cambridge University Press, Cambridge, UK, 1992.