## Abstract

This chapter provides a brief and coarse discussion on the theory of fractal interpolation functions and their recent developments including some of the research made by the authors. It focuses on fractal interpolation as well as on recurrent fractal interpolation in one and two dimensions. The resulting self‐affine or self‐similar graphs, which usually have non‐integral dimension, were generated through a family of (discrete) dynamic systems, the iterated function system, by using affine transformations. Specifically, the fractal interpolation surfaces presented here were constructed over triangular as well as over polygonal lattices with triangular subdomains. A further purpose of this chapter is the exploration of the existent breakthroughs and their application to a flexible and integrated software that constructs and visualises the above‐mentioned models. We intent to supply both a panoramic view of interpolating functions and a useful source of links to assist a novice as well as an expert in fractals. The ideas or findings contained in this paper are not claimed to be exhaustive, but are intended to be read before, or in parallel with, technical papers available in the literature on this subject.

### Keywords

- approximation
- attractor
- fractal
- interpolating function
- iterated function system
- recurrent
- self‐affinity
- self‐similarity
- surface construction

## 1. Introduction

In the mathematical field of numerical analysis, *interpolation* is a method of constructing new data points within the range of a discrete set of known data points. Interpolation by fractal (graph of) functions, as defined in Refs. [1, 2], is based on the theory of *iterated function systems* and can be seen as an alternative to traditional interpolation techniques, aiming mainly at data which present detail at different scales, some degree of self‐similarity or self‐affinity. A *fractal interpolation function* can be seen as a continuous function whose graph is the attractor of an appropriately chosen iterated function system. This attractor is called a *fractal interpolation surface*, since the graph of a continuous function of two variables defined over a connected open subset of ^{2} is a *topological surface*, if the so formed graph belongs to the three‐dimensional space and has *Hausdorff‐Besicovitch dimension* between 2 and 3. It is called a *fractal interpolation volume* whenever the graph has dimension greater than three. The key difficulty in constructing fractal interpolation surfaces (or volumes) involves ensuring continuity. Another important element necessary in modelling complicated surfaces of this type is the existence of the *contractivity*, or *vertical scaling*, *factors*.

### 1.1. Historical background

Massopust introduced in Ref. [3] fractal surfaces constructed as attractors of iterated function systems. He considered the case of a triangular domain with coplanar boundary data. Later on, Geronimo and Hardin in Ref. [4] presented a slightly more general construction of such fractal surfaces. They examined the case when the domain is a polygonal region with arbitrary interpolation points but with identical contractivity factors. In Ref. [5], ^{m}‐valued multivariate fractal functions were investigated. The latter two constructions use the *recurrent* iterated‐function‐system formalism. The construction of Wittenbrink [6] either produces discontinuous surfaces (and volumes) or reduces to the case, where the contractivity factors must be constant. Zhao [7] allows the contractivity factors to become a continuous ‘contraction function’ and uses consistent triangulation in order to guarantee continuity. All of the previous constructions are based on triangular subdomains. As it is always possible to construct fractal surfaces as tensor products of univariate continuous fractal functions, Massopust in Ref. [8], Section 9.4 suggests a construction by taking the tensor product of two univariate fractal interpolation functions. The derived function is uniquely determined by its evaluation along a pair of adjacent sides of the rectangular domain.

Two piecewise self‐affine models for representing discrete image data on rectangular lattices by using fractal surfaces are proposed in Ref. [9]. In Ref. [10], the piecewise self‐affine IFS model is extended from ^{3} to ^{n} (*n* is an integer greater than 3), which is called the *multi‐dimensional piecewise self‐affine fractal interpolation model*. The same methodology is repeated in Ref. [11]. According to Ref. [12], the self‐affine iterated function system model is extended from ^{3} to ^{n} (*n* is an integer greater than 3), which is called the *multi‐dimensional self‐affine fractal interpolation model*. Vijender and Chand in Ref. [13] proposed a class of affine fractal interpolation surfaces that stitch a given set of surface data arranged on a rectangular grid. The created fractal interpolation surfaces are a blend of the affine fractal interpolation functions constructed along the grid lines of a given interpolation domain.

