Construction of Fractal Surfaces via Solutions of Partial Differential Equations

A new construction of fractal interpolation surfaces, using solutions of partial differential equations, is presented. We consider a set of interpolation points placed on a rectangular grid and a specific PDE, such that its Dirichlet’s problem is uniquely solvable inside any given orthogonal region. We solve the PDE, using numerical methods, for a number of regions, to construct two functions H and B, which are then used to produce the fractal surface, as the attractor of an appropriately chosen recurrent iterated function system.


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 (X, ρ), together with a finite set of continuous, contractive mappings w i : X → X, i = 1, 2, . . ., N (N ≥ 2).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][8][9][10][11][12][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 I = [0, 1] 2 , 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.

Iterated Function Systems.
As mentioned above, an Iterated Function System {X; w 1−N } is defined as a pair consisting of a complete metric space (X, ρ), together with a finite set of continuous, contractive mappings w i : X → X, with respective contraction factors s i , for i = 1, 2, . . ., N (N ≥ 2).The attractor of an IFS is the unique set E, for which E = lim k → ∞ W k (A 0 ) for every starting compact set A 0 , where and H (X) is the complete metric space of all nonempty compact subsets of X with respect to the Hausdorff metric h (for the definition of the Hausdorff metric and properties of H (X), h see [15]).A simple example of an IFS defined on R 2 is the one that produces the well-known Sierpinski's Triangle (see Figure 1(a)), which consists of the three mappings: ( The attractor of the following IFS looks like a natural fern (see Figure 1 ( 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 A 0 and produces the sets , and so forth.The sequence A 0 , A 1 , A 2 , . . .obviously converges to the attractor of the IFS.The RIA, on the other hand, starts with an arbitrary point x 0 and selects randomly (using a set of probabilities that sum to 1) one of the mappings of the IFS, to produce a new point x 1 = w i1 (x 0 ) (where w i1 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 A (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, p 1 = 0.01, p 2 = 0.69, p 3 = 0.15, p 4 = 0.15).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.

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 {X; w 1−N }, together with an irreducible row-stochastic matrix The recurrent structure is given by the (irreducible) connec- , which is defined by where ν, μ = 1, 2, . . ., N. The transition probability for a certain discrete time Markov process is p ν,μ , 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 W needs a little more effort.First, we define the mappings for all A ∈ H (X) and the metric space equipped with the metric Then H (X), h is a complete metric space.The map W is now defined by where Next, we construct a sequence of sets that converge to the attractor.Let A ∈ H (X). We define the sequences {A n } n∈N in H (X) and {A n } n∈N in H (X) as follows: 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): Both DA and RIA can be modified to construct the attractor of any RIFS.The modification of DA is described above, where the sequences {A n } n∈N are been defined.We note that in many cases, the functions w i of the RIFS are defined on a compact subset of X, which we call X i , instead of the space X itself, for i = 1, 2, . . .N. In this case the initial vector A 0 is defined as A 0 = (B 1 , B 2 , . . ., B N ) t , where B i ⊆ w i (X i ), i = 1, 2, . . .N. The modification of the RIA is much simpler.We start with an initial point x 0 and an initial state i 0 .Then we select arbitrarily, using the set of probabilities p i0,1 , p i0,2 , . . ., p i0,N , one of the mappings of the RIFS and apply it to x 0 to produce a new point x 1 .Suppose that w i1 is the selected map, then x 1 = w i1 (x 0 ) and the system transits to state i 1 .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 N initial points x 0,1 , x 0,2 , . . ., x 0,N , such that x 0,i ∈ w i (X i ), for i = 1, 2, . . ., N. For each point x 0,i , i = 1, 2, . . ., N, we select randomly one of the maps of the RIFS, using the set of probabilities {p i,1 , p i,2 , . . ., p i,N }, and apply it to x 0,i to produce the new point x 1,i , where j 1,i ∈ {1, 2, . . ., N} indicates that the selected map is w j1,i .After the completion of this step, the new set of points is {x Again, for each point x ( j1,i) 1,i , i = 1, 2, . . ., N, we select randomly one of the maps of the RIFS, using the set of probabilities {p j1,i,1 , p j1,i,2 , . . ., p j1,i,N }, and apply it to x 2,i , where j 2,i ∈ {1, 2, . . ., N} indicates that the selected map in the second step is w j2,i .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 A), if some of the transition probabilities are zero.

Fractal Interpolation Surfaces
In this section, we present a general construction of fractal surfaces, that interpolate points placed on a rectangular grid.
To complete the definition of the RIFS, we consider N 1 • N 2 mappings of the form with for all (x, y) ∈ I i, j , where Q i, j is a Lipschitz continuous function on [0, 1] n and s i, j ∈ (−1, 1), for all (i, j) ∈ C 1 .The parameters s i, j , (i, j) ∈ C 1 are called vertical scaling factors.
It is easy to show that there exists a metric ρ θ (equivalent with the Euclidean metric) such that w i, j is a contraction with respect to the metric h for all (i, j) ∈ C 1 .We confine the map w i, j , so that it maps the interpolation points that lie on the vertices of J J(i, j) to the interpolation points that lie on the vertices of I i, j .Hence, we obtain the following relations: for all δ 1 , δ 2 ∈ {0, 1}, (k, l) = J(i, j), and (i, j) ∈ C 1 .Using the first four ( 12) equations, we can compute a 1,i , a 2, j , b 1,i , b 2, j , in terms of the interpolation points and the vertical scaling factors, for all (i, j) ∈ C 1 .
Consequently, the RIFS {[0, 1] 2 × R, w i, j , P; (i, j) ∈ C 1 } has a unique attractor G.In general G is a compact subset of R 3 containing the points of Δ.The following Proposition gives conditions so that G is the graph of a continuous function f .These conditions involve points that lie on ∂I i, j × R, for all (i, j) ∈ C 1 , (where ∂I i, j is the boundary of I i, j ).

