---
title: Fractal Interpolation Surfaces (FISs)
url: https://www.emergentmind.com/topics/fractal-interpolation-surfaces-fiss
type: topic
---

# Fractal Interpolation Surfaces (FISs)

Searching arXiv for recent and foundational papers on fractal interpolation surfaces.
Fractal interpolation surfaces (FISs) are graphs of continuous bivariate functions constructed to interpolate prescribed data while exhibiting recursive geometric structure. In the literature represented here, an FIS is typically realized as the attractor of an iterated function system (IFS) or recurrent iterated function system (RIFS) on a subset of $\mathbb{R}^3$, with the defining function satisfying a Read–Bajraktarević-type self-referential equation. The subject encompasses rectangular-grid constructions, recurrent and generalized affine variants, operator-theoretic formulations, calculi for derivatives and fractional integrals, dimension theory for box and Hausdorff dimensions, and specialized geometric settings such as triangular-domain “wedding cake surfaces” [1810.09701].

## 1. Definitions and core construction paradigms

On rectangular grids, the basic interpolation problem is to find a continuous function $f:I\times J\to\mathbb{R}$ satisfying $f(x_i,y_j)=Z_{i,j}$ for prescribed data $\{(x_i,y_j,Z_{i,j})\}$, where $I=[x_0,x_N]$ and $J=[y_0,y_M]$ are compact intervals with ordered partitions. A standard construction chooses contracting horizontal maps $u_i:I\to I_i=[x_{i-1},x_i]$ and $v_j:J\to J_j=[y_{j-1},y_j]$, together with continuous vertical maps $F_{i,j}:I\times J\times\mathbb{R}\to\mathbb{R}$ satisfying interpolation boundary conditions and a uniform contraction condition in the last variable. The maps
\[
W_{i,j}(x,y,z)=(u_i(x),v_j(y),F_{i,j}(x,y,z))
\]
then define an IFS on $K=I\times J\times\mathbb{R}$, and Hutchinson’s theorem yields a unique attractor that is the graph of a continuous interpolation function [1810.09701].

Equivalently, the interpolant is the unique fixed point of a Read–Bajraktarević operator acting on a suitable space of continuous functions with prescribed nodal values. In the rectangular-grid setting this fixed-point form is commonly written piecewise as
\[
(Tg)(x,y)=F_{i,j}(u_i^{-1}(x),v_j^{-1}(y),g(u_i^{-1}(x),v_j^{-1}(y))),
\]
for $(x,y)\in I_i\times J_j$. The resulting function satisfies a self-referential equation on each cell; this formulation is central both for existence theory and for later operator-theoretic and approximation analyses [2005.06241].

Several variants of the basic scheme appear in the literature. Recurrent fractal interpolation surfaces (RFISs) replace the finite independent-cell IFS by a recurrent system encoded by a row-stochastic matrix and domain-to-region assignments. In that setting, each map may originate from an expanded domain $D'_{ij}$ rather than the whole rectangle, and the graph of the interpolant is the unique invariant set of the recurrent system [1902.01165]. Another generalization treats countably infinite data on a rectangular grid, producing countable fractal interpolation surfaces (CFISs) by a countable IFS that remains hyperbolic in an equivalent metric [2010.05467].

A different constructive direction builds surfaces from one-dimensional recurrent fractal interpolation curves combined with Lipschitz mixing functions. In that framework, one defines
\[
F(x,y)=a(x,y)\,r_x(x)+b(x,y)\,r_y(y)+c(x,y),
\]
or more generally
\[
F(x,y)=\sum_{i=1}^N \lambda_i(x,y)\,f_i(x)+\sum_{j=1}^M \mu_j(x,y)\,g_j(y),
\]
thereby obtaining a broad class of continuous fractal surfaces without using a full two-dimensional IFS as the primary construction mechanism [1303.0615].

## 2. Self-reference, scaling structure, and interpolation constraints

The defining mechanism of an FIS is the coupling of horizontal contraction with vertical scaling. In generalized affine rectangular constructions, one chooses a continuous vertical-scaling function $S:D\to(-1,1)$, continuous functions $g,h:D\to\mathbb{R}$ satisfying interpolation constraints, and defines
\[
F_\omega(x,y,z)=S(L_\omega(x,y))\,(z-g(x,y))+h(L_\omega(x,y)),
\]
with $\Psi_\omega(x,y,z)=(u_i(x),v_j(y),F_\omega(x,y,z))$. The attractor is then the graph $\Gamma f$ of the unique continuous function interpolating the data [2312.15192].

