Metric-Aware Processing of Spherical Imagery
Michael Kazhdan
Johns Hopkins University
Hugues Hoppe
Microsoft Research
Abstract
Processing spherical images is challenging. Because no spherical
parameterization is globally uniform, an accurate solver must ac-
count for the spatially varying metric. We present the first efficient
metric-aware solver for Laplacian processing of spherical data. Our
approach builds on the commonly used equirectangular parameter-
ization, which provides differentiability, axial symmetry, and grid
sampling. Crucially, axial symmetry lets us discretize the Laplacian
operator just once per grid row. One difficulty is that anisotropy
near the poles leads to a poorly conditioned system. Our solution is
to construct an adapted hierarchy of finite elements, adjusted at the
poles to maintain derivative continuity, and selectively coarsened to
bound element anisotropy. The resulting elements are nested both
within and across resolution levels. A streaming multigrid solver
over this hierarchy achieves excellent convergence rate and scales
to huge images. We demonstrate applications in reaction-diffusion
texture synthesis and panorama stitching and sharpening.
Keywords: Equirectangular, panoramas, Laplace-Beltrami.
1 Introduction
Spherical images are often used to represent photographed or ren-
dered environments. Many common formats are based on the
equirectangular map, which parameterizes ray directions by zenith
and azimuth angles (Figure 1). This simple representation offers
a number of advantages. Because the map is cylindrical, lines of
constant inclination in the spherical scene are mapped to horizon-
tal lines in the parametric domain, i.e. the map has axial symmetry.
In addition, the map has a simple analytic description that is glob-
ally defined and infinitely differentiable (except at the poles). And
unlike the conformal Mercator projection, the equirectangular map
has a compact domain, which can be discretized using a regular grid
of sample points or finite elements.
The equirectangular representation does have drawbacks. One lim-
itation is that a uniform grid in the parametric domain oversamples
the sphere near the poles. Fortunately, the situation is not so dire
because no region is undersampled; the nonuniformity only con-
tributes an over-abundance of samples near the poles. Thus, as long
as the spherical signal is properly bandlimited during analysis and
discretization, the reconstruction can have uniform appearance and
be free of aliasing artifacts.
The more significant obstacle is that the distortion near the poles
complicates signal processing. Indeed, taking the metric into ac-
Figure 1: Equirectangular representation of a spherical panorama.
Figure 2: Our adapted grid hierarchy is regularly spaced along
meridians but adaptively spaced across circles of latitude. The mul-
tiresolution refinement (left-to-right) supports a set of nested finite
elements, which enable an efficient, metric-aware multigrid solver.
count when constructing the Laplace-Beltrami operator (the exten-
sion of the Laplacian to non-Euclidean domains) results in a poorly
conditioned system. For applications that require the solution of a
Poisson-like system (e.g. diffusion, wave-propagation, gradient-
domain image processing, etc.), this leads to unacceptably slow
convergence, even with conventional multigrid techniques.
One can try to address this convergence problem by using a spher-
ical parameterization with bounded distortion, but this presents ad-
ditional difficulties. There are two strategies. (1) The first is to
assume that the distortion is small enough, and define a linear sys-
tem that is metric-unaware by assigning equal weight to each el-
ement of the parameterization; unfortunately the inaccuracies lead
to undesirable artifacts for common signal-processing applications.
(2) The other strategy is to define a metric-aware system, relying
on the bounded distortion to provide good conditioning. However,
a critical aspect is that the geometry of element neighborhoods is
spatially varying, so it is necessary to compute and store the coef-
ficient masks of the discretized Laplace-Beltrami operator at each
grid node. Performing this over high-resolution grids is unwieldy.
In this work we introduce a new representation that builds on the
simple analytic description of the equirectangular map, thereby in-
heriting the advantages of differentiability, axial symmetry, and reg-
ularity. To address the poor conditioning due to the parameteric
distortion at the poles, we adaptively refine the domain grid, gener-
ating spherical elements with bounded anisotropy (Figure 2).
By retaining the underlying grid structure, we obtain finite-element
basis functions that nest within and across the multiresolution lev-
els. Thus, the representation is ideally suited for streaming multi-
grid processing. Additionally, thanks to the axial symmetry of the
adapted equirectangular grid, we need only compute and store the
mask coefficients once per row of the domain grid — a big savings
in both space and time. Combining these benefits, we present the
first spherical solver that is both efficient and metric-aware.
We demonstrate the new representation and multigrid solver using
two quite different applications: reaction-diffusion simulation (for
texture synthesis) and spherical gradient-domain image processing
(for panorama stitching and enhancement). We envision that the ef-
ficient spherical solver may also provide a contribution to other dis-
ciplines, where the accurate solution to a spherical Poisson equation
allows for the modeling of phenomena such as weather patterns,
electromagnetic potentials, and plate tectonics.
2 Related Work
Spherical Parameterization The challenge of effectively repre-
senting functions over a spherical domain has been well studied
in computer graphics. Many approaches define a partition of the
sphere by recursively subdividing the faces of platonic solids, e.g.,
cubes/octahedra [Miller and Hoffman 1984; G
´
orski et al. 2005;
Wan et al. 2007; Fu et al. 2009] and icosahedra [Fekete 1990;
Schr
¨
oder and Sweldens 1995a]. The use of a platonic solid as a base
mesh lets the final spherical partition inherit the associated symme-
tries, and the subdivision of the base faces ensures the topological
regularity of the final facets.
Spherical Processing Metric-aware spherical processing is gen-
erally difficult. G
´
orski et al. [2005] use the isolatitude property
of HEALPix samples (also present in the isocube map [Wan et al.
2007]) to achieve efficient integration over the sphere, thus enabling
spherical harmonic transforms. Schr
¨
oder and Sweldens [1995b]
demonstrate smoothing and sharpening of environment maps using
spherical wavelets; they compute local integrals using numerical
quadrature. Laplacian processing requires the evaluation of deriva-
tives over neighborhoods of elements, and the geometry of these
neighborhoods generally varies. Kazhdan et al. [2010] adapt a so-
called TOAST parameterization to perform Laplacian processing of
spherical images, but resort to a metric-unaware solver.
Metric-Aware Signal Processing Turk [1991] computes
reaction-diffusion over an irregular sampling of a triangle mesh,
adjusting diffusion based on induced Voronoi regions. Witkin
and Kass [1991] simulate reaction-diffusion over a surface patch
by accounting for the parameterization Jacobian in the diffusion
algorithm. Similarly, Stam [2003] simulates fluid flow over an
arbitrary subdivision surface by computing the spatially varying
metric on the subdivided mesh. Texture synthesis schemes such
as [Ying et al. 2001; Lefebvre and Hoppe 2006] use the varying
metric to create predistorted texture that looks correct when
mapped on the surface.
Multigrid Solvers It is well known that anisotropy leads to poor
multigrid convergence [Briggs et al. 2000]. If the anisotropy is
axis-aligned and spatially homogeneous, an effective solution is
semicoarsening, whereby the grid is only coarsened in one direc-
tion [Schaffer 1998]. Our adapted equirectangular grid is related to
semicoarsening in that we coarsen the rows independently. How-
ever, our approach is more general because we apply different
coarsening to different rows, resulting in a domain that is not re-
stricted to have a tensor-product structure.
Algebraic multigrid [St
¨
uben 2001] is a more general framework for
handling anisotropy. Its adaptive clustering leads to irregular grids
and thus requires storing matrices explicitly. In comparison, our
multigrid algorithm maintains the regularity of the system and does
not require constructing the system matrices.
3 Motivation and Approach
Spherical Poisson Problem Many interesting image process-
ing operations can be expressed in the gradient domain, includ-
ing image stitching, dynamic range reduction, removal of shadows
or reflections, and gradient-based sharpening [Agrawal and Raskar
2007]. The basic strategy is to extract gradient fields from one or
more images, modify the gradient values, and find the new image
that best approximates the desired gradients.
Using the notation of [Kazhdan et al. 2010], in the spherical setting
the user provides a vector field
˜
V : S
2
→T S
2
specifying the desired
spherical gradients of an image, a scalar field
˜
W : S
2
→R specifying
the desired values, and interpolation weights α,β , and the system
finds the best-fitting function
˜
U : S
2
→ R, minimizing:
min
˜
U
Z
S
2
α k∇
˜
U −
˜
V k
2
+ β (
˜
U −
˜
W )
2
. (1)
The function
˜
U is obtained by solving a screened Poisson equa-
tion [Bhat et al. 2008]:
(α ∆ −β)
˜
U = α ∇ ·
˜
V − β
˜
W . (2)
Here ∆ is the Laplace-Beltrami operator which accounts for the ge-
ometry of the sphere. This continuous formulation is discretized by
introducing a set of finite element basis functions and applying the
Galerkin method (Section 4). The result is a linear system
(α L −βD) u = f , (3)
where matrix L is the discretized Laplace-Beltrami operator, D is
the mass matrix, u is the vector of coefficients associated with the
finite elements, and f is the vector giving the inner products of the
constraint function with the finite elements.
Motivation for a Metric-Aware Solution One could use any
of the spherical parameterizations introduced in the previous sec-
tion to perform Laplacian processing on the sphere. Because
these parameterizations bound the distortion, a metric-unaware
formulation might seem adequate. (For example, in performing
gradient-domain processing over an adapted TOAST parameteri-
zation, Kazhdan et al. [2010] show that results obtained using a
metric-unaware solution are very similar to results obtained when
incorporating the metric.) However, even for parameterizations like
HEALPix [G
´
orski et al. 2005] and Rhombic Dodecahedron [Fu
et al. 2009], which are specifically designed to have low distortion,
the spatial nonuniformity of the finite differences/elements intro-
duces significant errors.
Figure 3 provides a visualization of this error in simple cases where
the analytic solution is known. It considers solutions
˜
U : S
2
→R to
the Poisson equation
∆
˜
U = −l ·(l + 1)
˜
Y
m
l
, (4)
where the constraint is a scaled version of a spherical harmonic
˜
Y
m
l
of degree l. Because these spherical harmonics are eigenvectors of
the Laplace-Beltrami operator, with eigenvalues −l ·(l + 1), we ex-
pect the solutions of (4) to be the original spherical harmonics
˜
Y
m
l
.
However, as Figure 3 indicates, relying on the mesh connectivity
alone to define a metric-unaware finite element system results in
significant error. It is only when derivatives and integrals are com-
puted with respect to the sphere metric that the solution accurately
reproduces the original spherical harmonics.
Resolution
N = 32 N = 64 N = 128 N = 256 N = 512 N = 1024
TOAST 6.8·10
−1
6.8·10
−1
6.8·10
−1
6.8·10
−1
6.8·10
−1
6.8·10
−1
Cube 3.4·10
−1
3.4·10
−1
3.4·10
−1
3.4·10
−1
3.4·10
−1
3.4·10
−1
HEALPix 2.7·10
−1
2.7·10
−1
2.7·10
−1
2.7·10
−1
2.7·10
−1
2.7·10
−1
Rhombic 1.2·10
−1
1.2·10
−1
1.2·10
−1
1.2·10
−1
1.2·10
−1
1.2·10
−1
Adaptive 1.4·10
−3
4.0·10
−4
1.2·10
−4
3.5·10
−5
9.7·10
−6
2.6·10
−6
Figure 3: Solutions to the Poisson equation (4) obtained using
a metric-unaware solver with various spherical parameterizations,
compared to the expected harmonic solutions (degree 1,2,3) in the
top row, and our metric-aware solver in the bottom row. (Functions
are visualized by scaling points on the unit sphere in proportion to
the absolute function value and coloring the points with positive
value in red and negative value in blue.) The table shows the RMS
errors for the degree-3 harmonic at solution convergence.
Our Approach In the discretization of the metric-aware finite el-
ement system, the coefficients of matrices L and D in (3) are expen-
sive to compute because they require inner products of functions
defined on the sphere. To avoid a situation in which the construc-
tion of the linear system is computationally more expensive than
its solution, we seek to exploit symmetries. In particular, the ax-
ial symmetry of the equirectangular parameterization lets us amor-
tize the cost of computing the coefficients across all finite elements
at a fixed latitude. Thus, if we think of a sphere discretization as
having complexity O(N
2
), computing the system coefficients using
our partition requires only O(N) work. In contrast, using a par-
tition based on subdivision of platonic solids, we could get up to
48-fold (cube/octahedron) or 120-fold (dodecahedron/icosahedron)
symmetry, still requiring O(N
2
) computation to set up the system.
(For further analysis, please see the discussion in Appendix A.)
It is important to note that although our discretization of the sphere
is not isotropic, the anisotropy is accounted for by defining all sys-
tem coefficients in terms of integrals over the sphere. Because the
solver is metric-aware, it converges to the correct solution as the
resolution N is increased, independent of the axis of the equirectan-
gular parameterization.
4 Geometry-Aware Finite Element Solution
In this section we describe the basic approach using a regular tes-
sellation of the equirectangular domain. Although this tessellation
results in a poorly conditioned solver, the design of the system in-
corporates many of the basic building blocks that we leverage in
our well-conditioned adapted discretization, described in Section 5.
The main tasks in defining a metric-aware spherical solver are to
choose the set of spherical finite elements that span the (approxi-
mating) subspace of functions, and to compute the coefficients of
the linear system using the sphere’s metric. We first show how
an equirectangular parameterization can be used to transform B-
splines defined on a planar domain into B-spline-like functions on
the sphere. Using these functions, we derive the coefficients of the
linear system. Finally, we describe how the nesting of B-splines
can be leveraged to define a multigrid solver for the spherical sys-
tem and briefly describe the implications for solving large systems.
4.1 Choosing the Finite Elements
The equirectangular parameterization Φ : [0, π]×[0,2π] →S
2
maps
a rectangle domain to the sphere as
Φ(θ, ϕ) = (sinθ cos ϕ,sin θ sin ϕ,cos θ).
To establish the subspace of spherical functions, we first create a set
of finite element basis functions over the domain and then use the
parameterization inverse to pull-back the functions to the sphere, as
shown in Figure 4.
Figure 4: Using the equirectangular parameterization (top), we
define our spherical finite elements (right) as the pull-back of pla-
nar B-splines (left) through the parameterization inverse.
We partition the planar (θ ,ϕ) domain into a grid of cells, such that
cells [0,N −1] ×[0,2N −1] cover the desired range [0,π] ×[0,2π].
Let B
N
i, j
(θ, ϕ) be the biquadratic B-spline centered on cell (i, j).
Because B-splines form a partition of unity, their sum (over the
slightly extended grid [−1, N] ×[−1,2N]) is the constant function
with value one:
∑
i∈[−1,N]
∑
j∈[−1,2N]
B
N
i, j
(θ, ϕ) = 1 for θ ∈ [0, π],ϕ ∈ [0,2π].
Second-order B-splines have many benefits, including a tensor-
product structure that simplifies computation, continuous first
derivatives, a small support that results in sparse linear systems,
and a nesting structure compatible with a multigrid solver.
In the domain interior, we set the test function to equal the B-spline,
F
i, j
= B
i, j
, and define the spherical function as its pull-back
˜
F
i, j
=
F
i, j
◦Φ
−1
, illustrated by functions F
A
and
˜
F
A
in Figure 4.
Next we create test functions near the domain boundaries in accor-
dance with the spherical topology. For the spherical function to
be 2π-periodic in the ϕ direction, we define boundary functions as
F
N
i,0
= B
N
i,0
+ B
N
i,2N
and F
N
i,2N−1
= B
N
i,2N−1
+ B
N
i,−1
, illustrated by F
B
and
˜
F
B
in Figure 4. And, for the spherical function to have con-
stant value and a well-defined (zero) derivative at each of the two
poles, we create special polar functions as follows. At the north
pole (θ = 0), we sum the B-spline functions along the two grid
rows i = −1 and i = 0, thus forming
F
N
0
(θ) = B
N
−1
(θ) + B
N
0
(θ) with B
N
i
(θ) =
∑
j∈[−1,2N]
B
N
i, j
(θ, ϕ),
illustrated by F
C
and
˜
F
C
in Figure 4. The function at the south pole
is defined similarly as F
N
N−1
(θ) = B
N
N−1
(θ) + B
N
N
(θ). Thus, the
dimension of our discretization at resolution N is (N −2) ×2N + 2.
In defining these functions, our approach is inspired by previous
work in polar subdivision (e.g. [Kar
ˇ
ciauskas and Peters 2006])
which supports the processing of quad-dominant meshes by using
polar caps to handle extraordinary vertices.
4.2 Defining the Linear System
To form the discretized linear system (3), for each pair of spherical
elements
˜
F
I
(p) = F
I
◦Φ
−1
(p) and
˜
F
J
(p) = F
J
◦Φ
−1
(p) we must
compute the coefficients:
D
I,J
=
˜
F
I
,
˜
F
J
=
Z
S
2
˜
F
I
(p) ·
˜
F
J
(p)d p
L
I,J
=
∆
˜
F
I
,
˜
F
J
=
Z
S
2
∆
˜
F
I
(p) ·
˜
F
J
(p)d p.
These integrals can be expressed in terms of the functions F
I
and F
J
defined over the parameterization domain:
D
I,J
=
Z
π
0
Z
2π
0
F
I
(θ, ϕ)F
J
(θ, ϕ)sin θ dϕ dθ
L
I,J
= −
Z
π
0
Z
2π
0
sinθ
∂ F
J
∂ θ
∂ F
I
∂ θ
+
1
sinθ
∂ F
J
∂ ϕ
∂ F
I
∂ ϕ
dϕ dθ . (5)
Because the functions F
I
are tensor products of B-splines, com-
puting the system coefficients reduces to integrating products of
polynomials with trigonometric functions. For integrals of the form
R
x
n
sin(x)dx we can obtain a closed-form value using the identities:
Z
x
n
sinx dx = −x
n
cosx + n
Z
x
n−1
cosx dx
Z
x
n
cosx dx = x
n
sinx −n
Z
x
n−1
sinx dx.
For the cosecant term arising in the Laplace-Beltrami coefficients,
R
x
n
/sin(x)dx, we can compute the integral by evaluating the poly-
logarithm function (e.g. using the Cephes Library [Cephes 1995]).
4.3 Nesting and Multigrid
The B-spline basis has a nesting property that enables efficient
multigrid solutions. In this section, we derive the prolongation ma-
trix to transition between levels of a multigrid hierarchy.
In 1D, a B-spline at resolution N can be expressed as the linear
combination of B-splines at resolution 2N. For second-order B-
splines, the nesting relationship is:
B
N
i
(x) = η
0
B
2N
2i−1
(x) + η
1
B
2N
2i
(x) + η
2
B
2N
2i+1
(x) + η
3
B
2N
2i+2
(x),
with interpolation weights η
0
= η
3
= 0.25 and η
1
= η
2
= 0.75.
For 2D tensor-product B-splines, the 2D interpolation weights are
simply the products of the 1D interpolation weights:
B
N
i, j
(x,y) =
3
∑
k,l=0
η
k
η
l
B
2N
2i−1+k,2 j−1+l
(x,y).
In the context of our spherical elements, the linearity of the pull-
back operator implies that we can continue to use the nesting struc-
ture of B-splines to express coarser spherical elements
˜
F
N
as the
linear combination of finer spherical elements
˜
F
2N
. One important
aspect is to show that this nesting property is satisfied for the spe-
cial polar elements. The finite element F
N
0
(θ) associated with the
north pole can be expressed as
F
N
0
(θ) = B
N
−1
(θ) + B
0
(θ) =
η
0
B
2N
−3
(θ) +
η
1
B
2N
−2
(θ) +
η
1
B
2N
−1
(θ) + η
0
B
2N
−1
(θ)+
η
0
B
2N
0
(θ) + η
1
B
2N
0
(θ)+
η
1
B
2N
1
(θ)+
η
0
B
2N
2
(θ) .
Using the fact that θ ∈ [0,π] and that the second-order B-spline
B
i
(x) is supported in the range [i −1,i + 2], the expression for the
north-pole element becomes:
F
N
0
(θ) = (η
0
+ η
1
)
B
2N
−1
(θ) + B
2N
0
(θ)
+ η
1
B
2N
1
(θ) + η
0
B
2N
2
(θ)
= F
2N
0
(θ) + η
1
B
2N
1
(θ) + η
0
B
2N
2
(θ).
Finally, since the B-splines form a partition of unity, we get:
F
N
0
(θ) = F
2N
0
(θ) +
η
1
B
2N
1
(θ) + η
0
B
2N
2
(θ)
4N
∑
j=−1
B
4N
j
(ϕ)
= F
2N
0
(θ) +
4N−1
∑
j=0
η
1
F
2N
1, j
(θ, ϕ) + η
0
F
2N
2, j
(θ, ϕ)
.
Thus, the element associated with the north pole at resolution N can
be expressed as a linear combination of elements at resolution 2N.
An analogous argument holds at the south pole.
Gathering the interpolation weights into the matrix P
2N
N
we obtain
the prolongation operator that can be used to transition between the
different resolutions of a multigrid solver.
Note that by construction, if L
N
(resp. D
N
) denotes the Laplace-
Beltrami matrix (resp. mass matrix) at resolution N, we have:
L
N
=
P
2N
N
T
L
2N
P
2N
N
and D
N
=
P
2N
N
T
D
2N
P
2N
N
.
Thus, the system matrices defined in Equation (5) automatically sat-
isfy the Galerkin condition and we can compute the matrices for the
coarser resolutions directly, without having to perform the sparse
matrix multiplications required by an algebraic multigrid solver.
4.4 Efficient System Solution
Creating the linear system using an equirectangular parameteriza-
tion provides two efficiency advantages.
First, due to rotational symmetry, for any pair of spherical elements
˜
F
i, j
and
˜
F
i
0
, j
0
, and any longitudinal offset δ , we have:
h
˜
F
i, j
,
˜
F
i
0
, j
0
i = h
˜
F
i, j+δ
,
˜
F
i
0
, j
0
+δ
i
h∆
˜
F
i, j
,
˜
F
i
0
, j
0
i = h∆
˜
F
i, j+δ
,
˜
F
i
0
, j
0
+δ
i.
Thus, we can re-use the computation performed in defining the row
of the system matrix corresponding to element F
i,0
to set the values
for all other elements F
i, j
. This is crucial for large datasets, where
the system matrix becomes too large to store and computation of the
system coefficients at high precision would overwhelm the solver
time if it had to be performed for every element independently.
The second advantage is that the resulting matrix L is banded, be-
cause elements in one image row only overlap elements in nearby
rows. Specifically, for any |i −i
0
| > 2, elements F
i, j
and F
i
0
, j
0
are
always disjoint. Thus, a θ-major data ordering lets us implement
a streaming solver akin to that described in [Kazhdan and Hoppe
2008]. Such a solver requires only a few rows of the grid to be res-
ident in memory at any given time — providing a way to solve a
system of complexity O(N
2
) with only O(N) working memory.
V-Cycles
0 1 2 3 4 5
Regular 1.0·10
0
8.8·10
−2
7.1·10
−2
6.2·10
−2
5.7·10
−2
5.4·10
−2
Adapted 1.0·10
0
4.1·10
−3
3.4·10
−5
1.4·10
−7
2.8·10
−10
7.2·10
−13
Residual (Degree 2)
Figure 5: Visualizations (top) and residual norms (bottom) for the
multigrid solution to the Poisson equation with degree-2 spherical
harmonic constraints, using a regular equirectangular parameter-
ization (first row) and our adapted equirectangular elements (sec-
ond row). The visualizations show the parameterization alignment
(left), the solution (center), and isophotes highlighting the solver
imprecision near the north pole (right).
5 Adapted Elements for Better Conditioning
While the finite element system developed in the previous section
offers a multiresolution hierarchy, it is difficult to directly leverage
this hierarchy in a classical multigrid solver. The challenge is that
the severe anisotropy in elements near the poles hinders the conver-
gence of traditional relaxation techniques (e.g. Gauss-Seidel).
Figure 5 demonstrates this phenomenon, providing both a visual-
ization and residual magnitudes for the multigrid solutions to the
Poisson equation, obtained when the input constraint is the degree-
2 spherical harmonic from Figure 3. As these results indicate, the
solver defined by the regular equirectangular parameterization fails
to converge to the correct solution, and barely improves with in-
creased number of V-cycles.
To overcome this, we use the equirectangular map, but adaptively
discretize its domain to bound the anisotropy of the resulting ele-
ments (Figure 2). Specifically, in discretizing the sphere at resolu-
tion N we still use N circles of latitude, but partition the k-th circle
of latitude into N(k) intervals (rather than always 2N), choosing
N(k) so that it is a power of two and the lengths of the longitudinal
and latitudinal arcs are within a factor of
√
2. Using the fact that the
circle at latitude θ has length 2π sinθ and that the longitudinal arcs
have length π/N, we set N(k) to be the power of two that satisfies:
π
N
1
√
2
≤
2π sin(θ
k
)
N(k)
<
π
N
√
2 with θ
k
=
k + 0.5
N
π.
The motivation is similar to that in the isocube map [Wan et al.
2007]. What distinguishes our method is that we restrict N(k) to
be a power of two. As a result, the geometric structure of samples
across adjacent rows remains regular, allowing us to leverage the
axial symmetry in efficiently defining a geometry-aware system.
Note that, as in [Grinspun et al. 2002], our implementation has the
B-splines centered on the facets of the parameterization and defined
independently of the size and shape of adjacent faces. Thus, the
function-space does not exhibit any T-junction discontinuities.
Defining the Linear System Thanks to the nesting of B-splines,
our adapted discretization defines a space of functions that is a sub-
space of the original space of functions. Thus, the matrix and con-
straint coefficients defined by the adapted system can also be ex-
pressed in terms of the coefficients of the non-adapted system.
The key to establishing the transition between the adapted and non-
adapted system is the definition of a refinement operator R that ex-
presses elements in the adapted system as linear combinations of
elements in the non-adapted system from Section 4. Using R, the
non-adapted system (αL −β D)u = f can be transformed into the
adapted system (α
¯
L −β
¯
D) ¯u =
¯
f , with:
¯
L = R
T
L R,
¯
D = R
T
DR,
¯
f = R
T
f and u = R ¯u.
We define the operator R per row of the domain grid. To do this,
we express an element defined by sampling a circle at resolution
N(k) = 2
c
in terms of the elements defined by sampling the circle at
resolution 2N = 2
C
. Assuming that C ≥c, this can be accomplished
in an iterative manner:
• If c = C then the refinement operator is the identity and we
are done with the current row.
• Otherwise, we refine along the circle, using the weights
{η
0
,η
1
,η
2
,η
3
} = {0.25, 0.75,0.75,0.25} to express a uni-
variate B-spline sampled at 2
c
locations in terms of univariate
B-splines sampled at 2
c+1
locations and repeat the refinement
of the current row with the exponent c replaced by c + 1.