ACommon approach to simulating fractal surfaces is viathe sample paths of fractional Brownian motions and their multidimensional extensions to (e.g. [1], [2]). These random fields are self-similar in distribution in that when sampled at various scales the distribution of the sample is the same up to a constant scaling factor. However, they are indexed by . Suppose instead one wanted to simulate a fractal surface with some more complex geometry or topology, e.g., a fractal cylinder. The analogous approach to simulating such an object would require a random field indexed by a manifold, such as a surface in , that possessed the same properties as fractional Brownian motion, i.e., self-similarity, almost surely H¨older continuous sample paths, and stationary increments.

Constructing such fields has been a subject of ongoing research for a number of years and until recently, extensions were known only for limited classes of surfaces, such as spheres or hyperboloids, and only for certain ranges of self-similarity, i.e., only with Hurst index . Then in [3] self-similar Gaussian random fields were constructed over wide classes of surfaces, including arbitrary compact manifolds, both with and without boundary. These extensions were shown to have H¨older continuous sample paths, be self-similar in distribution, and have stationary increments, thus possessing all the basic properties one would expect of an extension of fractional Brownian motion. In this paper we describe the simulation of these fields for surfaces in .

The simulation of Gaussian random fields presents non-trivial challenges, in particular if one considers index sets other than . One common approach is the use of Fourier or wavelet techniques (e.g. [1], ch. 2, [2]), but if one is working on a general surface the Fourier transform is no longer available. On the other hand, one can view the Fourier transform as the spectral decomposition of the Laplacian, the sine and cosine functions being the eigenfunctions. These we do have on a surface, and so the analogous approach is to build and simulate random fields over manifolds and surfaces using Fourier series of eigenfunctions of the Laplace-Beltrami operator. One novel feature of this approach is that using the Dirichlet Laplacian for a surface with boundary allows us to produce fields with almost sure boundary values. We begin the next section with a derivation in the simple case of the interval [0, 1], before moving on to surfaces.

A. A Basic Example

As motivation and a base case for our construction, consider Brownian motion on [0, 1]. Being a Gaussian process, is determined by its covariance function,

Now notice that is the Green’s function of the Laplacian on , i.e., for

Of course we must specify boundary conditions to have a well defined Laplacian, and we see that these come from : If

then F(0) = 0 and . Thus we can say that is the unique (up to equality in distribution) Gaussian process on [0, 1] whose covariance is the Green’s function of acting on smooth functions such that 0.

then we can write (1) as

and in this way we see that starting from , we uniquely determine as the inverse of the integral operator K defined by the covariance of . If and are the eigenvalues and eigenfunctions respectively of , with the above boundary conditions, then a calculation shows

The Spectral Theorem and the functional calculus associated with it then yield

Fig. 1. fractional Gaussian bridge with from top to bottom.

and

denoting the inner product,. If we now write

it is easily seen that the series defining k(x, y) converges absolutely and uniformly on [0, 1] and moreover it must be that . Then if are i.i.d. standard normal random variables, an application of Fubini’s Theorem yields

and we thus arrive at the well known Fourier series expansion of ,

Now let W be the white noise (or isonormal process, cf [?]) on and denote its action a function by

Then is a set of i.i.d. standard normal random variables, denoted again by , and if

then

and

Thus we can write and arrive at the stochastic steady-state equation

Notice that working backwards we can define to be the unique Gaussian process satisfying (4). This is similar, but not the same as, the classical equation

To see that (4) is in fact different, note that is self adjoint on , while d/dt is not (integrate by parts).

Now suppose we replace the mixed Laplacian above by the Dirichlet Laplacian in (4), acting on functions such that f(0) = f(1) = 0. The Dirichlet Laplacian on [0, 1] has eigenvalues and eigenfunctions given by

respectively. Then we have as above a process defined by

If we now let (the case above being ) we can extend (4) to

and write

which is the homogeneous Riesz field constructed in [3]. This is a self-similar Gaussian process with almost surely Holder continuous sample paths and such that almost surely, i.e., it is a Gaussian bridge.

