Warren J., Weimer H. Subdivision Methods for Geometric Design. A Constructive Approach

Подождите немного. Документ загружается.

94 CHAPTER 4 A Differential Approach to Uniform Subdivision

THEOREM

4.1

Let p[x] be a function whose integral I [x] p[x] is a continuous function. Then

D [x] I [x] p[x] == p[x].

Proof By definition, I [x] p[x] is the integral



∞

p[x − t] dt. By a simple change of

variables, the integral can be rewritten as



−∞

p[t] dt. Taking the derivative

of this integral with respect to

x corresponds to taking the limit of the

expression







−∞

p[t] dt −



x−

−∞

p[t] dt



 → 0. Combining the two integrals yields a new expression of the form





x−

p[t] dt. Because the value of this integral is a continuous function (by

assumption), the limit as

 →∞is exactly p[x ].

(A similar theorem holds in the case when p[x] is the Dirac delta function and the

integral

I[x] p[x] is the step function c

[x].) Because D[x] is the left inverse of I[x],

multiplying both sides of equation 4.1 by

D [x]

yields a differential equation for

the truncated powers of the form

D [x]

[x] == δ[x]. (4.2)

Thus, the truncated powers

[x]

are the solutions to the differential equation

(m)

[x] == δ[x]. Solutions to differential equations whose right-hand side is the Dirac

delta

δ[x] are often referred to as a Green’s function, associated with the differential

operator. Note that the Green’s functions for this differential operator are unique

only up to addition of any polynomial of degree

m − 1.

Having derived a differential equation for the truncated powers, we next con-

struct the differential equation for a basis function

[x] using Theorem 2.4. In

particular, if

[x] is the generating function (1 − x)

, the B-spline basis function

[x] has the form



i =0

[[ i ]] c

[x − i], where d

[[ i ]] is the coefficient of x

in d

[x].

Substituting this relation into equation 4.2 reveals that

[x] satisfies the differential

equation

D [x]

[x] ==



i =0

[[ i ]] δ[x − i]. (4.3)

4.1 Subdivision for B-splines 95

An equivalent method of arriving at this differential equation starts from the

integral definition of B-spline basis functions expressed in terms of the integral

operator

I [x]:

[x] = I [x](n

m−1

[x] − n

m−1

[x − 1]).

(4.4)

Multiplying both sides of this equation by D[x] yields a differential recurrence for

the B-spline basis functions of the form

D [x]n

[x] == n

m−1

[x] − n

m−1

[x − 1].

Iterating this recurrence

m times and substituting δ[x] for n

[x]

in the base case yields

equation 4.3 as desired.

To derive the differential equation for an arbitrary B-spline

p[x], we recall

that

p[x] can be written as a linear combination of translated basis functions:



[[ i ]] n

[x − i], where p

is the initial vector of control points. Substituting this

definition into the previous equation reveals that

p[x] satisfies a differential equation

of the form

D [x]

p[x] ==



[[ i ]] δ[x − i], (4.5)

where dp

is the vector of coefficients of the generating function d

[x]p

[x]. Note

that these coefficients are simply the

mth differences of successive control points

on the grid Z. Because the right-hand side of this equation is zero between

integer knots, the segments of

p[x] are polynomials of degree m − 1. Due to the

placement of the Dirac delta functions at the integer knots, the function

p[x] has at

least

m−2

continuity at these knots.

4.1.2 A Finite Difference Equation for B-splines

Our next task is to develop a finite difference equation that models equation 4.5.

Although finite difference equations and their associated finite difference schemes

are the subject of entire books (see [96, 144] for examples), the basic idea behind

finite differencing is relatively simple: develop a set of discrete linear equations

that approximates the behavior of a continuous differential equation at a set of

sample points. If these sample points are chosen to lie on a uniform grid, these

linear equations can be expressed succinctly using generating functions. The most

important generating function used in this approach is the finite difference mask

1 − x. The coefficients of this mask are the discrete analogs of the left-hand side of

the differential equation (i.e., of the differential operator).

96 CHAPTER 4 A Differential Approach to Uniform Subdivision

