Symmetric Box-splines on the A* Lattice
Minho Kim
School of Computer Science,
University of Seoul,
Siripdae- 13, Dongdaemun-gu, Seoul, 130-743, 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 non-Cartesian root lattices,
such as the BCC (Body-Centered 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 non-Cartesian root lattice
An* by a symmetric box-spline 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 Cartesian-shifted box-spline family Mr, document
the main properties of M,* and the partition induced by its knot planes and construct, in n variables, the
optimal quasi-interpolant of M'.
1. Introduction
Box-splines shifted on the Cartesian lattice are a useful generalization of uniform B-splines to several vari-
ables. In particular, a family Mr of n-variate box-splines is justly popular due to their linear independence
and approximation properties (Section 3.3). Members of Mr are defined by r-fold convolution, in the n
directions of the Cartesian grid plus a diagonal, so that the footprint of these box-splines is asymmetrically
distorted in the diagonal direction. To make reconstruction of vector fields less biased, convolution and
shifts on 2- and 3-dimensional non-Cartesian lattices have recently been advocated [33, 16, 17, 18, 15, 24].
For example, Kim et al. [24] show that reconstruction by a trivariate C1 6-direction box-spline of data on
the FCC (Face-Centered Cubic) lattice both is more time-efficient by ;.'. and results in less aliasing of level
sets than the standard C1 tri-quadratic B-spline for the same number of samples on the Cartesian grid.
Entezari et al. [16] show that the quality of reconstruction of the C2 tri-cubic B-spline on the Cartesian
grid is matched by reconstruction on the BCC lattice with the 8-direction C2 box-spline, 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 tensor-product B-spline of
the same smoothness and approximation order.
In this paper, we generalize the bivariate box-splines on the hexagonal lattice and the trivariate box-splines
on the BCC lattice to symmetric n-variate box-splines 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 well-known in crystallography
and discrete geometry. There it occurs (and is therefore defined as) a lattice embedded in an n-dimensional
hyperplane of R"1. By contrast to this standard formulation, we re-define 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 box-spline Mj on the A* lattice simplifies to the classical construction of n-variate box-splines by
projection: The shifts of the symmetric linear box-spline 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 (non-trivially) in the service
of bringing together ideas from signal processing and spline theory to show that the best reconstruction
lattices have an associated symmetric box-spline 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 quasi-interpolant 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 quasi-interpolant 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 box-splines 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 box-splines on non-Cartesian lattices to box splines on the Cartesian grid and shows how
quasi-interpolation can be inherited by change-of-variables. Section 5 shows that the lattice A* allows for a
symmetric box-spline 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 box-spline family I/, on A* .
2. Related Work
Piecewise linear hat functions, in particular the shifts of the bivariate 3-direction linear box-spline and of
the trivariate 4-direction linear box-splines 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 box-splines, obtained by convolution along the mesh directions, require structured meshes. For a
small sample of the literature on the bivariate 3-direction box-spline 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(Bernstein-Bezier)-
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 8-direction box-spline, M2.
On the n-dimensional 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 well-known 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 bottom-up 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 hex-splines on the hexagonal lattice which share many properties with
the box-splines on the hexagonal lattice. Similarly, the BCC lattice is the optimal 3D sampling lattice for
functions with isotropic and band-limited 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 8-direction
box-splines 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 k-th unit vector,
In the n x n identity matrix,
0:= [0 ... 0] the zero vector,
j : [1 ... 1] t the 'diagonal vector',
H' the n-dimensional 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 multi-set (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), box-spline 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 j-th 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 n-dimensional 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 n-dimensional 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( G-tk), (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 band-limited 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 low-pass 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 multi-dimensional signal processing. Assuming the frequency of the input
signal is isotropic and band-limited, we can reconstruct the original signal using a sphere-shaped 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 band-limited frequencies, the n-dimensional optimal
sampling lattice is the dual of the n-dimensional 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 (n-dimensional) 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 (Face-Centered Cubic) lattice.
Aj Dj is equivalent to the BCC (Body-Centered Cubic) lattice.
An* is the optimal sampling lattice in 2- and in 3-dimensions [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 n-dimensional equilateral simplex.
Lemma 2 (Classic geometric construction of An and A* in R"+1).
(i) Let a0 be an equilateral n-dimensional 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 non-invertible 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
-In-tl E nxn hence U1
-j
jt -1 1 O
-In-1 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 box-splines M,
We briefly review the later-referenced facts about box-splines following de Boor et al. [13] and introduce the
box-splines Mr.
A box-spline ML is defined by its matrix (multi-set) 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 box-spline
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 m-dimensional half-open 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)]
E-l{} = 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 box-spline 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 box-spline 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 box-spline Ms provides a fast way of approximating a function f with a spline
Q=f E S= [13]. Here we focus on the quasi-interpolant 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 multi-index 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 a-power function
"la :x a"/a?! f x(v)av
l a( (V)!*
The Box-Spline Mr. Box-splines defined by possibly repeated (n + 1) distinct convolution directions are
also called box-splines 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 box-spline with multiplicity r in each direction is defined by the n x r(n+ 1) matrix of directions T, [13,
page 80] with the multi-set
Tr : T and we abbreviate Mr := MTT. (25)
j-1
As pointed out in Section 2, this family of box-splines has been widely used. Since T1 [1 -1] in the
univariate case, M, can be viewed as a generalization of the uniform B-splines of odd degree to arbitrary
dimensions.
4. Box-splines on Non-Cartesian Lattices
By (12), given a square generator matrix G, any weighted sum of the shifts of the (scaled) box-spline
Ms := detGIMGs (26)
on the (possibly non-Cartesian) 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 (Quasi-interpolant). 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 quasi-interpolant 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(G-x k)As ((fo G) (. + k))
kcZ"
SI det G|MG(x j)A4 ((f oG) ( + G-j)),
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 non-standard representation of the efficient reconstruction
lattice A*,. This representation, AZ", is introduced in Section 5.2 and Section 5.3 defines the family of
box-splines M, := M, o A 1 on A*.
5.1. Bias of box-splines M1
The box-spline 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 box-spline ,/ /
domain lattice is A, / ,/
domain is R" / /
Table 2: Box-spline 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:
Ap-tA = I J-- / I,. (32)
n+l
Therefore M1 is a biased reconstruction filter.
By contrast, the domain lattice of the box-spline 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
(non-parallel) 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 box-splines 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 box-spline 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 box-splines and (bottom) shifts of (the support) of linear bivariate
box-splines.
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 n-dimensional simplex conv({Aij : 1 < j <
n} U {0}) is equilateral, we verify that the vectors ij satisfy
|Aii 1|2 -= 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):
IAj|l= 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
box-spline 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 Box-spline 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 box-spline 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 non-parallel 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 non-parallel 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)
(n-2) 2 2
non-parallel 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 axis-aligned, 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) (v-j) = -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 Zj-1 + 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 (Quasi-interpolant for Mf). The quasi-interpolant 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 quasi-interpolant 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. Ia|c= 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)
J-2
( 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/31-1 /31-2
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 non-standard 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 well-known box-spline 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
quasi-interpolant construction for M5.
8. Acknowledgment
This work was supported in part by NSF Grant CCF-0728797.
References
[1] E. Arge and M. Daehlen. Grid point interpolation on finite regions using C1 box splines. SIAM Journal of Numerical
Analysis, 29(4):1136-1153, 1992.
[2] A. Bejancu. Semi-cardinal interpolation and difference equations: From cubic B-splines to a three-direction box-spline
construction. Journal of Computational and Applied Mathematics, 197(1):62-77, December 2006.
[3] M. W. Bern, D. Eppstein, and J. G. Erickson. Flipping cubical meshes. Engineering with Computers, 18(3):173-187,
October 2002.
[4] G. Casciola, E. Franchini, and L. Romani. The mixed directional difference-summation algorithm for generating the B6zier
net of a trivariate four-direction box-spline. Numerical Algorithms, 43(1):1017-1398, 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 226-233, New York, NY, USA, 2002. ACM.
[6] Y.-S. Chang and H. Qin. A unified subdivision approach for multi-dimensional non-manifold modeling. Computer-Aided
Design, 38(7):770-785, 2006.
[7] C. K. Chui, K. Jetter, and J. D. Ward. Cardinal interpolation by multivariate splines. Mathematics of Computation,
48(178):711-724, April 1987.
[8] C. K. Chui and M.-J. Lai. Algorithms for generating B-nets and graphically displaying spline surfaces on three- and
four-directional meshes. Computer Aided Geometric Design, 8(6):479-493, 1991.
[9] J. H. Conway and N. J. A. Sloane. Sphere Packings, Lattices and Groups. Springer-Verlag New York, Inc., New York,
NY, USA, 3rd edition, 1998.
[10] C. de Boor and K. Hollig. Approximation order from bivariate C1-cubics: A counterexample. I of the American
Mathematical Society, 87(4):649-655, 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):13-28, 1983.
[12] C. de Boor, K. Hollig, and S. Riemenschneider. Bivariate cardinal spline interpolation on a three direction mesh. Illinois
Journal of Mathematics, 29:533-566, 1985.
[13] C. de Boor, K. Hollig, and S. Riemenschneider. Box splines. Springer-Verlag New York, Inc., New York, NY, USA, 1993.
[14] D. E. Dudgeon and R. M. Mersereau. Multidimensional Digital Signal Processing. Prentice-Hall, 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 11-18, 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 227-255. 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):313-328, 2008.
[19] P. O. Frederickson. Quasi-interpolation, extrapolation and approximation on the plane. In Proc. Manitoba Conf. on
Numerical Mathematics, pages 159-176, Winnipeg, Canada, 1971.
[20] H. Freudenthal. Simplizialzerlegungen von beschrankter Flachheit. Annals of Mathematics, 43:580-582, 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):401-417, November 2005.
[22] K. Jetter and P. G. Binev. Cardinal interpolation with shifted 3-directional box splines. Proc. Royal Soc. Edinburgh, Sect.
A, 122:205-220, 1992.
[23] R.-Q. Jia. Approximation order from certain spaces of smooth bivariate splines on a three-direction mesh. Trans. Amer.
Math. Soc., 295:199-212, 1986.
[24] M. Kim, A. Entezari, and J. Peters. Box spline reconstruction on the face-centered cubic lattice. IEEE Transactions on
Visualization and Computer Graphics, 14(6):1523-1530, 2008.
[25] H. R. Kiinsch, E. Agrell, and F. A. Hamprecht. Optimal lattices for sampling. IEEE Transactions on Information Theory,
51(2):634-647, 2005.
[26] M.-J. Lai. Fortran subroutines for B-nets of box splines on three- and four-directional meshes. Numerical Algorithms,
2(1):33-38, 1992.
[27] A. Levin. Polynomial generation and quasi-interpolation in stationary non-uniform subdivision. Computer Aided Geo-
metric Design, 20(1):41-60, 2003.
[28] J. Martinet. Perfect Lattices in Euclidean Spaces. Springer-Verlag 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 265-272, New York, NY, USA,
2007. ACM.
[30] R. M. Mersereau. The processing of hexagonally sampled two-dimensional signals. I of IEEE, 67(6):930-949,
June 1979.
[31] X. Shi and R. Wang. Spline space and its B-splines on an n + 1 direction mesh in R". J Comput. Appl. Math.,
144(1-2):241-250, 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 91-98, 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. Hex-splines: A novel spline f ..... for
hexagonal lattices. IEEE Transactions on Image Processing, 13(6):758-772, June 2004.
|