The study of fractals is a field in science that unifies mathematics, theoretical physics, art and computer science. Therefore, it is not difficult to find applications of fractal interpolation functions in almost every scientific field wherein information available in finite number of grid points has to be modelled with a continuous function. Applications of this theory include geometric design, data visualisation, reverse engineering, physics, geology, image encoding and compression (see Ref. [14]), signal processing and wavelet theory. Fractal interpolation also provides a good representation of economic time series such as the stock market fluctuation and weather data. The reason for this variety of applications lies in the underlying complicated mathematical structure of fractal (graph of) functions, produced with simple recursive construction. It has been noted that, for certain problems, they provide better approximants than their classical non‐recursive counterparts.

### 1.2. Clarifications and interpretations

Although, in most of the cases, the words *function*, *map* and *mapping* can be used as nouns interchangeably in several parts of mathematics, they generally do not have the same meaning. Originally, the word ‘map’ comes from the medieval Latin *Mappa mundi*, wherein *mappa* meant napkin or cloth and *mundi* the world. ‘Mapping’ as a noun is the process of making maps. In mathematics, we think of a *map* as a way of sending elements from one set to elements of another set. A *mapping* shows how the elements are paired. A *function* is a way or rule of associating elements from one set, the *domain*, to elements of another set, the *codomain*. So, a domain is where a function maps from, and a *range*, as a subset of the codomain, is where it maps to. The definition of a *function* requires that each element in the domain corresponds to one and only one element of the codomain. Moreover, a function is commonly used as a special type of mapping, namely it is used as a mapping from a set into the set of numbers, that is, into *transformation*.

An *affine map* (from the Latin, *affinis*, “connected with”) between two (vector) spaces consists of a *linear map* followed by a translation. Similarly, one can define an *affine transformation*. In a geometric setting, these are precisely the functions that map straight lines to straight lines. Be aware that the term *linear function* refers to two distinct but related notions. It may be a linear map or a polynomial of degree one or less, including the zero polynomial, because its graph, when there is only one independent variable, is a non‐vertical straight line.

Another common error involves the incorrect use of the notions ‘self‐affine function’ and ‘self‐similar function’. There is no ‘self‐affine’ or ‘self‐similar’ function but a function with a self‐affine or self‐similar graph. Therefore, an object or a set, but not a function, can be self‐affine or self‐similar. Each part of a *self‐affine* object is an image of the whole object (either strictly or in a statistical sense) scaled differently in different directions. Self‐affine sets form an important class of sets, which include self‐similar sets as a particular case. A *self‐similar* object is exactly or approximately similar to a part of itself (i.e. the whole has the same shape as one or more of the parts).

The aim of our paper is to review the usage of fractal interpolation functions in order to construct self‐affine graphs generated by iterated function systems. Furthermore, we compare and contrast several constructions presented in Refs. [3–8] by pointing out some of their ambiguities, limitations and restrictions. Particularly, in Section 2, we briefly review the theory of iterated function systems. In Section 3, we revisit the 1D fractal interpolation theory and state the prerequisites of the constructions. Necessary conditions for the attractor of an iterated function system to be the graph of a continuous function interpolating a given set of data are also given. In Section 4, we revisit the two‐dimensional fractal interpolation theory. A comparison to existing methods and some examples of fractal interpolation surfaces constructed by them are also presented. The corresponding algorithms for constructing these surfaces are developed and illustrated through several graphic examples. Finally, Section 5 summarises our conclusions and points out areas of future work.

## 2. Fractal image generation

In mathematics, an *iterated map* is a map composed with itself, possibly *ad infinitum*, in a process called iteration. *Iteration* means the act of repeating a process with the aim of approaching a desired goal, target or result. More formally, let *X* be a set and *f*: *X*→*X* be a map. Define *fk* as the *k*‐th iterate of *f*, where *k* is a non‐negative integer, by *f*^{0} = id_{X} and *fk*^{+1} = *f* ◦ *fk*, where id_{X} is the identity map on *X* and *f* ◦ *g* denotes map composition.

