Wavelets for Computer graphics: A Primer Part 2 EricJ. StolInitz Tony D Derose David H salesin University of Washington 1 Introduction It is often convenient to put the different scaling functions d(x)for a given level j together into a single row matrix, Wavelets are a mathematical tool (x)]. verall shape plus detai ess of whether the func a curve or a where nr is the dimension of y. We can do the same for the surface, wavelets provid levels of detail present y,(x)I In Part I of this primer we discu simple case of Haar wavelets in one and two dimensions, and showed how they can be If is the dimension of w. Because wd is the orthogonal com- used for image compression. In Part 2, we present the mathematical theory of multiresolution analysis, then develop bounded-interval nt of in pr, the dimensions of these spaces satisfy m spline wavelets and describe their use in multiresolution curve and surface editing The condition that the subs ted is equivalent to requir ing that the scaling function That is, forall=1, 2 there must exist a matrix 2 Multiresolution analysis The Haar wavelets we discussed in Part I are just one of many bases that can be used to treat functions in a hierarchical fashion In this each scaling function at levelj-l must l ection, we develop a mathematical framework known as multires combination of“ finer” scaling function olution analysis for studying wavelets [2, 11]. Our examples will w and J-I have dimensions m and m- 学 continue to focus on the Haar basis, but the more general mathe- tively, P is an mr x m-matrix(taller than it is wide) matical notation used here will come in handy for discussing other wavelet bases in later sections ince the wavelet space w- is by definition also a subspace ofpd, we can write the wavelets w-(x)as linear combinations of the scal Multiresolution analysis relies on many results from linear algebra. ing functions (x). This means there is an m x natrix of con- Some readers may wish to consult the appendix in Part 1 for a brief stants g satisfying v(x)=()Q As discussed in Part l, the starting point for multiresolution analysis is a nested set of vector spaces In the Haar basis, at a particular level aling functions and n=2 wavelets I cvcv Asj increases, the resolution of functions in y increases. The basis the four scaling functions in I functions for the space V are known as scaling funct The next step in multiresolution analysis is to definewavelet space 10 For each j, we define Wa as the orthogonal complementof /o in //+ This means that wo includes all the functions in pi that are orthog- al to all those in w under some chosen inner product. The func tions we choose as a basis for wo are called wavelets of ne d In the case of wavelets constructed on the real line. the columns of p/ are shifted 2.1 A matrix formulation for refinement nother, as are the columns of @. One column The rest of our discussion of multiresolution analysis will focus on characterizes each matrix, so P and g wavelets defined on a bounded domain, although we will also refer pletely determined by se (,p-1,P0,p1,…)and wavelets on the unbounded real line wherever appropriate. In the (.. 9-1, 9o, 91, .. ) which also do not depend on j. Equa- bounded case, each space W has a finite basis, allowing us to use ma- pressions of the form trix notation in much of what follows, as did lounsbery et al.[101 and Quak and Weyrich [13] ox)=∑p2x-0 ter graphics: A primer, part 2, IEEE Computer Graphics and Applica- ions,15(4):75-85,July199
Wavelets for Computer Graphics: A Primer Part 2y Eric J. Stollnitz Tony D. DeRose David H. Salesin University of Washington 1 Introduction Wavelets are a mathematical tool for hierarchically decomposing functions. They allow a function to be described in terms of a coarse overall shape, plus details that range from broad to narrow. Regardless of whether the function of interest is an image, a curve, or a surface, wavelets provide an elegant technique for representing the levels of detail present. In Part 1 of this primer we discussed the simple case of Haar wavelets in one and two dimensions, and showed how they can be used for image compression. In Part 2, we present the mathematical theory of multiresolution analysis, then develop bounded-interval spline wavelets and describe their use in multiresolution curve and surface editing. 2 Multiresolution analysis The Haar wavelets we discussed in Part 1 are just one of many bases that can be used to treat functions in a hierarchical fashion. In this section, we develop a mathematical framework known asmultiresolution analysis for studying wavelets [2, 11]. Our examples will continue to focus on the Haar basis, but the more general mathematical notation used here will come in handy for discussing other wavelet bases in later sections. Multiresolution analysis relies on many results from linear algebra. Some readers may wish to consult the appendix in Part 1 for a brief review. As discussed in Part 1, the starting point for multiresolution analysis is a nested set of vector spaces V0 V1 V2 As j increases, the resolution of functions inVj increases. The basis functions for the space Vj are known as scaling functions. The next step in multiresolution analysis is to definewavelet spaces. For each j, we define Wj as the orthogonal complement of Vj in Vj+1. This means that Wj includes all the functions inVj+1 that are orthogonal to all those in Vj under some chosen inner product. The functions we choose as a basis for Wj are called wavelets. 2.1 A matrix formulation for refinement The rest of our discussion of multiresolution analysis will focus on wavelets defined on a bounded domain, although we will also refer to wavelets on the unbounded real line wherever appropriate. In the bounded case, each spaceVj has a finite basis, allowing us to use matrix notation in much of what follows, as did Lounsbery et al. [10] and Quak and Weyrich [13]. y Eric J. Stollnitz, Tony D. DeRose, and David H. Salesin. Wavelets for computer graphics: A primer, part 2. IEEE Computer Graphics and Applications, 15(4):75–85, July 1995. It is often convenient to put the different scaling functions j i (x) for a given level j together into a single row matrix, j (x) := [ j 0(x) j mj1 (x)], where mj is the dimension of Vj . We can do the same for the wavelets: j (x) := [ j 0(x) j nj1(x)], where nj is the dimension of Wj . Because Wj is the orthogonal complement of Vj in Vj+1, the dimensions of these spaces satisfy mj+1 = mj + nj . The condition that the subspacesVj be nested is equivalent to requiring that the scaling functions berefinable. That is, for all j = 1, 2, : : : there must exist a matrix of constants Pj such that j1 (x) = j (x) Pj . (1) In other words, each scaling function at levelj 1 must be expressible as a linear combination of “finer” scaling functions at level j. Note that since Vj and Vj1 have dimensions mj and mj1 , respectively, Pj is an mj mj1 matrix (taller than it is wide). Since the wavelet space Wj1 is by definition also a subspace ofVj , we can write the wavelets j1 (x) as linear combinations of the scaling functions j (x). This means there is an mj nj1 matrix of constants Qj satisfying j1 (x) = j (x) Qj . (2) Example: In the Haar basis, at a particular level j there are mj = 2j scaling functions and nj = 2j wavelets. Thus, there must be refinement matrices describing how the two scaling functions in V1 and the two wavelets in W1 can be made from the four scaling functions inV2 : P2 = 2 6 4 1 0 1 0 0 1 0 1 3 7 5 and Q2 = 2 6 4 1 0 1 0 0 1 0 1 3 7 5 Remark: In the case of wavelets constructed on the unbounded real line, the columns of Pj are shifted versions of one another, as are the columns of Qj . One column therefore characterizes each matrix, so Pj and Qj are completely determined by sequences (: : : , p1, p0, p1, : : :) and (: : : , q1, q0, q1, : : :), which also do not depend on j. Equations (1) and (2) therefore often appear in the literature as expressions of the form (x) = X i pi (2x i) (x) = X i qi (2x i).
These equations are referred to as nwo-scale relations for notation can also be used for the decomposition process outlined in Section 2.1 of Part I Note that equations(1)and(2)can be expressed as a single equation Consider a function in some approximation space / Let's assume sing block-matrix notation we have the coefficients of this function in terms of some scaling function basis. We can write these coefficients as a column matrix [型-1|吵-1]=型[P|]」 The coefficients c could. for ex- ample, be thought of as pixel colors, or alternatively, as the xor Example: Substituting the matrices from the previous ex y-coordinates of a curve's control points in R- ple into Equation(3) along with the appropriate basis Suppose we wish to create a low-resolution version C-l of C wit a smaller number of coefficients mr. The standard approach for the m-I values of g-l is to use some form of linear filter 1010 down-sampling on then entries of C. This process can be t=[的或的 10-10 expressed as a matrix equation where A is an m-I x n matrix of constants(wider than it is tall) ions and their refinement matrices P, the wavelet matrices o are Since C- contains fewer entries than C, this filtering process somewhat constrained( though not completely determined ). In fact, since all functions in /-(x)must be orthogonal to all functions possible to capture the lost detail as another column matrix D-I in (x), we know(h lv)=0 for all k and e comp To deal with all these inner products simultaneously, lets define some new notation for a matrix of inner products. We will denote where B is an -xm matrix of constants related to 4. The pair of K4-4)] the matrix whose(k, e) entry is(d-v- Armed with this notation, we can rewrite the orthogonality condi- the coefficients C into a low-resolution version C-I and detail D tion on the wavelets as is called analysis or decomposition y1y-)]=0 If A and B are chosen appropriately, then the original coefficients C can be recovered from C- and D- by using the matrices P Substituting Equation(2)into Equation(4)yield and o from the previous section: ]=0 C=PC-I+OD A matrix equation with a right-hand side of zero like this one is Recovering C from C-and D- is called synthesis or reconstruc- known as a homogeneous system of equations. The set of all pos- tion In this context, P and g are called synthesis filters sible solutions is called the null space of [(-)], and the columns ofg must form a basis for this space. There are a multitude In the unnormalized haar basis. the matrices a of bases for the null space of a matrix, implying that there are different wavelet bases for a given wavelet space w. Ordinarily, we 1100 uely determine the o matrices by imposing further constraints in addition to the orthogonality requirement given above. For exam- 0011 ple, the Haar wavelet matrices can be found by requiring the least number of consecutive nonzero entries in each column The literature on wavelets includes various terminologies for or Ithors refer to a collection of functions that orthogonal to scaling functions but not to each other aspre-wavelets, reserving the term "wavelets" for functions that are orthogonal to ch other as well. Another common approach is to differentiate be- Remark: Once again, the matrices for wavelets constructed tween an orthogonal wavelet basis, in which all functions are m on the unbounded real line have a simple structure: The i tually orthogonal, and a semi-orthogonal wavelet basis, in which A are shifted versions of each other. as are the rows of B the wavelets are orthogonal to the scaling functions but not to each The analysis Equations(6)and(7)often appear in the litera- other. The Haar basis is an example of an orthogonal wavelet basis, while the spline wavelets we will describe in Section 3 are example Finally, it is sometimes desirable to define wavelets that are not qu orthogonal to scaling functions in order to have wavelets with small supports. This last alternative might be termed a non-orthogonal wavelet basis, and we will mention an example when we describe where the sequences multiresolution surfaces in Section 4.3 (.. b-1, bo, b1,.. )are the entries espectively. Similarly, Equation(8) 2.2 The filter bank The previous section showed how scaling functions and wavelets could be related by matrices. In this section, we show how matrix
These equations are referred to as two-scale relations for scaling functions and wavelets, respectively. Note that equations (1) and (2) can be expressed as a single equation using block-matrix notation: j1 j1 = j Pj Qj . (3) Example: Substituting the matrices from the previous example into Equation (3) along with the appropriate basis functions gives [ 1 0 1 1 1 0 1 1 ] = [ 2 0 2 1 2 2 2 3] 2 6 4 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 3 7 5 It is important to realize that once we have chosen scaling functions and their refinement matrices Pj , the wavelet matrices Qj are somewhat constrained (though not completely determined). In fact, since all functions in j1 (x) must be orthogonal to all functions in j1 (x), we know h j1 k j j1 ` i = 0 for all k and `. To deal with all these inner products simultaneously, let’s define some new notation for a matrix of inner products. We will denote by [hj1 j j1i] the matrix whose (k, `) entry is h j1 k j j1 ` i. Armed with this notation, we can rewrite the orthogonality condition on the wavelets as [hj1 j j1i] = 0. (4) Substituting Equation (2) into Equation (4) yields [hj1 j j i] Qj = 0. (5) A matrix equation with a right-hand side of zero like this one is known as a homogeneous system of equations. The set of all possible solutions is called the null space of [hj1 j j i], and the columns ofQj must form a basis for this space. There are a multitude of bases for the null space of a matrix, implying that there are many different wavelet bases for a given wavelet spaceWj . Ordinarily, we uniquely determine the Qj matrices by imposing further constraints in addition to the orthogonality requirement given above. For example, the Haar wavelet matrices can be found by requiring the least number of consecutive nonzero entries in each column. The literature on wavelets includes various terminologies for orthogonality. Some authors refer to a collection of functions that are orthogonal to scaling functions but not to each other aspre-wavelets, reserving the term “wavelets” for functions that are orthogonal to each other as well. Another common approach is to differentiate between an orthogonal wavelet basis, in which all functions are mutually orthogonal, and a semi-orthogonal wavelet basis, in which the wavelets are orthogonal to the scaling functions but not to each other. The Haar basis is an example of an orthogonal wavelet basis, while the spline wavelets we will describe in Section 3 are examples of semi-orthogonal wavelet bases. Finally, it is sometimes desirable to define wavelets that are not quite orthogonal to scaling functions in order to have wavelets with small supports. This last alternative might be termed a non-orthogonal wavelet basis, and we will mention an example when we describe multiresolution surfaces in Section 4.3. 2.2 The filter bank The previous section showed how scaling functions and wavelets could be related by matrices. In this section, we show how matrix notation can also be used for the decomposition process outlined in Section 2.1 of Part 1. Consider a function in some approximation spaceVj . Let’s assume we have the coefficients of this function in terms of some scaling function basis. We can write these coefficients as a column matrix of values Cj = [c j 0 c j mj1] T . The coefficients c j i could, for example, be thought of as pixel colors, or alternatively, as the x- or y-coordinates of a curve’s control points in IR2 . Suppose we wish to create a low-resolution versionCj1 of Cj with a smaller number of coefficients mj1 . The standard approach for creating the mj1 values of Cj1 is to use some form of linear filtering and down-sampling on themj entries of Cj . This process can be expressed as a matrix equation Cj1 = Aj Cj (6) where Aj is an mj1 mj matrix of constants (wider than it is tall). Since Cj1 contains fewer entries than Cj , this filtering process clearly loses some amount of detail. For many choices of Aj , it is possible to capture the lost detail as another column matrix Dj1 , computed by Dj1 = Bj Cj (7) where Bj is an nj1mj matrix of constants related toAj . The pair of matrices Aj and Bj are called analysis filters. The process of splitting the coefficientsCj into a low-resolution versionCj1 and detail Dj1 is called analysis or decomposition. If Aj and Bj are chosen appropriately, then the original coefficients Cj can be recovered from Cj1 and Dj1 by using the matrices Pj and Qj from the previous section: Cj = Pj Cj1 + Qj Dj1 . (8) Recovering Cj from Cj1 and Dj1 is called synthesis orreconstruction. In this context, Pj and Qj are called synthesis filters. Example: In the unnormalized Haar basis, the matrices A2 and B2 are given by: A2 = 1 2 1 1 0 0 0 0 1 1 B2 = 1 2 1 1 0 0 0 0 1 1 These matrices represent the averaging and differencing operations described in Section 2.1 of Part 1. Remark: Once again, the matrices for wavelets constructed on the unbounded real line have a simple structure: The rows of Aj are shifted versions of each other, as are the rows ofBj . The analysis Equations (6) and (7) often appear in the literature as c j1 k = X ` a`2k c j` dj1 k = X ` b`2k c j` where the sequences (: : : , a1, a0, a1, : : :) and (: : : , b1, b0, b1, : : :) are the entries in a row of Aj and Bj , respectively. Similarly, Equation (8) for reconstruction often appears as c j k = X ` pk2` c j1 ` + qk2` dj1 ` . 2
must be made. In the case of the spline wavelets presented in the next section, the wavelets are constructed to have minimal support, but they are not orthogonal to one another(except for the piecewise D constant case). Wavelets that are both locally supported and mutu- ally orthogonal (other than Haar wavelets) were thought to be im- Figure 1 The filter bank. possible until Daubechies'ground-breaking work showing that cer- tain families of spaces W actually do admit mutually orthogonal Note that the procedure for splitting C into a low-resolution wavelets of small support [5] low-resolution version C-. Thus, the original coefficients can be 3 Splin expressed as a hierarchy of lower-resolution versions C,..., C- and details D, .. D-, as shown in Figure 1. This recursive pr Until now, the only specific wavelet basis we have considered is the cess is known as a filter bank. Haar basis. Haar basis functions have a number of advantages. in- cluding transform of the original coefficients, known as wavelet transform Note that the total size of the transform Co. Do. DI..D-I is the orthogonality, same as that of the original version C, so no extra storage is re- quired. (However, the wavelet coefficients may require more bits very small supports, nonoverlapping scaling functions(at a given level), and In general, the analysis filters A and B are not necessarily trans- nonoverlapping wavelets(at a given level), posed multiples of the synthesis filters, as was the case for the Haar basis. Rather, A and B are formed by the matrices satisfying the which make them useful in many applications. However, despite ese advantages, the Haar basis is a poor choice for applications such editing 8 and animation [9 because of its lack of B There are a variety of ways to construct wavelets with k continu ous derivatives One such class of wavelets can be constructed from Note that y are both square matrices. Thus, piecewise-polynomial splines. These spline wavelets have been de veloped to a large extent by Chui and colleagues 3, 4]. The Haar ba- combining equations(3)and(9)gives sis is in fact the simplest instance of spline wavelets, resulting when [P g (10) In the following, we briefly sketch the ideas behind the construc- B on of endpoint-interpolating B-spline wavelets. Finkelstein and Although we have not yet gotten specific about how to choose ma. Salesin [8] developed a collection of wavelets for the cubic case, trices A, B, P, and @, it should be clear from Equation (10)that Although the derivations for arbitrary degree are too involved to the two matrices in that equation Haar), linear, quadratic, and cubic cases in Appendix A. The next 2.3 Designing a multiresolution analysis three sections parallel the three steps described in Section 2.3 for The multiresolution analysis framework presente very designing a multiresolution analysis olution analysis specifically suited to a particule . s general. In practice you often have the freedom to litres- 3.1 B-spline scaling Our first step is to define the scaling functions for a nested set of 1. Select the scaling functions (x)for each=0, function spaces. We'll start with the general definition of B-splines, This choice determines the nested approximation spaces v,the then specify how to make uniformly spaced, endpoint-interpolating synthesis filters P and the smoothness-that is, the number of continuous derivatives--of the analy B-splines from these. (More detailed derivations of these and other splines appear in a number of standard texts [1, 7]) k. with k>d an This choice determines the L- norm and the orthogonal comple- non-decreasing values xo, .. xhd+i called knots, the nonuniform B ment spaces Wo. Although the standard inner product is the com spline basis functions of degree d are defined recursively as follows mon choice, in general the inner product should be chosen to cap- For i=0..k, and for r=1..d let hat is meaningful in the context of the Icat M(x):= x< x< xa+ 3. Select a set of wavelets w(x)that span W for each=0, This choice determines the synthesis filters. Together, the syn thesis filters P and g determine the analysis filters A and B by N(x) =-N-(x) ri+r+1-x xit+I --xH+ onal basis for W and to have small support(the support of a f tion(r)is the set of points x wheref(x)#0). However, orthogonal- ity often comes at the expense of increased supports, so a tradeoff The endpoint-interpolating B-splines of degree d on [0, 1] result
Cj HH Bj HHjH Dj1 A -j Cj1 HH Bj1HHHj Dj2 A - j1 Cj2 C1 HH B1HHHj D0 A -1 C0 Figure 1 The filter bank. Note that the procedure for splitting Cj into a low-resolution part Cj1 and a detail part Dj1 can be applied recursively to the low-resolution version Cj1 . Thus, the original coefficients can be expressed as a hierarchy of lower-resolution versionsC0 , : : : ,Cj1 and details D0 , : : : , Dj1 , as shown in Figure 1. This recursive process is known as a filter bank. Since the original coefficients Cj can be recovered from the sequence C0 , D0 , D1 , : : :, Dj1 , we can think of this sequence as a transform of the original coefficients, known as awavelet transform. Note that the total size of the transform C0 , D0 , D1 , : : :, Dj1 is the same as that of the original version Cj , so no extra storage is required. (However, the wavelet coefficients may require more bits to retain the accuracy of the original values.) In general, the analysis filters Aj and Bj are not necessarily transposed multiples of the synthesis filters, as was the case for the Haar basis. Rather, Aj and Bj are formed by the matrices satisfying the relation j1 j1 Aj Bj = j . (9) Note that Pj Qj and Aj Bj are both square matrices. Thus, combining Equations (3) and (9) gives Aj Bj = Pj Qj 1 (10) Although we have not yet gotten specific about how to choose matrices Aj , Bj , Pj , and Qj , it should be clear from Equation (10) that the two matrices in that equation must at least be invertible. 2.3 Designing a multiresolution analysis The multiresolution analysis framework presented above is very general. In practice you often have the freedom to design a multiresolution analysis specifically suited to a particular application. The steps involved are 1. Select the scaling functions j (x) for each j = 0, 1, : : : . This choice determines the nested approximation spacesVj , the synthesis filters Pj , and the smoothness—that is, the number of continuous derivatives—of the analysis. 2. Select an inner product defined on the functions in V0 , V1 , : : : . This choice determines the L2 norm and the orthogonal complement spaces Wj . Although the standard inner product is the common choice, in general the inner product should be chosen to capture a measure of error that is meaningful in the context of the application. 3. Select a set of wavelets j (x) that span Wj for each j = 0, 1, : : : . This choice determines the synthesis filtersQj . Together, the synthesis filters Pj and Qj determine the analysis filters Aj and Bj by Equation (10). It is generally desirable to construct the wavelets to form an orthogonal basis for Wj and to have small support (the support of a function f(x) is the set of points x where f(x) 6= 0). However, orthogonality often comes at the expense of increased supports, so a tradeoff must be made. In the case of the spline wavelets presented in the next section, the wavelets are constructed to have minimal support, but they are not orthogonal to one another (except for the piecewiseconstant case). Wavelets that are both locally supported and mutually orthogonal (other than Haar wavelets) were thought to be impossible until Daubechies’ ground-breaking work showing that certain families of spaces Vj actually do admit mutually orthogonal wavelets of small support [5]. 3 Spline wavelets Until now, the only specific wavelet basis we have considered is the Haar basis. Haar basis functions have a number of advantages, including simplicity, orthogonality, very small supports, nonoverlapping scaling functions (at a given level), and nonoverlapping wavelets (at a given level), which make them useful in many applications. However, despite these advantages, the Haar basis is a poor choice for applications such as curve editing [8] and animation [9] because of its lack of continuity. There are a variety of ways to construct wavelets with k continuous derivatives. One such class of wavelets can be constructed from piecewise-polynomial splines. These spline wavelets have been developed to a large extent by Chui and colleagues [3, 4]. The Haar basis is in fact the simplest instance of spline wavelets, resulting when the polynomial degree is set to zero. In the following, we briefly sketch the ideas behind the construction of endpoint-interpolating B-spline wavelets. Finkelstein and Salesin [8] developed a collection of wavelets for the cubic case, and Chui and Quak [4] presented constructions for arbitrary degree. Although the derivations for arbitrary degree are too involved to present here, we give the synthesis filters for the piecewise-constant (Haar), linear, quadratic, and cubic cases in Appendix A. The next three sections parallel the three steps described in Section 2.3 for designing a multiresolution analysis. 3.1 B-spline scaling functions Our first step is to define the scaling functions for a nested set of function spaces. We’ll start with the general definition of B-splines, then specify how to make uniformly spaced, endpoint-interpolating B-splines from these. (More detailed derivations of these and other splines appear in a number of standard texts [1, 7].) Given positive integers d and k, with k d, and a collection of non-decreasing valuesx0, : : : , xk+d+1 called knots, the nonuniform Bspline basis functions of degree d are defined recursively as follows. For i = 0, : : : , k, and for r = 1, : : : , d, let N0 i (x) := 1 if xi x < xi+1 0 otherwise Nr i(x) := x xi xi+r xi Nr1 i (x) + xi+r+1 x xi+r+1 xi+1 Nr1 i+1 (x). (Note: The fractions in these equations are taken to be 0 when their denominators are 0.) The endpoint-interpolating B-splines of degree d on [0, 1] result when the first and last d + 1 knots are set to 0 and 1, respectively. In 3
A The second step of designing a multiresolution analysis is the choice N of an inner product. We'll simply use the standard inner product M N fx)g(x)dx N 3.3 B-spline wavelets To complete our development of a multiresolution analysis based on degree 2 B-spl need to find basis functions for th wa th the wavelets are determined by matrices O satisfying Equation(5) Figure 2 B-spline scaling functions for v(d) with degree d which we repeat here for convenience: l型)Q=0 (11) this case, the functions No(x).., NR(x) form a basis for the space of piecewise-polynomials of degree d with d-I continuous deriva- Since this is a homogeneous system of linear equations, there tives and breakpoints at the interior knots xl <xd2<--<xk not a unique solution. We must therefore impose additional conc. p To make uniformly spaced B-splines, we choose =2+d-1 tions. To get wavelets with small supports, for example, we require and xd+l, . .. xk to produce 2 equally spaced interior intervals. This each column of o to have a minimal number of consecutive non- zeros. This constraint imposes a banded structure ono similar to onstruction gives 2+ d B-spline basis functions for a particular that of Pl. For each column g of @, Equation(11)leads to a small egree d and level j. We will use these functions as the endpoint interpolating B-spline scaling functions. Figure 2 shows examples homogeneous system that we solve for the non-zero entries in q f these functions at level= I(two interior intervals) for various The matrices that result and the corresponding B-spline wavelets are degrees d. Note that the basis functions defined here are not normal- shown in Appendix A ized in the l norm Finkelstein and Salesin 8] took this approach to construct cubic B- If v(d) denotes the space spanned by the B-spline scaling functions spline wavelets. Chui and Quak [4] derived slightly different spline of degree d with 2 ul the spaces/(),v(d), als, it is not difficult to show that Note that both approaches result in semi-orthogonal wavelet bases sted as required by multiresolution The wavelets are orthogonal to scaling functions at the same level, but not to each other, except in the piecewise-constant case The rich theory of B-splines can be used to develop expressions for the entries of the refinement matrix P(see Chui and Quak [4]or 3.4 B-spline filter bank Quak and Weyrich[13] for details ) The columns of P are sparse, reflecting the fact that the B-spline basis functions are locally sup- At this point, we have completed the steps in designing a multires- ported. The first and last d columns of p are relatively compli olution analysis. However, to use spline wavelets, we must im ated, but the remaining (interior) columns are shifted versions of plement a filter bank procedure incorporating the analysis filters.A column d+ 1. Moreover. the entries of these interior columns are and B. These matrices allow us to determine C-I and D-I from Cy up to a common factor of 1/2, given by binomial coefficients. using matrix multiplication as in Equations(6)and cussed earlier in Section 2, the analysis filters are uniquely deter- Example: In the case of cubic splines(d= 3), the matrix Pr mined by the inverse relation in Equation(10) forj≥3 has the form B」=[Py 44 However, as Quak and Weyrich [13] point out, when implementing the filter bank procedure for spline wavelets, it is generally not a good idea to form the filters A and B explicitly. Although P and g are sparse, having only O(d)entries per column, t and B are in gen eral dense, so matrix-vector multiplication would require quadratic nstead of linear time Fortunately, there is a better approach. The idea is to compute C- and D- from C by solving the sparse linear system In order to solve this system for where blank entries are taken to be zero. and the dots indi- D cate that the previous column ated, shifted down by trix [p gl into a banded matrix simply by interspersing the columns of P and g. The resulting banded system can then be
N0 0 N0 1 degree 0 N1 0 N1 1 N1 2 degree 1 N2 0 N2 1 N2 2 N2 3 degree 2 N3 0 N3 1 N3 2 N3 3 N3 4 degree 3 Figure 2 B-spline scaling functions for V1(d) with degree d = 0, 1, 2, and 3. this case, the functions Nd 0 (x), : : : , Nd k (x) form a basis for the space of piecewise-polynomials of degree d with d1 continuous derivatives and breakpoints at the interior knots xd+1 < xd+2 < < xk. To make uniformly spaced B-splines, we choose k = 2j + d 1 and xd+1, : : : , xk to produce 2j equally spaced interior intervals. This construction gives 2j + d B-spline basis functions for a particular degree d and level j. We will use these functions as the endpointinterpolating B-spline scaling functions. Figure 2 shows examples of these functions at level j = 1 (two interior intervals) for various degrees d. Note that the basis functions defined here are not normalized in the L2 norm. If Vj (d) denotes the space spanned by the B-spline scaling functions of degree d with 2j uniform intervals, it is not difficult to show that the spaces V0 (d), V1 (d), : : : are nested as required by multiresolution analysis. The rich theory of B-splines can be used to develop expressions for the entries of the refinement matrix Pj (see Chui and Quak [4] or Quak and Weyrich [13] for details). The columns ofPj are sparse, reflecting the fact that the B-spline basis functions are locally supported. The first and last d columns of Pj are relatively complicated, but the remaining (interior) columns are shifted versions of column d + 1. Moreover, the entries of these interior columns are, up to a common factor of 1=2d , given by binomial coefficients. Example: In the case of cubic splines (d = 3), the matrix Pj for j 3 has the form Pj = 1 8 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 8 4 4 6 2 3 2 11 2 1 4 4 1 6 1 4 4 1 6 4 1 1 4 6 1 4 4 1 11 2 3 2 2 6 4 4 8 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 where blank entries are taken to be zero, and the dots indicate that the previous column is repeated, shifted down by two rows each time. 3.2 Inner product The second step of designing a multiresolution analysis is the choice of an inner product. We’ll simply use the standard inner product here: hf j gi := Z 1 0 f(x) g(x) dx. 3.3 B-spline wavelets To complete our development of a multiresolution analysis based on B-splines, we need to find basis functions for the spacesWj that are orthogonal complements to the spacesVj . As shown in Section 2.1, the wavelets are determined by matricesQj satisfying Equation (5), which we repeat here for convenience: [hj1 j j i] Qj = 0. (11) Since this is a homogeneous system of linear equations, there is not a unique solution. We must therefore impose additional conditions. To get wavelets with small supports, for example, we require each column of Qj to have a minimal number of consecutive nonzeros. This constraint imposes a banded structure on Qj similar to that of Pj . For each column q of Qj , Equation (11) leads to a small homogeneous system that we solve for the non-zero entries in q. The matrices that result and the corresponding B-spline wavelets are shown in Appendix A Finkelstein and Salesin [8] took this approach to construct cubic Bspline wavelets. Chui and Quak [4] derived slightly different spline wavelets using derivative and interpolation properties of B-splines. Note that both approaches result in semi-orthogonal wavelet bases: The wavelets are orthogonal to scaling functions at the same level, but not to each other, except in the piecewise-constant case. 3.4 B-spline filter bank At this point, we have completed the steps in designing a multiresolution analysis. However, to use spline wavelets, we must implement a filter bank procedure incorporating the analysis filtersAj and Bj . These matrices allow us to determineCj1 and Dj1 from Cj using matrix multiplication as in Equations (6) and (7). As discussed earlier in Section 2, the analysis filters are uniquely determined by the inverse relation in Equation (10): Aj Bj = Pj Qj 1 However, as Quak and Weyrich [13] point out, when implementing the filter bank procedure for spline wavelets, it is generally not a good idea to form the filtersAj and Bj explicitly. Although Pj and Qj are sparse, having onlyO(d) entries per column, Aj and Bj are in general dense, so matrix–vector multiplication would require quadratic instead of linear time. Fortunately, there is a better approach. The idea is to computeCj1 and Dj1 from Cj by solving the sparse linear system Pj Qj Cj1 Dj1 = Cj . In order to solve this system for Cj1 Dj1 , we first make the matrix Pj Qj into a banded matrix simply by interspersing the columns of Pj and Qj . The resulting banded system can then be 4
Figure 3 Changing a curves overall sweep without affecting its character. Given the original curve(a), the system extracts the over- all sweep(b). If the user modifies the sweep(c), the system can I pply the detail (d) Figure 5 Changing the character of a curve without affecting its sweep through synthesis +AC C+pp-1. P+I AC Note that editing the sweep of the curve at lower levels of The middle of the dark curve is pulled, using editing at var- els of smoothing. A change in a control point incl ing affects larger portions of the high-resolution curve(o lowest level, when j=0, the entire curve is affected. At the d effect, while a change in a control point row effect level, when j=J, only the narrow portion influenced by one control point is affected. The kind of flexibility that this multireso- solved in linear time using LU decomposition [12]. Thus we can lution editing allows is suggested in Figures 3 and 4 compute the entire filter bank operation without ever forming and using A or B explicitly 4. 2 Editing the character of the curve Multiresolution curves also naturally support changes in the char 4 Application II: Multiresolution curves and surfaces acter of a curve. without affecting its overall Let c be control points of a curve, and letC, Do,..., D-I denote its wavelet transform. Editing the character of the curve is simply a matter of re- We presented two applications of wavelets in Part I: compression of placing the existing set of detail coefficients,.,D-Iwithsome ages. Our third application of wavelets in computer graphics is new set and reconstructing. To avoid coo curve design and editing, as described in detail by Finkelstein and system artifacts. all detail coefficients are expressed in terms of the Salesin (8J. Their multiresolution curves are built from a wavelet curve s local tangent and normal, rather than thex and y directions cussed in the previous section s cubic B-splines, which we dis- Figure 5 demonstrates how the character of curves in an illustration can be modified with various detail styles. ( The interactive illus- Multiresolution curves conveniently support a variety of opera- tration system used to create this figure was described by Salisbury eta.[l4].) changing a curve's overall"sweep"while maintaining its fine de 4.3 Multiresolution surfaces tails, or"character"(Figures 3 and 4); changing a curve's"character without affecting its overall Multiresolution editing can be extended to surfaces by using ten- sor products of B-spline scaling functions and wavelets. Either the standard construction or the nonstandard construction described g a curve at any continuous level of detail, allowing an ar- in Part 1 for Haar basis functions can be used to form a two- ary portion of the curve to be affected through direct manipu- dimensional basis from a one-dimensional B-spline basis. We can then edit surfaces using the same operations described for curves. For example, Figure 6 shows a bicubic tensor-product B-spline sur- smoothing at continuous levels to remove undesirable features face after altering its at different levels of detail from a curve approximating or"fitting"a curve within a guaranteed maximum We can further generalize multiresolution analysis to surfaces of error tolerance, for scan conversion and other applications arbitrary topology by defining wavelets based on subdivision sur- faces, as described by Lounsbery et al. [10]. Their nonorthogonal Here we'll describe briefly just the first two of these operations, wavelet basis, in combination with the work of Eck et al. 6], al which fall out quite naturally from the multiresolution representa lows any polyhedral object to be decomposed into scaling function and wavelet coefficients. Then a compression scheme similar to the one presented for images in Section 3. 3 of Part 1 can be used to dis- play the object at various levels of detail simply by leaving out sma 4.1 Editing the sweep of the curve wavelet coefficients during reconstruction. An example of this tech- nique is shown in Figure 7. Editing the sweep of a curve at an level of the wavelet transform is simple. Let C be the points of the original curve o, let C be a low-resolution version of C, and let c be 5 Conclusion an edited version of C, given by C-C+ AC. The edited ver- Our primer has only touched on a few of the many uses for wavelets sion of the highest resolution curve C in computer graphics. We hope this introduction to the topic has ex-
(a) (b) (c) (d) Figure 3 Changing a curve’s overall sweep without affecting its character. Given the original curve (a), the system extracts the overall sweep (b). If the user modifies the sweep (c), the system can reapply the detail (d). 1 2 3 4 Figure 4 The middle of the dark curve is pulled, using editing at various levels of smoothing j. A change in a control point inC1 has a very broad effect, while a change in a control point inC4 has a narrow effect. solved in linear time using LU decomposition [12]. Thus we can compute the entire filter bank operation without ever forming and using Aj or Bj explicitly. 4 Application III: Multiresolution curves and surfaces We presented two applications of wavelets in Part 1: compression of one-dimensional signals and compression of two-dimensional images. Our third application of wavelets in computer graphics is curve design and editing, as described in detail by Finkelstein and Salesin [8]. Their multiresolution curves are built from a wavelet basis for endpoint-interpolating cubic B-splines, which we discussed in the previous section. Multiresolution curves conveniently support a variety of operations: changing a curve’s overall “sweep” while maintaining its fine details, or “character” (Figures 3 and 4); changing a curve’s “character” without affecting its overall “sweep” (Figure 5); editing a curve at any continuous level of detail, allowing an arbitrary portion of the curve to be affected through direct manipulation; smoothing at continuous levels to remove undesirable features from a curve; approximating or “fitting” a curve within a guaranteed maximum error tolerance, for scan conversion and other applications. Here we’ll describe briefly just the first two of these operations, which fall out quite naturally from the multiresolution representation. 4.1 Editing the sweep of the curve Editing the sweep of a curve at an integer level of the wavelet transform is simple. Let CJ be the control points of the original curve f J (t), let Cj be a low-resolution version of CJ , and let Cˆj be an edited version of Cj , given by Cˆj = Cj + Cj . The edited version of the highest resolution curveCˆ J = CJ +CJ can be computed Figure 5 Changing the character of a curve without affecting its sweep. through synthesis: Cˆ J = CJ + CJ = CJ + PJ PJ1 Pj+1 Cj . Note that editing the sweep of the curve at lower levels of smoothing j affects larger portions of the high-resolution curvef J (t). At the lowest level, when j = 0, the entire curve is affected. At the highest level, when j = J, only the narrow portion influenced by one original control point is affected. The kind of flexibility that this multiresolution editing allows is suggested in Figures 3 and 4. 4.2 Editing the character of the curve Multiresolution curves also naturally support changes in the character of a curve, without affecting its overall sweep. Let CJ be the control points of a curve, and letC0 , D0 , : : :, DJ1 denote its wavelet transform. Editing the character of the curve is simply a matter of replacing the existing set of detail coefficientsDj , : : : , DJ1 with some new set Dˆj , : : : , Dˆ J1 , and reconstructing. To avoid coordinatesystem artifacts, all detail coefficients are expressed in terms of the curve’s local tangent and normal, rather than the x and y directions. Figure 5 demonstrates how the character of curves in an illustration can be modified with various detail styles. (The interactive illustration system used to create this figure was described by Salisbury et al. [14].) 4.3 Multiresolution surfaces Multiresolution editing can be extended to surfaces by using tensor products of B-spline scaling functions and wavelets. Either the standard construction or the nonstandard construction described in Part 1 for Haar basis functions can be used to form a twodimensional basis from a one-dimensional B-spline basis. We can then edit surfaces using the same operations described for curves. For example, Figure 6 shows a bicubic tensor-product B-spline surface after altering its sweep at different levels of detail. We can further generalize multiresolution analysis to surfaces of arbitrary topology by defining wavelets based on subdivision surfaces, as described by Lounsbery et al. [10]. Their nonorthogonal wavelet basis, in combination with the work of Eck et al. [6], allows any polyhedral object to be decomposed into scaling function and wavelet coefficients. Then a compression scheme similar to the one presented for images in Section 3.3 of Part 1 can be used to display the object at various levels of detail simply by leaving out small wavelet coefficients during reconstruction. An example of this technique is shown in Figure 7. 5 Conclusion Our primer has only touched on a few of the many uses for wavelets in computer graphics. We hope this introduction to the topic has ex- 5