Proposition 1.
Let h ∈ C( (i, j)∈C1 ∂I i, j ) be a Lipschitz function, that interpolates the points of Δ (i.e., h(x i , y j ) = z i, j , (i, j) ∈ C 1 ).If the RIFS defined above satisfies the conditions where (k, l) = J(i, j), for all (x, y) ∈ J J(i, j) , (i, j) ∈ C 1 , then its attractor G is the graph of a continuous function f that interpolates the data points.
The proof of the above proposition, together with several methods that one may use to select suitable functions Q i, j , 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 h ∈ C( (i, j)∈C1 ∂I i, j ) be a Lipschitz function that interpolates the points of Δ.Consider the case where where H and B are Lipschitz functions defined on [0, 1] 2 , such that The unique attractor The functions H and B 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 s i, j = 0, for all i, j, the resulting attractor will be the graph of H. Otherwise, the attractor can be obtained as a sum of H and H − B multiplied by s i , s 2 i , and so forth.Examples of such constructions for the 1-dimensional case may be found in [15, page 218].

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): (ii) find a function B that solves the Laplace's equation with boundary conditions where {(x 0 , y 0 ), . . ., (x N , y N )} 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 H and B 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 C k Δ , which consists of continuous functions that are defined only on (i, j)∈C1 ∂I i, j , interpolate the points of Δ, and have continuous partial derivatives up to kth order in (i, j)∈C1 ∂I i, j .Let R([0, 1] 2 ) be the set containing all the subsets of [0, 1] 2 of the form [a 1 , b 1 ] × [a 2 , b 2 ] and a kth order Partial Differential Equation (PDE) defined on the set R ∈ R([0, 1] 2 ), such that its solution g satisfies where v ∈ C k (∂R).The problem of finding a function, which solves a specified partial differential equation (PDE) in the interior of a given region (I i, j ), 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 R ∈ R([0, 1] 2 ) and any v ∈ C k (∂R) and that the solution is a Lipschitz function, we consider the operator that assigns any function v defined on ∂R to the solution g of the corresponding PDE with boundary conditions as in (28).Now, consider a function h ∈ C k Δ .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 Q i, j (x, y) = H(T 1,i (x), T 2, j (y)) − s i, j • B(x, y), for all (x, y) ∈ J J(i, j) , (i, j) ∈ C 1 , as in Corollary 1.We compute H and B as follows: H| Ii,j = P Ii,j h| ∂Ii,j , for all (i, j) ∈ C 1 .In this case, the conditions of Corollary 1 are satisfied for the RIFS {[0, 1] 2 × R, w i, j , P; (i, j) ∈ C 1 }, for any selection of the stochastic matrix P. Thence, the unique attractor G of the corresponding RIFS is the graph of a continuous function f , 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 N 1 • N 2 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 N 1 • N 2 regions I i, j , for all (i, j) ∈ C 1 , while the points of Δ form the M 1 • M 2 domains J k,l .The points of Δ are chosen so that inside each domain lies at least one region.Additionally, we select a function J that associates each region to a specific domain.From this function we compute the connection vector V .
(4) We solve (numerically) the partial differential equation inside the region I i j , so that its solution g satisfies the boundary conditions g| ∂Iij = h| ∂Iij , for all (i, j) ∈ C 1 (Dirichlet's problem).Hence, we obtain the function H (see Figure 4(c)).
(5) We solve (numerically) the partial differential equation inside J kl , so that its solution g satisfies the boundary conditions g| ∂Jkl = h| ∂Jkl , for all (k, l) ∈ D 1 .Hence, we obtain the function B (see Figure 4(d)).

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 R 2 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 α ∈ N. The algorithm produces ((N 1 − 1) • α n + 1) × ((N 2 − 1) • α n + 1) points (n is the number of iterations), so that their projections on R 2 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 N 1 − 1 must be a multiple of α (to be more specific: N 1 − 1 = α • (M 1 − 1)).
(1) INPUT: The N 1 × N 2 interpolation points, the number α, that defines which points belong in Δ and which points belong in Δ, a set of N 1 • N 2 vertical scaling factors, the connection vector V , the number of iterations n, and the PDE.
(2) The initial set G is set to contain all the interpolation points (N 1 × N 2 points).
(3) The projections of the interpolation points on R 2 form N 2 horizontal lines and N 1 vertical lines.For each one of the N 2 horizontal lines we compute (N 1 − 1) • α n + 1 points of a function that interpolates the corresponding points.Similarly, for each one of the N 1 vertical lines we compute (N 2 −1)•α n +1 points of a function that interpolates the corresponding points.
(4) Form the function H.
For i := 1 to N 1 DO: (a) For j := 1 to N 2 DO: (i) We solve numerically, inside the region I i, j , the partial differential equation with Dirichlet's boundary conditions defined by the computed points of the interpolation functions.

Figure 2 :
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.

Figure 5 :
Figure 5: An example of the algorithm, where N 1 = 3, N 2 = 3, α = 2, and n = 2.In this case, there are four regions and one domain.Thus, the RIFS is actually an IFS.The connection vector is V = (1, 1, 1, 1).The algorithm is expected to produce 81 points for the attractor.In all the figures the squares represent the projection of the points on R 2 .(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.

Figure 7 :
Figure 7: Some examples of the proposed construction.At the 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.