The same structural idea appears in earlier rectangular-grid formulations using a vertical contraction factor function $s:E\to(0,1)$ and maps
\[
F_{ij}(x,y,z)=s(L_{ij}(x,y))\,(z-g(x,y))+h(L_{ij}(x,y)).
\]
Here $g$ and $h$ are continuous Lipschitz functions, with $g$ matching corner data and $h$ matching all grid data. The attractor is again the graph of a continuous interpolation function [1208.2081].

Function-valued scaling factors were subsequently emphasized in constructions where the vertical scaling vanishes on the cell boundaries. In that setting, each map has the form
\[
w_{ij}(x,y,z)=\bigl(L_{x_i}(x),L_{y_j}(y),S_{ij}(L_{x_i}(x),L_{y_j}(y))\,z+D_{ij}(x,y)\bigr),
\]
with $S_{ij}:E_{ij}\to(-1,1)$ satisfying
\[
S_{ij}(x_{i-1},v)=S_{ij}(x_i,v)=S_{ij}(u,y_{j-1})=S_{ij}(u,y_j)=0.
\]
Because $S_{ij}=0$ on each edge of $E_{ij}$, the prescribed boundary curves and grid data are preserved without tearing, and the attractor is a graph [1404.1300]. This boundary-vanishing condition is one of the clearest mechanisms for enforcing continuity across cell interfaces.

In the bivariate $\alpha$-fractal operator framework, the vertical scaling is encoded by a scale function $\alpha\in C(I\times J)$ with $\|\alpha\|_\infty<1$, and for a seed function $f$ and a bounded linear operator $L$ preserving grid values one defines
\[
F_{i,j}(x,y,z)=\alpha(u_i(x),v_j(y))\,z+f(u_i(x),v_j(y))-\alpha(u_i(x),v_j(y))\,(Lf)(x,y).
\]
The corresponding fractal function $f_{\alpha,L}$ satisfies
\[
f_{\alpha,L}(x,y)=f(x,y)+\alpha(x,y)\,[\,f_{\alpha,L}(u_i^{-1}(x),v_j^{-1}(y))-(Lf)(u_i^{-1}(x),v_j^{-1}(y))\,]
\]
on each cell. This recasts fractal interpolation as an operator-induced perturbation of a prescribed function rather than solely as a data-interpolation scheme [1810.09701].

A related nonlinear formulation for countable grids replaces the linear operator $L$ by a linear or nonlinear operator reproducing the boundary data and uses a family of scale functions $\{\varphi_{ij}\}$ with $\sup|\varphi_{ij}|<1$. The resulting map $f\mapsto f^{\varphi,L}$ generates a parameterized family of bivariate fractal functions that interpolate the seed function at the grid points while deviating from it in a controlled self-referential way [2010.05467].

## 3. Operator-theoretic formulations and approximation theory

The bivariate fractal operator formalism places FISs within functional analysis. For $f\in C(I\times J)$, the operator
\[
F_{\alpha,L}:C(I\times J)\to C(I\times J),\qquad f\mapsto f_{\alpha,L}
\]
is called the bivariate $\alpha$-fractal operator. It is a bounded linear operator on $C(I\times J)$, and several norm estimates are available. Writing $A=\|\alpha\|_\infty<1$ and $M=1+\|Id-L\|$, one has
\[
\|f_{\alpha,L}-f\|_\infty\le \frac{A}{1-A}\,\|Id-L\|\,\|f\|_\infty,
\]
and
\[
\|F_{\alpha,L}\|\le 1+\frac{A\|Id-L\|}{1-A}.
\]
Under the small-scale bound $A<1/\|Id-L\|$, the operator is one-to-one and has closed range; under the stronger condition $A<1/(1+\|Id-L\|)$, it is a topological automorphism of $C(I\times J)$ [1810.09701].

The same work extends the operator to $L^p(I\times J)$ for $1\le p<\infty$ by complexification and density, giving a bounded linear operator with perturbation estimate
\[
\|f_{\alpha,L}-f\|_p\le \frac{A}{1-A}\,\|Id-L\|\,\|f\|_p.
\]
If $A\cdot M^2<1$, then $F_{\alpha,L}$ is invertible on $L^p$ and hence Fredholm of index $0$ [1810.09701].