To discretize equation 4.5, our first task is to construct a difference mask d

[x]

that models the action of the differential operator D[x]

on the

-integer grid.

Luckily, we can rely on our observation from equation 3.7:

(1)

[x] = lim

k→∞



p[x] − p



x −



(4.6)

Given a generating function p

[x] whose coefficients are samples of a continuous

function

p[x] on the knot sequence

Z, we can compute a new generating function

whose coefficients approximate

D [x] p[x] on

Z by multiplying p

[x] with the gen-

erating function

(1 − x). The coefficients of 2

(1 − x)p

[x] can be simply viewed

as discrete divided differences of

taken on the grid

Z. Following a similar line

of reasoning, the discrete analog of the

mth derivative operator D[x]

Z is a

generating function

[x] of the form

[x] = (2

(1 − x))

== 2

(1 − x)

Note that multiplying the difference mask d

[x] by a generating function p[x] yields

a generating function whose coefficients are the

mth differences of the coefficients

p[x]. For example, the generating function d

[x]p[x] has coefficients that are

second differences of the coefficients

p[[ i ]] of p[x ] and that have the form p[[ i ]] −

2p[[i − 1]] + p[[ i − 2]]

. The corresponding coefficients of d

[x]p[x] are mth-divided

differences of the coefficients of

p[x] with respect to the grid

The difference mask

[x] has the property that if the coefficients of p

[x]

are samples of a polynomial of degree m − 1 on the grid

Z then d

[x]p

[x] is

exactly zero. This observation is simply a rephrasing of the fact that the

mth divided

difference of the sample of a polynomial function of degree

m − 1 on a uniform grid

is exactly zero. Because

m remains constant through the entire construction that

follows, we drop the superscript

m in d

[x] and instead use d

[x].

To complete our discretization of equation 4.5, we must settle on a discrete

analog of the Dirac delta

δ[x]. Recall that by definition this function is zero every-

where but the origin; there it has an impulse with unit integral. The discrete analog

δ[x] on the grid

Z is a vector that is zero everywhere except at the origin, where

the vector has value

. (The factor of 2

is chosen so that this discrete approxima-

tion has unit integral with respect to the grid

Z.) The coefficient of this vector

can be represented as the generating function

because all terms of this generating

function of the form

have coefficient zero for i = 0.

Finally, integer translates of the Dirac delta

δ[x − i] can be modeled on the

grid

Z by multiplying the generating function for δ[x], 2

, by the monomial x

4.1 Subdivision for B-splines 97

Absorbing these powers into the generating function dp[x] == d

[x]p

[x] yields a

finite difference equation of the form

[x]p

[x] == 2

. (4.7)

Given p

[x], this finite difference equation would ideally have a single solution

[x]. Unfortunately, the solution p

[x] to this equation is not unique because the

difference equation

[x]p

[x] == 0 has m solutions corresponding to the samples

of polynomials of degree

m − 1. However, there does exist one sequence of the

solutions

[x] that converges to the B-spline

p[x]. The next section demonstrates a

simple method of constructing this sequence of solutions.

4.1.3 The Associated Subdivision Scheme

For each value of k, equation 4.7 characterizes a solution p

on the grid

Z that is

independent of the solution

k−1

on the coarser grid

k−1

Z. However, the beauty of

this particular finite difference equation is that there exists a simple recurrence for

constructing an entire sequence of solutions

to equation 4.7. In particular, these

solutions can be constructed using a subdivision scheme.

THEOREM

4.2

Given an initial set of control points

, consider the function

[x]

satis-

fying the recurrence relation

[x] = s

k−1

[x]p

k−1

], where

k−1

[x] =

k−1

]

[x]

(4.8)

for all k > 0. Then p

[x] satisfies the finite difference relation of equa-

tion 4.7 for all

k ≥ 0.

Proof The proof is by induction on k.Fork == 0, equation 4.7 follows reflexively.

Next, assume that

k−1

[x] satisfies equation 4.7; that is,

k−1

[x]p

k−1

[x] == 2

k−1

Replacing x with x

in these generating functions and multiplying each side

2 yields a new equation of the form

