The PDE surface method in higher dimensions
A. Woodland1, H. Ugail2, F. Labrosse1
1Department of Computer Science, University of Wales, Aberystwyth
2School of Informatics, University of Bradford
June 6, 2008
Abstract
This paper presents a method to extend PDE surfaces to high dimen-
sional spaces. We review a common existing analytic solution, and show
how it can be used straightforwardly to increase the dimension of the
space the surface is embedded within. We then further develop a numeri-
cal scheme suitable for increasing the number of variables that parametrise
the surface, and investigate some of the properties of this solution with a
view to future work.
1 Introduction
The concept of image manifolds is one that has become increasingly common
in vision [7, 1, 5], and recently to some extent in graphics literature. Images,
which can be considered points in a high-dimensional image space, of an object
or environment undergoing a series of transformations (which are typically pa-
rameterised in some way) are not randomly placed in image space. Instead under
certain circumstances these can be approximated as n-manifolds, where n is the
number of independent variables controlling the image transformations. This
mapping between the high-dimensional image space, and a lower dimensional
space has been particularly useful in the area of vision [8], however developing
efficient in-memory representations of image manifolds has proven hard, and
many researchers have used techniques like PCA to simplify the representa-
tions. This is acceptable if the intended use is recognition of images, however
for the purpose of synthesising new images it proves to lose too many of the
details that make images appear realistic to a human observer. A potential so-
lution to this problem is to develop exisiting surface representation techniques
to model these structures in image space directly, thereby representing all of the
images that make up the manifold.
The method of representing surfaces using partial differential equations, re-
ferred to as the PDE surface method henceforth, has been proposed [2, 11] as
an alternative to methods such as NURBS [9] in a typical CAD application. In
this paper we look at extending this method to model surfaces (or more gen-
erally manifolds) which are embedded in much higher dimensional spaces than
required for a CAD application, making it suitable for modelling image man-
ifolds. Additionally we look at increasing the number of parameters used to
control the manifold, which can be more than the typical (u,v) parametrisation
1
of a surface used in a CAD application. One such existing example of using the
PDE surface method to model a volume is presented in [6].
This paper is split into five sections. In Section 2 we provide an overview of
the problem, and a common analytic solution which is then straightforwardly
extended to spaces of much higher dimension. Section 3 discusses approaches
to increasing the dimension of the parameter space, which requires a numerical
solution to the resulting PDE. This work is conducted with a view to the specific
application of image manifolds, and some early results are presented and briefly
discussed in Sections 4 and 5 respectively, along with desription of some future
work we intend to undertake as a result of this work.
2 Increasing the embedding space
There have been many different surface modelling techniques proposed in the
past, in particular parametric surface modelling.
A parametric surface X(u,v) can be expressed
X(u,v) = (x(u,v),y(u,v),z(u,v)) x,y,z ?R (1)
We now assume our surface to be periodic in v, and restrict it to the finite
domain ? = {u,v : 0 ? u ? 1;0 ? v ? 2?}. We therefore have that
parenleftbigg ?2
?u2 +?
2 ?2
?v2
parenrightbigg2
X(u,v) = 0, (2)
where ? is a ?smoothing parameter?, which controls the length over which the
boundary conditions influence the interior of the surface. This equation is often
referred to as the biharmonic
Existing work using the PDE surface method [11] has used the method
of separation of variables [4] to find an analytic solution to the PDE above,
for which a computationally efficient implementation is also possible. We also
assume now that sufficient boundary conditions have been imposed upon the
surface. The basic solution to this is now outlined.
Giventhat Eqn. (2) is elliptic, theboundaryconditions continuous andclosed
we use the following general solution:
X(u,v) = A0(u)+
?summationdisplay
n=1
[An(u)cos(nv)+Bn(u)sin(nv)], (3)
where A0, An and Bn are vector valued functions as follows:
A0 = a00 +a01u+a02u2 +a03u3 (4)
An = an1e?nu +an2ue?nu +an3e??nu +an4ue??nu (5)
Bn = bn1e?nu +bn2ue?nu +bn3e??nu +bn4ue??nu (6)
By using Fourier analysis the boundary conditions can be written in a form
whereby the values of the constants in the general solution can be calculated.
In practise not all boundary conditions have a Fourier series as simple as that
of a circle or ellipse, and for that reason the Fourier series is commonly truncated
at some given value of n, typically ([11]) N = 6. A remainder term can be used
2
to approximate the higher order components of the boundary condition that
would otherwise be lost:
X(u,v) = A0(u)+
Nsummationdisplay
n=1
[An(u)cos(nv)+Bn(u)sin(nv)]+R(u,v), (7)
where R(u,v) is usually chosen to be of the form
R(u,v) = r0e?u +r1e??u +r2ue?u +r3ue??u, (8)
where R0(v),...,R3(v) are calculated as the difference between the real bound-
ary conditions, and the approximation of them. Choosing ? = ?(N +1) gives a
decay rate similar to the decay rate of the actual solution given that the mode
(N +1) is the dominant mode of the modes that have been truncated.
Given the solution outlined above it is straight forward to increase the di-
mensionality of the space in which the surface is embedded; one simply has to
increase the dimensionality of the vector-valued functions A0, An and Bn. This
in turn simply requires one to specify boundary conditions for the PDE in each
additional dimension, and solve the corresponding linear system of equations.
3 The parameter space
Extending the surface to be controlled by more than 2 parameters (an n-
manifold) is a more challenging problem ? the solution outlined in Section 2 is
no longer applicable, and we must turn now to a numerical solution [10].
In order to solve for the case with more than two parameters we first derive
a numerical scheme for the solution of the lower order case, namely Laplace?s
equation (Eqn. (9)) which is simpler than we previously solved analytically. A
similar scheme for Eqn. (2) is straightforward to derive using the same method.
parenleftbigg ?2
?u2 +
?2
?v2
parenrightbigg
X(u,v) = 0 (9)
?
?
?
f(u,v)
Figure 1: Discretisation of the domain ?
Firstly we discretise our domain ? and produce a finite grid (Fig. 1) that
represents it. Then we seek a discrete approximation in terms of the grid of the
various terms involved, firstly wrt. u. We assume here a uniform grid to simplify
3
the problem, although this need not always be the case. In this case we have
used a centred difference scheme, because it is a more accurate approximation
of the derivatives.
?
?u ?
1
?
bracketleftbigg
f(u+ ?2 ,v)?f(u? ?2 ,v)
bracketrightbigg
(10)
A discrete approximation of ??v is produced similarly. We can now proceed
further to produce discrete approximations on the gird of higher order terms:
?2
?u2 ?
1
?2 [f(u+?,v)?2f(u,v)+f(u??,v)], (11)
and similarly
?2
?v2 ?
1
?2 [f(u,v +?)?2f(u,v)+f(u,v ??)], (12)
By substituting Eqn. (11) and (12) into Eqn. (9) and simplifying slightly we
get:
1
?2 [f(u+?,v)+f(u??,v)f(u,v +?)+f(u,v ??)?4f(u,v)] = 0. (13)
The resulting schemes for 2-D and 3-D versions of 2nd, 4th and 6th order
equations are illustrated in Figure 3. Similar schemes were produced for 4 and
5 dimensional versions, but for obvious practical reasons they are not illustrated.
0 ????2? ...
Figure 2: At the edges of the grid ?imaginary? points are calculated to fill in for
the higher order schemes
We then apply this scheme to the entire grid, filling in the specified boundary
values at the boundaries. This gives a system of linear equations.
Ax = B. (14)
4
Notice that in the case of the higher order schemes near the edge of our grid
we will need to compute values for points outside of the domain itself; in order
to achieve the desired continuity right up to the boundary. This is because
the higher order schemes express a given point in terms of more than just its
immediate neighbours, yet we seek a solution to these equations right up to the
edge of the domain where our positional boundary conditions have been defined.
These points are often referred to as ghost points, illustrated in Fig. 2.
This scheme produces a system of equations of the form:
?
??
??
??
??
??
??
??
??
??
??
??
...
...
...
...
1 ... 1 ?4 1 ... 1
...
...
...
...
?
??
??
??
??
??
??
??
??
??
??
??
?
??
??
??
??
??
??
??
??
??
??
?
x1
x2
...
...
...
...
...
...
xn2?1
xn2
?
??
??
??
??
??
??
??
??
??
??
?
=
?
??
??
??
??
??
??
??
??
??
??
?
b1
b2
...
...
...
...
...
...
bn2?1
bn2
?
??
??
??
??
??
??
??
??
??
??
?
(15)
Having produced a suitable linear system of equations that approximates the
solution to our given PDE we now seek an efficient solution of that system. In
all but the most trivial of cases any direct solution (e.g. inversion, or via some
decomposition) is not feasible. Therefore we turn to an iterative approach, for
which the reader is referred to the texts [12, 3].
4 Results
For the purposes of testing, a simple example application, Figure. 4 was pro-
duced to demonstrate the two different solutions running side by side for the
biharmonic (Eqn. 2). In both cases the boundary conditions specified were
identical. The implementation used here utilises a simple multigrid solver, with
symmetric Gauss-Seidel relaxations performed at each stage.
Some example computation times with the schemes we produced are pre-
sented in Table 1. Note that in this table the 4th order 3-D scheme we use is
the same as in [6]. This implementation has not been particularly aggressively
optimised yet, and in practice significant performance gains can most probably
be made.
In Figure 5 the computational cost of increasing the dimension of the space in
which the surface is embedded is investigated. The time is clearly proportional
to the dimensionality of the space.
5 Conclusion
The original motivation for this work was that of modelling image manifolds.
The results presented here demonstrate the feasibility of this idea from a purely
5
(a) 2nd order 2-D scheme (b) 2nd order 3-D scheme
(c) 4th order 2-D scheme (d) 4th order 3-D scheme
(e) 6th order 2-D scheme
Figure 3: Visual illustrations of several of the schemes we used.
6
(a) Biharmonic example (analytic) (b) Biharmonic example (numeric)
Figure 4: Side by side visual comparison of the results of the two different
methods for a simple 2-D biharmonic. The numerical solution has a relatively
high maximum residual error (?max = 0.01), but only required a small number of
iterations. If we continued until a smaller maximum residual error is obtained,
then the numeric solution does converge to be the closer to the analytic one.
Method/scheme Computation time
per dimension
Max Residual error
(Iterations)
2nd Order 2-D numeric 0.010s 0.00854 (13)
2nd Order 3-D numeric 0.320s 0.00949 (15)
2nd Order 4-D numeric 29.250s 0.0118 (17)
4th Order 2-D numeric 0.020s 0.0098 (90)
4th Order 3-D numeric 1.180s 0.00998 (137)
4th Order 4-D numeric 197.430s 0.00985 (180)
6th Order 2-D numeric 0.090s 0.00999 (544)
6th Order 3-D numeric 7.030s 0.0118 (1000)
6th Order 4-D numeric 1034.530s 0.0442 (1000)
6th Order 6-D numeric dnf dnf
Table 1: Computational results from numerical schemes. In all cases a 20d grid
was used. The final scheme failed to run because of memory requirements.
computational point of view, and have opened a lot of potential avenues for
further work. Work is therefore now progressing on a more detailed study of
the application of the techniques developed here, along with several other surface
modelling techinques, to the specific problem of modelling image manifolds.
The multigrid solver we used for the results presented here was relatively
simplistic, and does not exploit the potential for speed increases by parallel
solving of the system of equations. In the future however with large image
manifolds it is anticipated that it will be necessary to use parallel solvers on
large clusters, and there is a great potential for further speed increases in other
areas.
Another interesting avenue for further work is the possibility of using the
multigrid solutions to form the basis of a dynamic level of detail mechanism
for the rendering of PDE surfaces analagous to current dynamic level of detail
techniques used when rendering triangle meshes. This could also potentially be
7
0
50
100
150
200
250
300
350
400
450
500
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 100000
Time (sec)
Dimension
Analytic Computation
Figure 5: Computational results from analytic solution. As we increase the di-
mension of the space the surface is embeded in the computational time increases
proportionally.
applied when rendering images using image manifolds.
Acknowledgements
The authors would like to thank the VVG network for the grant which made
this work possible.
References
[1] M. Bichsel and A. P. Pentland. Human face recognition and the face image
set?s topology. CVGIP: Image Underst., 59(2):254?261, 1994.
[2] M. I. G. Bloor and M. J. Wilson. Using partial differential equations to
generate free-form surfaces. Comput. Aided Des., 22(4):202?212, 1990.
[3] William L. Briggs, Van Emden Henson, and Steve F. McCormick. A Multi-
grid Tutorial. Other Titles in Applied Mathematics. SIAM, second edition,
2000.
[4] Ruel V. Churchill and James Ward Brown. Fourier Series and Boundary
Value Problems. McGraw-Hill, Inc., 1978.
[5] David L. Donoho and Carrie Grimes. Image manifolds which are isometric
to Euclidean space. Journal of Mathematical Imaging and Vision, 23(1):5?
24, 2005.
8
0
0.5
1
1.5
2
2.5
3
3.5
0 5 10 15 20 25 30 35 40
Maximum residual error
Number of iterations performed
1st Dimension
2nd Dimension
3rd Dimension
Figure 6: Maximum residual error vs. iterations with the numerical biharmonic
scheme as used in Figure 4(b).
[6] Haixia Du and Hong Qin. Free-form geometric modeling by integrating
parametric and implicit PDEs. IEEE Transactions on Visualization and
Computer Graphics, 13(3):549?561, 2007.
[7] Haw-minn Lu, Yeshaiahu Fainman, and Robert Hecht-Nielsen. Image man-
ifolds. In Nasser M. Nasrabadi and Aggelos K. Katsaggelos, editors, Appli-
cations of Artificial Neural Networks in Image Processing III, volume 3307,
pages 52?63. SPIE, 1998.
[8] S. Naya, S. Nene, and H. Murase. Subspace methods for robot vision. IEEE
Transactions on Robotics and Automation, 12:750?758, 1996.
[9] Les Piegl and Wayne Tiller. The NURBS book. Springer, 2nd edition, 1997.
[10] J.W. Thomas. Numerical Partial Differential Equations: Finite Difference
Methods, volume 22 of Texts in Applied Mathematics. Springer, second
edition, 1995.
[11] Hassan Ugail, Malcolm I. G. Bloor, and Michael J. Wilson. Techniques for
interactive design using the PDE method. ACM Trans. Graph., 18(2):195?
212, 1999.
[12] R. S. Varga. Matrix Iterative Analysis. Englewood Cliffs, NJ, USA, 1962.
9