A multivariate extension on $I^k$ yields a related operator $\mathcal F^\alpha_{\Delta,D}$ associated with a partition net $\Delta$, a bounded operator $D:C(I^k)\to C(I^k)$, and a scale function $\alpha\in C(I^k)$ with $\|\alpha\|_\infty<1$. On each sub-rectangle, the corresponding Read–Bajraktarević operator is a contraction in both the uniform norm and the $L^p$ norm, with contraction constant $\|\alpha\|_\infty$. The associated approximation estimate is
\[
\|f_{\Delta,D}^\alpha-f\|_\infty\le \frac{\|\alpha\|_\infty}{1-\|\alpha\|_\infty}\,\|f-Df\|_\infty.
\]
This shows that the fractal perturbation remains quantitatively controlled by the approximation quality of the chosen base operator [2310.12276].

Approximation theory also appears through finite-dimensional classes of fractal polynomials. If $P_{m,n}$ denotes the space of bivariate polynomials of partial degrees $\le m,n$, then
\[
P_{m,n}^\alpha=F_{\alpha,L}(P_{m,n})
\]
is finite-dimensional of the same dimension as $P_{m,n}$ when $\alpha$ is small enough. Every $f\in C(I\times J)$ admits a best approximant in $P_{m,n}^\alpha$, and the union of these classes is dense in $C(I\times J)$ under suitable smallness assumptions on $\|\alpha\|_\infty$ [1810.09701].

These operator frameworks suggest that FISs are not only interpolation objects but also structured perturbations of classical approximation schemes. A plausible implication is that the choice of $L$ or $D$ determines how much of the classical approximation architecture is retained after fractalization.

## 4. Regularity, calculus, and constrained geometric design

The calculus of bivariate fractal interpolation functions has been developed for partial derivatives, partial integrals, mixed Riemann–Liouville fractional integrals, mixed Riemann–Liouville fractional derivatives, and integral transforms. For a rectangular FIS defined by
\[
f(x,y)=\alpha_{ij}\,f(u_i^{-1}(x),v_j^{-1}(y))+q_{ij}(u_i^{-1}(x),v_j^{-1}(y)),
\]
differentiation on a cell yields
\[
\partial_x f(x,y)=\frac{\alpha_{ij}}{a_i}\,\partial_x f(u_i^{-1}(x),v_j^{-1}(y))+\frac{1}{a_i}\,\partial_x q_{ij}(u_i^{-1}(x),v_j^{-1}(y)),
\]
and similarly for $\partial_y f$. Under the contraction condition $|\alpha_{ij}|<\min(a_i,c_j)$ and appropriate smooth matching of the $q_{ij}$ on cell boundaries, one obtains $f\in C^1$ or $C^2$ according to the smoothness of the seed data [2005.06241].

Partial integrals preserve the self-referential structure. If
\[
I_x f(x,y)=\int_a^x f(s,y)\,ds,
\]
then
\[
I_x f(u_i(x),v_j(y))=\alpha_{ij}\,a_i\,I_x f(x,y)+\widehat q_{ij}(x,y),
\]
so $I_x f$ is again a fractal interpolation function, now with vertical scaling $\alpha_{ij}a_i$. The mixed Riemann–Liouville fractional integral of order $\gamma=(p,q)$ satisfies
\[
I^{(p,q)}f(u_i(x),v_j(y))=\alpha_{ij}\,a_i^p c_j^q\,I^{(p,q)}f(x,y)+\widehat q_{ij}^{(p,q)}(x,y),
\]
which shows that mixed fractional integrals again belong to the same general recursive class [2005.06241].

A more design-oriented development is the bicubic partially blended rational quartic surface. There, rational quartic fractal interpolation functions are first constructed along horizontal and vertical grid lines, using scaling factors $\alpha_{i,j}$, $\beta_{i,j}$ and shape parameters. On each rectangular cell, these boundary curves are blended by a bicubic Coons-type patch. The full surface $S$ is assembled cellwise, and by construction $S\in C^1(K)$ because each boundary FIF is $C^1$ and the blending functions enforce first-derivative matching across cell boundaries [1910.09822].

That construction also includes explicit positivity and monotonicity constraints. The scaling factors must satisfy $|\alpha_n|<a_n<1$, and additional inequalities guarantee that a univariate rational quartic fractal interpolant lies in a prescribed vertical strip or stays above a line. In the surface setting, analogous inequalities are imposed separately on the horizontal and vertical boundary FIFs so that they lie above a plane of the form
\[
t=c\Bigl[1-\frac{x-x_1}{x_m-x_1}-\frac{y-y_1}{y_n-y_1}\Bigr].
\]
The resulting surface interpolates nodal values and prescribed first partial derivatives at the grid nodes [1910.09822].