k−1

] == 2

98 CHAPTER 4 A Differential Approach to Uniform Subdivision

Substituting the right-hand side of equation 4.8 into p

[x] == s

k−1

[x]p

k−1

]

and multiplying by d

[x] yields a finite difference equation of the form

[x]p

[x] == 2d

k−1

Combining the last two equations completes the induction and proves the

theorem.

For polynomial splines, the difference mask d

[x] in equation 4.8 has the form

(1 − x))

. Therefore, the subdivision mask s

k−1

[x] has the form

k−1

[x] = 2

k−1

]

[x]

== 2

m(k−1)

(1 − x

)

(1 − x)

== 2



1 + x



This mask, s

k−1

[x], is exactly the Lane-Riesenfeld mask for B-splines derived in

Chapter 2. However, independent of our previous knowledge concerning B-splines,

we can still deduce several important properties of this scheme directly from the

structure of the subdivision mask

k−1

[x]. For example, the convergence of the co-

efficients of the solution

[x] to a limit function p

∞

[x] can be proved directly from

the structure of the subdivision mask

k−1

[x] using the techniques of Chapter 3.

Further, the function

∞

[x] possesses the following properties:

■

∞

[x] is a C

m−2

piecewise polynomial function of degree m − 1, with knots at

the integers

Z. This observation follows from the fact that p

∞

[x]

is a linear

combination of integer translates of the truncated power

[x].

■

The associated basis function n[x] for this scheme is non-negative and is

supported on the interval

[0, m]. This observation follows from the fact that

the subdivision mask

k−1

[x] has only m +1 non-zero coefficients, all of which

are positive.

■

The basis function n[x] has unit integral. This observation follows from the

fact that the coefficients of the generating function

i =1

i −1

k−i

] are con-

verging to the basis function

n[x] when plotted on the grid

Z. Because the

masks

k−1

[x] satisfy s

k−1

[1] == 2 for all k > 0, the sum of the coefficients in

this product of generating functions is

. Therefore, the integrals of these

bounded approximations to

n[x] are converging to one, and thus the integral

n[x] is also one.

The beauty of this differential approach is that all that is needed to construct

the subdivision mask

k−1

[x] is the original differential operator governing the spline.

4.2 Subdivision for Box Splines 99

Given this operator, we simply discretized it over the grid

Z and constructed the

discrete difference mask

[x]. The corresponding subdivision mask s

k−1

[x] had the

form

k−1

]

[x ]

. The rest of this chapter and the next demonstrate that this approach

can be applied to a wide range of differential operators, including those in more

than one variable.

4.2 Subdivision for Box Splines

Chapter 2 constructed box-spline scaling functions in terms of an integral recur-

rence on a set of direction vectors. In this section, we recast box splines from a

differential viewpoint. As in the univariate case, our approach is to construct ex-

plicitly the integral operator used in defining the box-spline scaling functions. Using

the differential operator that is the inverse of the integral operator, we next show

that the cone splines are Green’s functions associated with this operator. Because

the box splines can be written as linear combinations of cone splines, the partial

differential equation for box splines follows immediately. We conclude the section

by deriving the finite difference equation associated with this differential equation

and an associated subdivision mask that relates solutions to this difference equation

on successive grids.

4.2.1 A Differential Equation for Box Splines

The integral operator I [x] was the crux of our initial integral construction for uni-

variate B-splines. Using the inverse differential operator

D [x], we were able to con-

struct the differential equation for B-splines from the original integral definition.

Two-dimensional versions of these operators play analogous roles for box splines.

Given a direction

{a, b}, the integral operator I [x

, y

] computes the directional in-

tegral of a function

p[x, y] via the definition

I [x

, y

] p[x, y] =



∞

p[x − at, y − bt] dt.

One simple application of this operator is the recursive definition for cone

splines given in equation 2.21. If the set of direction vectors



 can be decomposed

into

 ∪{{a, b}}, the cone spline c



[x, y] can be expressed in terms of the cone spline



[x, y] via



[x, y] =



∞



[x − at, y − bt] dt

= I [x

, y

] c



[x, y].