A *contraction mapping*, or *contraction*, on a metric space (*X*, *ρ*) is a function *f* from *Χ* to itself, that is, a *transformation*, with the property that there is a non‐negative real number *s* < 1 such that for all *x* and *y* in *X*, *ρ*(*f*(*x*), *f*(*y*)) ≤ *s*⋅*ρ*(*x*, *y*), where *ρ* is a *distance function* between elements of *X*. The smallest such value of *s* is called the *Lipschitz constant* or *contractivity factor* of *f*. If the above condition is satisfied for *s* ≤ 1, then the mapping is said to be *non‐expansive*. A contraction mapping has at most one *fixed point*, that is, a point *x** in *X* such that *f*(*x**) = *x**. Moreover, the *Banach fixed point theorem*, also known as the *contraction mapping theorem* or *contraction mapping principle*, states that every contraction mapping on a non‐empty, complete metric space has a unique fixed point, and that for any *x* in *X* the iterated function sequence *x*, *f*(*x*), *f*(*f*(*x*)), *f*(*f*(*f*(*x*))), … converges to this fixed point. Furthermore, this fixed point can be found as follows: Start with an arbitrary element *x*_{0} in *X* and define an iterative sequence by *xn* = *f*(*xn*_{−1}) for *n* = 1, 2, 3, …. This sequence converges and its limit is *x**.

### 2.1. Iterated function systems

An *iterated function system*, or *IFS* for short, is defined as a collection of a complete metric space (*X*, *ρ*), for example, (^{n}, ‖⋅‖) or a subset, together with a finite set of continuous transformations {*w*_{i}: *X*→*X*, *i* = 1, 2, …, *M*}. It is often convenient to write an IFS formally as {*X*; *w*_{1}, *w*_{2}, …, *wM*} or, somewhat more briefly, as {*X*; *w*_{1−M}}. If *wi* are contractions with respective contractivity factors *si*, *i* = 1, 2, …, *M*, the IFS is termed *hyperbolic*.

Hutchinson in Ref. [15] showed that, for the metric space ^{n}, such a (hyperbolic) system of functions has a unique *compact* (closed and bounded) *fixed set S*. One way for constructing a fixed set is to start with an initial point or set *S*_{0} and iterate the actions of the *wi*, taking *Sn*_{+1} to be the union of the image of *Sn* under the *wi*; then taking *S* to be the closure of the union of the *Sn*. Symbolically, the unique fixed (non‐empty compact) set *S* ⊂

The set *S* is thus the fixed set of the *Hutchinson operator*

where *A* is any subset of ^{n}. The operator *H* itself is a contraction with contractivity factor *s* = max{*s*_{1}, *s*_{2}, …, *sM*} (Ref. [2], Theorem 7.1, p. 81 or Ref. [16]). The existence and uniqueness of *S*, which is called the *attractor*, or *invariant set*, of the IFS, are a consequence of the contraction mapping principle as is the fact that *A* in ^{n}), where *X*) is the metric space of all non‐empty, compact subsets of *X* with respect to some metric, for example, the *Hausdorff metric*. The operator *H* is also called the *collage map* to alert us to the fact that *H*(*A*) is formed as a union or ‘collage’ of sets. If *X* is a Euclidean space and the *wi* are *similitudes*, that is, *similarity transformations*, then the attractor is called a *self‐similar set*. These sets are usually fractals.

A fractal derived by an IFS is made up of the union of several copies of itself, each copy being transformed by a function (hence ‘function system’). A canonical example is the Sierpiński gasket; see Figure 1. The functions are normally contractions which means they bring points closer together and make shapes smaller. Hence, such a shape is made up of several possibly overlapping smaller copies of itself, each of which is also made up of copies of itself, ad infinitum. This is the source of its self‐similar nature. Note that this infinite process is not dependent upon the starting shape being a triangle—it is just clearer that way. The first few steps starting, for example, from a square also tend towards a Sierpiński gasket; see Figure 2.

Sometimes each function *wi* is required to be a linear, or more generally an affine transformation, and hence represented by a matrix. Formally, a transformation *w* is *affine*, if it may be represented by a matrix **A** and translation **t** as *w*(**x**) = **Ax** + **t**, or, if *X* = ^{2},

The *code* of *w* is the 6‐tuple (*a*, *b*, *c*, *s*, *d*, *e*) and the *code of an IFS* is a table whose rows are the codes of *w*_{1}, *w*_{2}, …, *wM*. For the three‐dimensional case, this becomes

However, IFS’s may also be built from non‐linear functions, including *projections* and *Möbius transformations*.

### 2.2. Recurrent iterated function systems

