Construction of Fractal Surfaces via Solutions of Partial Differential Equations

1. Introduction

Fractals are known to produce extremely complicated and impressive shapes, which many times resemble objects of the physical world. They have been used extensively in various scientific areas such as computer graphics, image compression and processing, biology and metallurgy. They usually emerge as attractors of Iterated Function Systems, a notion introduced by Barnsley and Demko in [1]. An Iterated Function System (IFS) is defined as a pair consisting of a complete metric space , together with a finite set of continuous, contractive mappings . IFSs are able to produce extremely complicated sets using only a handful of mappings. A more general definition is that of the Recurrent Iterated Function System (RIFS), which can also be used to construct fractal sets.

In this paper, we will deal with the construction of fractal surfaces, a subject firstly addressed by Massopust in [2] (see also [3]), using mainly the notion of fractal interpolation. Fractal Interpolation is an alternative to traditional interpolation techniques, introduced by Barnsley in [4], that gives a broader set of interpolants. Fractal Interpolation Functions are constructed as attractors of specific IFSs or RIFSs. Using this method, one can construct not only interpolants with nonintegral dimension but also smooth, nonpolynomial interpolants, or even splines and Hermite functions (see [5, 6]). Fractal interpolation functions are highly irregular and cannot be described using elementary functions such as polynomials (excluding the trivial cases where the fractal function is actually a polygonal line, a spline, or some other ordinary interpolants).

In the years that followed, the problem of the construction of fractal surfaces has been dealt by many other authors (see [7–13]), in the case where the interpolation points or the vertical scaling factors are confined in some manner. In some cases, the construction uses interpolation points, that are restricted to be collinear in the borders of , and in some others, it uses maps with equal vertical scaling factors. Recently, in [14], a general solution for arbitrary data was given. This construction uses the values of the interpolation points, as all the previous attempts did, as well as the values at the borders of each specified region (which are chosen arbitrarily). The fractal surface emerges as the attractor of an appropriately chosen RIFS. The construction presented here uses the results of [14] and extends a new method, which uses the solutions of certain Partial Differential Equations, that are then used to construct fractal interpolation surfaces on arbitrary data, placed on rectangular grids.

2. Mathematical Background

2.1. Iterated Function Systems

As mentioned above, an Iterated Function System   is defined as a pair consisting of a complete metric space , together with a finite set of continuous, contractive mappings , with respective contraction factors , for (). The attractor of an IFS is the unique set , for which for every starting compact set , where (1) and is the complete metric space of all nonempty compact subsets of with respect to the Hausdorff metric (for the definition of the Hausdorff metric and properties of see [15]). A simple example of an IFS defined on is the one that produces the well-known Sierpinski's Triangle (see Figure 1(a)), which consists of the three mappings: (2) The attractor of the following IFS looks like a natural fern (see Figure 1(b)):

Figure 1: Two known attractors that arise from IFS. (a) Sierpinski's Triangle, (b) a fern.

(3)

There are two known algorithms, that are used to construct the attractor of an IFS, the Deterministic Algorithm (DA) and the Random Iteration Algorithm (RIA). The first starts with an arbitrary compact set and produces the sets , and so forth. The sequence obviously converges to the attractor of the IFS. The RIA, on the other hand, starts with an arbitrary point and selects randomly (using a set of probabilities that sum to 1) one of the mappings of the IFS, to produce a new point (where is the selected map). The operation is continued indefinitely. The produced points trace out the attractor. While the first algorithm constructs the attractor itself, it uses a great amount of memory, especially if one applies many iterations. Thus, many times, the RIA is preferred. However, one should choose carefully the set of probabilities; otherwise some sections of the attractor will not be traced out as clear as others. In fact, the set of probabilities determine a unique “invariant measure” that is supported on (the attractor of the IFS). If each pixel in a computer screen is colored according to the visitation frequency of the random iterates of RIA, then one obtains a pictorial representation of the measure . The interested reader may find more about this subject (including some interesting photos) in [15] (mainly chapter 9) and [16]. More about the DA and RIA can also be found in [15]. In Figure 1, the attractor of the Sierpinski's triangle (a) was constructed using DA, while the fern-like attractor (b) was constructed using RIA with a suitable selection for the probabilities (in particular, ). We should, also, note that in both algorithms it is preferred to start with points which we know they belong in the attractor. Thus, we will always get subsets of the attractor itself, instead of points that lie near it.

2.2. Recurrent Iterated Function Systems

