Scale-Space Derived From B-Splines.
AbstractIt is well-known that the linear scale-space theory in computer vision is mainly based on the Gaussian kernel. The purpose of the paper is to propose a scale-space theory based on B-spline kernels. Our aim is twofold. On one hand, we present a general framework and show how B-splines provide a flexible tool to design various scale-space representations: continuous scale-space, dyadic scale-space frame, and compact scale-space representation. In particular, we focus on the design of continuous scale-space and dyadic scale-space frame representation. A general algorithm is presented for fast implementation of continuous scale-space at rational scales. In the dyadic case, efficient frame algorithms are derived using B-spline techniques to analyze the geometry of an image. Moreover, the image can be synthesized from its multiscale local partial derivatives. Also, the relationship between several scale-space approaches is explored. In particular, the evolution of wavelet theory from traditional scale-space filtering can be well understood in terms of B-splines. On the other hand, the behavior of edge models, the properties of completeness, causality, and other properties in such a scale-space representation are examined in the framework of B-splines. It is shown that, besides the good properties inherited from the Gaussian kernel, the B-spline derived scale-space exhibits many advantages for modeling visual mechanism with regard to the efficiency, compactness, orientation feature, and parallel structure.
Introduction
Scale is a fundamental aspect of early image representation. Koenderink [1] emphasized that the
problem of scale must be faced in any imaging situation. A multiscale representation is of crucial
importance if one aims at describing the structure of the world. Both the psychophysical and
physiological experiments have confirmed that multiscale transformed information appears in the
visual cortex of mammals. This leads to the motivation for the interpretation of image structures in
terms of spatial scale in computer vision. Some researchers such as Burt and Adelson [6], Koenderink
[1], Marr [10], Witkin [11], Rosenfeld [2] had exposed the necessity and advantages of using operators
of different sizes for extracting multiscale information in an image. For a more detailed review, see
[3].
The Gaussian scale-space approach of a signal as introduced by Witkin [11], is an embedding
of the original signal into a one-parameter family of derived signal constructed by convolution with
Gaussian kernels of increasing width. One reason the traditional scale-space is mainly based on
the Gaussian kernel is that the Gaussian function is the unique kernel which satisfy the causality
property as guaranteed by the scaling theorem [14], [15], [16], [19]; it states that no new feature points
are created with increasing scale. Another reason is that the response of human retina resembles a
Gaussian function. Neurophysiological research by Young [22] has shown that there are receptive
field profiles in the mammalian retina and visual cortex, whose measured response profiles can be
well-modeled by superposition of Gaussian derivatives. Therefore, the Gaussian function is suitable
for modeling human visual system.
In practice, since the computational load becomes extremely heavy when the scale gets larger,
many techniques are proposed for efficient implementation of scale-space filtering. Among them
B-splines or binomials have been widely used to approximate the Gaussian kernel. Such examples
include, Wells [5], Ferrari et al. [8], [9], Poggio et al. [20], Unser, Aldroubi and Eden [25], [26], etc.
For a more compact representation, pyramid technique is another widely used representation that
combine the sub-sampling operation with a smoothing step. Historically they have yielded important
steps towards the scale-space theory. The low pass pyramid representation proposed by Burt [6], [7]
is a famous example, which is also closely related to B-spline techniques.
The general idea of representing a signal at multiple scales is not new to us. It is through wavelet
theory that these early ideas have been well formulated and refined. In fact, this is largely due
to the contribution of B-spline techniques. As will be shown later, the orthogonal multiresolution
pyramid originally proposed by S. Mallat [40] and the biorthogonal pyramid [28], [46], [47] in wavelet
theory can all be derived from B-splines [30], [33], [36]. Other types of wavelets such as the wavelets
on a interval [48], the periodic wavelets [37] and the cardinal spline wavelets [29] are all related to
B-splines.
Motivated by these observations, the purpose of this paper is to build a more general framework
of scale-space representation in the context of B-splines as an improvement of the traditional scale-space
theory. By and large, the paper is divided into two parts. Firstly, we present a systematic
development of the scale-space representation in the framework of B-splines. In particular, we
focus on two classes of scale-space design. It is shown that if an image is represented as a B-spline
efficient subdivision algorithm can be designed to give a geometric description at
varying degrees of detail. The various scale-space representations are derived from B-splines in
different forms. For continuous scale-space representation, a general algorithm is derived for fast
and parallel implementation at rational scales. Several classical fast algorithms are shown to be
special cases under some conditions. Differential operators have been used for multiscale geometric
description of an image. However, it is not clear whether an image can be synthesized from these
differential descriptors. Using B-spline techniques frame algorithms are designed to express the
image as combinations of multiscale local partial derivatives. These operators include the gradient
operator, second directional operator, Laplacian operator and oriented operators. At the same time
the intrinsic relationship between wavelet theory and the traditional linear scale-space approach is
exhibited. Although B-spline has been used in practice in place of Gaussian, there is little effort
to consider its scale-space behavior directly. Therefore the second part of the paper is devoted
to examining the properties of B-spline derived scale-space. In particular, the advantages of such
a scale-space representation are highlighted. It was shown that the B-spline derived scale-space
inherits most of the nice properties of the Gaussian derived scale-space. Nevertheless, the B-spline
kernel outperforms the Gaussian kernel in that it can provide more meaningful, efficient and flexible
description of image information for multiscale feature extraction.
The organization of the paper is as follows. In Section 2, some fundamental properties of B-splines
are reviewed, which also explain why a B-spline kernel is a good kernel for scale-space design.
Following this section, we categorize three types of scale-space representation. In particular, we focus
on the design of continuous scale-space and dyadic scale-space frame representation. The relation
with compact scale-space is discussed. From the viewpoint of B-splines the evolution of wavelets
from classical continuous scale-space is well understood. Moreover, the equivalence of several famous
scale-space methods is explored. In Section 4 a general procedure is presented to study the edge
models in the B-spline scale-space. In Section 5 some basic properties of the B-spline derived scale-space
are investigated in parallel with that of the Gaussian derived scale-space. These include the
the completeness property, causality property, orientation feature and so on. Conclusions are given
in Section 6. Finally, Appendix A, B are given in Sections 7, 8 respectively in order to make the
paper mathematically complete.
2.1 Notations and definitions
We adopt the convention of [25]. Let L 2 (R) be the Hilbert space of measurable, square integrable
functions on R and l 2 (Z) be the vector space of square summable sequences. We denote the central
continuous B-spline of order n by fi n (x), which can be generated by repeated convolution of a
B-spline of order 0,
z -
where the 0th-order B-spline fi 0 (x) is the pulse function with support [
]. The Fourier transform
of fi n (x) is
For an integer m - 1, fi n
m (x) is defined as the nth-order continuous B-spline dilated by a scale
factor m, i.e.,
The discrete B-spline of order n at scale m is defined as
z -
is a normalized sampled pulse of width m.
The discrete sampled B-spline b n
m (k) of order n and integer coarseness m - 1 is obtained by
directly sampling the nth-order continuous B-spline at the scale m:
We write b n
. Consequently, the frequency response of the directly sampled B-spline is an
aliased version of the frequency response of the continuous B-spline since
which is due to the Possion's summation formula.
The discrete convolution between two sequences a and b in l 2 (Z) is the sequence b   a:
Under this definition, the convolution is commutative. The convolution inverse (b) \Gamma1 of a sequence
b is defined by
where ffi(k) is the unit impulse whose value is 1 at 0 and 0 elsewhere.
The decimation operation [b] #m down-samples the sequence b by the integer factor m, i.e.,
Conversely, the operator [b] "m up-samples by the integer factor m, i.e., it takes a discrete signal b
and expands it by padding consecutive samples:
0; otherwise:
2.2 The similarity between Gaussian and B-spline
B-splines are good approximations of the Gaussian kernel which is commonly used in computer
vision. This is the consequences of the central limit theorem. For a review of the B-spline window
and other famous filters, see Torachi et al. [39]. In [28] Unser, Aldroubi and Eden have presented a
more general proof that B-splines converge to the Gaussian function in L p (R); 8p 2 [2; +1) as the
order of the spline n tends to infinity. Since the variance of the n-th order of B-spline is n+1, the
approximation relation is as follows:
Furthermore, by numerical computation [28] it was shown that the cubic B-spline is already near
optimal in terms of time/frequency localization in the sense that its variance product is within 2% of
the limit specified by the uncertainty principle. A graphical comparison between the Gaussian kernel
and the cubic B-spline is given in Fig. 1. Moreover, both the physiological and biological experiments
[22] have shown that the human visual system can be modeled with the Gaussian kernel. Therefore,
B-splines are also suitable for modeling biological vision due to their close approximation to the
Gaussian kernel.
2.3 Stable hierarchical representation of a signal by B-splines
Another significant property of the B-spline of a given order n is that it is the unique compactly
supported refinable spline function of order n which can provide a stable hierarchical representation
of a signal at different scales. It has been proven [38] that a compactly supported spline is m-refinable
and stable if and only if it is a shifted B-spline. Let h ? 0 and define the polynomial spline space S n
consisting of the dilated and shifted B-splines of order n (n is odd, which we will assume throughout
the paper) by
Then
and
The embedding property (13) follows from the fact that the B-spline fi n (x) is m-refinable, i.e., it
satisfies the following m-scale relation,m fi n ( x
The m-refinability of the B-splines can be easily verified [29], [38]. It also establishes the intrinsic
relationship between the continuous B-spline and discrete B-spline. If we take it is just the
commonly used two-scale relation and B n
2 (k) is the discrete binomial. H. Olkkonen [53] has used the
binomial kernels for designing multiresolution wavelet bases.
From the m-scale relation (15), we can also establish the relationship between the discrete sampled
B-spline and the discrete B-spline:
Since B-splines provide a stable multiresolution representation of a signal at multiple scales, it
is preferable to select B-splines as smoothing kernels to extract multiscale information inherent in
an image. Therefore, it is not surprising that many vision models [40], [29], [30] are derived from
B-splines.
One can refer to [25], [26], [33], [47] for a more complete and extensive exposition of the B-spline
methods. For example, its minimal support and m-refinability properties have led to fast
implementation of the scale-space algorithms [27], [5], [8], [9] in computer vision. In the following
section, we will classify scale-space representations into three types and show how B-splines are used
as a flexible tool for designing an efficient visual model according to the above requirements.
3 Scale-space representations designed from B-splines
In this section we focus on the fast implementation of continuous scale-space filtering and the design of
dyadic scale-space frame representation. Their relations with the compact scale-space representation
or compact wavelet models are indicated.
3.1 Implementation of continuous scale-space filtering using B-splines
3.1.1 Discrete signal approximation using B-spline bases
In practice, a discrete sets of points are given. Because spline spaces S n
provide close and stable
approximations of L 2 , it is reasonable to parameterize the discrete signal or image using B-spline
bases. We use the translated B-splines of order n 1 as bases to approximate the signal,
i.e., the signal f(x) 2 L 2 (R) is projected into the spline space S n 1
1 at resolution 1. In (17) we have
assumed that the sampling rate is 1 for convenience. We can call this procedure the generalized
sampling of the original data, where the sampling basis function is taken as the B-spline. There are
different types of approximation (see [25]). A common approach ([25], [26]) is the direct B-spline
transform where the exact or reversible representation of a discrete signal f(k) in the space of B-splines
is obtained by imposing the interpolation condition: 8k 2 Z; ~
f(k). Thus, the
coefficient c(k) can be computed as,
where (b n
denotes the inverse of the discrete sampled B-spline which can be computed recursively.
3.1.2 Fast algorithm for continuous scale-space filtering at rational scales
In this section we use the B-splines to derive a filter bank algorithm for fast implementation of
continuous scale-space filtering.
The linear scale-space representation is to make a map of a signal at multiple scales by changing
the scale parameter continuously. In the language of wavelet transform the traditional scale-space
approach can be regarded as a continuous wavelet transform of the signal f 2 L 2 ,
Z
is the scaled wavelet. Because the geometric features of an image are
characterized using differential descriptors, /(x) is often taken as certain derivative of a smooth
kernel or has certain order of vanishing moments. Under different physical meanings various linear
scale-space representations are proposed [4] where different kernels / are assumed. Here we also use
the B-spline of order n 2 to approximate the wavelet /(x):
In the classical scale-space theory, two frequently used multiscale edge detection filters are the
famous Marr-Hildreth operator [10] and Canny operator [12], which are obtained by taking the first
and second derivative of the Gaussian kernel respectively. Since B-splines are good approximations
of the Gaussian kernel, we shall use the derivatives of B-splines instead. In such cases, the coefficients
in (20) are the coefficients of first and second order difference operators respectively, i.e.
Such spline wavelets are shown in Fig. 2. Using these spline wavelets we obtain the approximate
Marr-Hildreth operator [10] and Canny operator [12] respectively. Higher order of derivatives of B-splines
can detect edges with higher singularity [42] and the coefficients g are the binomial-Hermite
sequences. Explicitly, in the Fourier domain the rth-order difference of the B-spline of order n can
be written as
which is also the rth-order derivative of the B-spline of order n + r. We remark that B-splines can
also efficiently approximate other kinds of wavelets such as the generalized edge detectors in the
-space representation [4] and the coefficients g can be computed numerically.
Since a real number can be approximated arbitrarily close by a rational number
Z, we take rational scale and derive a general filter bank implementation of the scale-space representation
of (19) using the m-refinable relation (15). The cascaded implementation of (19) with /
given in (20) and f approximated by (17) is,
The derivation of this algorithm is given in the Appendix A. This algorithm extend that in [27]. If
the scale is taken as an integer, i.e., when then the resulting formula is similar to that in
[27]:
c   g "m 1
The implementation of the above algorithm is illustrated in the block diagram in Fig. 3.
In the filter bank implementation of (23), we can interchange the order of convolution. The result
is then equivalent to the discrete B-spline filtering of the difference of the discrete sampled signal
(with a down-sampling)
is the signal sampling followed by difference operation
We note that the computational complexity is mainly due to the discrete B-spline filtering (25)
which can be implemented efficiently. By (4), it turns into the cascaded convolution with the 0th-
order of discrete B-spline, and can be implemented via running average sum technique. If we define
such a running average operation as
R
then it can be realized using the following iterative strategy
Therefore, starting from the initial coefficient R 0 using only the
addition. Then after a down-sampling with a factor m 2 , we obtain the scale-space filtering at the
rational scale m1
m2 . Suppose m 2 is fixed, as usually the case in practice, the computational cost is
independent of the scale. Fig. 4 shows an example of the scale-space filtering of a simulated signal
using the above algorithm.
We record two important properties of this procedure.
ffl Efficiency: In practice, m 2 is usually fixed. The computational complexity at each scale m 1
is O(N ). The computational complexity is largely due to the convolution with the smoothing
B-spline kernel which can be reduced by running average sum technique. In contrast
with the existing procedure based on direct numerical integration or FFT-based scheme, the
computational complexity does not increase with the increasing number of values of the scale
parameter.
ffl Parallelism: The structure of the above algorithm is parallel and independent across scales.
This makes it inexpensive to run on arrays of simple parallel processors. In other words this
can be ideally suited for VLSI implementation.
One can recall that the above subdivision algorithm is also similar to the ' a trous algorithm [49],
[32] for fast computation of continuous wavelet transform except with some constraints on their
filters. However, the ' a trous algorithm can only compute the wavelet transform at dyadic scales.
The above algorithm can compute the wavelet transform efficiently at any scales. As will be shown in
Section 3.2, if the scale is restricted to dyadic, the above algorithm is similar to the ' a trous algorithm.
However, using B-splines more efficient scheme can be obtained which only need addition operation.
3.1.3 Extension to 2D images
Although the above algorithm is derived in the one-dimensional case, it can be easily extended to two
dimensions. The inner product B-spline fi n (x; used as a basis to parameterize
the image and approximate the two-dimensional wavelet kernels. For example, we can use the tensor
product B-splines to approximate the Maar-Hildreth's LoG operator [10]:
where
1). Since this two-dimensional
kernel is represented by the separable one-dimensional B-spline bases, by performing the above one-dimensional
fast algorithm along the horizontal and the vertical orientation respectively the LoG
operator can be computed efficiently. Fig. 5 shows such results for Lena image at 3 different scales.
3.2 Dyadic scale-space frame representation
3.2.1 One-dimensional signal representation by its local partial derivatives at dyadic
scales
The above continuous scale-space is too redundant for some applications. In addition, as stated by
Witkin [11], an initial representation ought to be as compact as possible, and its elements should
correspond as closely as possible to meaningful objects or events in the signal-forming process. A
description that characterizes a signal by its extrema and those of its first few derivatives is a
qualitative description to "sketch" a function. If we sample the scale of the above continuous scale-space
as dyadic while keeping the time variable continuous, we can obtain a more compact scale-space
representation. Moreover, such a representation is shift invariant and therefore is suitable for some
pattern recognition applications. In particular, using B-spline techniques efficient frame algorithms
can be designed to express the signal in terms of its local partial derivatives.
We now show the relationship between this type of transform and the above continuous scale-space
implementation.
If we use the approximation b n 2 +n 1 +1   c to replace the original signal f , and the width of the
B-spline m is restricted to be dyadic, say 2 m , we will get a recursive relation for the wavelet transform
between the dyadic scales. In this case, the formula (23) becomes
where
It is easy to derive the two-scale relations for
. From property (4) and the property of
z-transform,
is the z-transform of B n
We have the following relation
z \Gammaj 2
In the time domain, it becomes
By this relation, we can get a fast recursive implementation of dyadic-scale space filtering,
Or simply written as
is the binomial kernel.
The approximate Marr-Hildreth operator and Canny operator at dyadic scales can now be computed
as
"2
where is the first or second order of difference operator given in (21). The recursive
refinement relation (36) and (37) can be rewritten in the z-transform domain as
2: (39)
By requiring the reconstruction filter ~
G to satisfy the following perfect reconstruction condition
we can reconstruct the signal from its multiscale partial derivatives
where ~
are the time responses of ~
G k which are given explicitly in Appendix B. Since all of
these filters are linear combinations of binomial and divided by 1
using Pascal triangular algorithm
only addition operation and bit shift operation are needed. This is very suitable for hardware
implementation.
3.2.2 Image representation by its local directional derivatives
The tensor product B-spline basis is fi n (x; (y). It is interesting that we can still derive
an efficient frame algorithm to characterize an image from its local differential components. A fast
algorithm for the gradient case has been proposed [42] and further refined in [34]. Now we consider
the case of the second directional derivative.
For edge detection an approach is to detect the zero-crossings of the second directional derivative
of the smoothed image f   fi n
along the gradient orientation
@x
@y
We can still derive a subdivision algorithm to compute the three local partial derivative components
or wavelet transforms:
where the three directional wavelet components are defined in the Fourier domain as
where G (1) and G (2) are the transfer functions of the first and second order difference operator. From
these definitions we can obtain a recursive algorithm for the computation of the three local partial
derivative components:
where I   (h; g) "2 j \Gamma1 represents the separable convolution of the rows and columns of the image with
the one dimensional filters [h] respectively. The symbol d denotes the Dirac filter
whose impulse is 1 at the origin and 0 elsewhere. We can also reconstruct the image from these
dyadic wavelet transforms using the following formula
2, is the FIR of the transfer function (!). The
reconstruction formula (48) follows from the following perfect reconstruct identity:
G (2) (! x )U(! y
G (2) (! y )U(! x )
G (1) (! y
If we define the three corresponding reconstruction wavelets
G (2) (! x )H 2 (! y
G (2) (! y )H 2 (! x
G (1) (! x ) ~
G (1) (! y
it can be shown [34] that an image f(x; y) can be represented as
One can also notice that in the above decomposition and reconstruction formula all the filters are
binomial which require only the addition operation. For illustration Figure 6 shows the above decomposition
and reconstruction results of a square image at the scales 1, 2, 4. Like the compact wavelet
decomposition [40], the above algorithms also decompose an image into horizontal, vertical and diagonal
components. However, this transform has explicit physical meaning and is shift-invariant. This
can be useful for certain pattern recognition tasks.
3.2.3 Image representations by its isotropic and multi-orientational derivative component

One can obtain a more compact isotropic wavelet representation of an image that is complete and
efficient using a radial B-spline as the smoothing kernel in two dimensions. This representation is
important because it indicates that an image can be recovered from its multiscale LoG-like compo-
nents. The radial B-spline OE(x; y) is a non-separable function of two variables defined by its Fourier
where the radius ae = min(
and the wavelet is defined by
One may notice that such a wavelet is isotropic and LoG-like, which resembles the human visual
system. With these definitions we still have a filter bank implementation of the decomposition and
reconstruction. We omit the details, and just give the decomposition formulas,
r are the 2D non-separable radial filters corresponding to h and g respectively. In this
decomposition, two components are obtained at each resolution. By designing the filter ~
g r from the
relation (40) the corresponding 2D non-separable radial filter ~
r can be computed numerically via
its Fourier transform ~
G r . Then the reconstruction formula similar to (41) can be obtained. Also, it
is easy to check using the same arguments as the in 1D case that an image can be represented as [34]
is the 2D reconstruction wavelet defined by -
G r (ae) -
OE(ae).
One can build a wavelet representation having as many orientation tunings as desired by using
non-separable wavelet bases. A generalized Pythagorean theorem has been proved to decompose an
image into a finite number of equally spaced angles [34]:
(2m\Gamma1)!!n
1). If we multiply the above isotropic wavelet
(51) by the angular part H k Z, we can extract the orientational information
in the dyadic scale-space tuned to n orientations [34]:
Such wavelets can be called orientation tuned LoG-like filters. An image can be represented by its
multiscale and multi-orientational components,
is the oriented wavelet for reconstruction. Similarly, the pyramid-like filter bank implementation
of such a representation can be obtained as follows,
Therefore, through such an approach we can analyze the directional information of an image feature
at a certain angle in dyadic scale space. In Fig. 7 we show a multiscale orientational decomposition
and reconstruction, where the orientation number is chosen as 3.
3.2.4 Some comments on the application of dyadic scale-space representation
The dyadic scale-space frame representation from B-spline gives rise to many applications, since
it provides a invertible, translation-invariant, and pyramid-like compact representation of a signal.
One example is the fingerprint based compression [41] by combining with other techniques. Many
image features such as ridges, corners, blobs, junctions are usually characterized by local differential
descriptors [3]. It is usually enough to consider their behaviors only at dyadic scales. The proposed
algorithms provide efficient ways and are easy for hardware implementation. For example, in multi-scale
shape representation usually the computation of the curvature function is treated in continuous
scales [51]. In fact sometimes it is enough to consider its behavior only at the dyadic scales [50]. We
have used the above algorithms to efficiently compute the geometric descriptors for multiscale shape
analysis [35].
3.3 Compact scale-space representation
While the dyadic scale-space frame approach provides a more compact representation, it is still over-complete
for signal representation. For image compression applications, compact representation is
preferred. In order to give a complete picture we mention briefly the discrete wavelet transform.
While the scale-space technique has existed for a long time, it was the orthogonal multiresolution
representation proposed by S. Mallat [40] that makes the mathematical structure of the image more
explicit. This is an extension or refinement of traditional scale-space theory. This approach restricts
the scale to dyadic and samples the time variable. The starting point is to orthogonalize the B-spline
basis, and then decompose the signal approximated at a fine scale space S n
space S n
by imposing the orthogonal condition
The detail irregular information of the signal is contained in the subspace W n
This defines an
orthogonal multiscale representation. After the B-spline basis is converted into an orthogonal basis,
the two-scale relation still exists which results in an efficient pyramidal algorithm. The perfect reconstruction
condition (40) still exists. However, additional conditions on the filters H; G; ~
G are imposed
to ensure the orthogonality. There are several ways to achieve a compact multiresolution by imposing
the biorthogonal instead of the orthogonal condition (57). All these compact multiresolutions are
related to B-splines. A detailed study can be found in [30], [36].
From the above analysis, it is easy to see that the dyadic scale-space frame representation lies
between the continuous scale-space and the compact representation. Which kind of representation
to select depends on the problem at hand. For example, in multiscale feature extraction one may
compute the differential operation either at the continuous scales or only dyadic scales. Therefore,
the continuous or dyadic scale-space frame representation is more useful. However, for compression
applications, the compact multiresolution is the favorite.
3.4 Relations between the existing scale-space algorithms in computer vision
Before the appearance of wavelets, the B-spline technique has been widely used in computer vision.
Examples include Wells [5], Burt [6], [7], Ferrari et al. [8], [9]. We shall show that under certain
circumstances, they are either equivalent or the special cases of the general algorithm given in
Section 3.
3.4.1 Relation to Ferrari et al.'s method
Ferrari et al. [8], [9] have proposed B-spline functions to realize the 2D image filtering recursively.
Their idea is to use B-splines to as the filter kernel:
l
are the interpolation coefficients at the knots (kM; lN ), which can be computed
using the usual method. Then using the properties of discrete B-spline, they derive the recursive 2D
image filtering. In this way, the computational load can be greatly reduced. Therefore, this approach
is a special case of our proposed algorithm for continuous scale-space filtering.
3.4.2 Relation to Wells' method
In [5], Wells proposed an approach for efficient synthesis of Gaussian filters by cascaded uniform
filters. It is easy to show that his method is equivalent to using cascaded 0th discrete B-spline to
approximate the Gaussian kernel. By this approach, the cascaded convolution with a 0th-degree
B-spline (uniform filters) can be realized by running average sum technique as discussed above.
Obviously, his method is a special case of our recursive algorithm.
3.4.3 Relation to Burt's Laplacian pyramid algorithm
Burt [6], [7] has introduced the following low pass filter for the generation of Gaussian or Laplacian
pyramids
If the parameter a is taken as a=3/8, then w(j) can be re-written as
This is equivalent
to the special case of in the formula (35). In this case, the filter is also equivalent to the operator
used for the generation of the dual cubic spline pyramid representation as discussed in [29].
Edge patterns in the B-spline based scale-space
The study of edge patterns in scale-space is very important for many applications. Much work has
been done on the study of edge patterns in the Gaussian based scale-space. In this section, we want
to investigate the edge behaviors in the B-spline based scale-space.
It is well-known that there exists a similar "uncertainty principle" between good edge localization
and noise removal. At finer scale, better localization can be achieved with the cost of noise pollution,
and vice versa. Many researchers have studied the localization of operators based on Gaussian kernel
in scale-space. Berzins [24] studied the accuracy of Laplacian operator, M. Shah, A. Sood and Jain
[23] considered the localization of pulse and staircase edge models, I. J. Clark [21] investigated the
phantom edges in a scale-space. For corner detection, A. Huertas and H. Asada [51], R. T. Chin
[52] have also considered the behavior of edge models. We found that their derivations are all based
on Gaussian kernel. It is necessary to study the behavior of various edge models in B-spline based
scale-space which give some a priori knowledge of various patterns in an image.
Here we present a more general proof with the assumption that the primitive ' in the definition
dx '(x) is symmetric and with compact support [\Gammaw=2; w=2] and its derivatives have the
shape as in Fig. 2. In practice, the support of the Gaussian kernel is usually truncated to a finite
interval. Obviously, the truncated Gaussian and B-spline kernel are included in this assumption.
We shall show that these edge models have the same behaviors as that derived from the traditional
Gaussian kernel. First we adopt the following accurate definition of an edge [21].
Definition 1. A point x 0 is called an authentic edge of a signal f(x), if jW 1 f(s; x 0 )j is maxima,
will be a phantom edge.
We shall use W 1 f and W 2 f to denote two types of wavelet transform where the wavelets are the
first and second derivatives of ' respectively. I. J. Clark [21] analyzed these two types of detection.
Generally, zero-crossing detection is equivalent to the extrema detection. However, extrema detection
includes both maxima and minima detection. Only the edge point detected by local maxima belongs
to authentic edge and the edge point detected by local minima corresponds to phantom edge. It is
shown that zero-crossing edge detection algorithms can produce edges which do not correspond to
significant image intensity changes. Such edges are called phantom or spurious.
Now, as an example we study the behavior of staircase edge model in scale-space. The staircase
edge model can be represented as
amplitudes of the edge, d is the distance between the two abrupt changes
at is the step function,
whose derivative in the distributional sense is ffi(x). Hence, from (19),
d
s (x)dx
x
s
s
i.e., the wavelet transform is just the sum of two dilated smoothing functions. At the location
@
@x
s
s
s
s
s
Therefore, only at a small scale s ! 2d
w , the location of edge at can be detected exactly.
Similarly, at the location
@
@x
s
d
s
s
s
d
s
i.e., the location of edge at can only be detected exactly at a small scale s ! 2d
w . For large
scale, it will spread like a cone in scale-space and will be influenced by another cone at As a
consequence, the edge location will be mis-detected due to the superposition of two diffused cones.
One may deduce that there will be another point x 0 such that @
due to different
parities in the values of @
at the location of However, it is easy to
find such a point corresponding to a local minimum, which means it is a phantom edge point. This
cannot be distinguished from others by zero-crossing detection. Fig. 8 illustrates the behavior of this
type of edge in the scale-space. That was why local maxima is preferred for edge detection in [41],
[42].
Similarly, for zero-crossing detection of W 2 f(s; x), we can draw the same conclusion. In the
above analysis, we only consider one type of edge model. Other types of edges such as the step,
pulse, ramp, roof can be treated in a similar way [34]. Also, our derivation is based on a more
general assumption on the kernels which include both the truncated Gaussian and the B-spline.
5 Discussions on the properties of the B-spline based scale-space
We now discuss the advantages and properties of the B-spline based scale-space.
ffl Efficiency: It is a basic requirement that any algorithm should be able to capture and process
the meaningful information contained in the signal as fast as possible. Obviously, the major
weakness of the traditional Gaussian based scale-space is the lack of efficient algorithms. On
the contrary, B-spline techniques facilitate computational efficiency. The computational complexity
is scale independent. Moreover, in contrast with the scale-space based on the Gaussian
kernel, the B-spline representation of a signal is determined directly on an initial discrete set
of points, avoiding problems caused by discretization in continuous scale-space. B-splines also
have been used as smoothing windows for efficient computation of Gabor transform to extract
frequency information [31].
ffl Parallelism: Data parallelism is common in computer vision which arises from the nature of
an image. As Koendrink [1] pointed out, it seems clear enough that visual perception treats
images on several levels of resolution simultaneously and that this fact must be important for
the study of perception. In this paper, efficient parallel structure of an image is exhibited using
B-spline techniques. It may provide a good interpretation of human visual system which can
process the hierarchical information simultaneously. B-splines provide a flexible way to process
the multiscale information using either coarse to fine strategy or in parallel. This is also easy
for hardware implementation.
ffl Completeness and invertibility: We usually use the zero-crossings or the local extrema as
meaningful description of a signal. It is clearly important, therefore, to characterize in what
sense the information in an image or a signal is captured by these primal sketches uniquely. For a
Gaussian based scale space, the completeness property is guaranteed by the fingerprint theorem
[13]: the map of the zero crossing across scales determines the signal uniquely for almost all
signals in the absence of noise. Such results have theoretical interest in that they answer the
question of what information is conveyed by the zero and level crossings of multiscale Gaussian
filtered signals. Poggio and Yuille's proof is heavily dependent on the Gaussian kernel and they
conjectured that under certain conditions Gaussian kernel is necessary for fingerprints theorem
to be true. However, later Wu and Xie [17] gave a negative answer and present a more general
proof, which states that fingerprint theorem holds for any symmetry kernel. Therefore, the
fingerprint theorem is also true in the case of B-splines for continuous scale-space
representation.
Differential operators have also been widely used for multiscale geometric description of images,
but it has not been clear such representations are invertible. As shown in the paper, using B-
splines, efficient frame algorithms can be designed to express an image from its local derivatives
at dyadic scales.
ffl Compactness: For compression application, we require a representation to be as compact as
possible so that an image can be represented by the corresponding primitives using less storage.
In [13] Poggio and Yuille conjectured that the fingerprints are redundant and the appropriate
constraints derived from the process underlying signal generation should be used to characterize
how to collapse the fingerprints into a more compact representation. In the paper, the more
compact dyadic scale-space representations is proposed. We can use such representation for
compression applications by combining with other techniques.
ffl Causality: Since edge points are important features, it is natural to require that no new
features are created as the scale increases. A multiscale feature detection method that does
not introduce new features as the scale increases is said to possess the property of causality.
Causality is in fact equivalent to the maximum principle in the theory of parabolic differential
equations [18]. The Gaussian scale-space is governed by the heat diffusion equation and therefore
possesses the causality property. Such continuous causality property of the Gaussian kernel
is not shared by the B-spline. However in the discrete sense, M. Aissen, I. J. Schoenberg and
A. Whitney [54] proved that for a discrete scale-space kernel h, the number of local extrema or
zero-crossings in f out = h   f in does not exceed the number of local extrema or zero-crossings
in f in if and only if its generating function H(z) =P
h(n)z n is of the form
where
It is easy to verify that the discrete B-spline kernel in (25) satisfies such a condition. Therefore,
the causality property still holds for discrete B-spline filtering in the discrete sense.
The number of local extrema or zero-crossings of the derivative of the discrete signal does not
increase after running average sum. This justify the use of the discrete smoothing kernel in
practice.
ffl Orientation: Orientation analysis is an important task in early vision and image processing,
for example in texture analysis [45]. The Laplacian multiresolution [6] does not introduce
any spatial orientation selectivity into the decomposition process. Daugman [45] showed that
these impulse response can be approximated by Gaussian windows modulated with a wave.
It is meaningful to combine both the orientation analysis and the scale feature [43]. In this
paper, efficient pyramid-like algorithm is designed using B-spline technique to analyze and
synthesize an image from its multi-orientational information at any number of angles in the
dyadic scale-space. Note that the usual wavelet transforms can decompose an image in only
three orientations.
There are other advantages of the B-spline kernel. In the time-frequency analysis, the Gaussian
kernel is the optimal function that minimizes the uncertainty principle. The cubic B-spline is already
a good approximation to the Gaussian function [28], see also Fig. 1. As the order goes infinity, the
B-splines converge to the Gaussian in both the time and frequency domain. Moreover, a B-spline
resembles the response of receptive field [22] and is also suitable for modeling the human visual
system. Edge detection is an ill-posed problem. From the view point of regularization theory, cubic
spline is proved optimal. The connection between the regularized edge detection and the smoothing
spline problem proposed by Schoenberg, Reinsh in statistics is noted by Poggio et al. [20]. It was
shown that the cubic B-spline rather than the Gaussian kernel is optimal for edge detection.
B-splines are the shortest basis functions that provide a stable multiresolution analysis of a signal
[33]. This explains why many wavelet models of a vision [40], [28],[46], [47], [48], [37] are derived
from B-splines [33], [36]. For the derivative operations, B-spline approach is very intrinsic which
elucidates the relationship between derivative and difference which are usually characterized by the
two-scale difference relations. B-splines play an important role to bridge the traditional scale-space
theory, dyadic scale-space frame and compact multiresolution representation.
6 Conclusions
This paper describes a B-spline based visual model. For a long time, the Gaussian kernel has been
commonly used in computer vision. In this paper a general framework for scale-space representation
using B-splines is presented. In particular, the design of two types of scale-space representations
is given in detail. A fast algorithm for continuous scale-space filtering is proposed. In the case
of dyadic scale some efficient frame algorithms are designed to express the image from its local
differential descriptors. The intrinsic relationship with the compact wavelet models is also indicated.
Several algorithms are proved to be special cases of our proposed algorithm.
To our knowledge, the scale-space property based on B-splines has not been fully studied before.
We examine the property of B-spline based scale-space in parallel with Gaussian kernel. Our results
indicate that B-splines possess almost the same properties as the Gaussian kernel. Moreover, the
B-spline kernels outperform the Gaussian in many aspects, notably, computational efficiency.

Acknowledgment

The authors wish to thank the referees for their comments which greatly
improve the presentation of the paper. The first author wishes to express his thanks to Prof. Wu
for providing the reference [17].

Appendix

A: Derivation of fast implementation of continuous
wavelet transform at rational scales.
We use the m-scale relation (15) to derive the filter bank implementation of scale-space filtering at
rational scales. From (19), (20), (17),
l
l
Using the m-refinable relation (15),
Hence,
(j)B n1
(j)B n1
where the following property of B-spline is used [25]:
Substituting (70) into (69) gives
l
l
l
If we take then the above formula can be written as,
If we take
we get an interpolation formula, and the size of the transformed data m 2
is as large as that of the original sampling signal data,

Appendix

B: Derivation of the reconstruction filter responses
For the first order difference the perfect reconstruction condition (40) gives
~
and the corresponding FIR is
j+l
j+l
\Gamman - l -
j+l
For the second order difference
~
4z
and the corresponding FIR is



--R

The structures of images
Edge and curve detection for visual scene analysis


Efficient synthesis of Gaussian filters by cascaded uniform filters
The Laplacian pyramid as a compact image code
Fast hierarchical correlations with Gaussian-like kernels
Efficient two-dimensional filters using B-spline functions
Recursive algorithms for implementing digital image filters
Theory of edge detection

A computational approach to edge detection
Fingerprints theorems for zero-crossings
Scaling theorems for zero-crossings
Uniqueness of Gaussian kernel for scale-space filtering
Scaling theorem for zero-crossings

Reconstruction from zero-crossings in scale-space
Scaling theorems for zero-crossings of bandlimited signals
Computational vision and regularization theory
Authenticating edges produced by zero-crossing algorithms
The Gaussian derivative model for machine vision: Visual cortex simulation.
Pulse and staircase edge models
Accuracy of Laplacian edge detection

Part II-efficient design and applications
Fast implementation of continuous wavelet transforms with integer scale
On the asymptotic convergence of B-spline wavelets to Gabor functions
The polynomial spline pyramid
A family of polynomial spline wavelet transforms
Fast Gabor-like windowed Fourier and continuous wavelet transform
A practical guide to the implementation of the wavelet transform
Ten good reasons for using spline wavelets

Multiscale curvature based shape representation using B-spline wavelets

Periodic orthogonal splines and wavelets
Characterization of compactly supported refinable splines
Window functions represented by B-spline functions
A theory for multiresolution signal decomposition: wavelet representation
Characterization of signals from multiscale edges
IEEE Trans.
Wavelets for a vision
Subband image coding using watershed and watercourse lines of the wavelet transform
Complete discrete 2D Gabor transform by neural networks for image analysis and compression
Biorthogonal bases of compactly supported wavelets
On compactly supported spline wavelets and a duality principle
Wavelets on a bounded interval
The discrete wavelet transform: wedding the
Multiscale corner detection by using wavelet transform
The curvature primal sketch


On the generating functions of totally positive sequences I
--TR

--CTR
Anant Madabhushi , Jayaram K. Udupa , Andre Souza, Generalized scale: theory, algorithms, and application to image inhomogeneity correction, Computer Vision and Image Understanding, v.101 n.2, p.100-121, February 2006
Yu-Ping Wang , Ruibin Qu, Fast Implementation of Scale-Space by Interpolatory Subdivision Scheme, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.21 n.9, p.933-939, September 1999
Yu-Ping Wang , Jie Chen , Qiang Wu , Kenneth R. Castleman, Fast frequency estimation by zero crossings of differential spline wavelet transform, EURASIP Journal on Applied Signal Processing, v.2005 n.1, p.1251-1260, 1 January 2005
Huamin Feng , Wei Fang , Sen Liu , Yong Fang, A new general framework for shot boundary detection and key-frame extraction, Proceedings of the 7th ACM SIGMM international workshop on Multimedia information retrieval, November 10-11, 2005, Hilton, Singapore
Luca Lombardi , Alfredo Petrosino, Distributed recursive learning for shape recognition through multiscale trees, Image and Vision Computing, v.25 n.2, p.240-247, February, 2007
b-spline;image modeling;scaling theorem;wavelet;scale-space;fingerprint theorem