An *IFS with probabilities*, written formally as {*X*; *w*_{1}, *w*_{2}, …,*wM*; *p*_{1}, *p*_{2}, …, *pM*} or, somewhat more briefly, as {*X*; *w*_{1−M}; *p*_{1−M}}, gives to each transformation in *H* a probability or weight. If the weights of transformations differ, so do the measures on different parts of the attractor. A non‐self‐similar attractor, however, is more easily represented with a *recurrent iterated function system*, or *RIFS* for short. Each transformation has, instead of a single weight for the next iteration, a vector of weights for each transformation, {*X*; *w*_{1−M}; *pi*_{,j} ∈ [0, 1]; *i*, *j* = 1, 2, …, *M*}, so that the matrix of weights is a recurrent Markov operator for the Hutchinson operator’s transformation (see Ref. [17]). Therefore, the attractor of a RIFS need not exhibits the self‐similarity or self‐tiling properties characteristic of the attractor of an IFS, such as Eq. (1).

The most common algorithm to compute fractals derived by IFS is called the *chaos game* or *random iteration algorithm*. It consists of picking a random point in the plane, then iteratively applying one of the functions chosen at random from the function system and drawing the point. An alternative algorithm, the *deterministic iteration algorithm*, or *DIA* for short, is to generate each possible sequence of functions up to a given maximum length and then to plot the results of applying each of these sequences of functions to an initial point or shape.

## 3. Fractal interpolation functions in R

FIFs are suitable for data sets with points linearly ordered with respect to their abscissa. This is often sufficient, for example, when interpolating time series data. In practice, however, there are many cases where the data are suitable for fractal interpolation but define a curve rather than a function, for example, when modelling coastlines or plants. There exist methods for constructing *fractal interpolation curves* based on the theory of FIFs. These methods use various approaches, such as generalisations to higher dimensions, use of index coordinates or application of reversible transformations.

Based on a theorem of Hutchinson ([15], p. 731) and using IFS theory [16], Barnsley introduced a class of functions (see Ref. [1]) which he called *fractal interpolation functions*. He basically worked with affine fractal interpolation functions, in the sense that they are obtained using affine transformations. Their main difference from elementary functions is their fractal character. Since their graphs usually have non‐integral dimension, they can be used to approximate image components such as the profiles of mountain ranges, the tops of clouds and horizons over forests, to name but a few. For a short survey on fractal interpolation functions see Ref. [19].