A more general concept, that allows the construction of even more complicated sets, is that of the Recurrent Iterated Function System, or RIFS for short, which consists of an IFS , together with an irreducible row-stochastic matrix  , such that (4) The recurrent structure is given by the (irreducible) connection matrix   , which is defined by (5) where . The transition probability for a certain discrete time Markov process is , which gives the probability of transfer into state , given that the process is in state . Condition (4) says that whichever state the system is in (say ), a set of probabilities is available that sum to one and describe the possible states to which the system transits at the next step.

In this case, the construction of the contractive map needs a little more effort. First, we define the mappings (6) for all and the metric space (7) equipped with the metric (8) Then is a complete metric space. The map is now defined by (9) where , for . If are contractions, then is a contraction and there is such that and for . The set is called the attractor of the RIFS .

Next, we construct a sequence of sets that converge to the attractor. Let . We define the sequences in and in as follows: , and , for , where . Evidently (10) We emphasize that the attractor of an RIFS depends not only on the corresponding IFS but also on the stochastic matrix. For example, the following IFS equipped with different stochastic matrices produces different attractors (see Figure 2): (11)(12)

Figure 2: The attractors of the IFS defined by (11)-(12), equipped with various stochastic matrices. Note that the attractor shown in (a) is a solid rectangle. It is the unequal weights of P and the finite number of iterations of the RIA algorithm that give the resulting figure its nonuniform appearance. The 0 entries in the matrices for (b) and (c) give those sets a nonintegral dimension less than 2.

Both DA and RIA can be modified to construct the attractor of any RIFS. The modification of DA is described above, where the sequences are been defined. We note that in many cases, the functions of the RIFS are defined on a compact subset of , which we call , instead of the space itself, for . In this case the initial vector is defined as , where . The modification of the RIA is much simpler. We start with an initial point and an initial state . Then we select arbitrarily, using the set of probabilities , one of the mappings of the RIFS and apply it to to produce a new point . Suppose that is the selected map, then and the system transits to state . We continue this operation until enough points are produced. We note that if the connection matrix of the RIFS is reducible, this procedure will not work. The produced graph will depend on the selection of the initial state. Therefore, for RIFS with reducible connection matrices the following modification of the RIA is needed: we start with initial points , such that , for . For each point , we select randomly one of the maps of the RIFS, using the set of probabilities , and apply it to to produce the new point , where indicates that the selected map is . After the completion of this step, the new set of points is . Again, for each point , we select randomly one of the maps of the RIFS, using the set of probabilities , and apply it to to produce the new point , where indicates that the selected map in the second step is . The procedure continues indefinitely. The produced points trace out the attractor of the corresponding RIFS. Finally, we should note that, as was the case in the IFS, the probabilities of the RIFS determine a unique invariant measure , which can be visualized in a computer screen (Figure 2(a) is a close approximation of this procedure). It is possible (see Figures 2(b) and 2(c)) that there will be “holes” in the support of (or equivalently in the graph of the attractor ), if some of the transition probabilities are zero.

3. Fractal Interpolation Surfaces

In this section, we present a general construction of fractal surfaces, that interpolate points placed on a rectangular grid. The fractal surfaces emerge as attractors of specific RIFS. Consider a data set: (13)

such that ,  ,  , where , which contains in total points. Consider, also, a data set :

(14)

such that ,  , which contains in total points. To simplify the notation we set (15) The points of divide into regions: (16) for all , while the points of divide into domains: (17) for all . We make the additional assumption that for every there is at least one interpolation point, that lies in the interior of . Furthermore, we define a labeling map (18) and the 1-1 functions (an enumeration of the sets and , respectively) such that

(19) We define the stochastic matrix by (20) where is the number of non zero elements of its th row. Consequently, the connection matrix is defined as in Section 2.2 and the connection vector as follows: In Figure 3, we show an example of the association between the regions and the domains, which is implied by the labeling map , that defines the stochastic matrix, the connection matrix and the connection vector of the RIFS. The connection matrix and the connection vector, that corresponds to the shown association, are given as follows: (21)

Figure 3: The interpolation points divide in 16 regions and in 4 domains. The association between regions and domains is shown with the arrows. For example, ,  ,  ,  , and so forth. Thus, the connection vector is .

To complete the definition of the RIFS, we consider mappings of the form (22) with (23)

