Multiscale Curvature-Based Shape Representation Using ... - CiteSeerX

Multiscale Curvature-Based Shape Representation Using B-spline Wavelets

Contents 1 Introduction

2

2 Fast implementation of B-spline scale-space

4

2.1 De nition and properties of B-splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Parameterization and diusion of a curve using B-spline bases . . . . . . . . . . . . . . 2.3 Flexible and ecient implementation of B-spline scale-space using lter bank techniques 2.3.1 Implementation of scale-space at arbitrary scale . . . . . . . . . . . . . . . . . . 2.3.2 Implementation of scale-space at dyadic scales . . . . . . . . . . . . . . . . . . 2.3.3 B-spline wavelet transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 A theory of curve representation and evolution in B-spline scale space

4 5 6 6 6 7

7

3.1 Some properties of B-spline representation of a curve . . . . . . . . . . . . . . . . . . . 7 3.2 The properties of evolving curve in B-spline scale-space . . . . . . . . . . . . . . . . . 9 3.3 Evolution of random noise in B-spline scale-space . . . . . . . . . . . . . . . . . . . . . 12

4 The curvature scale space image and its properties

14

4.1 The computation of curvature function using B-spline wavelets . . . . . . . . . . . . . 14 4.2 Curvature scale-space image of the curve and its properties . . . . . . . . . . . . . . . 15

5 Application to data compression and implementation examples

16

5.1 An algorithm for curve data compression . . . . . . . . . . . . . . . . . . . . . . . . . . 16 5.2 Experimental examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

6 Conclusions and discussions

19

0

Multiscale Curvature-Based Shape Representation Using B-spline Wavelets



Wang Yuping, S. L. Leey, and K. Toraichi z August, 1997 Abstract This paper presents a new multiscale curvature-based shape representation technique with application to curve data compression using B-spline wavelets. The properties and behavior of evolving curves in B-spline scale-space are investigated. They are shown to be similar to that in Gaussian scale-space, but B-spline representation enjoys a number of advantages over Gaussian kernels, for instance, the availability of fast algorithms. As an application of the multiscale shape representation, we introduce an algorithm and present some illustrative examples on curve data compression. The B-spline wavelet transforms are used to estimate eciently the multiscale curvature functions. Based on the curvature scale-space image, we introduce a coarse-to- ne matching algorithm to automatically detect the dominant points and use them as control knots for curve interpolation.

Keywords: B-spline, diusion equation, curve evolution, curvature function, curve tting. EDICS category: 1.6, 1.8, 3.4

This work is supported by Wavelets Strategic Research Programme funded by the National Science and Technology Board and Ministry of Education under the Grant RP960 601/A. y The authors are with the Wavelets Strategic Research Programme, Department of Mathematics, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260. Phone: (65)8742773, Fax: (65)7795452, Email: fwyp, [email protected] z Institute of Information Science and Electronics, University of Tsukuba, Tsukuba, Ibaraki 305, Japan. 

1

1 Introduction Two-dimensional shape representation and analysis is an important issue in many applications of machine vision, such as character recognition, automatic visual inspection, ngerprint identi cation and contour matching for medical imaging ( see [13], [14], [18], [19]). More speci cally, the contour of a shape contains a substantial amount of information about the original shape [10], [25]. One of the most important parameters characterizing contours is the curvature at each of their constituent points, particularly at places characterized by in ection or curvature modulus maxima, which are usually called visual primitives [10]. The curvature function of the curve is invariant under ane geometric transformation, i.e., scaling, rotation and translation. There are psychological results showing that curvature plays a fundamental role in human shape perception and representation [24]. Since the visual primitives, usually characterized by high curvature of the curve, occur at dierent density, it is necessary that these features be described in a multiscale sense. Thus the curvature scale-space image was introduced by Mokhtarian and Mackworth [12] as a tool for representing planar curves. This representation is computed by convolving a path-based parametric representation of the curve with a Gaussian function with varying standard deviation. In other words, the evolution of the curve

C(u; s) : S  [0;  ] ?! R ; 1

2

(1)

where u 2 S 1 is the arc-length parameter of the curve and s is the scale level, is according to the following heat diusion equation,

8 @C @2C < @s = @u2 : C(u; 0) = C (u)

(2)

0

where C0 (u) = C(u; 0) is the initial curve. The solution vector of this classical heat equation, C(u; s) = [x(u; s); y(u; s)]T , can be obtained by convolution with the Gaussian. Most classical multiscale representations of planar curves ([10], [11], [12]) are based on Gaussian function because the evolution of curve in such a scale-space satis es the causality criteria, i.e., no new feature points are created when the scale increases ([8]). With the advent of wavelet transforms, several approaches on the representation and analysis of planar curves using this tool have been introduced. G. H. Chang, C. C. Jay, Kuo [17] have used orthogonal and biorthogonal wavelet expansion for multiscale representation of curves and studied the 2

