# Construction of Fractal Surfaces via Solutions of Partial Differential Equations

Authors: P. Bouboulis
Journal: SRX Mathematics
 (2010)

---

## Abstract

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 functionsH and B, which are then used to produce the fractal surface, as the attractor of an appropriately chosen recurrent iterated
function system.


---

## Body

## 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 (X,ρ), together with a finite set of continuous, contractive mappings wi: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–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.


## 2. Mathematical Background

### 2.1. Iterated Function Systems

As mentioned above, anIterated Function System  {X;w1-N} is defined as a pair consisting of a complete metric space (X,ρ), together with a finite set of continuous, contractive mappings wi:X→X, with respective contraction factors si, for i=1,2,…,N (N≥2).  The attractor of an IFS is the unique set E, for which E=limk→∞Wk(A0) for every starting compact set A0, where
(1)W(A)=⋃i=1Nwi(A)∀A∈ℋ(X),
and ℋ(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 〈ℋ(X),h〉 see [15]). A simple example of an IFS defined on ℝ2 is the one that produces the well-known Sierpinski's Triangle (see Figure 1(a)), which consists of the three mappings:
(2)w1(xy)=(0.5000.5)(xy),w2(xy)=(0.5000.5)(xy)+(050),w3(xy)=(0.5000.5)(xy)+(500).
The attractor of the following IFS looks like a natural fern (see Figure 1(b)):Two known attractors that arise from IFS. (a) Sierpinski's Triangle, (b) a fern.
(a)(b)(3)w1(xy)=(0000.16)(xy)+(00),w2(xy)=(0.850.04-0.040.85)(xy)+(01.6),w3(xy)=(0.2-0.260.230.22)(xy)+(01.6),w4(xy)=(-0.150.280.260.24)(xy)+(00.44).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 setA0 and produces the sets A1=W(A0)=⋃i=1Nwi(A0),A2=W(A1), and so forth. The sequence A0,A1,A2,… obviously converges to the attractor of the IFS. The RIA, on the other hand, starts with an arbitrary point x0 and selects randomly (using a set of probabilities that sum to 1) one of the mappings of the IFS, to produce a new point x1=wi1(x0) (where wi1 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, p1=0.01,p2=0.69,p3=0.15,p4=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.


### 2.2. Recurrent Iterated Function Systems

A more general concept, that allows the construction of even more complicated sets, is that of theRecurrent Iterated Function System, or RIFS for short, which consists of an IFS {X;w1-N}, together with an irreducible row-stochastic matrix  P=(pν,μ)N(pν,μ∈[0,1]:ν,μ=1,…,N), such that
(4)∑μ=1Npν,μ=1,ν=1,…,N.
The recurrent structure is given by the (irreducible) connection matrix   C=(cν,μ)N, which is defined by
(5)cν,μ={1,ifpμ,ν>0,0,ifpμ,ν=0,
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 mapW needs a little more effort. First, we define the mappings
(6)Wi,j:ℋ(X)→ℋ(X),withWi,j(A)={wi(A),pj,i>0,∅,pj,i=0
for all A∈ℋ(X) and the metric space
(7)ℋ̃(X)=ℋ(X)N=ℋ(X)×ℋ(X)×⋯×ℋ(X)
equipped with the metric
(8)h̃((A1A2⋮AN),(B1B2⋮BN))=max{h(Ai,Bi);i=1,2,…,N}.
Then 〈ℋ̃(X),h̃〉 is a complete metric space. The map W is now defined by
(9)W:ℋ̃(X)→ℋ̃(X):W(A1A2⋮AN)=(W11W12⋯W1NW21W22⋯W2N⋮⋮⋮WN1WN2⋯WNN)·(A1A2⋮AN)=(⋃j∈I(1)w1(Aj)⋃j∈I(2)w2(Aj)⋮⋃j∈I(N)wN(Aj)),
where I(i)={j:pj,i>0}, for i=1,2,…,N. If wi are contractions, then W is a contraction and there is E=(E1,E2,…,EN)t∈ℋ̃(X) such that W(E)=E and Ei=⋃j∈I(i)wi(Ej), for i=1,2,…,N. The set G=⋃i=1NEi is called the attractor of the RIFS {X;w1-N,P}.Next, we construct a sequence of sets that converge to the attractor. LetA∈ℋ(X). We define the sequences {An}n∈ℕ in ℋ̃(X) and {An}n∈ℕ in ℋ(X) as follows: A0=(A,A,…,A)t, An=W(An-1), and An=⋃i=1N(An)i, for n∈ℕ, where An=((An)1,(An)2,…,(An)N). Evidently
(10)G=limn→∞An.
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)w1(xy)=(0.5000.5)(xy)+(00),w2(xy)=(0.5000.5)(xy)+(050),(12)w3(xy)=(0.5000.5)(xy)+(500),w4(xy)=(0.5000.5)(xy)+(5050).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.
(a)(b)(c)Both DA and RIA can be modified to construct the attractor of any RIFS. The modification of DA is described above, where the sequences{An}n∈ℕ are been defined. We note that in many cases, the functions wi of the RIFS are defined on a compact subset of X, which we call Xi, instead of the space X itself, for i=1,2,…N. In this case the initial vector A0 is defined as A0=(B1,B2,…,BN)t, where Bi⊆wi(Xi),i=1,2,…N. The modification of the RIA is much simpler. We start with an initial point x0 and an initial state i0. Then we select arbitrarily, using the set of probabilities pi0,1,pi0,2,…,pi0,N, one of the mappings of the RIFS and apply it to x0 to produce a new point x1. Suppose that wi1 is the selected map, then x1=wi1(x0) and the system transits to state i1. 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 x0,1,x0,2,…,x0,N, such that x0,i∈wi(Xi), for i=1,2,…,N. For each point x0,i, i=1,2,…,N, we select randomly one of the maps of the RIFS, using the set of probabilities {pi,1,pi,2,…,pi,N}, and apply it to x0,i to produce the new point x1,i(j1,i), where j1,i∈{1,2,…,N} indicates that the selected map is wj1,i. After the completion of this step, the new set of points is {x1,1(j1,1),x1,2(j1,2),…,x1,N(j1,N)}. Again, for each point x1,i(j1,i), i=1,2,…,N, we select randomly one of the maps of the RIFS, using the set of probabilities {pj1,i,1,pj1,i,2,…,pj1,i,N}, and apply it to x1,i(j1,i) to produce the new point x2,i(j2,i), where j2,i∈{1,2,…,N} indicates that the selected map in the second step is wj2,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.


## 2.1. Iterated Function Systems

As mentioned above, anIterated Function System  {X;w1-N} is defined as a pair consisting of a complete metric space (X,ρ), together with a finite set of continuous, contractive mappings wi:X→X, with respective contraction factors si, for i=1,2,…,N (N≥2).  The attractor of an IFS is the unique set E, for which E=limk→∞Wk(A0) for every starting compact set A0, where
(1)W(A)=⋃i=1Nwi(A)∀A∈ℋ(X),
and ℋ(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 〈ℋ(X),h〉 see [15]). A simple example of an IFS defined on ℝ2 is the one that produces the well-known Sierpinski's Triangle (see Figure 1(a)), which consists of the three mappings:
(2)w1(xy)=(0.5000.5)(xy),w2(xy)=(0.5000.5)(xy)+(050),w3(xy)=(0.5000.5)(xy)+(500).
The attractor of the following IFS looks like a natural fern (see Figure 1(b)):Two known attractors that arise from IFS. (a) Sierpinski's Triangle, (b) a fern.
(a)(b)(3)w1(xy)=(0000.16)(xy)+(00),w2(xy)=(0.850.04-0.040.85)(xy)+(01.6),w3(xy)=(0.2-0.260.230.22)(xy)+(01.6),w4(xy)=(-0.150.280.260.24)(xy)+(00.44).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 setA0 and produces the sets A1=W(A0)=⋃i=1Nwi(A0),A2=W(A1), and so forth. The sequence A0,A1,A2,… obviously converges to the attractor of the IFS. The RIA, on the other hand, starts with an arbitrary point x0 and selects randomly (using a set of probabilities that sum to 1) one of the mappings of the IFS, to produce a new point x1=wi1(x0) (where wi1 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, p1=0.01,p2=0.69,p3=0.15,p4=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.


## 2.2. Recurrent Iterated Function Systems

A more general concept, that allows the construction of even more complicated sets, is that of theRecurrent Iterated Function System, or RIFS for short, which consists of an IFS {X;w1-N}, together with an irreducible row-stochastic matrix  P=(pν,μ)N(pν,μ∈[0,1]:ν,μ=1,…,N), such that
(4)∑μ=1Npν,μ=1,ν=1,…,N.
The recurrent structure is given by the (irreducible) connection matrix   C=(cν,μ)N, which is defined by
(5)cν,μ={1,ifpμ,ν>0,0,ifpμ,ν=0,
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 mapW needs a little more effort. First, we define the mappings
(6)Wi,j:ℋ(X)→ℋ(X),withWi,j(A)={wi(A),pj,i>0,∅,pj,i=0
for all A∈ℋ(X) and the metric space
(7)ℋ̃(X)=ℋ(X)N=ℋ(X)×ℋ(X)×⋯×ℋ(X)
equipped with the metric
(8)h̃((A1A2⋮AN),(B1B2⋮BN))=max{h(Ai,Bi);i=1,2,…,N}.
Then 〈ℋ̃(X),h̃〉 is a complete metric space. The map W is now defined by
(9)W:ℋ̃(X)→ℋ̃(X):W(A1A2⋮AN)=(W11W12⋯W1NW21W22⋯W2N⋮⋮⋮WN1WN2⋯WNN)·(A1A2⋮AN)=(⋃j∈I(1)w1(Aj)⋃j∈I(2)w2(Aj)⋮⋃j∈I(N)wN(Aj)),
where I(i)={j:pj,i>0}, for i=1,2,…,N. If wi are contractions, then W is a contraction and there is E=(E1,E2,…,EN)t∈ℋ̃(X) such that W(E)=E and Ei=⋃j∈I(i)wi(Ej), for i=1,2,…,N. The set G=⋃i=1NEi is called the attractor of the RIFS {X;w1-N,P}.Next, we construct a sequence of sets that converge to the attractor. LetA∈ℋ(X). We define the sequences {An}n∈ℕ in ℋ̃(X) and {An}n∈ℕ in ℋ(X) as follows: A0=(A,A,…,A)t, An=W(An-1), and An=⋃i=1N(An)i, for n∈ℕ, where An=((An)1,(An)2,…,(An)N). Evidently
(10)G=limn→∞An.
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)w1(xy)=(0.5000.5)(xy)+(00),w2(xy)=(0.5000.5)(xy)+(050),(12)w3(xy)=(0.5000.5)(xy)+(500),w4(xy)=(0.5000.5)(xy)+(5050).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.
(a)(b)(c)Both DA and RIA can be modified to construct the attractor of any RIFS. The modification of DA is described above, where the sequences{An}n∈ℕ are been defined. We note that in many cases, the functions wi of the RIFS are defined on a compact subset of X, which we call Xi, instead of the space X itself, for i=1,2,…N. In this case the initial vector A0 is defined as A0=(B1,B2,…,BN)t, where Bi⊆wi(Xi),i=1,2,…N. The modification of the RIA is much simpler. We start with an initial point x0 and an initial state i0. Then we select arbitrarily, using the set of probabilities pi0,1,pi0,2,…,pi0,N, one of the mappings of the RIFS and apply it to x0 to produce a new point x1. Suppose that wi1 is the selected map, then x1=wi1(x0) and the system transits to state i1. 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 x0,1,x0,2,…,x0,N, such that x0,i∈wi(Xi), for i=1,2,…,N. For each point x0,i, i=1,2,…,N, we select randomly one of the maps of the RIFS, using the set of probabilities {pi,1,pi,2,…,pi,N}, and apply it to x0,i to produce the new point x1,i(j1,i), where j1,i∈{1,2,…,N} indicates that the selected map is wj1,i. After the completion of this step, the new set of points is {x1,1(j1,1),x1,2(j1,2),…,x1,N(j1,N)}. Again, for each point x1,i(j1,i), i=1,2,…,N, we select randomly one of the maps of the RIFS, using the set of probabilities {pj1,i,1,pj1,i,2,…,pj1,i,N}, and apply it to x1,i(j1,i) to produce the new point x2,i(j2,i), where j2,i∈{1,2,…,N} indicates that the selected map in the second step is wj2,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.


## 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)Δ={(xi,yj,zi,j)∈I×ℝ;i=0,1,…,N1,j=0,1,…,N2(xi,yj,zi,j)}⊂I×ℝ,such that0=x0<x1<⋯<xN1=1,  0=y0<y1<⋯<yN2=1,  N1,N2∈ℕ, where I=[0,1]2, which contains in total (N1+1)·(N2+1) points. Consider, also, a data set Δ̂⊂Δ:(14)Δ̂={(x̂k,ŷl,ẑk,l)∈I×ℝ;k=0,1,…,M1,l=0,1,…,M2},such that0=x̂0<x̂1<⋯<x̂M1=1,  0=ŷ0<ŷ1<⋯<ŷM2=1M1,M2∈ℕ, which contains in total (M1+1)·(M2+1) points. To simplify the notation we set
(15)C0={(i,j);i=0,1,…,N1,j=0,1,…,N2},D0={(k,l);k=0,1,…,M1,l=0,1,…,M2}C1={(i,j);i=1,2,…,N1,j=1,2,…,N2},D1={(k,l);k=1,2,…,M1,l=1,2,…,M2}.
The points of Δ divide [0,1]2 into N1·N2 regions:
(16)Ii,j=[xi-1,xi]×[yj-1,yj]
for all (i,j)∈C1, while the points of Δ̂ divide [0,1]2 into M1·M2 domains:
(17)Jk,l=[x̂k-1,x̂k]×[ŷl-1,ŷl]
for all (k,l)∈D1. We make the additional assumption that for every (k,l)∈D1 there is at least one interpolation point, that lies in the interior of Jk,l×ℝ. Furthermore, we define a labeling map
(18)𝒥:C1→D1:𝒥(i,j)=(k,l)
and the 1-1 functions Φ, Φ̂ (an enumeration of the sets C1 and D1, respectively) such that(19)Φ:C1→{1,2,…,N1·N2}:Φ(i,j)=i+(j-1)N1,Φ̂:D1→{1,2,…,M1·M2}:Φ̂(k,l)=k+(l-1)M1.
We define the (N1·N2)×(N1·N2) stochastic matrix P by
(20)pν,μ={1γν,ifIΦ-1(ν)⊆J𝒥(Φ-1(μ)),0,otherwise,
where γν is the number of non zero elements of its νth row. Consequently, the connection matrix C is defined as in Section 2.2 and the connection vector V={v1,v2,…,vN1·N2} as follows: vν=Φ̂(𝒥(Φ-1(ν))),ν=1,2,…,N1·N2. 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)C=(00110011000000001100110000000000000000000011001100000000110011000000000000110011110011000000000000000000110011000011001100000000110011000000000000000000001100110000000011001100001100110000000000110011000000001100110000000000000000001100110000000000001100111111111111111111),V=(2,1,4,3,4,1,3,2,1,4,3,2,2,1,3,4).Figure 3
The interpolation points divide[0,1]2 in 16 regions and in 4 domains. The association between regions and domains is shown with the arrows. For example, 𝒥(1,1)=(1,2),  𝒥(1,2)=(1,1),  𝒥(1,3)=(2,2),  𝒥(1,4)=(2,1), and so forth. Thus, the connection vector is V=(2,1,4,3,4,1,3,2,1,4,3,2,2,1,3,4).To complete the definition of the RIFS, we considerN1·N2 mappings of the form
(22)wi,j:J𝒥(i,j)×ℝ→Ii,j×ℝ:wi,j(xyz)=(T1,i(x)T2,j(y)Fi,j(x,y,z))=(T1,i(x)T2,j(y)si,jz+Qi,j(x,y)),
with
(23)T1,i(x)=a1,ix+b1,i,T2,j(y)=a2,jy+b2,jfor all(x,y)∈Ii,j, where Qi,j is a Lipschitz continuous function on [0,1]n and si,j∈(-1,1), for all (i,j)∈C1. The parameters si,j,(i,j)∈C1 are called vertical scaling factors. It is easy to show that there exists a metric ρθ (equivalent with the Euclidean metric) such that wi,j is a contraction with respect to the metric h for all (i,j)∈C1. We confine the map wi,j, so that it maps the interpolation points that lie on the vertices of J𝒥(i,j) to the interpolation points that lie on the vertices of Ii,j. Hence, we obtain the following relations:
(24)T1,i(x̂k-1)=xi-1,T1,i(x̂k)=xi,T2,j(ŷl-1)=yj-1,T2,j(ŷl)=yj,Fi,j(x̂k-δ1,ŷl-δ2,ẑk-δ1,l-δ2)=zi-δ1,j-δ2
for all δ1,δ2∈{0,1},  (k,l)=𝒥(i,j), and  (i,j)∈C1. Using the first four (12) equations, we can compute a1,i,  a2,j,  b1,i,  b2,j, in terms of the interpolation points and the vertical scaling factors, for all (i,j)∈C1.Consequently, the RIFS{[0,1]2×ℝ,wi,j,P;(i,j)∈C1} has a unique attractor G. In general G is a compact subset of ℝ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 ∂Ii,j×ℝ, for all (i,j)∈C1, (where ∂Ii,j is the boundary of Ii,j).Proposition 1.
Leth∈C(⋃(i,j)∈C1∂Ii,j) be a Lipschitz function, that interpolates the points of Δ (i.e., h(xi,yj)=zi,j,(i,j)∈C1). If the RIFS defined above satisfies the conditions
(25)Fi,j(x̂k-1,y,h(x̂k-1,y))=h(xi-1,T2,j(y)),Fi,j(x̂k,y,h(x̂k,y))=h(xi,T2,j(y)),Fi,j(x,ŷl-1,h(x,ŷl-1))=h(T1,i(x),yj-1),Fi,j(x,ŷl,h(x,ŷl))=h(T1,i(x),yj),
where (k,l)=𝒥(i,j), for all (x,y)∈J𝒥(i,j),(i,j)∈C1, 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 functionsQi,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.
Leth∈C(⋃(i,j)∈C1∂Ii,j) be a Lipschitz function that interpolates the points of Δ. Consider the case where
(26)Qi,j(x,y)=H(T1,i(x),T2,j(y))-si,j·B(x,y),∀(x,y)∈J𝒥(i,j),(i,j)∈C1,
where H andB are Lipschitz functions defined on [0,1]2, such that
(27)H(x,y)=h(x,y),for(x,y)∈∂Ii,j,B(x,y)=h(x,y),for(x,y)∈∂J𝒥(i,j).
The unique attractor G of the corresponding RIFS {[0,1]2×ℝ,wi,j,P;(i,j)∈C1} is the graph of a continuous function f that interpolates the points of Δ.The functionsH 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 si,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 si, si2, 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 functionH|[xi-1,xi] that solves the Laplace's equation with boundary conditions H(xi)=yi,  H(xi-1)=yi-1, for i=1,…,N,(ii)
find a functionB that solves the Laplace's equation with boundary conditions B(x0)=y0,  B(xN)=yN,where{(x0,y0),…,(xN,yN)} 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 setC̃Δk, which consists of continuous functions that are defined only on ⋃(i,j)∈C1∂Ii,j, interpolate the points of Δ, and have continuous partial derivatives up to kth order in ⋃(i,j)∈C1∂Ii,j. Let ℛ([0,1]2) be the set containing all the subsets of [0,1]2 of the form [a1,b1]×[a2,b2] and a kth order Partial Differential Equation (PDE) defined on the set R∈ℛ([0,1]2), such that its solution g satisfies
(28)g|∂R=v,
where v∈Ck(∂R). The problem of finding a function, which solves a specified partial differential equation (PDE) in the interior of a given region (Ii,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∈ℛ([0,1]2) and any v∈Ck(∂R) and that the solution is a Lipschitz function, we consider the operator
(29)𝒫R:Ck(∂R)→Ck(R):𝒫R(v)=g,
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 functionh∈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 Qi,j(x,y)=H(T1,i(x),T2,j(y))-si,j·B(x,y), for all (x,y)∈J𝒥(i,j),(i,j)∈C1, as in Corollary 1. We compute H and B as follows:
(30)H|Ii,j=𝒫Ii,j(h|∂Ii,j),B|𝒥(i,j)=𝒫J𝒥(i,j)(h|∂J𝒥(i,j)),
for all (i,j)∈C1. In this case, the conditions of Corollary 1 are satisfied for the RIFS {[0,1]2×ℝ,wi,j,P;(i,j)∈C1}, 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 N1·N2 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 N1·N2 regions Ii,j, for all (i,j)∈C1, while the points of Δ̂ form the M1·M2 domains Jk,l. 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 V.(3)
We construct arbitrary(N1+1)·(N2+1) 1-dimensional functions that interpolate the points of Δ. This is the function h (see Figure 4(b)).(4)
We solve (numerically) the partial differential equation inside the regionIij, so that its solution g satisfies the boundary conditions g|∂Iij=h|∂Iij, for all (i,j)∈C1 (Dirichlet's problem). Hence, we obtain the function H (see Figure 4(c)).(5)
We solve (numerically) the partial differential equation insideJkl, so that its solution g satisfies the boundary conditions g|∂Jkl=h|∂Jkl, for all (k,l)∈D1. Hence, we obtain the function B (see Figure 4(d)).(6)
Finally, we construct the attractor of the corresponding RIFS (see Figure4(e)).An example of the construction. HereN1=N2=5, the regions have dimensions 0.25×0.25 and the domains 0.5×0.5. We used Laplace's partial differential equation (∂2g/∂x2+∂2g/∂y2=0) with Dirichlet's boundary conditions. (a) The set of interpolation points. (b) The function h. (c) The function H. (d) The function B. (e) The fractal interpolation surface.
(a)(b)(c)(d)(e)


## 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 functionH|[xi-1,xi] that solves the Laplace's equation with boundary conditions H(xi)=yi,  H(xi-1)=yi-1, for i=1,…,N,(ii)
find a functionB that solves the Laplace's equation with boundary conditions B(x0)=y0,  B(xN)=yN,where{(x0,y0),…,(xN,yN)} 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 setC̃Δk, which consists of continuous functions that are defined only on ⋃(i,j)∈C1∂Ii,j, interpolate the points of Δ, and have continuous partial derivatives up to kth order in ⋃(i,j)∈C1∂Ii,j. Let ℛ([0,1]2) be the set containing all the subsets of [0,1]2 of the form [a1,b1]×[a2,b2] and a kth order Partial Differential Equation (PDE) defined on the set R∈ℛ([0,1]2), such that its solution g satisfies
(28)g|∂R=v,
where v∈Ck(∂R). The problem of finding a function, which solves a specified partial differential equation (PDE) in the interior of a given region (Ii,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∈ℛ([0,1]2) and any v∈Ck(∂R) and that the solution is a Lipschitz function, we consider the operator
(29)𝒫R:Ck(∂R)→Ck(R):𝒫R(v)=g,
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 functionh∈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 Qi,j(x,y)=H(T1,i(x),T2,j(y))-si,j·B(x,y), for all (x,y)∈J𝒥(i,j),(i,j)∈C1, as in Corollary 1. We compute H and B as follows:
(30)H|Ii,j=𝒫Ii,j(h|∂Ii,j),B|𝒥(i,j)=𝒫J𝒥(i,j)(h|∂J𝒥(i,j)),
for all (i,j)∈C1. In this case, the conditions of Corollary 1 are satisfied for the RIFS {[0,1]2×ℝ,wi,j,P;(i,j)∈C1}, 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 N1·N2 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 N1·N2 regions Ii,j, for all (i,j)∈C1, while the points of Δ̂ form the M1·M2 domains Jk,l. 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 V.(3)
We construct arbitrary(N1+1)·(N2+1) 1-dimensional functions that interpolate the points of Δ. This is the function h (see Figure 4(b)).(4)
We solve (numerically) the partial differential equation inside the regionIij, so that its solution g satisfies the boundary conditions g|∂Iij=h|∂Iij, for all (i,j)∈C1 (Dirichlet's problem). Hence, we obtain the function H (see Figure 4(c)).(5)
We solve (numerically) the partial differential equation insideJkl, so that its solution g satisfies the boundary conditions g|∂Jkl=h|∂Jkl, for all (k,l)∈D1. Hence, we obtain the function B (see Figure 4(d)).(6)
Finally, we construct the attractor of the corresponding RIFS (see Figure4(e)).An example of the construction. HereN1=N2=5, the regions have dimensions 0.25×0.25 and the domains 0.5×0.5. We used Laplace's partial differential equation (∂2g/∂x2+∂2g/∂y2=0) with Dirichlet's boundary conditions. (a) The set of interpolation points. (b) The function h. (c) The function H. (d) The function B. (e) The fractal interpolation surface.
(a)(b)(c)(d)(e)


## 4. The Algorithm

In Section2.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 ℝ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 α∈ℕ. The algorithm produces ((N1-1)·αn+1)×((N2-1)·αn+1) points (n is the number of iterations), so that their projections on ℝ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 N1-1 must be a multiple of α (to be more specific: N1-1=α·(M1-1)).An example of the algorithm, whereN1=3,N2=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 ℝ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.
(a)(b)(c)(1)
INPUT: TheN1×N2 interpolation points, the number α, that defines which points belong in Δ and which points belong in Δ̂, a set of N1·N2 vertical scaling factors, the connection vector V, the number of iterations n, and the PDE.(2)
The initial setG is set to contain all the interpolation points (N1×N2 points).(3)
The projections of the interpolation points onℝ2 form N2 horizontal lines and N1 vertical lines. For each one of the N2 horizontal lines we compute (N1-1)·αn+1 points of a function that interpolates the corresponding points. Similarly, for each one of the N1 vertical lines we compute (N2-1)·αn+1 points of a function that interpolates the corresponding points.(4)
Form the function H.Fori:=1 to N1 DO:
(a)
Forj:=1 to N2 DO:
(i)
We solve numerically, inside the regionIi,j, the partial differential equation with Dirichlet's boundary conditions defined by the computed points of the interpolation functions.(5)
Form the function B. Fork:=1 to M1 DO:
(a)
Forl:=1 to M2 DO:
(i)
We solve numerically, inside the domainJk,l, 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. Forr:=1 to n DO
(a)
Fori:=1 to N1 DO:
(i)
Forj:=1 to N2 DO:
(A)
We deal with the points of the regionIi,j. First of all, we find the corresponding domain Jk,l, using the connection vector: (k,l)=Φ̂(vΦ(i,j)).(B)
As the parameters of the mapwi,j are computed by the functions H and B, we may use wi,j to map all the points of G, that lie inside Jk,l. Therefore, we obtain some new points inside the region Ii,j.(b)
The new set G contains all the points that were produced by the previous iterations (((N1-1)·αr+1)×((N2-1)·αr+1) points).(7)
OUTPUT: After the completion of the iterations, the final setG 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[0,1]2 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.The two fractal surfaces (shown above) were constructed using the same parameters, but different partial differential equations. In (a) Laplace's equation was used (∂2g/∂x2+∂2g/∂y2=0), while in (b) the corresponding PDE is: ∂2g/∂x2+∂2g/∂y2=80.
(a)(b)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.
(a)(b)(c)(d)



---

*Source: 432521-2010-01-13.xml*