Let the continuous function *f* be defined on a real closed interval *I* = [*x*_{0}, *xM*] and with codomain the metric space (*x*_{0} < *x*_{1} < … < *xM*. It is not assumed that these points are equidistant. The function *f* is called an *interpolation function* corresponding to the *generalised set of data* {(*xm*, *ym*) ∈ *K* : *m* = 0, 1, …, *M*}, if *f*(*xm*) = *ym* for all *m* = 0, 1, …, *M* and *K* = *I* × *xm*, *ym*) ∈ ^{2} are called the *interpolation points*. We say that the function *f interpolates* the data and that (the graph of) *f* passes through the interpolation points. The *graph* of *f* is the set of points *G*(*f*) = {(*x*, *f*(*x*) : *x* ∈ *X*}.

### 3.1. Affine fractal interpolation

Let us represent our, real valued, set of *data points* as {(*un*, *vn*) : *n* = 0, 1, …, *N*; *un* < *un*_{+1}} and the interpolation points as {(*xm*, *ym*) : *m* = 0, 1, …, *M*; *M* ≤ *N*}, where *un* is the sampled index and *vn* the value of the given point in *un*. Let {^{2}; *w*_{1−M}} be an IFS with affine transformations of the special form (see Eq. (2))

constrained to satisfy

for every *i* = 1, 2, …, *M*. Solving the above equations results in

i.e., the coefficients *ai*, *ci*, *di*, *ei* are completely determined by the interpolation points, while the *si* are free parameters satisfying |*si*| < 1 in order to guarantee that the IFS is hyperbolic with respect to an appropriate metric for every *i* = 1, 2, …, *M*. The transformations *wi* are *shear transformations*: line segments parallel to the *y*‐axis are mapped to line segments parallel to the *y*‐axis contracted by the factor |*si*|. For this reason, the *si* are called *vertical scaling* (or *contractivity*) *factors*.

The IFS {^{2}; *w*_{1−M}} has a unique attractor that is the graph of some continuous function which interpolates the data points. This function is called a *fractal interpolation function*, or FIF for short, because its graph usually has non‐integral dimension. A *section* is defined as the function values between interpolation points. It is a function with a self‐affine graph since each affine transformation *wi* maps the entire (graph of the) function to its section. The above function is known as *affine FIF*, or *AFIF* for short. For example, let {(0, 0), (0.4, 0.5), (0.7, 0.2), (1, 0)} be a given set of data points. Figure 3 shows the graph of an AFIF with *s*_{1} = 0.5, *s*_{2} = −0.2 and *s*_{3} = 0.4. The closeness of fit of a FIF is mainly influenced by the determination of its vertical scaling factors. No direct way to find the optimum values of these factors exists, and various approaches have been proposed in the literature.

### 3.2. Piecewise affine fractal interpolation

The *piecewise self‐affine* fractal model is a generalisation of the affine fractal interpolation model and has its mathematical roots embedded in RIFS theory. A pair of data points *addresses* or *address points*, is now associated with each *interpolation interval*. Each pair of addresses defines the *domain* or *address interval*. The constraints Eq. (4) become

subjected to *i* = 1, 2, …, *M*. Solving the above equations results in

for every *i* = 1, 2, …, *M*. The function constructed as the attractor of the above‐mentioned IFS is called *recurrent fractal interpolation function*, or *RFIF* shortly, corresponding to the interpolation points. A RFIF is a piecewise self‐affine function since each affine transformation *wi* maps the part of the (graph of the) function defined by the corresponding address interval to each section.

## 4. Fractal interpolation functions in R ^{2}

Let the discrete data {(*xi*, *yj*, *zi*_{,j} = *z*(*xi*, *yj*)) ∈ ^{3} : *i* = 0, 1, …, *N*; *j* = 0, 1, …, *M*} be known. Each affine mapping that comprises the hyperbolic IFS {^{3}; *w*_{1−N, 1−M}} is given by the following special form of Eq. (3)

with |*sn*_{,m}| < 1 for every *n* = 1, 2, …, *N* and *m* = 1, 2,…, *M*. The condition

ensures that

is a similitude and the transformed surface does not vanish or flip over.

A formal definition for the fractal interpolation surfaces, as presented in Ref. [4], with some generalisations is given below. Let *D* be a convex polygon in ^{2} with *ℓ* vertices and let *S* = {*q*_{0}, *q*_{1, …,} *qm*_{−1}} be *m*, with *m* > *ℓ*, distinct points in *D* such that *q*_{0}, *q*_{1,…,} *qℓ*_{−1} are the vertices of *D*. Given real numbers *z*_{0}, *z*_{1}, …, *zm*_{−1}, we wish to construct a function *f* such that *f*(*qj*) = *zj*, *j* = 0, 1, …, *m*−1 and whose graph is self‐affine or self‐similar. Let us denote by *C*(*D*) the linear space of all real‐valued continuous functions defined on *D*, that is,

A mapping Φ on *C*(*D*) which corresponds to piecing the *G*(Φ(*f*)) together using copies of parts of *G*(*f*) will be defined.

### 4.1. Affine fractal interpolation

The basic idea is to decompose *D* into *N* non‐degenerate subregions *Δ*_{1}, *Δ*_{2}, …, *Δ*_{N} with vertices chosen from *S* and define affine maps *Li*: *D* → *Δ*_{i} and *Fi*: *D* × *i* = 1, 2, …, *N* such that Φ defined by

for **x** ∈ *Δ*_{i} maps an appropriate subset of *C*(*D*) onto itself. The main work is in showing that Φ is well defined and contractive on some subset of *C*(*D*). If *Li* is invertible, *G*(*f*) is mapped onto *Li*, *Fi*). Henceforth, we assume that the set *S*. Let *k*: {1, 2, …, *N*} × {0, 1, …, *ℓ*−1} → {0, 1, …, *m*−1} be such that

Let *i* ∈ {1, 2, …, *N*}. Since *D* and *Δ*_{i} are non‐degenerate, there exists an invertible map satisfying

Also, given any necessary free parameters,