properties in wavelet representation. L. M. Reissell [16] constructed a family of compactly supported symmetric biorthogonal wavelets with interpolating scaling functions and used them for compression and curve-curve intersection. The advantage of these wavelet methods is that the approximation properties of wavelet multiresolution lters yield hierarchical curves that lie close to the original data, a property which is not shared by other lters [16]. However, because these methods do not rely on the curvature primitives, they may not satisfy the invariance criteria as required by a multiscale shape descriptor [12]. Lee J. S. et al. [9] employed the multiscale edge detection algorithm for corner detection. In fact, the smoothing kernels used in their methods are based on B-spline kernels which can lead to fast cascade lter bank implementation of the multiscale transform as indicated in our previous work [7]. This paper presents a curvature-based multiscale shape representation using B-spline wavelets. In a broad sense, we attempt to develop a B-spline based multiresolution shape representation theory and algorithm as an alternative to the Gaussian-based multiscale description. We note that the most of multiscale shape representation or analysis methods are based on Gaussian kernel. However as pointed out by Mokhtarian and Mackworth [12], the representation should satisfy the eciency and ease of implementation in order that it be useful for practical shape recognition tasks in computer vision. Obviously, Gaussian kernel does not satisfy such criteria since the computational burden is high at large scales. As our previous results have shown [6], [7], B-spline is a good approximation to Gaussian kernel. Besides inheriting almost all the good properties of Gaussian scale-space, B-spline derived scale-space exhibits many advantages. For example, it provides fast algorithms and parallel implementation at multiple scales. In a narrow sense, this paper introduces a multiscale curvaturebased algorithm using B-spline wavelets and apply it to curve data compression, which are dierent from that proposed in [17], [16]. B-spline methods have been widely used in curve design [1] and curve tting [19], [18]. We use its unique re nable property [4] as an alternative to Gaussian for multiscale shape analysis. Recently, G. Sapiro et al. [15] proposed an ane invariant multiscale representation of planar shapes using B-spline, but they select the scale as the order of the B-spline, an approach which is dierent from ours. We use dilated B-spline of xed order for scale-space ltering. The remainder of the paper is organized as follows. Section 2 reviews and presents fast lter bank implementation of scale-space ltering using B-spline technique. In Section 3, the geometric properties of evolving curve in B-spline scale space are presented. Section 4 studies the multiscale curvature image of a curve. Section 5 introduces an ecient scale-space algorithm for curve data compression. Some experiments are performed to illustrate the performance of the B-spline wavelets. 3

Finally, we give some concluding remarks on curve analysis based on B-spline scale-space.

2 Fast implementation of B-spline scale-space 2.1 De nition and properties of B-splines We now brie y describe the B-spline theory. For more details, see for example [1], [2], [6], [7]. The continuous central B-spline of order n is denoted by n (x), which can be generated by repeated n + 1 convolutions of B-spline of order 0,

z 0 0n}|+1 {0 n 0 n ? 1 (x) =  =  




(3)

where 0 (x) is the 0th-order B-spline, i.e., an impulse with support [ ?21 ; 21 ]. The Fourier transform of n (x) is

^n(!) = [ ^0 (!)]n+1 = (

sin !2

! ) 2

n+1

= sincn+1 !2 :

(4)

The discrete B-spline of order n at scale level m is de ned as

{ z 0 0n}|+1 n Bm = Bm  Bm 


 Bm0 ;

(5)

where Bm0 = m1 [1; 1;


; 1] is a normalized sampled pulse of width m. The continuous nth-order B-spline is related to the discrete nth-order B-spline through the following equation [4], [7] (suppose n is odd throughout the paper), 1 n( x ) =

m

m

X1 +

k=?1

Bmn (k) n (x ? k):

(6)

This property means that B-spline is m-re nable. When m = 2, it is just the famous two-scale relation for B-splines. The discrete sampled B-spline of order n at resolution m is denoted as bnm (k), which is obtained by directly sampling the nth-order continuous B-spline at the scale m

bnm (k) = m1 n( mk ) 8k 2 Z:

(7)

When m = 1, the discrete sampled B-spline of order n is just denoted as bn (k) = n (k). The inverse operator (bn )?1 is de ned as (bn )?1  bn (k) = (k): 4

(8)

From the m-re nable property (6) of the continuous B-spline, we can also get the m-re nable relation about the discrete sampled B-spline as follows,

bnm(k) = Bmn  bn(k):

(9)

2.2 Parameterization and diusion of a curve using B-spline bases Our strategy for multiresolution representation of a curve is to use the B-spline bases to rst parameterize the original curve and then diuse it using the B-spline kernel instead of the Gaussian kernel. Suppose that we have a discrete set of points Vi = [xi ; yi ]T (1  i  N ) on a curve. First we can obtain a parameterized curve using B-spline bases,

C (u) = 0

X k

c (k) l (u ? k) 0

(10)

where C0 (u) = [x(u; 0); y(u; 0)]T and the projection coecient vector c0 (k) are given by

c (k) = (bl )?  V(k): 1

0

(11)

Usually when we take the order of B-spline l=0, the projection coecients c0 (k) are just the original sampling points of the curve and often called control vertices. The polygon formed from the control vertices is called the control polygon. Unlike the Gaussian scale-space evolution, the evolution of curve at dierent instances or resolution levels s in B-spline scale-space is achieved by convolving the curve with the dilated B-spline kernel instead of the Gaussian kernel

C(u; s) = C  sn(u) 0

(12)

where sn (u) = 1s n ( us ) denotes the n-th order B-spline at the resolution s. Here for clarify, we use n to denote the dierent orders of the B-spline used for smoothing. Due to the central limit theorem, B-splines approximate the Gaussian kernel as the order tends to in nity. However, in practice, cubic B-spline already provides a good approximation of the Gaussian function [3]. The advantage of such an implementation is its eciency as detailed in the following section.

5

2.3 Flexible and ecient implementation of B-spline scale-space using lter bank techniques 2.3.1 Implementation of scale-space at arbitrary scale Because of the m-re nable property (6) of B-splines, we are able to obtain the evolving curve at a certain instant s in an ecient way. Since any real scale s can be represented or approximated by a rational number, we shall assume that s = mm21 . In [7], we derive a fast and parallel algorithm for the cascade implementation of the smoothing operation (12). Here we just present the formula. One can refer to [7] for more details. The evolving curve is

m ) = m (bl C(u; m 1

2

2

n

+ +1

 Bml 2  Bmn 1  c "m2  p"m1 )#m2 (u) 0

(13)

where c0"m2 represents the up-sampling operation applied to c0 (k) given in (11), i.e. taking the initial discrete curve approximation points and expanding them by adding m2 ? 1 zeros between consecutive points. The operation # m2 is the down-sampling operation by a factor of m2 , i.e. taking one sample among every m2 samples. If we take an integer scale , we obtain a similar formula as in [2]. Formula (13) can be realized by a moving average technique using the de nition of Bml 2 in (5). The complexity of the algorithm is O(N ), independent of the scale level, where N is the number of original curve sampling points.

2.3.2 Implementation of scale-space at dyadic scales If we restrict the scale level to the dyadics, mm12 = 2m , it is interesting to note that we can obtain the recursive smoothing between the levels 2m and 2m?1 . At level 2m , the formula (13) simpli es to

C(u; 2m ) = (bl

n

+ +1

 B nm  C )(u) = B nm  C (u) 0

2

0

2

(14)

where C 0 = bl+n+1  C0 is regarded as the original approximation of the curve. It is easy to establish the relationship of the evolving curve between two successive dyadic levels as follows ( [7], [6]):

C(u; 2m ) = 2n1

+1

0 1 n X @ n+1 A +1

i=0

j

C(u ? j 2m? ; 2m? ): 1

1

(15)

Therefore, curve evolving at dyadic levels can be implemented recursively. Since the convolution with the binomial kernel can be realized with the addition and bit shift operation [6], we can obtain an even more ecient scheme at dyadic scales. One can usually consider the evolving behavior of the curve only at dyadic scales and obtain a more ecient characterization as shown in [9]. 6

2.3.3 B-spline wavelet transform Similar to the Marr-Hilderth operator and Canny operator, if we de ne the wavelets as rst or second derivatives of nth-order B-spline, we will obtain the B-spline wavelet transform [6]. In this case, the B-spline kernel in the formula (12) is replaced by its rst or second derivative. It is easy to show that the wavelet transforms are obtained as

W1 C(u; s) =
u C(u; s) and W2C(u; s) =
2uC(u; s);

(16)

where
u f (u; s) = f (u + 21 ; s) ? f (u ? 12 ; s). Thus, if we take the rst or second dierence of the B-spline scale-space ltering (13) we will obtain the B-spline wavelet transform. These wavelet transforms are useful in the computation of curvature function in Section 4.

3 A theory of curve representation and evolution in B-spline scale space In this section we want to study the properties of multiscale B-spline representation and the behavior of the curve evolution in such a representation. Most existing multiscale shape analysis algorithms are derived from Gaussian function or in terms of heat diusion equation (2). We shall show that curve evolution using B-splines not only enjoys the similar properties, but is also computationally more ecient.

3.1 Some properties of B-spline representation of a curve Let us rst show the invariance of B-spline representation. Suppose the original curve V is geometrically transformed into V~ by rotation, translation and uniform scaling (A; T ),

V~ = AV + T

(17)

where V 2 R2 is a vector, A 2 GL+2(R) is the ane group, and T 2 R2 is the translation vector. Then

C~ (u) = AC (u) + T; 0

0

where C~ 0 (u) is the B-spline curve de ned by the points V~ . 7

(18)

Property 1 The B-spline scale-space representation is invariant under the ane geometric transform.

The invariance property is essential since it makes it possible to match a planar curve to another similar shape that has undergone a transformation consisting of arbitrary rotation, uniform scaling, and translation. It is also important to study the stability of B-spline scale-space representation. Stability means that if two curves have a small dierence their representations should also have a small dierence and if two representations have a small dierence, the curves they represent should also have a small shape dierence [12]. Since the shifted B-spline are unique spline that constitute a Riesz basis [4], it is not dicult to show that the B-spline representation is stable.

Property 2 The B-spline representation (10) provides stable representation of planar curves. The following property gives the constraint on the physical location of a planar curve as it evolves according to (12).

Property 3 Let X be a closed B-spline curve and let G be its convex hull. Then X remains inside G during evolution. This follows from the convex hull property of the B-spline representation (10). This property is useful for stereo matching application, where it is advantageous to carry out matching at coarser levels rst and then match at ner levels to increase accuracy. One of the requirement for scale-space representation is the causality which states that the number of in ection points or zero crossings of the curvature decreases while the scale increases. The following results has been proved by Goodman [5].

Property 4 If the control polygon is regular, then the number of in ection points in the B-spline curve does not exceed that of its control polygon.

There are other useful properties of B-spline representation. One is the local property which means changing one control point will aect only a small region of the curve about the point. Another is the variation diminishing property which states that the sign changes of the curve is not more than that of its control polygon. All these properties show that B-spline provides a good tool for geometric curve representation and description. 8

3.2 The properties of evolving curve in B-spline scale-space In order to study the evolving curve C(u; s) in B-spline scale-space, we use Fourier representation to characterize its evolving behavior. Suppose that the closed curve C(u; s) = [x(u; s); y(u; s)]T for s 2 [0;  ) is a period function of u with period T , the curve length of the initial curve C(u; 0). Then the curve can be expanded into the following Fourier series: 1 2 3 2 P a(l; s)ej 6 x ( u; s ) l ?1 4 5 = 64 P 1 =

y(u; s)

l=?1

=T )lu