Because the series converges in with probability 1, we can simulate these processes by taking partial sums (see fig. 1). Notice that as increases the sample paths become more regular. The reason can be seen from (6): Lower alpha emphasizes the larger eigenvalues and thus the higher frequency wave functions show more of a contribution. Similarly a larger suppresses their contribution and the lower frequency wave functions dominate, leading to less overall oscillation and a smoother sample path. We will see this same behavior in general below.

Fig. 2. Cylinders with Dirichlet boundary conditions for

B. Extension to Surfaces with Boundary

Suppose now D is a compact connected surface in

with smooth boundary. Then there is a canonical differential operator on corresponding to the Euclidean Laplace operator above, the Laplace-Beltrami operator, which we refer to also as simple the Laplacian of D (see e.g. [4]). As above we let and be the eigenvalues and eigenfunctions of the Dirichlet Laplacian of . Given a white noise W on D we start from

to obtain

where the change in the exponent from 1 to 2 ensures sample path regularity (see [3]). In this way we obtain a random field on D that is self-similar, almost-surely H¨older continuous of order for any and almost surely not H¨older continuous for , and has stationary increments. A few words are in order as to what we mean by self-similar: Because we are considering surfaces embedded in , that is, surfaces with Riemannian metrics induced by the Euclidean metric on , the usual definitions carry over, i.e., . For more general considerations of self-similarity on arbitrary manifolds see [3].

Now we have a generalization of the bridges above that possesses the properties we would like, however for a general surface D we no longer have simple analytical expressions for the functions . One way around this is to discretize and construct the analogous discretized version of (7). This requires some care, as there are a number of discretizations available for Laplace-Beltrami operators, some of which behave very differently (e.g. [5]).

Given a meshed surface, most discrete Laplacians are weighted graph Laplacians,

for some edge weights and where means x and y are adjacent in the mesh (i.e. they share an edge). Supposing we have a fairly uniform mesh, let us set . Denote this discrete Laplacian by L, and let its eigenvectors and eigenvalues be and respectively. Note that where the dimension of L (L is a matrix) is , and is a real valued function on our mesh. We then let

where we choose some and are again independent standard normals.

We would like to be able to take our mesh fine enough and large enough so that provides a good approximation of the full Riesz field . How good this approximation is, and in what sense, will depend completely on how well L approximates . For general surfaces, discrete Laplacians and their convergence to the continuum Laplacian is a very active area of current research, and we defer a more full discussion for the moment.

Some examples of planar domains are plotted in Figs. 3, 5, and 6. Note that with probability 1 the field is zero on the boundary. As an example of what is possible in three dimensions, consider the cylinder (Fig. 2). If we impose zero Dirichlet boundary conditions on the two boundary circles, we get analogues of the bridges in Fig. 1. Notice that again, controls the sample path (or surface) regularity and the fields are almost surely equal to zero on both boundary circles, i.e., the height of the field values plotted along normals is zero.

C. Surfaces Without Boundary

For surfaces without boundary there are no obvious boundary conditions, so we need to interpret our equations differ-

Fig. 3. Disk with 3 smaller disks removed,

ently. The basic issue is that on a compact manifold, e.g., the sphere , the lowest eigenvalue of the Laplacian is zero and is not defined. We can instead consider the series

where o is some fixed point serving as an “origin.” This series does converge and again defines a self-similar field with stationary increments and almost sure H¨older continuity as above, called the Riesz Field in [3]. Figure 4 shows spheres with increasing alpha, using the same discrete Laplacian L as before.

A. Convergence

For uniform or nearly uniform meshes, the above discretization produces good results. In general however, placing a uniform mesh on a surface is non-trivial (depending on the curvature and boundary), and in this case the above discrete Laplacian may be a poor approximation to the continuous Laplacian. In the case of highly curved surfaces, care must be taken in choosing an appropriate discretization.

Suppose we have a meshed surface and have chosen a discrete Laplacian, L. A simple condition for the almost sure approximation is as follows: Suppose that for any integer N > 0, for all we have

Fig. 4. Spheres with

where we can embed in denoting surface measure on D. Then for any we can choose such that with probability 1

This is a simple consequence of the fact that

converges in almost surely, for then we just bound the tail and the assumed convergence of the lower eigenfunctions and eigenvalues yields the claim.