## 5. Dimension theory on rectangular grids

A major line of work studies the box-counting dimension of FIS graphs. For rectangular-grid surfaces with vertical contraction factor function $s$, one defines
\[
S_{ij}=\max_{(x,y)\in E_{ij}}|s(x,y)|,\qquad s_{ij}=\min_{(x,y)\in E_{ij}}|s(x,y)|.
\]
From these values one forms nonnegative matrices whose spectral radii control the dimension estimates. If the boundary interpolation points are not all collinear and $\underline a>n$, then
\[
1+\frac{\log(\underline a)}{\log n}\le \dim_B A\le 1+\frac{\log(\bar a)}{\log n},
\]
while if $\bar a\le n$ then $\dim_B A=2$ [1208.2081].

In the constant-scale case $s(x,y)\equiv s$, one has $\bar a=\underline a=n^2 s$, and whenever $s>1/n$,
\[
\dim_B A=1+\frac{\ln(n^2 s)}{\ln n}=3+\frac{\ln s}{\ln n}.
\]
Thus the dimension ranges between $2$ and $3$ as $s$ varies from $1/n$ to values approaching $1$ [1208.2081].

For recurrent fractal interpolation surfaces on uniform grids, the dimension problem is formulated in terms of a compatible partition and uniform sums of the vertical scaling factors. One defines a nonnegative matrix
\[
G=(\gamma_{rt})_{r,t=1,\dots,m},
\]
decomposes it into irreducible components, and for each non-degenerate component $V$ sets
\[
d_V=\frac{\log \rho(G|_V)}{\log K}.
\]
With
\[
d^*=\max\{1,\max_{V\ \text{nondegenerate}} d_V\},
\]
the theorem states that
\[
\dim_B \Gamma f = 1+d^*.
\]
If $G$ is irreducible and nondegenerate, then $\dim_B\Gamma f=1+(\log\rho(G))/\log K$ when $\rho(G)>K$, and $\dim_B\Gamma f=2$ otherwise [1902.01165].

General recurrent constructions with function vertical scaling factors produce lower and upper bounds in terms of the spectral radii of $\underline S\cdot C$ and $\overline S\cdot C$, where $C$ is the connection matrix of the RIFS. After a bi-Lipschitz rescaling to an equal mesh, if $\underline\lambda>a$ then
\[
1+\log_a(\underline\lambda)\le \dim_B(\mathcal A)\le 1+\log_a(\overline\lambda),
\]
where $\underline\lambda$ and $\overline\lambda$ are the spectral radii of the two products; if $\overline\lambda\le a$, then $\dim_B(\mathcal A)=2$ [1307.3229].

The most systematic recent rectangular theory uses oscillation sums and vertical-scaling matrices. For a continuous function $\varphi$ on $D$, the level-$k$ oscillation sum is
\[
S_k(\varphi)=\sum_{\omega\in\Sigma^k} O(\varphi,D_\omega),
\]
and the graph dimension satisfies
\[
1+\liminf_{k\to\infty}\frac{\log(S_k(\varphi)+N^k)}{k\log N}\le \underline{\dim}_B\Gamma\varphi\le \overline{\dim}_B\Gamma\varphi\le 1+\limsup_{k\to\infty}\frac{\log(S_k(\varphi)+N^k)}{k\log N}.
\]
For generalized affine FISs, upper and lower vertical-scaling matrices $\overline V_n$ and $\underline V_n$ induce vector recurrences for oscillation vectors; in the sign-preserving case one obtains a common limit $\rho_S$ of their spectral radii and the exact formula
\[
\dim_B\Gamma f=
\begin{cases}
1+\dfrac{\log\rho_S}{\log N}, & \rho_S>N,\\[4pt]
2, & \rho_S\le N.
\end{cases}
\]
In the constant-scale case, this recovers $\rho_S=N^2|s|$ and hence
\[
\dim_B\Gamma f=1+\frac{\log(N^2|s|)}{\log N}=3+\frac{\log|s|}{\log N}
\]
when $|s|>1/N$ [2312.15192].

These results collectively show a recurring threshold phenomenon: when the effective cumulative vertical scaling exceeds the horizontal subdivision rate, the graph dimension rises above $2$; otherwise the graph has the “trivial” surface dimension $2$. That interpretation is directly supported by multiple formulations, though each paper encodes the threshold through a different matrix or spectral quantity.

