Fast and stable evaluation of boxsplines via the
Bizier form
Minho Kim and Jorg Peters
October 10, 2007
Abstract
To repeatedly evaluate linear combinations of boxsplines in a fast and
stable way, in particular along knot planes, we convert to and tabulate
the boxspline as piecewise polynomials in Bezier form. We show that
the Bezier coefficients can be stored as integers and a rational scale fac
tor and derive a hash table for efficiently accessing the B6zier pieces. The
preprocessing, the resulting evaluation algorithm and use in a widely avail
able raytracing package are illustrated for splines based on two trivariate
boxsplines, the 7direction boxspline on the Cartesian lattice and the
6direction boxspline on the FCC lattice.
1 Introduction
As a generalization of the univariate uniform Bspline to multivariate shift
invariant lattices, boxsplines are useful in many applications. For example,
they can be used to create a continuous field from data sampled at the lattice
points. This reconstruction typically uses quasiinterpolation or outright inter
polation and requires exact values of the box splines. But evaluation on lattice
edges and faces, a likely operation for a field so constructed, requires care. Al
ready de Boor [5] and later Kobbelt [12] observed a fundamental combinatorial
challenge due to the inclusion or exclusion of certain knot planes (Section 3.1.5)
and dealt with it in two different ways in their respective recursive boxspline
evaluation algorithms. Our interest was piked by the example of Figure la,
where the algorithm of [12] fails due to subtle numerical roundoff in the under
lying MATLAB routine.
An alternative approach to direct recursion is to evaluate after conversion to
the Bezier form. This approach was pioneered by Chui and Lai [2, 13] in two
variables, and recently by Casciola et al. [1] for a class of trivariate boxsplines.
Since any spline is a linear combination of shifts of one or more boxsplines,
we focus attention on the boxspline basis (or generator) functions. The point
of the conversion is that the total degree Bezier form of the polynomial pieces
(a) (b)
Figure 1: Isosurfaces for 10 1 (blue), 10 2 (green), 103 (red) and 10 10 (pur
ple) of the 6direction boxspline basis function (Section 6.2), (a) evaluated by
[12] and (b) the correct isosurface.
(Section 3.2.1) has a stable evaluation algorithm exactly on the knot planes
where the recursive algorithms encounter its challenges. In fact, de Casteljau's
algorithm is even of lower complexity along knot planes. The key challenge for
this approach are the derivation and exact representation of the change of basis.
Our first contribution, generalizing [13, 1], is Theorem 1:
(1) the Bezier coefficients, expressing a boxspline with integer directions, are
rational.
This allows us to use a simple interpolationbased approach for deriving a
changeofbasis matrix with exact integer entries, scaled by a rational number.
While both the recursive and the conversion approach benefit from localization,
i.e. from determining what boxsplines influence the evaluation, the conversion
approach must efficiently determine the simplex of Bezier evaluation. This
forces an understanding of the decomposition of the domain by the knot planes
implied by the boxspline directions. Our second contribution is
(2) an indexing strategy, based on the boxspline directions, for finding the
simplex of Bezier evaluation for a given parameter (Section 4.1).
The onetime preresolution of the combinatorics of the boxsplines required by
the indexing strategy is at the core of improved speed, by orders of magnitude
(Table 3), compared to evaluators that resolve the combinatorics at run time.
Conversion plus indexing, both precomputed, stored and quickly accessed, yield
an
(3) algorithm for fast evaluation of splines based on boxsplines that is stable,
in particular along knot planes (Algorithm 5.1).
2 Review of Existing Evaluation Techniques
Two different MATLAB packages for evaluating boxsplines, [5] and [12], are
based on the recursive formula (3). These packages, which are freely available
from the Internet, and immensely useful, because they accommodate arbitrary
matrices E and they are wellexplained in their companion papers. As the papers
point out, evaluators based on recursion face a key difficulty when evaluating a
combination of shifts of the characteristic function. Unless the combinatorics on
inclusion and exclusion of knot planes are correctly and consistently addressed,
evaluation along knot planes yield incorrect results. If the evaluation is correct,
it is called 'stable'.
De Boor [5] addresses the stability problem by perturbing input points that
are deemed too close to knot planes. Kobbelt [12] untangles the combinatorics
explicitly by deferring translation of input points until the base level of the re
cursion to avoid roundoff. This algorithm also precomputes the normals of the
knot planes in a deterministic way to avoid that a knot plane is doubly included
or completely excluded by two adjacent characteristic functions. We found that
[12] works well for bivariate boxsplines, but the released code fails in higher
dimensions as the example in Figure la illustrates. After analyzing the prob
lem in more detail than we had intended, we found the flaw in the application
of the MATLAB null function call. Generically, the null space is determined
by Singular Vector Decomposition [14]. But even minute SVD roundoff errors
create instability. Tellingly, we were often able to remove the instability, in the
trivariate cases we tested, by adding the 'r' parameter to the null function
call, i.e. by enforcing closetorational representation via Gaussian elimination.
The algorithms of Jetter and McCool [11, 15] approximate boxsplines by sam
pling in the Fourier domain and applying inverse FFT. This way, they leverage
the closed form of boxsplines in the Fourier domain.
Explicit formulas for the conversion of boxsplines to Bezier representation have
been derived by Chui and Lai [2, 13] in two variables and by Casciola et al. [1]
for a class of trivariate boxsplines. The approach of [2] generates the Bezier
form of bivariate boxsplines by comparing directional derivatives of boxsplines
with those of the Bezier form. Applying this approach to 3 and 4direction
bivariate boxsplines, [13] provides explicit Fortran codes and shows that the
Bezier coefficients are rational. Similarly, [1] converts the important class of
trivariate boxsplines spanned by four directions.
Condat and van de Ville [3] based the evaluation of bivariate 3direction box
splines on reduction of the boxsplines to cone splines, i.e. truncated powers.
Da(hlen [10] went one step further by converting the cone splines to simplex
splines of one dimension lower. The approach is shown to be efficient for bivari
ate boxsplines and, with explicit guidance along knot lines, for the 4direction
trivariate boxspline.
Except for [11], where the goal is interpolation, all the above aim at evaluating
individual boxsplines. Splines in boxspline form would off hand be evaluated
by evaluating shifts of the underlying boxsplines individually, and then adding
their contribution weighted by the coefficients.
3 Boxsplines and the Bezier form
In Section 3.1, we briefly review the basic definitions and properties of box
splines and in Section 3.2, we review the multivariate polynomials in Bezier
form. In Section 3.3, we prove that the Bezier coefficients of the polynomial
pieces of a boxspline with integer directions are rational.
3.1 Boxsplines
We use the notation made standard by the 'boxspline book' [6]. Here 'box
spline' denotes a smooth piecewise polynomial of finite support. A spline is a
linear combinations of shifts of a boxspline. If the shifts of a boxspline are
linearly independent, the boxspline is a basis function.
3.1.1 Definition
Geometrically, the value of a boxspline with direction matrix BE Rsxn at
x E ranE C R' is the shadowdensity [6] (13)
ME() := vollndimran3 ( 1 {} n ) // det I,
i.e. the scaled volume of the intersection of a cube L C R", n > s, with the
preimage E1 {} of x, an (n dimranE)dimensional subspace in R". The cube
or box gives the boxspline its name. In more detail,
C := [0..1)" is an ndimensional halfopen unit cube,
E is the s x n direction matrix, possibly with repeated columns, of the
boxspline M7 (Section 6.1 gives an example),
ranE is the subspace spanned by the column vectors {( : C E },
E l(x) is the preimage of x when viewing E as a linear transformation
E : RI R" and
vold (.) is the ddimensional volume of its argument.
In the following, we assume
dim ranE = s.
3.1.2 Degree and Continuity
A boxspline M7 is a piecewise polynomial on ranE. Its degree is less than or
equal to k := k(E) := #E s where #E denotes the number of columns in E.
The polynomial pieces join to form a function in C' 1(ranE) where [6] (p.9)
m := m(E) : min{#Z: Z A()} 1
and [6] (p.8) A() := {Z C : S\Z does not span}.
3.1.3 Spline Space
A spline fs spanned by a boxspline M7 is an infinite linear combination of the
shifts of the boxsplines on the integer grid [6] (p.33):
f := M(. j)a(j) (1)
where a : Z  R is a mesh function that returns the coefficient corresponding
to a mesh location.
3.1.4 Differentiation
With
Dz := IlczD( a composition of differential operators D( := C= ((1)D,
and
Vz := n czVC a composition of backward difference operators such that
vy: 0(. ),
a derivative of M7 in the directions Z C E equals the backward difference of
M \Z along them [6] (130):
DzMs = VzME\z.
3.1.5 Knot Planes
According to [6] (137), a boxspline is a piecewise polynomial with pieces delin
eated by a shiftinvariant mesh on Z~ generated by the collection of knot planes
(hyperplanes spanned by columns of E) H(E) [6] (p.16):
r(): U H + Z. (2)
The mesh F(E) decomposes RI into convex polytopes.
3.1.6 Recurrence Relation
As long as M7\ for E E is continuous at x = Et E R, t E R", the boxspline
M7 can be evaluated recursively with the recurrence [6] (143):
(n d)Ms(x) = EtM\M\(x) + (1 t) ME\((x E). (3)
Here t4 is the component of the vector t associated with the column E B by
X 7t.
3.1.7 Change of Variables
Using [6] (123), given an invertible linear transformation X E R 8X, one can
show that splines spanned by ME on the grid X 1Z~ can be expressed as splines
spanned by Mx7 on the integer grid Z":
E M (Ej)b(j)= det X Mx (X. j)b(X lj) (4)
jX1Zs jcZs
where b : X 1Z" R is a mesh function that returns the coefficient corre
sponding to a mesh location.
3.2 The B6zier Form of a Multivariate Polynomial
3.2.1 Definition
Let {vj E R' : j {1, ..., s + 1}} be a collection of vertices of a nondegenerate
simplex a. The map
/ R" R s+ : x 1 1 1x (5)
is called the barycentric coordinate function with respect to a [4]. The Bezier
form of an svariate polynomial of total degree d and coefficients c := {ca ER :
a d, a E Z28 } defined on a, is a polynomial
P(a) /. '"()
Ia=d
where
u is in barycentric coordinates w.r.t. a,
( d)= i ,
s 1 a(i)H
u" := HIl(u(i)a( and
{'.'(u) : (d)a : a = d} are the Bernstein basis polynomials of degree
d.
Denote the jth unit vector by ij. Then '. (ij) 1 if a dij and zero otherwise.
Therefore, Plj = P(ij) cdi. In other words, the Bezier form of a multivariate
polynomial interpolates its vertex coefficients {ca : a dij,j E {1, ..., s + 1}}.
3.2.2 Differentiation
The directional derivative of P along one of the edges of a is
D,vP =D, ij =d E (ca+, c, +i, .; (6)
a=d 3\ =dl
3.3 Rationality of boxsplines in B6zier Form
In [13], Lai proved for 3 and 4direction bivariate boxsplines and in [1] Casciola
et al. proved for 4direction trivariate boxsplines that the coefficients for those
boxsplines in Bezier form are rational numbers. In this section, we generalize
this observation.
We start by observing that rationality is preserved by boxsplines with rational
direction matrices.
Lemma 1. Let E Qsn and rank(E) = s. If x E Q8 then Ms(x) E Q.
Proof. We use induction on #E. Let E E QSXS. Then 1/detE E Q and
Ms(x) x= X /det  E Q where X7s is the characteristic function on EI, a
linear transformation of the unit cube O.
Now assume the claim holds for M \ That is M\ 4(x) E Q and M\ 4(x ) E
Q since x c E Q8. Since rank(E) s, there exists an invertible submatrix
E' QSXS of E. Since E' 1 E QSX, there also exists t E Q' so that x Et
Ec4 tE The recursion (3),
(n s)ME(x) M= Mt \(x) + (1 t)M\ (X E),
then implies that Ms(x) E Q. E
Now denote by F(B) the collection of all shifts of (sl)dimensional hyperplanes
spanned by columns of E (Equation (2)). Each H, E F(S) is defined by a
plane equation ( ji) = We denote by knotvertex the intersection x of s
hyperplanes H1, H, E F(E) whose normals span RS:
x = N where N : : and := (7)
7. Tj,
Lemma 2. Let EE Q"SX and rank(E) = s. Then the polynomial pieces of ME
can be represented in Bezier form with vertex coefficients in Q.
Proof. By Section 3.1.5, the piecewise polynomials of M7 can be expressed
over convex polytopes delineated by knot planes with rational normals. Since
E E QsXn, we have ni E Q8 and hence N 1 E QSS. Since the shifts are on
the integer grid, ji E Z~, and therefore all knotvertices are in Q8. Since any
sdimensional convex polytope can be decomposed into sdimensional simplices
without introducing any new vertex, the claim follows from Lemma 1. O
Lemma 2 can be extended to the main conclusion.
Theorem 1. Let EE ZSXn and rank(E) = s. Then the polynomial pieces of
M7 can be represented in Bezier form with coefficients in Q.
Proof. We use induction on #E. If E Z"8X, then the Bezier polynomial is
constant and equals the value at the vertex. By Lemma 2, this value is rational.
Since rank(E) s, for any w E Q8 there exists an y E Q" so that w = Ey
E1c y". By linearity of differentiation and [6] (130),
DwMs =M y Ms =M y V ME\.
By the induction hypothesis, M \ is in piecewise Bezier form with coefficients
in Q. And, since the knot planes are invariant under integer shifts and S E ZI,
V4ME\ is a difference of Bezier polynomials with coefficients in Q. Therefore
DMs is in piecewise Bezier form with coefficients in Q.
Now let vi and vj be any two knotvertices of a simplex (possibly obtained by
decomposition of an sdimensional convex polytope. Then w := vi vj E Q and
D, i ME has rational coefficients, i.e., in the notation of (6), (ca+i, ca+i ) E
Q~. Since the vertex coefficients are rational by Lemma 2, rationality of the
differences propagates rationality to all Bezier coefficients {ca}. D
4 Preprocessing boxsplines
We first discuss how to encode (or index) the total degree B6zier domain sim
plices for a particular direction matrix E and then how to find the change
ofbasis for conversion of boxsplines to the B6zier form. We note that the
combinatorial work of defining the partition into domain simplices of Bezier
pieces is done only once ever per boxspline generator or basis function since
the result is tabulated and quickly accessed by the following indexing strategy.
4.1 Indexing a Simplex
Unless the directions form a tensorproduct, the most convenient representation
of the polynomial pieces of a boxspline is the Bezier form on a simplex (Section
3.2.1). The challenge is to smartly index each simplex in suppME to derive, store
and efficiently access the Bezier coefficients. Our decomposition is inspired by
BSP (Binary Space Partitioning) trees (see e.g. [9, 16]). Let FL(E) C F(E)
be the set of knot planes of M7, each of which splits L into two sdimensional
subspaces. Then each path of the tree is converted into a pair of index vectors
(io, i) E Z' x {0, 1}q where q := q(E) := #FT(E) is the number of knot planes
in Fr(E).
To start with, we circumscribe the support of the boxspline. Since we assume
full rank of E, we may, possibly after change of variables (4), assume that F(E)
contains all axisaligned hyperplanes. Let I E Z" be the minimal set of grid
points such that
suppMs C Is + L,
i.e. I : {j E Z' : (j + LnsuppE) / 0}. These cubes are further partitioned
into convex polytopes by other knot planes in FL(E). If the convex polytope
is not already a simplex, we have several options, for example decomposition
into simplices without adding new vertices. A simple solution is to choose s
of the polytope vertices and define the Bezier polynomial with respect to this
(nondegenerate) simplex. de Casteljau's algorithm may then evaluate outside
the simplex. By shiftinvariance, each cube is partitioned alike. We now index
a simplex by
in E Is and
iA {0, 1}q.
In other words, in identifies a cube intersecting suppE and i, identifies a simplex
in that cube. The index vector i, is computed by membership test against all
the knot planes in FL(E):
ia(x) := U(NEx r1), U(t): 1, t > (8)
t < 0.
where Ns E Qqxs and Trs E Q define the knot planes in FL(E).
4.2 Computing the Change of Basis
An sdimensional Bezier polynomial of total degree d has (d+~) coefficients.
Since M7 is of degree n s, there are (n) coefficients. And since the Bezier
polynomials are linearly independent, we can compute the Bezier coefficients
by solving the (n) by (n) linear system obtained by evaluating (n) uniformly
distributed points in each simplex. We
(i) distribute the sample points uniformly so the system is wellformed,
(ii) avoid sampling near the knot planes so [12] or [5] evaluate stably, and
(iii) choose rational parameters so that Theorem 1 allows us to round to scaled
integers by leveraging rational computation in MATLAB.
Typically, boxsplines are symmetric and we can reduce the computation and
compress the output.
4.3 Barycentric Coordinates
We also precompute the matrices in Equation (5) that compute barycentric
coordinates u with respect to a simplex {vj} for a point x G R' in Cartesian
coordinates.
5 The Spline Evaluation Algorithm
Given the indexing (hash table), the table of Bezier coefficients and the pre
computed inverse matrices for computing barycentric coordinates, we can eval
uate splines, i.e. linear combinations of shifts of boxsplines, efficiently and sta
bly. The following algorithm evaluates a spline with boxspline directions E and
boxspline coefficients a at a point x. The steps are as follows. First, we find
the simplex index i, of the point x using the membership test (8). We then
compute the barycentric coordinates u of x with respect to the simplex. Shifts
of the index are used to pick up, via the hash table, all Bezier coefficients that
stem from boxspline shifts whose supports overlap x; and to form their linear
combination with weights from a. The coefficients of the resulting Bezier poly
nomial are stored in P. Finally, the algorithm evaluates this Bezier polynomial
with coefficient vector P. The following is the pseudocode for spline evaluation.
Algorithm 5.1: EVALUATESPLINE 3 (a,x)
i U(N(x [xj) 1)
u  ComputeBarycentric(i, x [xL)
P EiC. a(] io) CE (io, i,)
return EvaluateBezier(P, u)
Here the subscript 3, rather than an argument E, emphasizes that the algorithm
requires the preprocessing with respect to E according to Section 4,
a([x] in) G R is the boxspline coefficient with index [x] in 6 Z8,
x is the input point (in Cartesian coordinates) to be evaluated,
Ns and ril define the knot planes in FLr(), (Section 4.1)
ComputeBarycentric computes the barycentric coordinate using (5),
C s(io, i) is a vector of all coefficients of the Bezier polynomial pieces
with index (io, i), retrieved from the hash table, and
EvaluateBezier evaluates the Bezier polynomial with coefficients P at u.
6 Examples
We illustrate the initial conversion and generation of the index function for two
trivariate boxsplines that are useful for reconstructing volumetric data: the
7direction boxspline on the Cartesian lattice [17, 18, 8] and the 6direction
boxspline on the FCC lattice [7]. Both boxsplines are symmetric and of low
degree given their smoothness.
6.1 The 7direction boxspline
6.1.1 Definition
We consider the 7direction trivariate boxspline defined by the direction matrix
S1 0 0 1 1 1 1
r : 0 1 0 1 1 1 1 .
0 01 1 1 1 1
6.1.2 Degree and continuity
The boxspline is piecewise polynomial of degree k(B7) 7 3 4. It is re
markable, since ME1 6 C2 (R3), because at most 3 directions span a hyperplane,
m( 7) (7 3) 1 3.
6.1.3 Indexing tetrahedra
The cube 1 is decomposed into q = #FL(7) = 6 planes:
NJTx :
0 1 1
1 0 1
1 1 0
1 1 0
1 0 1
0 1 1
The knot planes in [FL(7) split L into the 24 tetrahedra listed in Table 1.
i, vertices i, vertices v, : (3, 5, 5)
010111 vc 1, I 110 111 110010 vc 10 100 101
110110 v, 1 ~i 101 100 100000 v', 0 001 000
110111 vc 1 5 111 101 100010 vc 0 101 001
010110 vc 1 ~i 100 110 110000 vc 0 000 100
001001 vc O 5 011 010 100111 vc 1 101 111
101000 Vc 05o 000 001 101011 Vc 1 011 001
101001 Vc 05o 001 011 101111 v c 1 111 011
001000 vc 05 010 000 100011 Vc 1 001 101
011111 v1 c 1 111 110 010100 vc 0 110 100
001101 vc 1 010 011 011000 v, o 000 010
001111 V c A1 011 111 011100 v6c o0 010 110
011101 vi v 1 110 010 010000 c 0 v o 100 000
Table 1: 7direction boxspline M7,. The knot planes in [FL(7) split L into 24
tetrahedra.
6.2 The 6direction boxspline on the FCC Lattice
6.2.1 Definition
A 6direction, trivariate boxspline can
matrix [7]:
1
0
1 1
1 0
0 1
be defined by the following direction
1 0
0 1
1 1
6.2.2 Spline space
The 6direction boxspline is associated with the FCC (FaceCentered Cubic)
lattice. To apply the indexing needed for Algorithm 5.1, we transform B6 ac
cording to (4):
1 0 0 1 0 1
6 : XG6 = 0 1 0 1 1 0 ,
0 1 1 0 0 1
where
where 1 1 1 0 1
XG:= 1 1 1 1 and X 1 1 0 .
1 1 1 0 1 1
X6 is a onetoone map between Z3 and the 3dimensional FCC lattice:
X. CZ3: fork 7E Xek =(ki+k2 k3,kl+k2+k3,1 kk2+
k3) = ((k + k2 + k3) (k + k2+k3) 2ki,(k + k2+ k3) :.) Z3
because kl + k2 + k3 is even;
Xg6Z3 C I7 forj e Z3, k1+k2 +k3 = 2(jl+j2 +j3) where k : X61j,
i.e. the sum of three components is always even.
6.2.3 Degree and continuity
The boxspline My and hence Ma, is piecewise polynomial of total degree
k(S6) 6 3 3. At most 3 directions in =6, and hence in E6, span a
hyperplane. Therefore m(B6) m(E6) (6 3) 1 2 and M7y, ME 6
C'(R3).
6.2.4 Indexing tetrahedra
We have q = #FU(6) = 5 planes defined by
[ 1 1 1 1 0 T '
N : 1 1 0 1 and : [1 2 11 1].
1 1 0 1 1
L is split into 10 tetrahedra as specified in Table 2.
7 Comparison and an Application
Table 3 illustrates the relative efficiency of [5], [12] and Algorithm 5.1. We
densely evaluated in an octant of the 6 and the 7direction trivariate box
splines. No linear combination of boxsplines is evaluated. All three MATLABe
i, vertices i, vertices
00000 000 100 010 001 11111 111 101 011 110
10010 vc 101 001 100 10011 vc 011 001 101
10001 vc 001 011 010 10000 vc 001 010 100
10110 vc 101 100 110 10111 vc 011 101 110
10100 vc 010 110 100 10101 vc 011 110 010
Table 2: The 6direction
into 10 tetrahedra.
Ia
boxspline M The knot planes in FC(E6) split L
thl
Ic
(d) (e)
Figure 2: Raytraced images of several level sets (10 1,10 2 and 10 3) of the
7direction box spline. In the bottom images, a random color is assigned to each
polynomial piece. The images are rendered by POVRay [19].
implementations are designed to handle vector input and avoid forloops. The
measurements used MATLAB on a Linux system with Intel Core"2 CPU 6400
@2.13GHz (2MB cache) and 2GB memory. The comparison shows Algorithm
5.1 to be faster by orders of magnitude. We explain the difference in speed
as the result of preresolution of boxspline combinatorics, as encoded in the
indexing and the tabulation of Bezier coefficients prior to running Algorithm
5.1. The other two packages are more general and resolve the combinatorics at
run time.
To compute highquality traces of level sets of a 3D field reconstructed by shifts
of a trivariate boxspline, we implemented a MATLAB script that exports the
Bezier form of a spline in boxspline form. Specifically, we output POVRay[19]
script format. POVRay is a popular and freely available raytracing engine.
The setup requires only adding one internal function to POVRay that evaluates
a Bezier polynomial, e.g. using de Casteljau's algorithm. Figures 2 and 3 show
examples.
(d) (e)
Figure 3: Raytraced images of several level sets (10 1,10 2 and 10 3) of the
6direction box spline. In the bottom images, a random color is assigned to each
polynomial piece. The images are rendered by POVRay [19].
time (x ratio to Alg. 5.1)
algorithm spline 213 313 413
[5] 7dir 20.273 (x144) 75.297 (x154) 187.716 (x153)
] 6dir 1.867 (x34) 7.088 (x39) 18.147 (x41)
[12] 7dir 52.728 (x375) 207.841 (x424) 550.423 (x450)
6dir 3.645 (x66) 14.035 (x78) 37.232 (x84)
7dir 0.141 0.490 1.223
Alg. 5.1
6dir 0.055 0.181 0.445
Table 3: Comparison of evaluation time of three packages for N3 points dis
tributed over: [0.5..3]3 for the 7direction boxspline and [1..3]3 for the 6
direction boxspline.
Such highquality traces is hard to obtain by subdivision, unless the raytracing
algorithm is carefully designed for recursion. Use of [5] or [12] is precluded by
their lack of speed.
Acknowledgments
References
[1] Giulio Casciola, Elena Franchini, and Lucia Romani, The mixed directional
differencesummation algorithm for generating the Bezier net of a trivariate
fourdirection boxspline, Numerical Algorithms 43 (2006), no. 1, 1017
1398.
(a) (b)
Figure 4: Level set of (a) 7direction and of the (b) 6direction boxspline, both
computed via the B6zier form.
[2] C'l ,i I. . K. Chui and MingJun Lai, Algorithms for generating Bnets and
graphically displaying spline surfaces on threeand fourdirectional meshes,
Computer Aided Geometric Design 8 (1991), no. 6, I7' 194.
[3] Laurent Condat and Dimitri Van De Ville, Threedirectional boxsplines:
Characterization and efficient evaluation, Signal Processing Letters, IEEE
13 (2006), no. 7, 417420.
[4] Carl de Boor, Bform basics, Geometric modeling, SIAM, Philadelphia,
PA, 1987, pp. 131148.
[5] Carl de Boor, On the evaluation of box splines, Numerical Algorithms 5
(1993), no. 14, 523.
[6] Carl de Boor, Klaus H6llig, and Sherman Riemenschneider, Box splines,
SpringerVerlag New York, Inc., New York, NY, USA, 1993.
[7] A. Entezari, May 2007, personal communication.
[8] Alireza Entezari and Torsten M6ller, Extensions of the ZwartPowell box
spline for volumetric data reconstruction on the Cartesian lattice, Visual
ization and Computer Graphics, IEEE Transactions on 12 (2006), no. 5,
13371344.
[9] Henry Fuchs, Zvi M. Kedem, and Bruce F. Naylor, On visible surface gen
eration by a prior tree structures, SIGGRAPH '80: Proceedings of the 7th
annual conference on Computer graphics and interactive techniques (New
York, NY, USA), AC': I Press, 1980, pp. 124133.
[10] Morten Daehlen, On the evaluation of box splines, Mathematical methods in
computer aided geometric design (San Diego, CA, USA), Academic Press
Professional, Inc., 1989, pp. 167179.
[11] Kurt Jetter and Joachim St6ckler, Algorithms for cardinal interpolation
using box splines and radial basis functions, Numerische Mathematik 60
(1991), no. 1, 97114.
[12] Leif Kobbelt, Stable evaluation of boxsplines, Numerical Algorithms 14
(1997), no. 4, 377382.
[13] MingJun Lai, Fortran subroutines for Bnets of box splines on three and
fourdirectional meshes, Numerical Algorithms 2 (1992), no. 1, 3338.
[14] MathWorks, Matlab 7 Function Reference: Volume 2 (FO), 2006,
mathworks. com/access/helpdesk/help/pdfdoc/matlab/refbook2.pdf.
[15] Michael D. McCool, Accelerated evaluation of box splines via a parallel
inverse FFT, Computer Graphics Forum 15 (1996), no. 1, 3545.
[16] Bruce F. Naylor, Binary space partitioning trees, Handbook of Data Struc
tures and Applications, Chapman & Hall/C'l; 2005 (',l. 11 i and Sahni,
eds.), 2005.
[17] J6rg Peters, C2 surfaces built from zero sets of the 7direction box spline.,
IMA Conference on the Mathematics of Surfaces (Glen Mullineux, ed.),
Clarendon Press, 1994, pp. 463474.
[18] J6rg Peters and Michael Wittman, Boxspline based CSG blends, Proceed
ings of the fourth AC'I' I symposium on Solid modeling and applications,
SIGGRAPH, AC', I Press, 1997, pp. 195205.
[19] The POV Team, POVRay: The Persistence of Vision Raytracer,
http://www.povray.org.