100 CHAPTER 4 A Differential Approach to Uniform Subdivision

As was the case for B-splines, the base case for this recurrence can be reduced to the

case when

 is empty. In this case, c

{}

[x, y] is taken to be the Dirac delta function in

two dimensions,

δ[x, y]. Unwinding the previous recurrence to this base case yields

an explicit integral expression for the cone spline



[x, y] of the form



[x, y] =





{a,b}∈

I [x

, y

]

δ[x, y]. (4.9)

Note that if the direction vector {a, b} appears in  with a particular multiplicity

the operator

I [x

, y

] appears in this product with the corresponding multiplicity.

Given the integral representation for cone splines, our next task is to construct

a differential equation that characterizes cone splines. As in the univariate case, the

key is to construct the inverse operator for

I [x

, y

]. In two dimensions, this operator

is simply the directional derivative. Given a direction vector

{a, b}, the directional

derivative of the function

p[x, y] along the direction {a, b} is given by

D [x

, y

] p[x, y] = lim

t→0

p[x, y] − p[x − at, y − bt]

Just as in the univariate case, the operators

I [x

, y

] and D[x

, y

] are inverses of each

other in Theorem 4.3.

THEOREM

4.3

Given a direction {a, b}, let

p[x, y] be a function such that its directional

integral

I [x

, y

] p[x, y] is continuous. Then

D [x

, y

] I [x

, y

] p[x, y] == p[x, y]. (4.10)

Proof The key observation is to restrict p[x, y] to lines parallel to the vector {a, b}.

For example, consider the line

{x = x

+ at, y = y

+ bt}. The restriction of

the function

p[x, y] to this line can be viewed as a univariate function:



p[t] = p[x

+ at, y

+ bt].

On this line, the left-hand side of equation 4.10 reduces to D [t]I [t]



p[t].

Because

I [x

, y

] p[x, y] is continuous, I [t]



p[t] must also be continuous.

Therefore, due to Theorem 4.1,

D [t] I [t]



p[t] is exactly



p[t]. Because this con-

struction holds independent of the choice of

, y

}, the theorem is proved.

4.2 Subdivision for Box Splines 101

Based on this theorem, we can derive the partial differential equation governing

cone splines. Multiplying both sides of equation 4.9 by the directional derivative

D [x

, y

] for each vector {a, b}∈ yields the differential equation for cone splines:





{a,b}∈

D [x

, y

]



[x, y] == δ[x, y]. (4.11)

(Again, the differential operator D[x

, y

] should appear with the correct multiplic-

ity.) Per Theorem 2.5, the box-spline scaling function



[x, y] can be expressed as

a linear combination of translates of the cone spline



[x, y], where the coefficients

of this linear combination are based on the discrete difference mask

d[x, y] of the

form

{a,b}∈

(1 −x

). Applying this theorem to equation 4.11 yields a differential

equation for the box-spline scaling functions:





{a,b}∈

D [x

, y

]



[x, y] ==



i, j

d[[i, j]] δ[x − i, y − j]. (4.12)

To complete our derivation, a box spline p[x] can be written as a linear combination

of translated scaling functions:



i, j

[[ i , j]]n



[x − i, y − j], where p

is the initial

vector of control points. Substituting this definition into equation 4.12 reveals that

p[x, y] satisfies a partial differential equation of the form





{a,b}∈

D [x

, y

]

p[x, y] ==



i, j

[[ i , j]] δ[x − i, y − j], (4.13)

where dp

is the vector of coefficients of the generating function d[x, y] p

[x, y].As

we shall see, these coefficients

are the discrete directional difference of the

control points

analogous to the directional derivative

{a,b}∈

D [x

, y

] p[x, y]

4.2.2 The Subdivision Scheme for Box Splines

In the univariate case, our approach is to derive a finite difference equation that

models the differential equation for B-splines. In the bivariate case, we follow a

similar approach. The key to constructing this finite difference scheme is building

the discrete analog of

D [x

, y

] on the grid

. If we replace t with

in the

definition of the directional derivative, we obtain

D [x

, y

] p[x, y] = lim

k→∞



p[x, y] − p



x −