With these definitions for *Li* and *Fi*, if *f* ∈ *C*(*D*) and *f*(*qj*) = *zj*, *j* = 0, 1, …, *ℓ*−1, then Φ_{Δi} ∈ *C*(*Δ*_{i}) and Φ(*f*)(*qk*_{(i, j)}) = *zk*_{(i, j)}, *j* = 0, 1, …, *ℓ*−1. If *Δ*_{i} and *Δ*_{i'} are adjacent polygons with common edge *f*) satisfies for all **x** ∈

When there is no proof that our construction always satisfies it, the surface may be not continuous and a geometrical and visual artefact, known as ‘discontinuity’, appears. When a surface suffers from discontinuities, a correct visualisation should render aligned horizontal holes over the surface. A surface with discontinuities should not be considered as a fractal interpolation surface; the function *f* is ambiguous on the common edge points and Φ(*f*) is not well defined.

Let the number of extreme points of the convex region *D* be 3. A triangular domain is formed, and the set *i* = 1, 2, …, *N* the invertible map *Li*: ^{2}→ ^{2} is defined as follows

and let the map *Fi*: ^{3}→

Then, the corresponding IFS is of the form {^{3}; *w*_{1−N}}, where *wi*(*x*, *y*, *z*) = (*Li*(*x*, *y*), *Fi*(*x*, *y*, *z*)), *i* = 1, 2, …, *N* and can be written in a matrix form

The real numbers *ai*, *bi*, *ci* and *di* for *i* = 1, 2, …, *N* are chosen to ensure that Condition (6) holds, that is, *Li*(*q*_{0}) = *qk*_{(i, 0)}, *Li*(*q*_{1}) = *qk*_{(i, 1)} and *Li*(*q*_{2}) = *qk*_{(i, 2)}. Thus, for *i* = 1, 2, …, *N*,

After selecting the scaling factors *si* with |*si*|< 1, the values *ei*, *fi* and *gi* are chosen to ensure that Condition (7) holds, that is, *Fi*(*q*_{0}, *z*_{0}) = *zk*_{(i, 0)}, *Fi*(*q*_{1}, *z*_{1}) = *zk*_{(i, 1)} and *Fi*(*q*_{2}, *z*_{2}) = *zk*_{(i, 2)}. That is, *si* ∈ (−1, 1) is chosen and then

for all *i* = 1, 2, …, *N.* We construct the IFS of the form {*S; w*_{1−N}}, with the intention to produce *G*(*f*) that is continuous and passes through the points (*qj*, *zj*), *qj* ∈ *S*, *j* = 0, 1, …, *m*−1 as the unique attractor *f* ∈ *C*(*D*), we must check if the ‘join‐up’ condition is satisfied for all possible starting points of the IFS. Figure 4 illustrates the surface graph (Level = 0) drawn with the original set of data {(0, 0, 0), (0.5, 0, 0.2), (1, 0, 0), (0, 0.5, 0.3), (0.5, 0.5, 0.5), (0, 1, 0)} where *s* = 0.6. See also Example 58 of Ref. [8].

#### 4.1.1. Coplanar boundary data, arbitrary contractivity factors

If *P* is a non‐vertical plane in ^{3}, let *Ĉ*(*D*) denote the collection of continuous functions *f*: *D* → *x*, *f*(*x*)) ∈ *P* for all *x* ∈ ∂*D*.

**Theorem 4.1.1.** Suppose the points {(*qj*, *zj*) :*qj* ∈ ∂*D*} are contained in a plane *P* ⊂ ^{3}. Let Φ be defined by (5), where *Li* and *Fi*, *i* = 1, 2, …, *N* are defined by Eqs. (6)–(8). Then, Φ: *Ĉ*(*D*) → *Ĉ*(*D*) is well defined and contractive in the sup‐norm ‖⋅‖_{∞}. Furthermore, for every *j* = 0, 1, …, *m*−1 and *f* ∈ *Ĉ*(*D*), Φ(*f*)(*qj*) = *zj*.

We call the unique attractor of the afore‐mentioned IFS a *self‐affine fractal interpolation surface*, or *SAFIS* for short, *with coplanar boundary data and arbitrary contractivity factors*. Figure 5 illustrates Example 1 constructed in Ref. [4].

#### 4.1.2. Arbitrary boundary data, identical contractivity factors

The *chromatic number* of a graph is the smallest number of colours needed to colour its vertices so that no two adjacent vertices share the same colour. Let *G* be the graph with *S* as its vertices and its edges correspond to the decomposition of *D* to *l = l*(*j*) ∈ {0, 1, 2}, *j* = 0, 1, …, *m*−1; see Figure 6.

For *i* = 1, 2, …, *N* and *j* = 0, 1, 2 let *k*(*i*, *j*) be determined by the condition *k*(*i*, *l*(*j′*)) = *j′* for all vertices *qj’* of Δ_{i}. Then, for each of the vertices *qj* of Δ_{i}, Eqs. (6) and (7) become

Let *Ĉ*(*D*) denote the collection of continuous functions *f* such that *f*(*qj*) = *zj*, *qj* ∈ ∂*D*.