(2

b(l; s)ej(2=T )lu

3 77 ; 5

(19)

where a(l; s) and b(l; s), the l-order frequency components of x(u; s) and y(u; s) at scale s, can be computed by

Z a(l; s) = T1 x(u; s)e?j(2=T )lu du; T

0

Z 1 b(l; s) = T y(u; s)e?j(2=T )lu du: T

(20)

0

If l = 0, the direct frequency components

Z 1 a(0; s) = T x(u; s)du T

0

Z b(0; s) = T1 y(u; s)du T

(21)

0

de ne the mass of the curve at scale s. The frequency components a(l; s), b(l; s) are generally complex and can be written as

a(l; s) = a0(l; s)e?j(l;s) b(l; s) = b0 (l; s)e?j (l;s) where a0 (l; s), b0 (l; s) are the amplitudes of a(l; s) and b(l; s) respectively, related to the size of the evolving curve, and (l; s), (l; s) are their corresponding complex angles which determine the orientation of the curve. Since x(u; s) and y(u; s) are real values, the amplitude a0 (l; s), b0 (l; s) are symmetric with respect to l and the angles (l; s), (l; s) are antisymmetric with respect to l. Therefore, (19) can be rewritten 9

as

3 1 2 3 2 a(0; s) + P 2 a ( l; s ) cos ((2 =T ) lu +  ( l; s )) 77 : l1 4 x(u; s) 5 = 664 P 5 y(u; s)

b(0; s) +

=1

l=1

0

2b0 (l; s)cos((2=T )lu + (l; s))

(22)

If the frequency components higher than zero order are zero, i.e., a0 (l; s) = 0 and b0 (l; s) = 0 for l  1, the curve is a point. If the frequency components higher than the rst order are zero, i.e., a0 (l; s) = 0 and b0 (l; s) = 0 for l  2, then the curve equation simpli es to

2 3 2 3 x ( u; s ) a (0 ; s ) + 2 a (1 ; s ) cos ((2 =T ) u +  (1 ; s )) 4 5=4 5; 0

y(u; s)

(23)

b(0; s) + 2b0 (1; s)cos((2=T )u + (1; s)) which is an ellipse equation. The center of this ellipse is (a(0; s); b(0; s)) and its two axes are 2a0 (1; s), 2b0 (1; s). The Fourier spectrum of (19) is

3 1 2 3 2 P a(l; s)(! ? j (2=T )l) 7 ?1 75 ; 4 x^(!; s) 5 = 664 l P 1 =

y^(!; s)

l=?1

b(l; s)(! ? j (2=T )l)

where x^(!; s) denotes the Fourier transform of x(u; s) with respect to the variable u. When s = 0

3 1 2 3 2 P a ( l; 0)  ( ! ? j (2 =T ) l ) 77 : ?1 4 x^(!; 0) 5 = 664 l P 1 5 =

y^(!; 0)

l=?1

b(l; 0)(! ? j (2=T )l)

(24)

where a(l; 0), b(l; 0) are the two frequency components of the initial curve C(u; 0). According to the evolution equation (12) and the convolution theorem of Fourier transform, the frequency spectrum of the curve at scale level s can be obtained from both the spectrum of the initial curve and that of the B-spline function. 3 1 2 3 2 P ^n (s!)(! ? j (2=T )l) a ( l; 0) 77 ?1 4 x^(!; s) 5 = 664 l=P 1 5 b(l; 0) ^n (s!)(! ? j (2=T )l) y^(!; s) ?1 3 2 lP 1  n ^ 66 l ?1 a(l; 0) (s T l)(! ? j (2=T )l) 77 : 1 5 4 P =

2

=

=

l=?1

b(l; 0) ^n (s 2T l)(! ? j (2=T )l)

(25)

Equation (25) shows the evolution of the frequency components, a(l; s) = a(l; 0) ^n (s 2 l) = a(l; 0)(sinc(  ls))n+1

T

T

T

T

b(l; s) = b(l; 0) ^n (s 2 l) = b(l; 0)(sinc(  ls))n+1 : 10

(26)

Equivalently, the amplitude and phases of the dierent frequency components of the evolving curve are related to those of the initial curve as follows:

a0 (l; s) = a0 (l; 0)sinc( T ls))n+1 ; b0 (l; s) = b0 (l; 0)(sinc( T ls))n+1