## 6. Geometric variants, coverings, and specialized surfaces

Not all FISs are based on rectangular grids. A recent development studies fractal interpolation surfaces over a triangular domain, called “wedding cake surfaces.” These surfaces are attractors of deterministic self-affine IFSs on $\mathbb{R}^3$ generated by a fractal interpolation algorithm. Their dimension theory differs from the rectangular box-dimension setting because the associated self-affine IFSs are not strongly irreducible, so the recent general theory for strongly irreducible and proximal self-affine systems does not apply directly. The main result is that the Hausdorff dimension of the self-affine set, or FIS, is the same as the affinity dimension outside a set of scaling parameters with zero Lebesgue measure; moreover, by computing the overlapping number for the associated Furstenberg IFS, the Hausdorff dimension is determined for every type of scaling parameter in a certain range of parameters [2510.04613].

This triangular-domain work indicates that FIS dimension theory is not exhausted by box-counting methods on Cartesian grids. A plausible implication is that the interpolation geometry of the domain may force fundamentally different irreducibility and overlap phenomena in the underlying linear cocycle.

Geometric localization has also been studied computationally. For a particular affine-plus-bilinear FIS framework, one can cover the attractor by a finite family of octahedra in the metric
\[
\rho((x,y,z),(x',y',z'))=|x-x'|+|y-y'|+\theta |z-z'|.
\]
If $\gamma_{k,l}$ denotes the fixed point of the $(k,l)$-th map and $C_{k,l}$ its Lipschitz constant, then radii $\rho_{k,l}$ solving a max-type system produce a cover
\[
G_f\subseteq \bigcup_{k=1}^n\bigcup_{l=1}^m O[\gamma_{k,l},\rho_{k,l}],
\]
where each $O[\gamma,r]$ is an octahedron. The method is described as computationally more efficient than a previous ball-cover procedure because it is based on finding the maximum of certain sets rather than using a sorting algorithm [2311.09512].

The covering construction does not change the underlying interpolation theory, but it clarifies that FIS research includes not only existence and dimension results but also localization algorithms for attractors. This suggests a computational strand of the subject concerned with certified enclosures and numerical geometry.

## 7. Relations to neighboring constructions and common misconceptions

A frequent misconception is that FISs are synonymous with bilinear fractal surfaces. In fact, bilinear examples are only one subclass. The literature includes generalized affine FISs with continuous vertical-scaling functions [2312.15192], recurrent surfaces with function vertical scaling factors [1307.3229], countable-grid constructions [2010.05467], surfaces assembled from one-dimensional recurrent fractal interpolation curves [1303.0615], bicubic partially blended rational quartic surfaces with $C^1$ regularity [1910.09822], and self-affine triangular-domain “wedding cake surfaces” studied via Hausdorff dimension rather than only box dimension [2510.04613].

Another misconception is that “fractal” necessarily implies graph dimension strictly larger than $2$. Several sources explicitly identify parameter regimes where the graph dimension is exactly $2$. For rectangular-grid models, this occurs when the relevant cumulative or spectral measure of vertical scaling does not exceed the subdivision threshold, such as $\bar a\le n$ [1208.2081], $\rho(G)\le K$ in the recurrent bilinear setting [1902.01165], or $\rho_S\le N$ in the vertical-scaling-matrix formulation [2312.15192].

It is also inaccurate to regard FISs as purely geometric constructions without analytic structure. Operator-theoretic studies show that the associated fractal operators can be bounded, injective, invertible, Fredholm of index $0$, extendable to $L^p$ spaces, and suitable for best approximation and Schauder-basis transport [1810.09701]. Likewise, calculus results show that partial integrals, partial derivatives, and mixed Riemann–Liouville fractional operators preserve a self-referential interpolation structure under explicit transformed scaling factors [2005.06241].

Finally, the role of the vertical scaling factor should not be reduced to a single roughness parameter. In some papers it is a constant, in others a continuous function, a cellwise bilinear function, a scale field interpolating nodal values, or a family of scaling functions attached to a recurrent system. Across these formulations, the scaling data determine contraction, continuity constraints, roughness, and dimension estimates; but the precise effect depends on the chosen construction. This suggests that “vertical scaling” is better understood as a structural component of the recursive model rather than as a one-parameter measure of irregularity alone.

Source: https://www.emergentmind.com/topics/fractal-interpolation-surfaces-fiss