**Theorem 4.1.2.** Suppose the graph associated with *Li* and *Fi*, *i* = 1, 2, …, *N* be determined by Eqs. (8) and (9) with *si* = *s* ∈ (−1, 1). Let Φ be defined by Eq. (5). Then, Φ: *Ĉ*(*D*) → *Ĉ*(*D*) is well defined and contractive in the sup‐norm ‖⋅‖_{∞}. Furthermore, for every *j* = 0, 1, …, *m*−1 and *f* ∈ *Ĉ*(*D*), Φ(*f*)(*qj*) = *zj*.

We call the unique attractor of the afore‐mentioned IFS a *self‐affine fractal interpolation surface*, or *SAFIS* for short, *with arbitrary boundary data and identical contractivity factors*. The colouring of Theorem 4.1.2 is also known as ‘consistent triangulation’ within the context of computational geometry. We will be using this term from now on, for simplicity. Figure 7 illustrates Example 2 constructed in Ref. [4].

#### 4.1.3. Arbitrary boundary data, arbitrary contractivity factors

Zhao in Ref. [7] deploys a piecewise linear function *γ* over the graph *G* of the IFS. Let *G* be the graph mentioned earlier and *S* its vertices. We assign a value *s* ∈ (−1, 1) to each vertex by defining for all *i* = 1, 2, …, *N* the piecewise linear function *γ* = *γi*: ^{2}→ (−1, 1) as follows

where *βs* is the barycentric interpolation function with the attribute *s* as coefficient

**Theorem 4.1.3.** Consider the IFS {*S*; *w*_{1−N}} and the corresponding graph *G* described above that also fulfils the conditions mentioned in Theorem 4.1.2, except for the usage of identical scaling factors over the whole surface. For each transformation *wi*, we substitute the scaling factors *si* with the function *γi* described above, so the map *Fi* becomes for all *i* = 1, 2, …, *N*

Then, the unique attractor of the IFS mentioned previously is the graph of a continuous function *f* that passes though the points *j* = 0, 1, …, *m*−1.

We call the unique attractor of the afore‐mentioned IFS a *self‐similar fractal interpolation surface*, or *SSFIS* for short, *with arbitrary boundary data and contractivity factors*. The assignment of the scaling factor value *sj*, *j* = 0, 1, …, *m*−1 to each vertex can be done either during the parameter identification process or can be inferred via a given set of vertical scaling factors {*s*_{1}, *s*_{2}, …, *sN*}, where each *si* corresponds to *Δ*_{i}, *i* = 1, 2, …, *N*. The technique is identical to the calculation of the vertex normal vectors on the Gouraud and Phong shading models, where the facets of the mesh are the sub‐regions *Δ*_{i}, *i* = 1, 2, …, *N*. Compare and contrast Figures 6 and 9 of Ref. [7] with Figure 8.

### 4.2. Recurrent affine fractal interpolation

Let *D* be a polygonal domain, *D* consisting of non‐degenerate triangles with non‐intersecting interiors whose union is *D* and *S* = {*q*_{0}, *q*_{1}, …, *qm*_{−1}} be the set of vertices of *M* triangles each of which is a union of some subset of *S* so that the first *n* points {*q*_{0}, *q*_{1}, …, *qn*_{−1}} ⊂ *S* are the vertices of

#### 4.2.1. Coplanar boundary data, arbitrary contractivity factors

We call *l*: {1, 2, …,*N*} → {1, 2, …, *M*} a *Δ‐labelling* which associates the vertices of *Δi* with the vertices of some *δk*. If *Pk* are non‐vertical planes in ^{3}, let *Ĉ*(*D*) denote the collection of continuous functions *fk* = *f*_{Δi}: *δk*→ *x*, *fk*(*x*)) ∈ *Pk* for all *x* ∈ ∂*δk*.

**Theorem 4.2.1.** Suppose the points {(*qj*, *zj*) : *qj* ∈ ∂*δk*} are contained in the planes *Pk* ⊂ ^{3} for every *k*. Let Φ be defined by Eq. (5), where *Li* and *Fi*, *i* = 1, 2, …, *N* are defined by Eqs. (6)–(8). Then, Φ: *Ĉ*(*D*) → *Ĉ*(*D*) is well defined and contractive in the sup‐norm ‖⋅‖_{∞}. Furthermore, for every *j* = 0, 1, …, *m*−1 and *f* ∈ *Ĉ*(*D*), Φ(*f*)(*qj*) = *zj*.