(27)

(l; s) = (l; 0); and (l; s) = (l; 0):

(28)

and

Equation (28) shows that phases of all the frequency components of a curve do not change when the curve evolves according to (12). Since the phase of the frequency component determines the orientation of the curve, we can draw the following conclusion.

Theorem 1 The orientation of a curve does not change when this curve evolves in the B-spline scale-space.

When l = 0, the direct frequency components become

a(0; s) = a(0; 0); and b(0; s) = b(0; 0):

(29)

This means that the direct frequency components remain unchanged across scales. Therefore, according to the de nition in (21) we have

Theorem 2 When a curve evolves in the B-spline scale-space, its center of mass does not move. This property shows that the center of a planar curve remains stationary as the curve evolves in B-spline scale-space. Since liml!1 (sinc( T ls))n+1 = 0, it follows from (26) that, when l gets larger, the high-frequency components are damped out due to smoothing and the evolving curve becomes increasingly smoother. Correspondingly the noise is removed. Moreover, the convergence is such that the shape of the curve approaches an ellipse. This is because the high frequency components of the Fourier spectrum die out much faster than the low frequency ones. Therefore, the limiting curve becomes approximately an ellipse, when only the zero and rst frequency components remain signi cant.

Theorem 3 When a curve evolves in the B-spline scale-space from s = 0 to s = 1, its shape becomes elliptical and nally shrinks to a point.

11

In [15], Sapiro et al. take the order of the B-spline as the scale parameter and when the order n ! 1 one also has limn!1(sinc( T ls))n+1 = 0. Therefore, the same conclusion about the shrinking property of the evolving curve also holds. If the curve evolves according to the Euclidean geometric heat equation, Gage and Himilton [20] proved that any smooth, embedded convex curve converges to a round point. Grayson [21] proved that any non-convex embedded plannar curve converges to a simple convex curve and hence any embedded curve converges to a round point. A similar conclusion is drawn in [13]. Theorem 3 shows that Gaussian kernel is not the only lter which shorten a curve. In the above analysis, if we use the Gaussian kernel instead of B-spline kernel, we will nd their dierence lies in the shrinking speed. For Gaussian kernel, the decreasing speed is exponential. While for B-spline of order n, the decreasing speed is polynomial O( sn1+1 ), a little slower. As an illustration, we show in Fig. 1 the evolution of a von Koch snow ake in B-spline scale-space. The von Koch curve is self similar at dierent resolution with fractal dimension 1.2618. It can be shown that, more and more details of the curve are smoothed out as the scale increases. At larger scale, the snow ake shrinks to a circle.

3.3 Evolution of random noise in B-spline scale-space In this section we study the in uence of random noise on B-spline scale-space representation of a curve. Two types of noise model are considered. One is Gaussian white noise f with the power spectrum

P (!) = 1:

(30)

Another is the fractional Brownian noise [23], which is self-similar in a stochastic sense with the power spectrum proportional to j!j21H +1 , for some exponent 1 > H > 0, i.e.,

P (!) =  j!j21H +1 ; 1 > H > 0:

(31)

where  is an absolute constant. The density of local extrema D for a stationary Gaussian process can be obtained according to the second and fourth order derivatives of the autocorrelation function R or equivalently according to the second and fourth order moments of its spectrum P [22]:

s

(4) (0) D = 21 ? R R(2) (0) ;

12

(32)

and the derivatives of the autocorrelation function at zero can be expressed in terms of the power spectrum by