for all , where is a Lipschitz continuous function on and , for all . The parameters are called vertical scaling factors. It is easy to show that there exists a metric (equivalent with the Euclidean metric) such that is a contraction with respect to the metric for all . We confine the map , so that it maps the interpolation points that lie on the vertices of to the interpolation points that lie on the vertices of . Hence, we obtain the following relations: (24) for all ,  , and  . Using the first four (12) equations, we can compute ,  ,  ,  , in terms of the interpolation points and the vertical scaling factors, for all .

Consequently, the RIFS has a unique attractor . In general is a compact subset of containing the points of . The following Proposition gives conditions so that is the graph of a continuous function . These conditions involve points that lie on , for all , (where is the boundary of ).

Proposition 1. Let be a Lipschitz function, that interpolates the points of (i.e., , ). If the RIFS defined above satisfies the conditions (25) where , for all , , then its attractor is the graph of a continuous function that interpolates the data points.

The proof of the above proposition, together with several methods that one may use to select suitable functions , so that the RIFS satisfies the conditions (25), can be found in [14]. Here, we are interested in the special case described in the following corollary.

Corollary 1. Let be a Lipschitz function that interpolates the points of . Consider the case where (26) where and are Lipschitz functions defined on , such that (27) The unique attractor of the corresponding RIFS is the graph of a continuous function that interpolates the points of .

The functions and mentioned in the corollary may be considered as some kind of building blocks for the constructed FIS. It is straightforward to see that in the trivial case where , for all , the resulting attractor will be the graph of . Otherwise, the attractor can be obtained as a sum of and multiplied by , and so forth. Examples of such constructions for the 1-dimensional case may be found in [15, page 218].

3.1. A Construction Using Partial Differential Equations

It is well known that in the case of the construction of affine FIFs of one variable, the attractor depends only on the choice of the interpolation points (and of course on the choice of the vertical scaling factors). It was shown in [14] that a similar property holds for multivariate FIFs of certain type (see Proposition 1). For example, in the case of fractal interpolation surfaces, one needs to know the interpolation points and the desired values of the constructed function on the borders of the interpolation grid. If all this information is available, then a unique FIS is determined. This fact closely resembles the case in many PDE-related problems with boundary conditions, that arise in physics and engineering (e.g., Laplace's equation with boundary conditions, Plateau's problem etc.). In fact, it is not hard to see that the problem of the well known self-affine fractal interpolation (see [15]) of one variable can be recast as follows (in the context of Corollary 1 for functions of one variable):

(i) find a function that solves the Laplace's equation with boundary conditions ,  , for , (ii) find a function that solves the Laplace's equation with boundary conditions ,  ,

where are the interpolation points. Similarly, the problem of bivariate fractal interpolation on collinear interpolation points on the boundary's grid (see [10]) can be recast as the problem of finding suitable and functions that solve the two-dimensional Laplace's equation on linear boundaries. Here we show that both cases are examples of a more general procedure.

Consider the set , which consists of continuous functions that are defined only on , interpolate the points of , and have continuous partial derivatives up to th order in . Let be the set containing all the subsets of of the form and a th order Partial Differential Equation (PDE) defined on the set , such that its solution satisfies (28) where . The problem of finding a function, which solves a specified partial differential equation (PDE) in the interior of a given region (), that takes prescribed values on the boundary of that region, is called a Dirichlet's problem. Assuming that the PDE is uniquely solvable for any and any and that the solution is a Lipschitz function, we consider the operator (29) that assigns any function defined on to the solution of the corresponding PDE with boundary conditions as in (28).

Now, consider a function . This function may consists of splines, Hermitian functions, or any other types of 1-dimensional interpolants, that have the desired properties. We will construct an RIFS, whose attractor will be the graph of a continuous function. As mentioned above, we study the case where , for all , as in Corollary 1. We compute and  as follows: (30) for all . In this case, the conditions of Corollary 1 are satisfied for the RIFS , for any selection of the stochastic matrix . Thence, the unique attractor of the corresponding RIFS is the graph of a continuous function , that interpolates the points of .

We summarize the procedure of the construction.

