Symmetric Boxsplines on the A* Lattice
Minho Kim
School of Computer Science,
University of Seoul,
Siripdae 13, Dongdaemungu, Seoul, 130743, Korea
Jorg Peters
Dept. of CISE,
University of Florida,
CSE Bldg., University of Florida, Gainesville, FL, 32611, USA
Abstract
Sampling and reconstruction of generic multivariate functions is more efficient on nonCartesian root lattices,
such as the BCC (BodyCentered Cubic) lattice, than on the Cartesian lattice. We introduce a new n x n
generator matrix A* that enables, in n variables, for efficient reconstruction on the nonCartesian root lattice
An* by a symmetric boxspline family Mr. A* is the hexagonal lattice and Aj is the BCC lattice. We point
out the similarities and differences of M* to the popular Cartesianshifted boxspline family Mr, document
the main properties of M,* and the partition induced by its knot planes and construct, in n variables, the
optimal quasiinterpolant of M'.
1. Introduction
Boxsplines shifted on the Cartesian lattice are a useful generalization of uniform Bsplines to several vari
ables. In particular, a family Mr of nvariate boxsplines is justly popular due to their linear independence
and approximation properties (Section 3.3). Members of Mr are defined by rfold convolution, in the n
directions of the Cartesian grid plus a diagonal, so that the footprint of these boxsplines is asymmetrically
distorted in the diagonal direction. To make reconstruction of vector fields less biased, convolution and
shifts on 2 and 3dimensional nonCartesian lattices have recently been advocated [33, 16, 17, 18, 15, 24].
For example, Kim et al. [24] show that reconstruction by a trivariate C1 6direction boxspline of data on
the FCC (FaceCentered Cubic) lattice both is more timeefficient by ;.'. and results in less aliasing of level
sets than the standard C1 triquadratic Bspline for the same number of samples on the Cartesian grid.
Entezari et al. [16] show that the quality of reconstruction of the C2 tricubic Bspline on the Cartesian
grid is matched by reconstruction on the BCC lattice with the 8direction C2 boxspline, but using only
71' of data. In both cases, concrete implementations have established a computational speed advantage
corresponding to the reduction of the number of convolution directions over the tensorproduct Bspline of
the same smoothness and approximation order.
In this paper, we generalize the bivariate boxsplines on the hexagonal lattice and the trivariate boxsplines
on the BCC lattice to symmetric nvariate boxsplines M* (Section 5.3) defined by convolving along the
Email addresses: minhokim@uos.ac.kr (Minho Kim), jorg@cise.ufl.edu (J6rg Peters)
Preprint submitted to Journal of Approximation Theory
August 12, 2009
(a) n 1 (b) n 2
Figure 1: Orthogonal projection of a slab of unit cubes along the diagonal direction for (a) n = 1 and (b)
n =2.
nearest neighbor directions of the A,* lattice (Section 3.2). The A,* lattice is wellknown in crystallography
and discrete geometry. There it occurs (and is therefore defined as) a lattice embedded in an ndimensional
hyperplane of R"1. By contrast to this standard formulation, we redefine the A,* lattice directly in R" by
introducing a new n x n generator matrix A. Then the geometric construction of the shifts of the symmetric
linear boxspline Mj on the A* lattice simplifies to the classical construction of nvariate boxsplines by
projection: The shifts of the symmetric linear boxspline on A, are the orthogonal projection of a slab of
thickness 1 decomposed into unit cubes along the diagonal of the cubes (Figure 1). By comparison, M1
(see Section 3.3) has the same preimage, but for n > 2, its support is distorted by its anisotropic direction
matrix (Figure 7 (d)). Nevertheless, we can take full advantage of the close relationship of M* and M1 to
analyze M>. That is, this paper can apply existing mathematical machinery (nontrivially) in the service
of bringing together ideas from signal processing and spline theory to show that the best reconstruction
lattices have an associated symmetric boxspline family, provided that the newly derived matrix A is used
to generate A* .
Specifically, this paper documents in any number of variables n, the support, its partition, the desirable
properties shared with M, and, for the important case r = 2, the quasiinterpolant construction associated
with M'. The four theorems of the paper summarize these results: Theorem 1 introduces the new square
generator matrix A, Theorem 2 lists the properties of M,*, Theorem 3 describes the support and partition
induced by M,*, and Theorem 4 presents the optimal quasiinterpolant of M'.
Overview of the paper. The paper combines ideas from signal processing and spline theory. So, after
a review of related work in Section 2, we recall the pertinent facts of both areas used in the later proofs.
Section 3 consists of subsection 3.1: lattice packing and optimal sampling, 3.2: the root lattices AS, and A*,
and their standard geometric construction as subsets of R" 3.3: the boxsplines Mr. Readers conversant
with optimal sampling lattices and box splines (in the notation of de Boor et al. [13]) might skip Section 3
after taking a look at our symbol glossary at its beginning. Sections 4 and 5 prepare for the main Section 6.
Section 4 relates boxsplines on nonCartesian lattices to box splines on the Cartesian grid and shows how
quasiinterpolation can be inherited by changeofvariables. Section 5 shows that the lattice A* allows for a
symmetric boxspline family M,* when represented in the form A*Z" where A is a square generator matrix
(Section 5.2) different from the standard geometric construction of A* embedded in R"+ Section 6 then
documents the properties of the symmetric boxspline family I/, on A* .
2. Related Work
Piecewise linear hat functions, in particular the shifts of the bivariate 3direction linear boxspline and of
the trivariate 4direction linear boxsplines are popular basis functions for the 2D and 3D Finite Element
Method, respectively. Linear hat functions apply to general triangular or tetrahedral meshes, but higher
degree boxsplines, obtained by convolution along the mesh directions, require structured meshes. For a
small sample of the literature on the bivariate 3direction boxspline see [11, 12, 23, 7, 22, 2]. Chui and Lai
[8] and Lai i;] derived efficient evaluation of convolutions of hat functions via the BB(BernsteinBezier)
form. Casciola et al. [4] extended this approach to three variables. Chang et al. [5, 6] proposed a volumetric
subdivision scheme based on the trivariate 8direction boxspline, M2.
On the ndimensional Cartesian grid, Arge and Daehlen [1] investigated interpolation by Mr, and Shi and
Wang [31] discussed the associated spline space. The literature refers to the space decomposition corre
sponding to the polynomial pieces of M, as (n + 1)directional mesh.
The root lattices AS, and A* are wellknown in crystallography, discrete geometry and related areas. Con
way and Sloane [9] provide a standard treatise of the subject. Here the lattices are embedded in R"+1
(Section 3.2). Hamitouche et al. [21] recognized the need for square generator matrices that embed the
AS, and A*, lattices in R". Their definition, in iterative bottomup fashion, is however unnecessarily more
complex and the resulting matrices are more complicated than the ones we will present in Section 5.2.
Frederickson [19] first discussed the (symmetric) bivariate splines on the hexagonal lattice. The hexagonal
lattice is known to be the optimal sampling lattice in two dimensions and is equivalent to the A* lattice.
Van De Ville et al. [33] proposed hexsplines on the hexagonal lattice which share many properties with
the boxsplines on the hexagonal lattice. Similarly, the BCC lattice is the optimal 3D sampling lattice for
functions with isotropic and bandlimited frequencies [17, 15] and is equivalent to the Aj* lattice [9]. Entezari
et al. [16, 17, 18] and Entezari [15] were the first to investigate the (symmetric) trivariate 4 and 8direction
boxsplines on the BCC lattice.
3. Notation and Background
The dimension of vectors and matrices is either explicitly given or is determined by context. Some of the
specific vectors and matrices are:
ie the kth unit vector,
In the n x n identity matrix,
0:= [0 ... 0] the zero vector,
j : [1 ... 1] t the 'diagonal vector',
H' the ndimensional hyperplane, embedded in R"+1, including 0 and with normal j,
J, : jjt the n x n matrix composed of Is only.
The dot product is defined as x y := xty E R.
Following the convention of de Boor et al. [13], an n x m matrix will be interpreted both as
a multiset (bag) of column vectors or
a linear transformation R"  R'.
For the matrices E and Z, E\Z: {C : C GE and C Z}.
Column vectors are used as either vectors or points depending on the context.
Linear transformations, e.g., P, (Section 3.2), boxspline matrices of directions (Section 3.3), e.g., E
and Tr, and lattice generator matrices (Section 3.1), such as G, A* and A, are typeset in upper bold.
Lattices are typeset in calligraphic upper case; e.g., G and A..
v(j) denotes the jth entry of the vector v.
X(i, j) denotes the (i, j)th entry of the matrix X.
conv(P) is the convex hull of the points in P.
A matrix B E Znxm, n < m, is unimodular [13, (11.57)] if
det Z = +1, VZ C B : Z is square and rankZ = n.
If n= m and B E Z~x is unimodular then B 1 E Z"n.
3.1. Lattice packing and optimal sampling
A lattice is a discrete subgroup of maximal rank in a Euclidean vector space I'] Given an m x n matrix
G with m > n and rankG n, all integer linear combinations of its columns, G' define (the points of)
an ndimensional lattice, say G, embedded in R":
,:= {Gj eT. :j e Z}.
G is called a generator matrix of ,, and we call the columns of G a basis of G. The choice of a generator
matrix for a lattice is not unique [28].
Lemma 1. IfU E Znx" is unimodular then G and GU generate the same lattice points: G = GUZ".
Proof. Since UZ" C Z", G(UZ") C G' Conversely, G' (GU)(U Z") C GUZ" since U 1Z C
Z". O
If a lattice can be obtained from another by rotation, reflection and uniform change of scale, we say they
are equivalent, written [9]. Any ndimensional lattice ,C has a dual lattice given by
:* {x : u E z, XVu E C}. (1)
If G is a square generator matrix of then Gt is a square generator matrix of * [9]. If G is the generator
matrix of C,, an orthogonal matrix B is in the symmetry group (or automorphism group) Aut(Q,), i.e. the set
of isometries with one invariant lattice point that transform C to itself, if and only if there is a unimodular
matrix U E Z"nx such that [9] GU = BG. Therefore, the order of Aut(,) tells how symmetric a lattice
is; C and * have the same symmetry group.
In many geometric problems related to lattices, root lattices defined via root systems [9] provide good
solutions due to their inherent symmetry. Symmetry also makes them good sampling lattices for the signals
with isotropic frequencies. In this paper we focus on the root lattices A., and An*.
Let IlG(x) := kEz" 6(x Gk) be the Dirac comb function that samples a function f on the lattice GC ,
G E ]Rnxn, and denote by f(w) = {f} (w) the Fourier transform of f. Since [14]
1
T {fJGJ} ( IW) deG f( Gtk), (2)
the Fourier transform of the sampling fLlG replicates f(w) on GtZ", the (scaled) dual lattice of G' .
The choice of G determines the 'packing density' as explained next.
J ) * * * .
(a) t) F0(00 * f (C) f S( ) (d) fr .
Zi A r D (n > 3
the definition of A and Aoo.o
The sphere packing problem, densely can we pack identical spheres in is one of the oldest
?r r x.
(a) E{fmci}() (b) fmcG1(x) (c) {fImcJ2} (w) (d) fmGc2(x)
Figure 2: Lattice packing (frequency domain) and efficient sampling (primal domain). The maroon, bold
star shapes in (a) and (c) represent the bandlimited Fourier transform u{f} of a given function f; the
gray replicas are the transforms I {fJUiG } and {fLUG2} of samples fIIIG1 and fLUG2 on lattices with
generator matrices, G1 and G2 respectively. From both transforms, the original signal can be reconstructed
by removing the replicas with a lowpass filter (thick circle) and applying the inverse Fourier transform.
But the denser packing of replicas in Figure (a) is more efficient since it corresponds to a sparser sampling
lattice in the primal space, Figure (b).
Z"n, o a l i t n o t (n> 3) (n 3)
2 2"/2 n+1)/2 nn/2 2 (n+2)/2 f 31.525 (n 3)
22(n l)( 1)/2 [2(" 1) (n > 3)
Table 1: Center density of several root lattices. TD, : {(il,..., i,) : k ik is even} [9]. See Section 3.2 for
the definition of A1, and A.
The sphere packing problem, I,,." densely can we pack identical spheres in I.' ', is one of the oldest
problems in geometry [9]. The lattice packing problem is to find the lattice that induces the densest sphere
packing when the spheres are located at the lattice points. The lattice packing problem is closely related to
the optimal sampling lattice for multidimensional signal processing. Assuming the frequency of the input
signal is isotropic and bandlimited, we can reconstruct the original signal using a sphereshaped filter in
the frequency domain (Figure 2(a) and 2(c)). Since the lattice in the frequency domain is the dual of
the sampling lattice, the more densely we can pack the spheres (reconstruction filters) in the frequency
domain, the sparser a sampling lattice we can choose in the space domain to reconstruct the original signal
(Figure 2). Therefore, for input signals with isotropic bandlimited frequencies, the ndimensional optimal
sampling lattice is the dual of the ndimensional optimal sphere packing lattice [17, 25, 15].
The density of a lattice packing is the proportion of the space occupied by the spheres when packed. The
center density of a lattice is the number of the lattice points per unit volume, which can be obtained by
dividing its density by the volume of the unit sphere [9]. Therefore, larger (center) density implies that its
dual is a more efficient sampling lattice. Table 1 and Figure 3 respectively show the center density and the
density of several important root lattices. Both imply poorer sampling efficiency of the Cartesian lattice Z"
compared to other root lattices.
hexagonal
 stn
FCC A
B IC D
0
1 2 3 4 5 6 7 8 9 10
dimension
Figure 3: Density of several root lattices up to dimension 10.
3.2. The root lattices A., and A* as subsets ofRn+1
The (ndimensional) lattice A., embedded in R"+1 has points {x Z"+1 :j 1 x = } Z"+ n HjI [28] and
can be generated by the (n + 1) x n matrix [9, page 109]
1
1 1
Ac 1 EZ(n+)x". (3)

1 1
We can easily check (see also [9, page 115]) that its dual A* can be generated by the (n + 1) x n matrix
1 * 1 n/(n+ 1)
1 1/(n+1) [ jt n/(n+1) 1
A : I j/(n + 1) eR(n+1)xn.
1 1/(n + 1) 0 1/(n+ 1)
1/(n + 1)
Some examples of A, and An* are:
A2 and A are equivalent to the hexagonal lattice.
As3 D3 is equivalent to the FCC (FaceCentered Cubic) lattice.
Aj Dj is equivalent to the BCC (BodyCentered Cubic) lattice.
An* is the optimal sampling lattice in 2 and in 3dimensions [30, 32, 17, 16, 25, 18, 29, 15]. In dimensions
higher than 3, Figure 3 shows that A, packs spheres better than the Cartesian lattice, making An* a better
sampling lattice than Z".
The basis of A, can be taken from an ndimensional equilateral simplex.
Lemma 2 (Classic geometric construction of An and A* in R"+1).
(i) Let a0 be an equilateral ndimensional simplex one of whose vertices is located at the origin. Then the
n edges of an emanating from the origin form a basis of a lattice equivalent to An.
(ii) A* can be generated by the noninvertible elementary matrix
Pn+1 : In+
1 j+1 e R(n+1)x(n+1)
n+1
the orthogonal projection of the (n+ 1)dimensional Cartesian lattice 7Zn+1 along the diagonal direction
j.
Proof. (i) Let
E Znxn, hence U 1
11
1
1 1
U :
1 1
By Lemma 1,
also generates An. Since
l vl 1 2 J
Iv2 
Sv2 Vv e AcU
v,I I1 2 Vvj,Vk E AcU, vj
the simplex conv({0} UUCeAcU{V}) is equilateral hence equivalent to any a.n
(ii) For
Intl E nxn hence U1
j
jt 1 1 O
In1 0 '
we verify that
A* := A*U
1 [ (n+l)In Jn ] R(n+l)xn
n + 1 j
where A* is the matrix of the first n columns of Pn+i. The last column of Pn
combination of the first n columns, A,. By Lemma 1, the claim follows.
(6)
1 is an integer linear
3.3. The boxsplines M,
We briefly review the laterreferenced facts about boxsplines following de Boor et al. [13] and introduce the
boxsplines Mr.
A boxspline ML is defined by its matrix (multiset) of directions E. Unless mentioned specifically, we
assume that E Gnxm (m > n) and ranE = R". Geometrically, the value at x E ranE of the boxspline
ML is defined as the (normalized) shadow density of the (m n)dimensional volume of the intersection
AcU [
In ] Z(n+l)xn
J
U 0 [ 0
U: 1
between the preimage of x and the mdimensional halfopen unit cube L := [0..1)m: [13, (1.3)] (see e.g.,
Figure 1)
Ms(x) : volm (E 1{x}) n ) / det E (7)
where E is viewed as a linear transformation E : R'" R" and the preimage of x is defined as [13, (1.7)]
El{} = E( t)l {} + kerE. (8)
Let H (E) be the collection of all the hyperplanes spanned by the columns of E. We call the shifts of all the
hyperplanes in H (E) knot planes: [13, page 16]
F (): U H + Z. (9)
HEM(B)
The boxspline Ms with E E Zn"X is a piecewise polynomial function on ranE. It is delineated by the
knot planes and is of degree less than or equal to [13, page 9]
k () := m dim rans. (10)
Specifically, k (E) = m n if ranE = R".
The centered boxspline M. of Ms is [13, (1.21)]
Ms :Ms(. + /2). (11)
Given an invertible linear map L on R", [13, (1.23)]
MS = det LIMLE o L. (12)
The Fourier transform of Ms is [13, (1.17)]
1 exp(i. w)
Ms(W) {Ms ) .i 1= (13)
If Ms is centered, i.e. if Ms= Mg, then [13, page 11]
Ms(~ ) = since w). (14)
By [13, page 9], the (closed) support of Ms consists of the set
suppM = E = { tc :0 < t1 < 1} (15)
where 0 := [0..1]m is the closed unit cube and tc is the element of t associated with C by Et. Assuming
ranE = R", the set of all bases of E is denoted [13, page 8]
3(E) : {ZC E : #Z rankZ n}. (16)
The support of Ms is composed of the parallelepipeds spanned by Z E 3(E): For ranE R" there exists
points cz Ef{0, 1}", Z G 3(E), so that ED is the essentially disjoint union of the sets [13, 1.53]
ZO+az, Z e (E). (17)
The cardinal spline space [13, (II.1)]
S := span (M( j))jc (18)
is the spline space spanned by the shifts of Ms on Z". The sequence (Ms( j))jc, is linearly inde
pendent if and only if 3 is unimodular [13, page 41].
The map
M=' f Ms( j) f (j) (19)
j CZ
reproduces the polynomials in HIIM : n S where H is the set of all the polynomials on R" [13,
page 52]. Specifically, Hm() C HMw where
HI is the set of polynomials of (total) degree up to a,
m(E) : min{#Z : Z C A(E)} 1 and
A(E) : {Z C : E\Z does not span}.
In other words, Ms*' can reproduce all the polynomials up to (total) degree m(E). The following quasi
interpolant Qg for a boxspline Ms provides a fast way of approximating a function f with a spline
Q=f E S= [13]. Here we focus on the quasiinterpolant that provides the maximal approximation order
m(E) + 1: [13, page 72]
(Qsf)(x) := Ms(x j)A (f(. +j)) (20)
where A. is the linear functional [13, (111.22)]
AS=f : g..(0) (D'f)(0) (21)
and a E Z' is a multiindex with Ia := E 1 a(v). The Appell sequence {ga} in (21) can be computed
either recursively as
go: ([13, (111.19)])
[ga : 0][ EoC,(s[] )(3sp
where
s(f) : =M=(j) f(j), (22)
j
or from the Fourier transform Ms: [13, (III.34)]
ga(0) ([ iD (1/Ms))(0). (23)
Note that U]" is the normalized apower function
"la :x a"/a?! f x(v)av
l a( (V)!*
The BoxSpline Mr. Boxsplines defined by possibly repeated (n + 1) distinct convolution directions are
also called boxsplines on the (n + 1)directional mesh [1]. Given the n x (n + 1) matrix of directions
T1: [ In j ] [ il .. in j ] 1), (24)
the boxspline with multiplicity r in each direction is defined by the n x r(n+ 1) matrix of directions T, [13,
page 80] with the multiset
Tr : T and we abbreviate Mr := MTT. (25)
j1
As pointed out in Section 2, this family of boxsplines has been widely used. Since T1 [1 1] in the
univariate case, M, can be viewed as a generalization of the uniform Bsplines of odd degree to arbitrary
dimensions.
4. Boxsplines on NonCartesian Lattices
By (12), given a square generator matrix G, any weighted sum of the shifts of the (scaled) boxspline
Ms := detGIMGs (26)
on the (possibly nonCartesian) lattice G can be expressed as a weighted sum of the shifts of Ms on the
Cartesian lattice Z" by change of variables:
E M(. j)a(j) = M(G k)a(Gk) (27)
jeGZ" kcZ"
where a : G R is the mesh function (spline coefficients) on GC In the bivariate setting, de Boor and
Hollig [10, page 650] already pointed to this relationship.
We denote the spline space spanned by the shifts of Ms on G_ by
S:= span Ms( j))
This notation becomes consistent with (18) by omitting G = In and defining
S=: SLn.
Lemma 3 (Quasiinterpolant). Let D' := neG Do' be the composition of directional derivatives D,
Enj1 v(j)Dj along the columns of G and {ga} the Appell sequence of AX (21). The quasiinterpolant Q!
for SS defined by the functional
A2 (f ( +j)) : AX ((fo G) (. + G j)) (28)
g .(O) (D0 f) (), j E GZ (29)
provides the same maximal approximation power as does Qg defined by As for S=.
Proof. If we define
(Q f) (x) : (Q G (f o G)) (G x)
then, since f f o G o G1,
(f Qf)(x) ((fo G) Q (fo G)) (G x) (i Qf) (x),
10
for f : f o G and x := G x, i.e., Qi has the same approximation power as Qg. Since
S ML(Gx k)As ((fo G) (. + k))
kcZ"
SI det GMG(x j)A4 ((f oG) ( + Gj)),
j GZ"
Dk If f (G +hGik)
D (f o G) = lim
hto h
f(G)
(DGif) o G,
the corresponding functional AX is for j E G _'
A2 (f (. +j)) = As ((fo G) (. + G j))
ga(0) (D' (fo G)) (G J)
a ga(O) (oDf)(j).
a
(by (21))
5. The Representation A*Z" of A*"
Next, in Section 5.1, we show the need for a nonstandard representation of the efficient reconstruction
lattice A*,. This representation, AZ", is introduced in Section 5.2 and Section 5.3 defines the family of
boxsplines M, := M, o A 1 on A*.
5.1. Bias of boxsplines M1
The boxspline family M, and the A* lattice have a close relationship that becomes apparent when we
compare the spline spaces
ASAP : span(M+( j))jeA),zn and ST := span(Mi(. j))jez,
where M := detA *Mp + 1 and A* was defined in (6). Since
AtAp = I, Ja/(n + 1) = I jj/(n + 1) (30)
and by Sylvester's determinant theorem,
1
det(I jjt/(n+ 1)) det(Ii n/(n+ 1)) = (31)
n+1
SA t,
Sdet A : det (A tA) 1/ /n 1. Since P,+ 1 I+ J /,+/(n + 1) A*T1, the two spaces are
related by
E M j)a(j)
j eAjZ,
SMi(A k)a(A4k)
kcZ"
where A* 1 is defined in the manner of (8). The equation is similar to (27) but A* is not a square matrix!
The spline space ST though widely used, corresponds to the Cartesian domain lattice that has poorer
sampling efficiency compared to other root lattices, as pointed out in Section 3.1. Moreover, while M+
(QE (f o G)) (G x)
spline space ST1 SP ST
M1 on Z" MI+ on ApZ" MI on AZ"
symmetric boxspline ,/ /
domain lattice is A, / ,/
domain is R" / /
Table 2: Boxspline spaces related by change of variables.
is symmetric, as shown below, M1 is not (Figure 4), since, according to (30), A* is not an orthonormal
transformation:
AptA = I J / I,. (32)
n+l
Therefore M1 is a biased reconstruction filter.
By contrast, the domain lattice of the boxspline M+ is the efficient sampling lattice A* and M+ is symmetric
since the directions, i.e. the columns of P,+1, are
isometric: they have the same lengths and
isotropic: the inner product (hence the angle) between any two directions is the same.
The support of M+ inherits the symmetry of A*Z (or An) since the directions in P,+1 are taken from the
(nonparallel) directions from the origin to the nearest lattice points (of which there are 2(n+ 1), the kissing
number of A* [9]).
The shifts of M+ are the boxsplines obtained by projecting a slab as shown in Figure 1. The lattice A*Z on
the hyperplane H' C R"+1 partitions the slab. The next lemma shows that this partition can serve as an
alternate preimage of (the shifts of) suppM1, besides the box 0 E R"+1 that defines it.
Lemma 4 (support of M1). Let : [0..1]"+1 and P,+1 := I+1 Jn+l/(n + 1). The preimage ofsuppM1
with respect to the map T1 decomposes into ker T1 span(j) and P,+ C0 Hj" C R"+ :
T1 l(suppM1) = P,+ 10span(j).
Therefore T10 suppM1A T1Pn+ 1.
Proof. Recall that Al is composed of the first n columns of P,+i. By (24) and (8),
T t(T1 tl) 1 (n + 1)In Jn ] A CR(n+1)xn
n + 1 j p
and therefore
T1 {x} Ax +span(j), x e Rn,j E R"+1. (33)
By (15), suppM1 T10 1 R" hence
T1 1(T10) = AT1 espan(j) Pn+10 span(j).
P,+ 1 is symmetric since the directions, i.e. the columns of P,+1, are isometric and isotropic. However,
the domain embedded in the hyperplane H!" makes Mp,+ difficult to use in applications. We therefore now
introduce a square generator matrix A* of A* .
A* Pn+I
* t
SI
* *
* *
** A
* *
= [0..1]"+
0 0
* 0
* 0
* 0
* 0
suppM1i T1i I R"
Pn+Figue n a o
Figure 4: Symmetry of the support of Mp,,^ and asymmetry of the support of Mi.
0 0
0 0
0 0 0 0 0
0 0 0 0 0
(a) Z2
(b) A+Z2 A2
Figure 5: Geometric construction of A., in R".
5.2. The new square generator matrices A and A
To obtain a boxspline with a symmetric footprint in R" (Figure 7(c) and 7(f)), we construct simple square
generator matrices for A., and A*. Consider a linear map that scales along the diagonal j by transforming
a point x E R" according to
Sx + (j x)j,
n
where c is the scaling factor.
Theorem 1 (Geometric construction of An, (Figure 5) and A* in R" (Figure 6) ).
(i) A, can be generated by
(ii) A* can be generated by
n
c*
:= I, + J
n
with c,
with c =
1+ n1.
1
17n
T1,
0 0 0
O O O
(a) Z2 (b) A*+Z2 A
Figure 6: Geometric construction of A,* in R".
(b) M+ on H1 C R2
(c) M* on A*Z C R
A
0* i
0 .
S 0 0 0 0 0 0
S* *
(d) M0 on Z2 (e) M3 on H2 C R3 (f) M 0 on 0A*Z2 C R2
Figure 7: (top) Shifts of linear univariate boxsplines and (bottom) shifts of (the support) of linear bivariate
boxsplines.
box splines .
(a) MI on
A
Proof. (i) Any vector ij ik for j / k is parallel to H 1 and hence its length remains 2, unchanged
by A and regardless of the dimension n. To show that the ndimensional simplex conv({Aij : 1 < j <
n} U {0}) is equilateral, we verify that the vectors ij satisfy
Aii 12 = Ai Ai + ) + ( 12 1) 2 _/2. (34)
The claim follows by Lemma 2. The two different choices of c, produce the equivalent result with
respect to Hj_ 1 because I, Jn/n projects i, on Hnl
(ii) Since
1 (/ ) (1) (c + 1) (c* + 1) = cc + c + c +1,
nc, + c, + c* = 0 and hence
AtA I, + (cc* + Cn + c*)J, n,.
Under the diagonal scaling A, the length of j becomes the same as those of the unit vectors (Figure 6):
IAjl= Ij, V1 << < n. (35)
As with A, two roots of c* result in equivalent transformations with respect to Hj For example, for
n 2,
S 1 1//3 1 1/3
2 1 1/4v 1 1/4v
and for n = 3, the BCC lattice, the two choices are
S5 1 1 1 1 1
A: 1 5 1 or 1 1 1 .
6 1 1 5 2 1 1 1
5.3. AZ" as the domain lattice of M,*
We now interpret the columns of the matrix T* : A* T as direction vectors in R".
Lemma 5. T* is isometric and isotropic.
Proof. Since (35) implies isometry, we need only verify isotropy,
1 1
(A (j)) (Ai,) n+ 1 Vij and (Aik) (Xij) 1, Vi ik.
Therefore MT* has the same symmetries as A* and A*Z" A* can serve as a domain lattice for the
boxspline family (Figure 7(c) and 7(f))
Jl/ : IdetA*IMT* Mr O A 1. (36)
Since, in contrast to (33), for M1
(A* l)t(AA* 1) In,
the symmetry of Pn,+ 1 is preserved when computing the preimage,
T*l{x} A* l{x} + kerT*. (37)
1 21 2 1
6. The Symmetric Boxspline Family M* on A*Z"
By (27), the weighted sum of the shifts of M* on A* Z A can be expressed as
j) (j)= M,(A*X .j)a(AXj).
jcA*Z"
Therefore M* inherits most of the properties of Mr. In particular, by (10), Mr, hence MA*, is a piecewise
polynomial of (total) degree less than or equal to (n + 1)r n. We now summarize its properties (and those
of its scaled copy MT*, cf. (36))
Theorem 2 (Properties of M,). The boxspline M, has the following properties:
(i) M, is centered.
(ii) MT* = M T*
(iii) M* = M () is an even function.
(iv) The sequence (A1 ( j)) eAz, is linearly independent.
(v) The map Mr*' (19) can reproduce all the polynomials of (total) degree up to 2r 1:
m(T*) = 2r 1.
Proof. (i) By (11),
'I "
since EcT,* = 0.
(ii) By (14),
F { MT*} (w)
i since (e w)
CT*
1i since (w ) H sinc(e w) i since( w)
.T* {T,* T,*
because since is an even function. The claim holds since the Fourier transform is invertible.
(iii) By (12) and (i),
(iv) Since any n directions in T1 span R",
detZ 1, VZE U Ti\{} B(T1) B (T,)
and the sequence (Mr( j))jz, is linearly independent, claim (iv) follows since by (38), the shifts
of M, on the integer grid and the shifts of M,* on AZ" are related by an invertible affine change of
variables.
F{MT*I}(W)
1)
+ 2
CT,)
MT; det () IlM T o (In) MT;(.).
Figure 8: Kuhn triangulation for n = 3.
(v) Due to (38), m(T*) = m(Tr). For M1, we have to remove at least 2 directions so that the remaining
directions in T1 no longer span R", hence
m(Ti) ((n + 1) (n 1)) 1 2 1 1.
In the same way, at most r(n 1) directions in T, span a hyperplane, therefore
m(Tr) = (r(n + 1) r(n 1)) 1 = 2r 1.
Note that m(T*) does not depend on the dimension n.
Next, we characterize the partition of R" induced by the knot planes in H (Tr). Since the knot planes
generated by T: are those of F (Tr) under invertible linear transformation, the mesh inherits the topology
of the (n + 1)directional mesh.
Lemma 6 (Partition by knot planes).
(i) There are n(n + 1)/2 nonparallel planes in H (Tr).
(ii) The knot planes in H (Tj) partition the unit cube 0 into n! simplices (i ..... 8)
n i
S: conv(V,), V, : {0} U U {i,()}, 7 E S (41)
i= l j
where S, be the set of all the permutations of {1, .. n}.
The partition {cr(a}Ts, is called Freudenthal triangulation [20] or Kuhn triangulation
Proof. (i) There are n planes generated by the n unit vectors in In and ( ",) additional nonparallel planes
are spanned by the diagonal direction j and n 2 additional unit vectors yielding a total of
n + )n n + 1=n(n 1) n(n +1)
(n2) 2 2
nonparallel planes in H (Tr).
(ii) Recall that T1 [I, j]. All planes with normal direction ij ie, j / k, intersect the interior of 0
and are generated by Ti\ {ij, ik} i.e., as knot planes of M1 generated by n 1 vectors including j. Unless
two vertices vj, vk are both in V, for some permutation i7, there exist indices a and 3 so that
Vj(a) = l,Vk(a) 0 and vj(B3) 0,Vk(f3) 1
and hence the knot plane with normal i, i, separates them,
(i ip) vj 1 >0 and (i, id) Vk = 1 < 0. (42)
(a) (b)
Figure 9: (a) Rhombic dodecahedron: support of M* for n = 3 and A* 11 .Figures (b)(c):
two decomposition of the support into parallelepipeds: (b) by (43) and (c) by (44).
Conversely, since knot planes excluding j are axisaligned, neither they nor their shifts on Z" intersect the
interior of the unit cube 0. It remains to show that no shifts of the knot planes with normal i, i3 separate
vertices of a simplex au for the same fixed permutation T7. Since shifts by j E Z" within the knot plane, i.e.
j (i i) j (a) j(3) 0, result in (i i) (x j) = 0, we can assume that j (i i) = j(a) j(3) > 0.
Then, for all v E {0, 1},
j(a) +j(3)+1 < 0 v(a) 1, v(3) 0
(i i) (vj) = j(a)+j(/) 1 < 1 v(a) 0,v(l) 1
Sj(a) + j () <0 v(a)= v ()
<0.
The case j (i, i3) < 0 corresponds to a flipped normal and yields (ia i3) (v j) > 0. E
With the help of Lemma 6, we can establish the structure of suppM* by first decomposing it into paral
lelepipeds. There are two decompositions (see Figure 9).
Theorem 3 (support of M*). The (closed) support of M* is the essentially disjoint union of the (n + 1)
parallelepipeds
{ZO: Z e B(Tj)} (43)
or, alternatively,
{ZD+(z : Z e B(T)} (44)
where (z := T*\Z. In either decomposition, all the parallelepipeds are congruent.
Proof. Due to the relation (38), we need only consider M1. Let Zj E B(T1) be a basis of T1 and
cj : T\z  .
For aj := az in (17), there are only two choices, aj E {0, j}, since for any C E Zj, 2C does not fit into
T10 (cf. Figure 4, right):
VC E Zj, + C E Zj1 + but C + T1n.
Now assume
aj = 0 and a0 = C for Zj,Zk E B(Ti), Zj / Zk.
This leads to a contradiction as we prove that the two parallelepipeds Zj +cj and Zk O+atk are not
essentially disjoint but rather share the point
p := Ck + .
Cez3nZk
To verify that p is in the interior of both parallelepipeds, let G denote the disjoint union and observe that
Zj = (Zj\Zk) (Zj n Za) so that for ca = 0,
P 2 1: C + 4 E: C + "a.
S CZy\Zk C(eZynZ
That is, p E Zj(0, 1)" + aj. (we use (0, 1)" rather than 0 to show essential disjointedness).
But also p E Za(0, 1)" + ak since
1 1 1
P= Ck2 4 C Z E (ECT+ CC O)
CeZ nZ3 CETI
c 2 4 4c+2
C Zj nZk C Zk
C CczZnZZC CcZk\z
This establishes that there are only the two listed alternatives.
Next, we prove that all parallelepipeds are congruent. To analyze the decomposition
{Z'*: Z* e (Ti)} (45)
we observe that
the matrices Xa, Ya E Rnxn
X, (j, k) :a and Y, (j, k) 1j:=k= a
X(j 0 otherwise 0 otherwise
satisfy the relations
XjJn Jn, YjJn = X, XjXt = J,, XjYj = Xj, YjXj Yj, X = Xj
and that
for Zj := I, X Yj,
Zj(I, + JJ )Z In + J and Z2 I,.
Then since A2 I, + Jn, for Z Z : A*Zj and Z% : A*Za,
Z* (AZjZkA 1)Z^ (46)
Secondly, we verify that
(A*ZjZA 1)(ArZZ )t= AZjZA ZA ZZ~A*~
A= (I + J)A*
(A = A)
In.
For AZj,A*Zk E B(TI), AZjZkA*1 is therefore an orthonormal (rigid) transformation. And hence, by
(46), all the parallelepipeds Z* o, for Z* E B(Ti) are congruent. The other decomposition is verified in the
same way. E
Lemma 3 is easily extended to T* since
T* { t : 0
{CT*
{T t:0
CTT
For Z E B(TT), the pair (Z, Cz) is a linear transformation of the pair (I,, j). Therefore ZO is decomposed
in the same way as the unit cube 0 is decomposed by the Kuhn triangulation and suppM* consists of (n+ 1)!
simplices. This count also agrees with the number of modular cells in the first neighbor polytope of A*, [21].
The two types of the decomposition of suppMt in Lemma 3 can be viewed as cubical meshes such that one
is the flip of the other [3] since each cubical mesh can be viewed as the projection of the (n + 1)dimensional
cube along one fixed diagonal in two opposite directions.
Next, we expand on Theorem 2(v), which showed that M' can reproduce all cubic polynomials. The
following lemma will simplify the proof.
Lemma 7. For an odd function f, fT2 f = 0.
Proof. By definition (22),
PbT2f M2(j) f (j)
ZM2( (j)f(j)
M2(j)f( j)
5M2(j)f(j).
(by Thm (2)(iii))
(f f( ))
(change of index)
Comparing the first to the fourth line, we see ~T2f = 0.
Theorem 4 (Quasiinterpolant for Mf). The quasiinterpolant of M2, defined by the functional
A (f( +j)) := A2 (f(. +j)) :
12
12S;
Dn (j),
j e A*Z
provides the maximal approximation power m(T) + 1 =4.
Proof. We derive the quasiinterpolant QT2 for ST2 defined by ATT (21). Then QT2 for ST2 defined by A
can be derived by (28).
Specifically, we compute ga(0) for each degree Ia.
1. Iac= 0
ga(0) go(0) 1 By [13, page 68].
2. Ia 1
By Lemma 7, PT2 "
0 and
a = E (P 2ii ) 9
/3#a
therefore g,(O) = 0.
3. I l = 2
By [13, page 11],
MT2 (W) : {MT2 (W)
JI sinc(l w) = SI sinc2( w).
(tT2 T1
Therefore, By (23), for j / k,
(DJD k1) (w)
( 1( 1
STi since( w) DjD sinc2(j w)sinc2(wj)sinc2(we)
Since since ) with the help of MAPLE, we can compute
Since sinc(0) 1, with the help of MAPLE, we can compute
(DDk
2) (0)
J2
( ksinc2 (wj)sinc2 (wj)sinc2 (w + Wk) (0)
Also,
D 2 (w)
eTHsinc2(. W)
Ssinc w)sin2
sinC2j c)sin1c2 )
Again, with the help of MAPLE, we can compute
D 2 (0)
j AT2 )
S 1
S111C4(
By (23), for j / k,
gij+ik (0) ([ iD 'i +i k 1 (0)
S1AYT2 )
ir' ) (0)
2 1
2 T2
g2ij (0
HI 0 (wT2a 0) 90
4. lal 3
By [13, (III.19)],
ro I[F S (i 3I[ ").g
g39a
= na (/T2[a) 90 + (/T ) )g 93+ S ((/1T2 f "3) 9
\ 1/311 /312
hence g(0) = 0 because
fT2 l0 = 0 by Lemma 7,
p = 0 for I~ 1 and
0T2 [a O 0for / 2
hence a /3 = 1 by Lemma 7.
Summing up,
AT2f = E g(O) (Df)(0)
a 
f() E (Dof)(O)
=(o) 1 2(Dff)(o)
cT1
Now, by (29),
Af= f (O) (D f)(O).
12
For discrete input f : A*Z * R, we approximate the directional derivative along ( E R" by finite differences,
e.g.,
D f f(. + () + f( ) 2f. (49)
Therefore
Af f (0) (f() + f() 2f(0))
( 1 + n f () ) 1 (f (() + f())
{CT,
When specialized to two variables, this agrees with Levin's formula [27].
7. Conclusion
We introduced a nonstandard representation A*Z" of the efficient reconstruction lattice A* that is based on
a new family of square generator matrices A. In this representation, A* naturally admits a symmetric box
spline family M,. We then documented, in any number of variables n, the support, the induced partition
of R" and the desirable properties shared with the wellknown boxspline family Mr. For the important
case r = 2 that provides a smooth field of low degree, we derived in any number of variables an optimal
quasiinterpolant construction for M5.
8. Acknowledgment
This work was supported in part by NSF Grant CCF0728797.
References
[1] E. Arge and M. Daehlen. Grid point interpolation on finite regions using C1 box splines. SIAM Journal of Numerical
Analysis, 29(4):11361153, 1992.
[2] A. Bejancu. Semicardinal interpolation and difference equations: From cubic Bsplines to a threedirection boxspline
construction. Journal of Computational and Applied Mathematics, 197(1):6277, December 2006.
[3] M. W. Bern, D. Eppstein, and J. G. Erickson. Flipping cubical meshes. Engineering with Computers, 18(3):173187,
October 2002.
[4] G. Casciola, E. Franchini, and L. Romani. The mixed directional differencesummation algorithm for generating the B6zier
net of a trivariate fourdirection boxspline. Numerical Algorithms, 43(1):10171398, September 2006.
[5] Y.S. C ... T. McDonnell, and H. Qin. A new solid subdivision scheme based on box splines. In SMA '02: 1
of the seventh ACM symposium on Solid modeling and applications, pages 226233, New York, NY, USA, 2002. ACM.
[6] Y.S. Chang and H. Qin. A unified subdivision approach for multidimensional nonmanifold modeling. ComputerAided
Design, 38(7):770785, 2006.
[7] C. K. Chui, K. Jetter, and J. D. Ward. Cardinal interpolation by multivariate splines. Mathematics of Computation,
48(178):711724, April 1987.
[8] C. K. Chui and M.J. Lai. Algorithms for generating Bnets and graphically displaying spline surfaces on three and
fourdirectional meshes. Computer Aided Geometric Design, 8(6):479493, 1991.
[9] J. H. Conway and N. J. A. Sloane. Sphere Packings, Lattices and Groups. SpringerVerlag New York, Inc., New York,
NY, USA, 3rd edition, 1998.
[10] C. de Boor and K. Hollig. Approximation order from bivariate C1cubics: A counterexample. I of the American
Mathematical Society, 87(4):649655, April 1983.
[11] C. de Boor and K. Hollig. Bivariate box splines and smooth pp functions on a three direction mesh. Journal of Compu
tational and Applied Mathematics, 9(1):1328, 1983.
[12] C. de Boor, K. Hollig, and S. Riemenschneider. Bivariate cardinal spline interpolation on a three direction mesh. Illinois
Journal of Mathematics, 29:533566, 1985.
[13] C. de Boor, K. Hollig, and S. Riemenschneider. Box splines. SpringerVerlag New York, Inc., New York, NY, USA, 1993.
[14] D. E. Dudgeon and R. M. Mersereau. Multidimensional Digital Signal Processing. PrenticeHall, Inc., Englewood Cliffs,
NJ, 1984.
[15] A. Entezari. Optimal Sampling Lattices and Trivariate Box Splines. PhD thesis, Simon Fraser University, 2007.
[16] A. Entezari, R. Dyer, and T. Mo11er. Linear and cubic box splines for the body centered cubic lattice. In VIS '04:
I of the conference on Visualization '04, pages 1118, Washington, DC, USA, 2004. IEEE Computer Society.
[17] A. Entezari, R. Dyer, and T. Mo11er. From sphere packing to the theory of optimal lattice sampling. In Mathematics and
Visualization, pages 227255. Springer Berlin Heidelberg, 2009.
[18] A. Entezari, D. Van De Ville, and T. Mo11er. Practical box splines for reconstruction on the body centered cubic lattice.
IEEE Trans. Vis. Comput. Graph, 14(2):313328, 2008.
[19] P. O. Frederickson. Quasiinterpolation, extrapolation and approximation on the plane. In Proc. Manitoba Conf. on
Numerical Mathematics, pages 159176, Winnipeg, Canada, 1971.
[20] H. Freudenthal. Simplizialzerlegungen von beschrankter Flachheit. Annals of Mathematics, 43:580582, 1942.
[21] C. Hamitouche, L. Ibafiez, and C. Roux. Discrete topology of A* optimal sampling grids. Interest in image processing
and visualization. Journal of Mathematical Imaging and Vision, 23(3):401417, November 2005.
[22] K. Jetter and P. G. Binev. Cardinal interpolation with shifted 3directional box splines. Proc. Royal Soc. Edinburgh, Sect.
A, 122:205220, 1992.
[23] R.Q. Jia. Approximation order from certain spaces of smooth bivariate splines on a threedirection mesh. Trans. Amer.
Math. Soc., 295:199212, 1986.
[24] M. Kim, A. Entezari, and J. Peters. Box spline reconstruction on the facecentered cubic lattice. IEEE Transactions on
Visualization and Computer Graphics, 14(6):15231530, 2008.
[25] H. R. Kiinsch, E. Agrell, and F. A. Hamprecht. Optimal lattices for sampling. IEEE Transactions on Information Theory,
51(2):634647, 2005.
[26] M.J. Lai. Fortran subroutines for Bnets of box splines on three and fourdirectional meshes. Numerical Algorithms,
2(1):3338, 1992.
[27] A. Levin. Polynomial generation and quasiinterpolation in stationary nonuniform subdivision. Computer Aided Geo
metric Design, 20(1):4160, 2003.
[28] J. Martinet. Perfect Lattices in Euclidean Spaces. SpringerVerlag Berlin Hidelberg, 2003.
[29] T. Meng, B. Smith, A. Entezari, A. E. Kirkpatrick, D. Weiskopf, L. Kalantari, and T. Moller. On visual quality of optimal
3D sampling and reconstruction. In GI '07: I 1 of Graphics Interface 2007, pages 265272, New York, NY, USA,
2007. ACM.
[30] R. M. Mersereau. The processing of hexagonally sampled twodimensional signals. I of IEEE, 67(6):930949,
June 1979.
[31] X. Shi and R. Wang. Spline space and its Bsplines on an n + 1 direction mesh in R". J Comput. Appl. Math.,
144(12):241250, 2002.
[32] T. Theufil, T. Mo11er, and M. E. Groller. Optimal regular volume sampling. In VIS '01: 1 of the conference on
Visualization '01, pages 9198, Washington, DC, USA, 2001. IEEE Computer Society.
[33] D. Van De Ville, T. Blu, M. Unser, W. Philips, I. Lemahieu, and R. V. de Walle. Hexsplines: A novel spline f ..... for
hexagonal lattices. IEEE Transactions on Image Processing, 13(6):758772, June 2004.