We call the unique attractor of the afore‐mentioned RIFS a *recurrent self‐affine fractal interpolation surface*, or *RSAFIS* for short, *with coplanar boundary data and arbitrary contractivity factors*.

#### 4.2.2. Arbitrary boundary data, identical contractivity factors

We call *l*: {0, 1, …, *m*−1} → {0, 1, …, *n*−1} a *Δ‐labelling* associated with *ql*_{(j)}, *ql*_{(j')}, *ql*_{(j”)}} are the vertices of some *δk* whenever {*qj*, *qj'*, *qj”*}are the vertices of some *Δi*; see Figure 10.

**Theorem 4.2.2.** Let *l* be a *Δ*‐labelling associated with the triangulations *Li*: ^{2} → ^{2} and *Fi*: ^{3}→*si* = *s* (|*s*| < 1). Let Φ be defined by Eq. (5). Then, Φ: *Ĉ*(*D*)→*Ĉ*(*D*) is well defined and contractive in the sup‐norm ‖⋅‖_{∞}. Furthermore, for every *j* = 0, 1, …, *m*−1 and *f* ∈ *Ĉ*(*D*), Φ(*f*)(*qj*) = *zj*.

Geronimo and Hardin in Ref. [4] have showcased a two‐dimensional multi‐resolution analysis. Based on their work, we have created acceptable colourings of graphs with any desirable density, by solving the problem on a small instance that we call *pattern graph*. If the solution has the identity of self‐similarity, the density of the graph is then enhanced by interpolating the pattern with a RIFS defined on the metric space of the undirected graphs. With the help of a geometric predicate, we unify the resulting points into a single graph that now has the desired density. The results are illustrated in Figure 11.

We call the unique attractor of the afore‐mentioned RIFS a *recurrent self‐affine fractal interpolation surface*, or *RSAFIS* for short, *with arbitrary boundary data and identical contractivity factors*. Figure 12 illustrates Example 3 constructed in Ref. [4].

#### 4.2.3. Arbitrary boundary data, arbitrary contractivity factors

The map *Fi* becomes for all *i* = 1, 2, …, *N*,

**Theorem 4.2.3**. Suppose the graph associated with *Li* and *Fi*, *i* = 1, 2, …, *N* are defined by Eqs. (6)–(8) and Φ be defined by Eq. (5). Then, Φ: *Ĉ*(*D*)→*Ĉ*(*D*) is well defined and contractive in the sup‐norm ‖⋅‖_{∞}. Furthermore, Φ(*f*)(*qj*) = *zj*, *j* = 0, 1, …, *m*−1 and *f* ∈ *Ĉ*(*D*).

We call the unique attractor of the afore‐mentioned RIFS a *recurrent self‐similar fractal interpolation surface*, or *RSSFIS* for short, *with arbitrary boundary data and contractivity factors*. We illustrate this construction in Figure 13.

## 5. Conclusions, extensions and future work

We have presented an overview of affine interpolation as well as of recurrent affine interpolation using fractal functions in one and two dimensions. The effectiveness of a self‐affine fractal model is limited to the types of data: only those data that are self‐affine, or nearly so, are well represented. Since most data are not self‐affine, the piecewise or recurrent self‐affine fractal model may be a suitable alternative tool.

Nevertheless, our future work should focus on the parameter identification of fractal interpolation surfaces. More methods must be explored, including the ones for self‐affine FIS’s and the vertex‐based ones, while a better sampling technique for the heights is needed. For the recurrent FIS’s, it is utmost important to connect the domains with the subdomains that they resemble the most instead of determining the connections through an arbitrary colouring scheme. Therefore, a premise for future work is to find the optimal connections between domains and subdomains through a statistical analysis and then, starting from a feasible solution, to perform a heuristic local search with the intention to minimise the distance between current connections, implied by the current colouring, and those of the optimal solution. Moreover, methods to construct FIS, other than those based on the iterated function systems, exist, including the ones that use wavelets and tensor products, in which they should be also studied.

Many areas of fractal functions and their applications are yet to be explored, for instance, calculating the Hausdorff dimension of a general FIF is a challenging open problem. By believing that the field of fractal functions has bright future, the reader, in order to be able to contribute to it, should leave the idea of staying in the comfort of well‐known classical approximation theory and start enjoying the benefits of the more versatile fractal approximants, to supplement the former if not to replace it.