(1) We take (as input) a set of interpolation points placed on a rectangular grid, a set of vertical scaling factors, and a specific partial differential equation (see Figure 4(a)). (2) We select a set (a subset of ) of points that are also placed on a rectangular grid. The points of form the regions , for all , while the points of form the domains . The points of are chosen so that inside each domain lies at least one region. Additionally, we select a function that associates each region to a specific domain. From this function we compute the connection vector . (3) We construct arbitrary 1-dimensional functions that interpolate the points of . This is the function (see Figure 4(b)). (4) We solve (numerically) the partial differential equation inside the region , so that its solution satisfies the boundary conditions , for all (Dirichlet's problem). Hence, we obtain the function (see Figure 4(c)). (5) We solve (numerically) the partial differential equation inside , so that its solution satisfies the boundary conditions , for all . Hence, we obtain the function (see Figure 4(d)). (6) Finally, we construct the attractor of the corresponding RIFS (see Figure 4(e)).

Figure 4: An example of the construction. Here , the regions have dimensions and the domains . We used Laplace's partial differential equation () with Dirichlet's boundary conditions. (a) The set of interpolation points. (b) The function . (c) The function . (d) The function . (e) The fractal interpolation surface.

4. The Algorithm

In Section 2.2, we described how RIA and DA can be modified to construct the attractor of any RIFS. As any fractal interpolation surface is the attractor of a specific RIFS, we should expect that those algorithms may be used to construct it. This is true; however in this case, the projections of the produced points on will form an nonuniform grid. This means that there will be “dense” and “sparse” areas, with a lot of produced points, and very few points respectively. Here, we describe an algorithm that solves this problem. For simplicity, we assume that the regions and the domains are squares with dimensions and , respectively. We choose and so that , where . The algorithm produces points ( is the number of iterations), so that their projections on form a square of uniformly distributed points. An example of the operation of the algorithm is given in Figure 5. We note that the number must be a multiple of (to be more specific: ).

Figure 5: An example of the algorithm, where and . In this case, there are four regions and one domain. Thus, the RIFS is actually an IFS. The connection vector is . The algorithm is expected to produce 81 points for the attractor. In all the figures the squares represent the projection of the points on . (a) The initial interpolation points are shown with black. The points of are the four corners. (b) The points produced after one iteration (c) The points produced after 2 iterations.

(1) INPUT: The interpolation points, the number , that defines which points belong in and which points belong in , a set of vertical scaling factors, the connection vector , the number of iterations , and the PDE. (2) The initial set is set to contain all the interpolation points ( points). (3) The projections of the interpolation points on form horizontal lines and vertical lines. For each one of the horizontal lines we compute points of a function that interpolates the corresponding points. Similarly, for each one of the vertical lines we compute points of a function that interpolates the corresponding points. (4) Form the function H.For to DO: (a)For to DO: (i)We solve numerically, inside the region , the partial differential equation with Dirichlet's boundary conditions defined by the computed points of the interpolation functions. (5) Form the function B. For to DO: (a)For to DO: (i)We solve numerically, inside the domain , the partial differential equation with Dirichlet's boundary conditions defined by the computed points of the interpolation functions. (6) Form the attractor of the RIFS. For to n DO (a)For to DO: (i)For to DO: (A)We deal with the points of the region . First of all, we find the corresponding domain , using the connection vector: . (B)As the parameters of the map are computed by the functions and , we may use to map all the points of G, that lie inside . Therefore, we obtain some new points inside the region . (b)The new set G contains all the points that were produced by the previous iterations ( points). (7) OUTPUT: After the completion of the iterations, the final set contains points that belong to the attractor of the RIFS.

As noted, the output of the algorithm is a set of points, that their projections on form a square of uniformly distributed points, placed on an orthogonal grid. We can use these points to produce a 3D surface using any 3D graphics suit. The produced surface depends greatly on the selection of the partial differential equation (see Figure 6). In fact, one can use different PDEs in each region/domain. In Figure 7, we give some examples of the application of the algorithm on arbitrary data. We note that in most of the examples we used the Laplace's PDE. For future research, it would be interesting to investigate more thoroughly how the choice of the PDE affects the resulting fractal surface.

Figure 6: The two fractal surfaces (shown above) were constructed using the same parameters, but different partial differential equations. In (a) Laplace's equation was used (), while in (b) the corresponding PDE is: .

Figure 7: Some examples of the proposed construction. At the left side, the interpolation points and the chosen borders are shown. At the right side, one can see the constructed fractal surfaces. We used Laplace's PDE in both examples.

Acknowledgments

The author would like to thank Professor M. Drakopoulos (Mathematics Department of the University of Athens) for his valuable help in producing the algorithms that solve Dirichlet's problem for Lapalce's PDE. This reserch is partially supported by the Special Account for Research Grants of the University of Athens no. 70/4/5626.