?R (0) =

Z1

(2)

?1 Z1

R(4) (0) =

?1

!2 P (!)d!; !4 P (!)d!:

(33)

If the noisy curve evolves as in (12), it is easy to verify that the power spectrum Ps of the curve at scale s is related to the power spectrum P0 of the curve at scale 0 by

Ps (!) = j ^n (s!)j2 P0 (!):

(34)

Since the B-spline of order n is asymptotically equal to Gaussian function by the following relation [7]:

n(x) 

s

6 exp(? 6x2 ); (n + 1) n+1

(35)

their Fourier transforms are related by + 1 !2 ): ^n (!)  21 exp(? n 24

(36)

From this relation and the following formula

Z1

?1

tme?t2 dt =

?( m2+1 )

(37)

 m2+1

we can easily estimate the moments of the power spectrum (33). For Gaussian white noise,

?R (0) 

Z1

(2)

?1 Z1

R(4) (0) 

?1

+ 1 !2 )d! = !2 41 exp(? n 12

?( 32 ) ; ( n12+1 s2 ) 32

+ 1 !2 )d! = !4 14 exp(? n 12

?( 52 ) : ( n12+1 s2 ) 52

(38)

Therefore, the density of local extrema of a Gaussian white noise at scale s is p 3 Ds = 22 1s p 1 : (39) n+1 Similarly for fractional Brownian noise, the second and fourth moments of its power spectrum are

Z1

+ 1 !2 )!?(2H +1) d! =  ?(1 ? H ) ; !2 14 exp(? n 12 ( n12+1 s2 )1?H ?1

?R (0)   (2)

Z1

+ 1 !2 )!?(2H +1) d! =  ?(2 ? H ) : !4 14 exp(? n 12 ( n12+1 s2 )2?H ?1

R (0)   (4)

13

(40)

and the density of local extrema of the fractional Brownian noise at scale s is

p

p 1 ? H; 0 < H < 1: Ds = 23 1s p 1 n+1

(41)

From (39), (41) we can arrive at the following conclusion about the in uence of scale and order of B-spline on the evolution of random noise.

Theorem 4 If a curve generated by either the Gaussian white noise or fractional Brownian motion,

the density of local extrema in B-spline scale-space decreases with scale as s?1 or decreases with the 1 order of B-spline as n? 2 .

Theorem 4 indicates that the in uence of the random noise in a curve can be weakened by increasing the scale or the order of B-spline.

4 The curvature scale space image and its properties 4.1 The computation of curvature function using B-spline wavelets The most important and signi cant description of a planar curve is its curvature function. The curvature function of the curve C0 (u) = [x(u); y(u)]T is de ned as

(u) = ((x_x(_u(u)y))(2u)+?(y_y_((uu)))x2()u3=)2

(42)

where x_ (u), x(u) represent the rst and second derivative of x(u), and similarly for y_ (u), y(u). If u is the normalized arc length parameter, then (42) can be rewritten as

(u) = x_ (u)y(u) ? y_ (u)x(u):

(43)

From Property 1 it is easy to show that the curvature function satis es the invariance property which is essential in computer vision. For example, it makes it possible to match a planar curve to another of similar shape that undergone an ane transformation consisting of arbitrary amounts of rotation, scaling and translation. For the evolving curve according to (12), the curve coordinates x(u; s), y(u; s) are given explicitly by

x(u; s) = x(u)  sn (u) y(u; s) = y(u)  sn(u): 14

(44)

Therefore,

@ (x(u)  n (u)) = x(u)  @ ( n (u)) x_ (u; s) = @u s @u s = x(u)  ( sn?1 (u + 1=2) ? sn?1 (u ? 1=2))

(45)

@ 2 (x(u)  n(u)) = x(u)  @ 2 ( n(u)) x(u; s) = @u s 2 @u2 s = x(u)  ( sn?2 (u + 1) ? 2 sn?2 (u) + sn?2 (u ? 1))

(46)

and

Hence, the computation of x_ (u) and x(u) are just the B-spline wavelet transform de ned in (16) and similarly for y_ (u) and y(u). This computation can be implemented using ecient lter bank algorithm or subdivision scheme as discussed in Section 2. According to the de nition (42), we can now express the curvature function at the scale s as follows:

s)W2 y(u; s) ? W1 y(u; s)W2 x(u; s) : (u; s) = W1 x(u;(W x(u; s)2 + W y(u; s)2 )3=2 1

1

(47)

4.2 Curvature scale-space image of the curve and its properties For dierent scale level, the function de ned implicitly by

(u; s) = 0

(48)

are called curvature scale-space image of a curve by Mokhtarian and Mackworth [12], or the curvature primal sketch by Asada and Brady [10]. Many tasks in computer vision such as curve compression, automatic object recognition can be performed by identifying signi cant points in the shape contour, which are generally points with high absolute curvature values. High curvature points also play a fundamental role in human vision system, both from the psychological and neurological points of view [24], [25]. Hierarchical or multiscale representation provides a robust way for the curvature estimation. The mathematical analysis of curvature changes in scale-space gives us an a priori knowledge for curvature identi cation. Asada and Brady [10] modeled two curvature changes as corners and smooth joints and gave an explicit mathematical characterization. They also studied the set of compound primitives such as End, Crank, Bump or Dent. In [11] Rattarangsi and Chin also presented a careful study of dierent types of shape models. All their approaches are based on Gaussian kernel. 15

In this paper we use B-spline kernels. Using a little more general mathematical analysis, it can be shown that the curvature edge changes in B-spline scale-space have the same behavior as that of Gaussian scale-space. Furthermore, in a discrete sense the extrema or zero-crossing number of the curve coordinates does not increase as the scale varies from low to high level. In other words, the scaling or causality property for Gaussian kernel also holds in the case of B-splines. Here we omit the detail which can be found in [7]. In fact, for discrete implementation Gaussian function is usually truncated, thus there is little dierence between Gaussian function and B-splines in practice.

5 Application to data compression and implementation examples The curvature scale-space image of a curve is equivalently characterized by its curvature extrema at dierent scales, which provide signi cant information about the object contour. The curvature extrema are commonly called dominant points. As Attneve pointed out [24], dominant points concentrate information about the shape of curves, which means they are very important in describing the shape. They can be used as a compact and eective representation for shape analysis and pattern recognition, for example, see [25]. In this section, we introduce an algorithm for curve data compression application making use of curvature image of B-spline scale-space.

5.1 An algorithm for curve data compression Generally our strategy for data compression is to estimate the curvature function which corresponds to the dominant points or corners of the curve, and then use these dominant points to reconstruct the original curve by interpolation or tting technique. We present a detailed description of our algorithm in four steps. 1. Obtaining the curvature image of B-spline scale-space. At rst we use the fast algorithm introduced in Section 2 to obtain the evolving curve C(ui ; sk ), 1  i  N; smin  k  smax at dierent scales, where N is the number of original curve sampling points and smin , smax are the minimal and maximal scales. Then by taking the rst and second dierence of the curve C(ui ; sk ) we obtain the wavelet transforms of the curve coordinates. Using the formula (47), we obtain the curvature estimation at dierent scales. At this step, one has the freedom to choose the discrete scale parameters using the algorithm given in Section 2. 16

2. Coarse to ne matching of the dominant points. The locations of peaks in the curvature scale-space image correspond to the corner points of the curve. It turns out that by looking at the movement of the peaks over several scales we can localize dierent types of corner points [10]. At ner scales we can obtain a good localization of corner points, but due to the in uence of noise or high frequency details, more corner points are detected as indicated in Theorem 4. At coarse scale, we will obtain the overall picture of the curve, but the locations of corner points are not accurate due to the smoothing procedure. Therefore, we have to combine the multiscale information to localize the dominant points of the curve. After extracting the local peaks among scales, we then match these ngerprints using the following coarse-to- ne strategy. (a) The local maxima points are classi ed into local positive maxima and local negative minima points. The polarization does not change with the scale. From the coarse to ne tracking, these two classes of points are matched across scales respectively. (b) A local maxima point P i+1 at level i +1 is said to correspond to the another local maxima point P i at the ner scale i, if

d(Pki+1 ; Pli)  "

(49)

where " is the predetermined threshold and the distance d(pik+1 ; pil ) can be measured in the sense of either l2 or l1 . The matching between the two points P i+1 and P i is searched in the neighborhood of P i+1 . If the two points match each other, then the location of P i+1 is corrected to be the same as that of P i . This matching procedure is implemented from large scale s2 to small scale s1 and repeated several times until the locations of the corner points are found. 3. Post-processing of the corner candidates. In the second phase, the candidate corners are detected and those with signi cant changes in direction are identi ed as corners. At ner scale, there may be more corner candidates whose locations are very close and are not necessary for restoration. In order to alleviate such a problem, we introduce the following geometric criteria to reduce the corner candidates. Our idea is, if the three consecutive corner candidates Pki?1 , Pki , Pii+1 are on the same straight line we regard that Pki is a phantom corner and remove it out. This is illustrated in Fig. 2. In implementation we use the following formula to judge Pki : 17

Pki?1 Pki + Pki Pki+1  Pki?1 Pki+1

(50)

where Pki?1 Pki denotes the length of the line segment between the two points Pki?1 and Pki . The parameter  is usually taken as 0:9  1. One may choose other criteria to re ne the detection of corner points. For example, if the two corner candidates are within a predetermined threshold, then one of them will be deleted. In this way, the data points needed to be store can be greatly reduced. 4. Curve reconstruction using interpolation or tting method. Once the above procedure is completed, we select the dominant points at a certain scale and store their positions. Thus a compressed form of curve is obtained. In order to recover the original curve, we use these corner points as control knots and selecting some interpolation methods such spline and polynomial interpolation to get the coordinates of the curve at other locations. One may obtain a better compression ratio at large scale since there are few corner points, but the recovered curve may also be poor. On the contrary, good reconstruction quality can be obtained with the cost of lower compression ratio at small scale. So a proper scale should be chosen depending on the needs at hand.

5.2 Experimental examples In this section we present one example to implement the above algorithm. As illustrated in Fig. 3, a typical terrain map containing dierent features at dierent locations are shown. In some areas, for example, around the x-coordinate 138.7 in the gure, there are abrupt changes of small size, which need ner scale representation, whereas other parts of the curve segment are very smooth that can be represented by fewer corner points. For example, the segment from the x-coordinate 138.7 to 138.85 is very smooth. Thus only few corner points are required for its description. The need for multiscale wavelet representation for adaptive processing of the curve is essential in this example. We rst obtain the curvature image at dierent scales. According to the algorithm presented in Section. 2, we can choose the scale level freely. In practise we usually choose the scale level as 1,3,5,7,9, ... The estimated curvature function at several scales are given in Fig. 4. We then do the coarse to ne matching across these scales. The detected corner points at level 1 are drawn in Fig. 5 as crosses. For illustration we draw the detected corner points in a small portion of the curve in 18

Fig. 7. It can be seen around the x-coordinate 138.7, more corner points have to be assigned in order to record the ner change of the curve orientation. Between the x-coordinate 138.7 and 138.8, only one point is necessary for the representation of this smooth segment. In Fig. 6 the reconstructed curve are drawn in dotted line, which gives almost the same curve as the original one. In the paper we select the spline tting for reconstruction. The compression ratio for this curve is de ned as the ratio between the number of the original curve points over that of the compressed data points. The compression ratio in this case is 9.2670, i.e., only 221 points are detected to represent the original curve with 2048 points. The reconstruction error is de ned as r Pn 2 i=1 (di ) error =

n

where di is the distance from reconstruction data points to original data points and n is the number of data points. For the scales 3, 5, 7 the detected corner points and recovered curve are drawn in Figure 8, 9 and 10. We can see that there is a tradeo in selecting the width or scale of the B-spline lter, a larger width will remove small details of the boundary curvature, a small width will permit false concavities and convexities. Therefore, at the ner scale the reconstruction error is lower with comparatively low compression ratio. At the coarser scale the reconstruction error is large with high compression ratio. The compression ratio is also related to the type of curve data at hand. An intermediate scale level is recommended in practice.

6 Conclusions and discussions In this paper we present a general theory of multiscale shape representation based on B-spline wavelets with applications to curve data compression. We now brie y summarize the main contributions of the proposed method. 1. For a long time, scale-space theory is mainly based on Gaussian function. So multiscale curve representation and analysis is mainly based in such a representation [10], [11]. According to the criteria of Mokhtarian and Mackworth [12] for multiscale curve representation, one requirement is eciency. Obviously, most of these approaches lack eciency. As an alternative, we found that B-spline kernels have many advantages in deriving scale-space representation [7]. Using B-spline wavelets [6], we can obtain ecient and parallel implementation of scale-space ltering, which are very suitable for multiscale curve representation. There has been a continuous eort for ecient multiscale shape representation. In [13] Li proposed the exponential kernel to approximate Gaussian function. 19

2. We have studied the properties of multiscale B-spline curve representation. Our results discussed in Section 2 shows that B-spline representation satis es very nice properties such as the invariance, stability and so on, as required by a shape descriptor. Furthermore, it inherits most of the curve evolving behavior from the Gaussian scale-space. Therefore, this result provides some a priori knowledge of scale-space matching and a sound theoretical foundation for shape descriptor based on B-spline wavelets. Moreover, in contrast with the classical scale-spaces, B-spline representation is de ned directly on an initial discrete set of samples, and the outputs of B-spline wavelet transform are directly used for curvature computation as shown in (47), thus avoiding possible problems caused by discretization of continuous Gaussian scale-space. In other words, the proposed approach is easy to implement. Furthermore, since B-spline representation is a piecewise polynomial representation, it is also easy to determine the properties of the shape of the curve. All of these indicate that the multiscale B-spline representation satis es the criteria as described in [12]. 3. We have used continuous B-spline wavelet transform or dyadic wavelet frame transform for the multiscale curvature estimation in the paper. This approach is dierent from [17] and [16]. We make use of wavelets for the characterization of the geometric invariance and the wavelet transforms are translation invariant. In [16], the curve coordinates are expanded in the biorthogonal wavelet bases and the geometric feature of the curve was not easily studied. The wavelet transform in our paper is in fact a generalization of that used in [9] and we present a systematic theory of this wavelet representation for shape analysis. 4. B-spline technique has been widely used in curve design and analysis. Examples include [19], [18], [1]. In this paper, we propose it for multiscale analysis. The fast wavelet transform is in fact one type of subdivision scheme. In [15], Sapiro et al. have applied B-splines for ane scale-space representation and they choose the order of the B-splines as the scale parameter which is dierent from ours. 5. As a nal remark, only the problem of representing the shape of a planar curve has been addressed in this paper. In the paper, we only consider the example of application to data compression. The method developed can also be used for other pattern recognition tasks such as shape information extraction for object recognition, noise elimination, human shape recognition, curvature estimation, and so on. We hope that the algorithm in this work can be further improved by developing more ecient coarse to ne matching algorithm, data tting algorithm and other methods. Acknowledgment: S. L. Lee acknowledges with thanks the data on the map from Ms. Hiyama of Tokyo University. Professor T. N. T. Goodman at University of Dundee provides us some useful 20

information on the property of B-spline during his visit of NUS.

References [1] R. H. Bartels, J. C. Beatty, B. A. Barsky, An introduction to splines for use in computer graphics and geometric modeling, Los Altos, Calif. : M. Kaufmann Publishers, 1987. [2] M. Unser, A. Aldroubi, S. J. Schi, Fast implementation of continuous wavelet transforms with integer scale, IEEE Trans. SP., vol. 42, no. 12, pp. 3519-3523, 1994. [3] M. Unser, A. Aldroubi, M. Eden, On the asymptotic convergence of B-spline wavelets to Gabor functions, IEEE Trans. IT., vol. 38, no. 2, pp. 864-872, 1992. [4] W. Lawton, S. L. Lee, Z. Shen, Characterization of compactly supported re nable splines, Advances in Computational Mathematics, vol. 3, pp. 137-145, 1995. [5] T. N. T. Goodman, In ections on curves in two and three dimensions, Computer Aided Geometric Design, vol. 8, pp. 37-50, 1991. [6] Wang Yuping, Design of multiscale edge detection lters: a B-spline based approach, preprint, Wavelets Strategic Research Programme, National University of Singapore. [7] Wang Yuping, S. L. Lee, Scale-space derived from B-splines, preprint, Wavelets Strategic Research Programme, National University of Singapore. [8] J. Babaud, A. P. Witkin, M.Baudin, R. O. Duda, Uniqueness of Gaussian kernel for scale-space ltering, IEEE Trans. PAMI., vol. 8, pp. 26-33, 1986. [9] J.-S. Lee, Y.-N. Sun, and C.-H. Chen, Multiscale corner detection by using wavelet transform, IEEE Trans. IP., vol. 4, no. 1, pp. 100-104, 1995. [10] A. Asada, M. Brady, The curvature primal sketch, IEEE Trans. PAMI., vol. 8, no. 1, pp. 2-3, 1986. [11] A. Rattarangsi, R. T. Chin, Scale-based detection of corners of planar curves, IEEE Trans. PAMI., vol. 14, no. 4, pp. 430-449, 1992. [12] F. Mokhtarian, A. K. Mackworth, A theory of multiscale, curvature-based shape representation for planar curves, IEEE Trans. PAMI., vol. 14, no. 8, pp. 789-805, 1992. 21

[13] Bingcheng, Li, Repeatedly smoothing, discrete scale-space evolution and dominant point detection, Pattern Recognition, vol. 29, no. 6, pp. 1049-1059, 1996. [14] M. S., Chong, R. K. L. Gay, H. N. Tan, J. Liu, Automatic representation of ngerprints for data compression by B-spline functions, Pattern Recognition, vol. 25, no. 10, pp. 1199-1210, 1992. [15] G. Sapiro, A. Cohen, A. M. Bruckstein, A subdivision scheme for continuous-scale B-splines and ane-invariant progressive smoothing, Journal of Mathematical Imaging and Vision, vol. 7, pp. 23-40, 1997. [16] L. M., Reissell, Wavelet multiresolution representation of curves and surfaces, Graphical Models and Image Processing, vol. 58, no. 3, pp. 198-217, 1996. [17] Gene C. H. Chang, C. C. Jay, Kuo, Wavelet descriptor of planar curves: theory and applications, IEEE Trans. IP., vol. 5, no. 1, pp. 50-70, 1996. [18] P. Saint-Marc, H. Rom, G. Medioni, B-Spline contour representation and symmetry detection, IEEE Trans. on PAMI, vol. 15, no. 11, pp. 1191-1196, 1993. [19] F. Cohen, J. Y., Wang, Modeling image curves using invariant 3-D object curve models-a path to 3-D recognition and shape estimation from image contours, IEEE Trans. on PAMI, vol. 16, no. 1, 1994. [20] M. Gage, R. Hamilton, The heat equation shrinking convex plane curves, J. Dierent. Geom., vol. 23, pp. 69-96, 1986. [21] M. Grayson, The heat equation shrinks embedded plane curves to round points, J. Dierent. Geom., vol. 26, pp. 285-314, 1987. [22] A. Papoulis, Probability, Random Variables and Stochastic Processes, McGraw-Hill, 1972. [23] G.W. Wornell, Wavelet-based representations for the 1/f family of fractal processes, Proceedings of the IEEE, Vol. 81, No. 10, October 1993. [24] F. Attneave, Some informational aspects of visual perception, Psychol. Rev. 61, pp. 183-193, 1954. [25] W. Richards, B. Dawson, D. Whittington, Encoding contour shape by curvature extrema, J. Opt. Soc. Am. A, vol. 3, no. 9, pp.1483-1491, 1986. 22

original data

shrinkage data

200

100

0

−100

−200 −200

250

200

200

150

150

100

100

50

50

0

0

−50

−50

−100

−100

−150 200 −100

0

shrinkage data

250

shrinkage data

0

100

200

−150 −100

shrinkage data 250

250

200

200

200

150

150

150

100

100

100

50

50

50

0

0

0

−50

−50

−50

−150 −100

−100 0

100

200

−150 −100

100

200

shrinkage data

250

−100

0

−100 0

100

200

−150 −100

0

100

200

Figure 1: Shrinkage of the Koch snow ake in the B-spline scale-space. As the scale level gets larger, the Koch snow ake becomes more and more smooth and shrinks to a circle. The evolution of a curve has a similar behavior as that of Gaussian scale-space except that the shrinking speed is dierent. ... .. ... .. .. .. ... .. ... ... .. ... ... ... .. .. .... ....... . . ... ... . . .. ... .. . . . .. ... ... ... ... ... .. ... . . ... ... . ... . ... ... . . . ... ... . .. ... .. . . ... ... . . .. ... . .. . .. ... . . . .. . ... . . .. . ... . . . . . . ... . . . . . . ... . . . . . . ... . . . .. . ... . . . . . .. ... . . . . . ... .. . . . . ... ... ... ... ... ... ... ... ... .. ... . . . . .. .. ... ... ... ... ... .. .... ... ... ... . . . . .. .. ... ... ... ... ... .. ........ .... ...... .... . ... . . . . . . ..... .. ... ...... ... ... ... ...... ... ...... ... ... ...... ... ... ...... ... . . . . . . . . . . ... ...... .... ... ...... ... ... ...... .... ... ..... ... ... .... ...... ..... .... ... . . . . . . . . . . ... ... ..... ..... .... ... ..... ..... ... ..... ..... ... ..... ..... ..... ... ...... ..... ... . . . . . . . . . . . . . ... ... ..... ....... ..... ... ...... ...... ... ........ ..... ... ......... ...... ... ........ ...... ... .......... . . . . . . . . . . . . . . . . . .. .. ..... .......... ... ...... ............. .......... .................. ........................... ............... ............... . . . . . . . . . . . . . . .......... ............... ............... ............... ............... ...............

P

i k

+1

P

i k

P ?1 i k

Figure 2: Post-processing of the corner candidates. If the corner candidates Pki?1 , Pki and Pki+1 are located in the same straight line, then the point Pki is removed out. 23

Original data 35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

Figure 3: Original terrain map curve with 2048 sampling points. Curvature function at multiple scales 1000

0

−1000

0

500

1000

1500

2000

2500

0

500

1000

1500

2000

2500

0

500

1000

1500

2000

2500

0

500

1000

1500

2000

2500

500

0

−500 100

0

−100

off

100 50 0 −50

Figure 4: Curvature function at scale 1, 3, 5, 7. The curvature scale-space image is a hierarchical organization of the in ection points of the curve as indicated by the zero-crossings the curvature function at multiple scales. In practice, only a discrete set of scales are sucient for application. Using B-spline technology, we can compute the scale-space image eciently at any rational scales. 24

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

Figure 5: Detected corner points in scale 1, which are denoted as cross \+ \ in the curve. The number of detected corner points is 221.

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

Figure 6: Recovered curve using spline tting. The recovered curve is drawn in dashed dot line. Compression ratio= 9.2670, error=0.0035. 25

35.2

35.15

35.1

35.05

35

34.95 138.55

138.6

138.65

138.7

138.75

138.8

138.85

138.9

138.95

139

Figure 7: Zoom in a small portion of the curve to show the multiscale property of the shape descriptor.

26

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

Figure 8: Curve data compression at scale 3. (a) detected corner points. number=210. (b) recovered curve. Compression ratio=9.7524, error=0.0031.

27

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

Figure 9: Curve data compression at scale 5. (a). detected corner points. number=160. (b) recovered curve. Compression ratio=12.8, error=0.0045.

28

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

138.6

138.7

138.8

138.9

139

139.1

139.2

139.3

139.4

139.5

35.4

35.3

35.2

35.1

35

34.9

34.8

34.7

34.6

34.5 138.5

Figure 10: Curve data compression at scale 7. (a) detected corner points. number=135. (b) recovered curve. Compression ratio=15.1704, error=0.0057.

29