a, y −



. (4.14)

102 CHAPTER 4 A Differential Approach to Uniform Subdivision

Now, consider a vector p

whose entries are samples of a function p[x, y] on the

grid

Z. The discrete analog of the directional derivative D[x

, y

] p[x, y] is a new

vector whose entries are coefficients of the generating function

(1 − x

[x, y].

In particular, the coefficients of the generating function

(1 −x

) exactly encode

the differences required by equation 4.14. To model the product of several dif-

ferential operators, we simply take the product of their corresponding difference

masks. Specifically, the discrete analog of the differential operator

{a,b}∈

D [x

, y

]

has the form

[x, y] =



{a,b}∈

(1 − x

). (4.15)

To complete our discretization of equation 4.13, we note that the analog of

δ[x, y] on the grid

is the constant 4

. Here, the factor of four arises from the move

to two dimensions. (In

n dimensions, the constant is (2

)

.) Given these discretiza-

tions, the finite difference equation corresponding to equation 4.13 for polynomials

has the form

[x, y]p

[x, y] == 4

, y

Just as in the univariate case, this equation has many solutions. However, there

exists a sequence of solutions

[x] for this difference equation that are related

by a simple subdivision scheme. Given an initial set of control points

, these

successive solutions

[x, y] can be constructed via the subdivision relation p

[x, y] =

k−1

[x, y]p

k−1

, y

], where

k−1

[x, y] =

k−1

, y

]

[x, y]

As in the univariate case, the generating function d

[x, y] divides d

k−1

, y

], yield-

ing a finite subdivision mask. In particular, the difference mask

(1 − x

) divides

k−1

(1−x

), yielding

(1+x

). Based on this observation, the subdivision mask

k−1

[x, y] has the form

k−1

[x, y] = 4



{a,b}∈



1 + x



. (4.16)

This subdivision mask s

k−1

[x, y] is exactly the subdivision mask for box splines de-

rived in Chapter 2. Therefore, the discrete solutions

produced by this scheme

converge to the box spline

p[x, y] that satisfies equation 4.13.

4.3 Subdivision for Exponential B-splines 103

4.3 Subdivision for Exponential B-splines

The first section in this chapter outlined a method of deriving the subdivision

scheme for B-splines directly from the homogeneous differential equation

(m)

[x] == 0. The crux of this construction is that the subdivision mask s

k−1

[x] for

B-splines has the form

k−1

]

[x ]

, where d

[x] is the discretization of the differen-

tial operator

D [x]

on the grid

Z. Remarkably, this construction can be extended

to inhomogeneous differential equations of the form



i =0

( i )

[x] == 0, (4.17)

where the β

are real constants. In this case, the resulting subdivision mask s

k−1

[x]

defines

m−2

splines whose pieces are exponential functions (the solutions to the

differential equation 4.17). These splines, the exponential B-splines, have been stud-

ied in a number of papers (e.g., by Dyn and Ron [57], Koch and Lyche [88], and

Zhang [168]).

In this section, we apply our differential approach to the problem of con-

structing the subdivision scheme for exponential splines. Instead of attempting

to construct the differential equation for exponential B-splines from an integral

definition, we simply construct a finite difference scheme of the form

[x]p

[x] == 2

, (4.18)

where d

[x] is a difference mask that annihilates samples of the appropriate expo-

nential functions. Given

[x], successive solutions p

[x] to this difference equation

are related by a subdivision mask

k−1

[x] of the form

k−1

]

[x ]

. The coefficients of

these solutions

[x] converge to exponential splines. The remarkable feature of

this construction is its simplicity. All that needs to be done is to develop a discrete

version

[x] of the differential operator in equation 4.17. Given this difference

mask, the subdivision mask for the scheme follows immediately.

4.3.1 Discretization of the Differential Equation

Before discretizing equation 4.17, we first rewrite this equation in terms of the

differential operator

D [x]. Due to the linearity of differentiation, linear combinations

of various powers of

D [x] can be used to model differential operators of higher order.

In particular, equation 4.17 can be rewritten as

(



i =0

D [x]

)p[x] == 0. If we assume,