Thus if one has a meshed surface, the Laplacian chosen should satisfy the above spectral convergence condition. Such conditions are known to hold for certain Laplacians, but not all [6], [7]. Of course, this is an active area of research and so we expect further results on spectral approximation of the Laplace-Beltrami operator to become applicable to the spectral synthesis we have used here.

B. Efficiency and Self-Similarity

There are two parameters which determine the computation time in simulating for a given surface: The number n of grid points in the mesh (i.e. how fine the mesh is) and the number N of eigenvectors appearing in

The dimension of the discrete Laplacian will be and in general N should be taken much smaller than , for Weyl’s asymptotic formula for the eigenvalues of (e.g. [4]) tells us that for a two dimensional surface as we are dealing with here. Thus

and so if our discrete Laplacian L satisfies the spectral convergence above, we see that the contribution, in mean square, of the Nth term is of the order which may rapidly become negligible depending on the fineness of the mesh.

On the other hand, for finer meshes, the computation of even the first few hundred eigenvectors and eigenvalues of the matrix L can become long. However, once computed, this spectral data can be stored for repeated use. All that is needed is the random numbers , which are not expensive. Moreover, the self-similarity of the Riesz fields allows one to compute the spectra of a surface of any size, and simply rescale the resulting field without computing the spectra for a rescaled surface. For example, once one has simulated the Riesz field on a sphere of radius 1, one can simulate it on a sphere of radius c simply by plotting for the appropriate value of . In this way, for very detailed images, the computational time required may be somewhat high, but it is a one time cost in the above sense.

We have restricted our attention to embedded surfaces in , however the theory in [3] applies to non-embedded manifolds, e.g., hyperbolic space or flat tori, as well as n-dimensional surfaces in . In this case the approximation takes the form

Of course working with a surface whose metric is not induced by the ambient Euclidean metric increases the difficulty in

Fig. 5. Fractal surface on the unit disk with Dirichlet boundary conditions () with negative values set to zero.

Fig. 6. Fractal surface over a teardrop shaped region with

computing accurate discretizations of the Laplacian. Also, we have here chosen to be fixed, but one could also allow to vary over the surface and thus obtain multifractal fields. We expect, similar to what is known on , that such fields with varying fractal index would still be continuous.

The first author thanks his advisor, Harold Parks, for his time and encouragement. Figure 4 was produced using the Matlab Toolboxes “Toolbox Graph” and “Toolbox Wavelets on Meshes,” written by Gabriel Pyre.

[1] M. F. Barnsley, R. L. Devaney, B. B. Mandelbrot, H.-O. Peitgen, D. Saupe, and R. F. Voss, The science of fractal images. New York: Springer-Verlag, 1988, with contributions by Yuval Fisher and Michael McGuire. [Online]. Available: http://dx.doi.org/10.1007/ 978-1-4612-3784-6

[2] G. Franceschetti and D. Riccio, Scattering, Natural Surfaces, and Fractals. Elsevier Science, 2006. [Online]. Available: http://books. google.com/books?id=v7pxQiUv9BgC

[3] Z. Gelbaum, “Fractional brownian fields over manifolds,” Trans. Amer. Math. Soc., accepted for publication, 2013.

[4] I. Chavel, Eigenvalues in Riemannian geometry, ser. Pure and Applied Mathematics. Orlando, FL: Academic Press Inc., 1984, vol. 115, including a chapter by Burton Randol, With an appendix by Jozef Dodziuk.

[5] M. Wardetsky, S. Mathur, F. Kalberer, and E. Grinspun, “Discrete laplace operators: No free lunch,” in Proc. Symp. Geom. Proc. 2007, 2007.

[6] T. K. Dey, P. Ranjan, and Y. Wang, “Convergence, stability, and discrete approximation of Laplace spectra,” in Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms. Philadelphia, PA: SIAM, 2010, pp. 650–663.

[7] M. Reuter, S. Biasotti, D. Giorgi, G. Patane, and M. Spagnuolo, “Discrete Laplace-Beltrami operators for shape analysis and segmentation,” COMPUTERS & GRAPHICS-UK, vol. 33, no. 3, SI, pp. 381–390, JUN 2009, IEEE International Conference on Shape Modeling and Applications, Tsinghua Univ, Beijing, PEOPLES R CHINA, JUN 26-28, 2009.