\section{Introduction}

Over the last decades, fractal interpolation and approximation have been extensively research. This research originated with \cite{B2} where a special set-up was used to define so-called affine fractal interpolation functions. The graphs of these affine fractal interpolation functions are the attractors of a class of iterated function systems and thus geometrically motivated. An analytic construction of general fractal functions originated in \cite{bedford,dubuc1,dubuc2} where the concept of a Read-Bajractarevi\'c operator is first encountered. Numerous constructions of fractal functions based on Read-Bajractarevi\'c operators satisfying given interpolation and approximation conditions followed. Some of these constructions are introduced and summarized in \cite{massopust,massopust1}. The number of publications in fractal interpolation theory is enormous and the interested reader may want to search for fractal functions using terms such as ``hidden variable,'' ``$V$-variable,'' ``coalescent,'' ``super,'' and ``$\alpha$-fractal functions'' to name just a few.

The intend of this chapter is to introduce the reader to the concept of fractal interpolation and its extensions from a global to a local setting, then to non-stationarity and finally to a quaternionic setting. In a certain sense, these are the main set-ups with the possible exclusion of \emph{unbounded} fractal interpolation \cite{m7}. It is understandable that such an endeavor must necessarily restrict itself to the main points of each construction and setting. However, the exposition will give the reader an overall perspective of the issues involved and the techniques used, and can be used as a starting point for a deeper investigation into each the topics.

The outline of this chapter is as follows. Section 2 introduces the global setting of fractal interpolation and exhibits a relationship between the Read-Bajractarevi\'c operator and the solution to a canonically associated system of functional equations. In the next section, we extend global interpolation to a local set-up giving more flexibility to the construction. This type of interpolation found deep applications to fractal imaging and fractal compression \cite{BH}. In Section 4, the recently introduced concept of non-stationary fractal interpolation is presented and it is shown that backward trajectories allow distinct features to be delineated at different interpolation scales. The final Section 5, describes the novel setting of fractal interpolation in the theory of quaternions and shows that the non-commutative character of quaternions introduces even more intricate fractal patterns.

\section{Global Fractal Interpolation}\label{GFI}

The purpose of the current section is to introduce global fractal interpolation and to relate the global fractal interpolant to the solution of a system of functional equations. We see that this system of functional equations defines in a canonical way a Read-Bajractarevi\'c (RB) operator and vice versa. This relationship will be encountered several times in the subsequent sections as well.

In the following, $(\mathsf{E},d_\mathsf{E})$ denotes a normed space and $(\mathsf{F},d_\mathsf{F})$ a Banach space. For $n\in \mathbb{N}$, we write $\mathbb{N}_n :=\{1,\ldots, n\}$ for the initial segment of the natural numbers $\mathbb{N}$ of length $n$.

For a given normed space $(\mathsf{E}, \n{\cdot}_\mathsf{E})$ and a map $f: \mathsf{E} \to\mathsf{E}$, we define the Lipschitz constant associated with $f$ by
\[
\Lip (f) := \sup_{x,y \in \mathsf{E}, x \neq y} \frac{\n{f(x)-f(y)}_\mathsf{E}}{\n{x-y}_\mathsf{E}}.
\]
The map $f$ is called Lipschitz if $\Lip (f) < + \infty$ and a contraction (on $\mathsf{E}$) if $\Lip (f) < 1$.

\subsection{Bounded Solutions}
Let $\mathsf{X}$ be a nonempty bounded subset of $\mathsf{E}$. Suppose we are given a finite family $\{l_i\}_{i = 1}^{n}$ of injective contractions $\mathsf{X}\to \mathsf{X}$ generating a partition of $\mathsf{X}$ in the sense that
\begin{align}
&\forall\;i, j\in \mathbb{N}_n, i\neq j: l_i(\mathsf{X})\cap l_j(\mathsf{X}) = \emptyset;\label{c1}\\
&\mathsf{X} = \bigcup_{i=1}^n l_i(\mathsf{X}).\label{c2}
\end{align}
For simplicity, we write $\mathsf{X}_i := l_i(\mathsf{X})$.

Given the above set-up, we are looking for a global function $\psi:\mathsf{X} = \bigcup\limits_{i=1}^n \mathsf{X}_i\to\mathsf{F}$ satisfying $n$ functional equations of the form
\begin{equation}\label{psieq}
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad\text{on $\mathsf{X}$ and for $i\in \mathbb{N}_n$},
\end{equation}
where for each $i\in\mathbb{N}_n$, $s_i$ is a given bounded function $\mathsf{X}\to\mathbb{R}$ and $q_i$ a bounded function $\mathsf{X}\to\mathsf{F}$. Recall that a function $f:\mathsf{X}\to\mathsf{F}$ is called \emph{bounded} if there exists a finite $M>0$ such that $\n{f(x)}_\mathsf{F} \leq M$, for all $x\in \mathsf{X}$.

The idea is to consider \eqref{psieq} as the fixed point equation for an associated affine operator acting on an appropriately defined function space.

To this end, let $\mathcal{B} (\mathsf{X},\mathsf{F}) := \{f:\mathsf{X}\to\mathsf{F} : \text{$f$ is bounded}\}$ denote the the Banach space of bounded functions equipped with the supremums norm $\n{f} := \sup\limits_{x\in \mathsf{X}} \n{f(x)}_\mathsf{F}$.

On the Banach space $\mathcal{B} (\mathsf{X},\mathsf{F})$, we define an affine operator $T: \mathcal{B} (\mathsf{X},\mathsf{F})\to \mathcal{B} (\mathsf{X},\mathsf{F})$, called a Read-Bajractarevi\'c (RB) operator, by
\begin{equation}\label{eq3.17}
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x),
\end{equation}
for $x\in \mathsf{X}_i$ and $i\in \mathbb{N}_n$, or, equivalently, by
\begin{align*}
T f (x) &= \sum_{i=1}^n (q_i\circ l_i^{-1})(x)\, \chi_{\mathsf{X}_i}(x) + \sum_{i=1}^n (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x)\, \chi_{\mathsf{X}_i}(x)\\
&= T(0) + \sum_{i=1}^n (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x)\, \chi_{\mathsf{X}_i}(x),\quad x\in \mathsf{X},
\end{align*}
where $\chi_S$ denotes the characteristic function of a set $S$: $\chi_S(x) = 1$, if $x\in S$, and $\chi_S(x) = 0$, otherwise.

The following result is well-known (see, for instance, \cite{B2,massopust1}) but for the sake of completeness we reproduce the proof. We also refer the interested reader to \cite{SB} where a similar set-up is considered.

\begin{theorem}\label{sol}
The system of functional equations \eqref{psieq} has a unique bounded solution $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item conditions \eqref{c1} and \eqref{c2} are satisfied, and
\item $s:= \max\limits_{i\in \mathbb{N}_n} \sup\limits_{x\in \mathsf{X}} |s_i(x)| < 1$.
\end{enumerate}
\end{theorem}

\begin{proof}
First note that, as the mappings $l_i$ are injective, the right-hand side of \eqref{psieq} can be written as the right-hand side of \eqref{eq3.17}.

As the functions $l_i$, $q_i$, and $s_i$ are all assumed to be bounded, $T$ maps $\mathcal{B} (\mathsf{X},\mathsf{F})$ into itself. For all $f,g\in \mathcal{B} (\mathsf{X},\mathsf{F})$, we have that
\begin{align*}
\sup_{x\in\mathsf{X}}\n{Tf(x) - Tg(x)}_\mathsf{F} &= \max_{i\in \mathbb{N}_n}\sup_{x\in\mathsf{X}_i} \n{(s_i\circ l_i^{-1})(x)\cdot (f-g)\circ l_i^{-1}(x)}_\mathsf{F}\\
&= \max_{i\in \mathbb{N}_n}\sup_{\xi\in\mathsf{X}} \n{s_i(\xi)\cdot (f-g)(\xi)}_\mathsf{F}\\
&\leq \max_{i\in \mathbb{N}_n}\sup\limits_{x\in \mathsf{X}} |s_i(x)| \sup_{x\in\mathsf{X}}\n{(f-g)(x)}_\mathsf{F},
\end{align*}
from which it follows that
\[
\n{Tf-Tg} \leq s \n{f-g}.
\]
Hence, $T$ is contractive on the Banach space $\mathcal{B} (\mathsf{X},\mathsf{F})$ and therefore, by the Banach Fixed Point Theorem, has a unique fixed point $\psi\in\mathcal{B} (\mathsf{X},\mathsf{F})$. This fixed point solves the functional equations \eqref{psieq}.
\end{proof}

\begin{remarks}\hfill
\begin{enumerate}
\item The fixed point $\psi\in\mathcal{B} (\mathsf{X},\mathsf{F})$ of the RB operator $T$ is also called a \emph{bounded fractal function}. In this context, Eqn. \eqref{psieq} is also referred to as a \emph{self-referential equation} for $\psi$.
\item The self-referential equation $T\psi = \psi$ expresses the fractal nature of the $\gr\psi$: It is made up of a finite number of copies of itself with each copy being supported on the partitioning sets $\mathsf{X}_i$. Hence, the terminology \emph{fractal function} for $\psi$.
\item The proof of Banach's Fixed Point Theorem also provides an algorithm for the construction of $\psi$: Choose \emph{any} function $\psi_0\in \mathcal{B} (\mathsf{X},\mathsf{F})$ and iteratively define the following sequence of functions:
$$
\psi_k := T \psi_{k-1}, \quad k\in \mathbb{N}.
$$
Then, $\psi$ is given by $\psi = \lim\limits_{k\to\infty}\psi_k$ where the limit is taking with respect to the norm $\n{\cdot}$ on  $\mathcal{B} (\mathsf{X},\mathsf{F})$.
\item The afore-mentioned algorithm for the construction of $\psi$ together with the proof of the Banach Fixed Point Theorem gives an error estimate as well, namely,
\[
\n{\psi - \psi_{k}} \leq \frac{s^k}{1-s} \n{\psi_1 - \psi_0}, \quad k\in \mathbb{N}.
\]
\item The fixed point $\psi$ depends on $n$, the partition $(\mathsf{X}_i : i\in\mathbb{N}_n)$, and the functions $s_i$ and $q_i$ with different choices yielding different fractal functions.
\item Emphasizing the dependence of $\psi$ on the functions $s_i$, the expression $s$-fractal function can be found in the literature. (See, for instance, \cite{N}.) In this context, one considers a fractal function as the image under an operator $\mathcal{F}^s$ associating with a given (non-fractal) function its fractal analogue.
\item Functional equations such as \eqref{psieq} exhibit connections to so-called \emph{fractels} \cite{bhm2,massopust1} and also to the approximation of rough functions \cite{BHVV}.
\end{enumerate}
\end{remarks}

Conditions \eqref{c1} and \eqref{c2} cannot be relaxed without adding some compatibility conditions to guarantee that the RB operator $T$ has the form given by Eqn. \eqref{eq3.17}. Should Eqn. \eqref{c1} not be satisfied, one would have to impose in our current setting the following compatibility conditions:
\begin{align}\label{wd}
\forall x_1,&x_2\in\mathsf{X}:\nonumber\\
& \; l_i (x_1) = l_j (x_2)\;\;\Longrightarrow\;\; q_i(x_1) + s_i (x_1) \psi (x_1) = q_j(x_2) + s_ j(x_2) \psi (x_2).
\end{align}
We refer to \cite{SB,SB1} for more details regarding this issue.\\

As an application of the above approach to obtain solutions to functional equations of the form \eqref{psieq} or equivalently finding the unique fixed points of the associated RB operator \eqref{eq3.17}, we provide the following example.
\begin{example}\label{ex1}
Let $\mathsf{E} := \mathbb{R} =: \mathsf{F}$ together with the Euclidean norm $\abs{\cdot}$. Further. let $\mathsf{X} := [0,1)\subset\mathsf{E}$. Assume that we are given two injective contractions $l_i:[0,1)\to [0,1)$, $i=1,2$, with $l_1(x) := \frac13 x$ and $l_2(x) := \frac23 x + \frac13$. Hence, $\mathsf{X}_1 = [0,\frac13)$ and $\mathsf{X}_2 = [\frac13,1)$. Clearly, $\mathsf{X} = \mathsf{X}_1\cup\mathsf{X}_2$ and $\mathsf{X}_1\cap\mathsf{X}_2 = \emptyset$.

Now choose $q_1(x) = -1$, $q_2(x) = {x}$, $s_1(x) = \frac12 \sin(x)$, and $s_2(x):=-\frac23 \cos(x)$. The system of functional equations and the associated RB operator read then
\[
\psi(\tfrac13 x) = -1 + \tfrac12 \sin(x)\psi(x)\quad\text{and}\quad\psi(\tfrac23 x + \tfrac13) = {x} - \tfrac23\cos(x)\psi(x),
\]
and
\[
Tf(x) = \begin{cases} -1 + \frac12\sin(3x) f(3x), & 0\leq x < \tfrac13;\\
{3x-1} - \frac23\cos(\tfrac12(3x-1))f(\tfrac12(3x-1)), & \tfrac13 \leq x < 1,
\end{cases}
\]
respectively.

As $s = \frac23 < 1$, $T$ is contractive. A plot of the solution, respectively, fixed point $\psi$, is shown in Fig. 1.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=7cm, height= 4cm]{ex1.pdf}
\caption{The solution/fixed point $\psi$.}\label{fig1}
\end{center}
\end{figure}
\end{example}
\subsection{$L^p$ solutions}
In the following, we set $\mathsf{X}\subset\mathsf{E}:=\mathbb{R}^m$ and $\mathsf{Y}:=\mathbb{R}^k$ where the Euclidean spaces $\mathbb{R}^m$ and $\mathbb{R}^k$ are endowed with the corresponding canonical Euclidean norms.

Recall that the (real) Lebesgue spaces $L^p (\mathsf{X}, \mathbb{R}^k)$, where $\mathsf{X}\subset \mathbb{R}^m$ is nonempty, are defined as consisting of (equivalence classes of) functions $f:\mathsf{X}\to\mathbb{R}^k$ for which
\[
\n{f}_p :=\begin{cases}\left(\displaystyle{\int_\mathsf{X} \n{f(x)}^p dx} \right)^{1/p}, & 1\leq p < \infty;\\ \\
\esssup_{x\in\mathsf{X}} \n{f(x)}, & p = \infty.
\end{cases}
\]
is finite. Here, $\n{f(x)} := \sqrt{\abs{f_1(x)}^2 + \cdots + \abs{f_k(x)}^2}$ with $f := (f_1, \ldots, f_k)$.

We ask under what conditions on the functions $q_i$ and $s_i$ the solution $\psi$ is an element of $L^p(\mathsf{X}, \mathbb{R}^k)$, for $1\leq p < \infty$ and a bounded nonempty set $\mathsf{X}\subset\mathbb{R}^m$.

To this end, note that in order for $\psi$ to be in $L^p(\mathsf{X}, \mathbb{R}^k)$, the RB operator $T$ must map $L^p(\mathsf{X}, \mathbb{R}^k)$ into itself. Therefore, the functions $q_i$ and $s_i$ must also be in $L^p(\mathsf{X}, \mathbb{R}^k)$. Moreover, $s_i$ needs to be in $L^\infty(\mathsf{X}, \mathbb{R}^k)$ for the product $s_i\cdot f$ to be in $L^p(\mathsf{X}, \mathbb{R}^k)$. Thus, as $\mathsf{X}$ is bounded it has finite measure and therefore $s_i\in L^\infty(\mathsf{X}, \mathbb{R}^k)$ implies that $s_i\in L^p(\mathsf{X}, \mathbb{R}^k)$ for all $1\leq p \leq \infty$.

Now it remains to be shown that the RB operator $T$ is contractive on $L^p(\mathsf{X}, \mathbb{R}^k)$. For this purpose, let $f,g\in L^p(\mathsf{X}, \mathbb{R}^k)$. Then, with $\mathsf{X}_i:= l_i(\mathsf{X})$,
\begin{align*}
\n{Tf - Tg}_p^p &= \int_\mathsf{X} \n{Tf (x) - Tg(x)}^p dx\\
&= \int_\mathsf{X} \n{\sum_{i=1}^n s_i(l_i^{-1}(x))\cdot (f-g)(l_i^{-1}(x)) \chi_{\mathsf{X}_i}(x)}^p dx\\
&\leq \sum_{i=1}^n \int_{\mathsf{X}_i} \abs{s_i(l_i^{-1}(x))}^p \n{(f-g)(l_i^{-1}(x))}^p dx\\
&= \sum_{i=1}^n \int_{\mathsf{X}} \abs{(l_i^{-1})'(x)} \abs{s_i(x)}^p \n{(f-g)((x)}^p dx.
\end{align*}

If, for all $i\in \mathbb{N}_n$, $\n{l_i^{-1})'}_\infty =: \lambda_i<\infty$ and $\n{s_i}_\infty = :s_i<\infty$, then
\[
\n{Tf - Tg}_p^p \leq \left(\sum_{i=1}^n \lambda_i s_i^p\right) \n{(f-g)((x)}_p^p,
\]
and $T$ is contractive provided that
\[
\sum_{i=1}^n \lambda_i s_i^p < 1.
\]
Hence, we arrived at the following result
\begin{theorem}\label{th3.2}
The system of functional equations
\[
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad\text{on $\mathsf{X}\subset\mathbb{R}^m$ and for $i\in \mathbb{N}_n$},
\]
has a unique solution $\psi\in L^p(\mathsf{X}, \mathbb{R}^k)$, $1\leq p \leq \infty$, respectively, the RB operator
\[
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f\circ l_i^{-1})(x),\quad{x\in\mathsf{X}_i},\;i\in \mathbb{N}_n,
\]
a unique fixed point $\psi\in L^p(\mathsf{X}, \mathbb{R}^k)$ provided that
\begin{enumerate}
\item $q_i\in L^p(\mathsf{X}, \mathbb{R}^k)$, $s_i\in L^\infty (\mathsf{X}, \mathbb{R}^k)$ and
\item $\sum\limits_{i=1}^n \lambda_i s_i^p < 1$, where $\lambda_i = \n{(l_i^{-1})'}_\infty $ and $s_i = \n{s_i}_\infty$.
\end{enumerate}
\end{theorem}

\begin{remark}
In a similar fashion, one can derive conditions such that the unique solutions/fixed points $\psi$ are elements of H\"older or Sobolev spaces. See, for instance, \cite{PRM,massopust,massopust1}.
\end{remark}
\subsection{Continuous Solutions}
So far, we only considered bounded solution/fixed points $\psi$ for \eqref{psieq} and \eqref{eq3.17}. However, in some instances, a continuous or even differentiable solution is required. We only present a result for continuous $\psi$ and make some remarks about how to obtain differential solutions.

\begin{theorem}\label{th3.3}
The system of functional equations \eqref{psieq} has a unique continuous solution $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item $\mathsf{X} = \bigcup\limits_{i=1}^n l_i(\mathsf{X})$,
\item the functions $l_i$, $q_i$, and $s_i$ are continuous,
\item and for all $i,j\in \mathbb{N}_n$ and $x_1, x_2\in X$:
\begin{equation}\label{joinupc}
\lim_{x\to x_1} f_j(x) = f_i(x_2)\;\;\Longrightarrow\;\; \lim_{x\to x_1} q_j(x) + s_j(x) \psi(x) = q_i(x_2) + s_i(x_2) \psi(x_2).
\end{equation}
\end{enumerate}
\end{theorem}
\begin{proof}
We refer the interested reader to \cite{SB} or \cite{massopust}. In the former reference, the proof follows the functional equation setting and in the latter the RB operator setting.
\end{proof}

\begin{example}

\end{example}
We connect up with the previous Example 1 but choose as $\mathsf{X} := [0,1]$. We modify the functions $q_i$ to be $q_1(x) := x$ and $q_2(x) := 1-x$. but keep $s_1$ and $s_2$ unchanged. Note that here we have $\mathsf{X}_1 \cap \mathsf{X}_2 = \{\frac13\}$ and we need to ensure that conditions \eqref{wd} are satisfied. In particular, as we have $l_1(1) = \frac13 = l_2(0)$, the following equality has to hold:
\begin{equation}\label{3.7}
q_1(1) + s_1(1)\psi(1) = q_2(0) + s_2(0)\psi(0).
\end{equation}
The functional equations \eqref{psieq} imply for $x\in \{0,1\}$
\[
\psi(0) = q_1 (0) + s_1(0) \psi(0)\quad\text{and}\quad \psi(1) = q_2 (1) + s_2(1)\psi(0),
\]
which gives the values of $\psi$ at the endpoints of $\mathsf{X}$:
\[
\psi(0) = \frac{q_1(0)}{1-s_1(0)}\quad\text{and}\quad \psi(1) = \frac{q_2(1)}{1-s_2(1)}.
\]
The validity of \eqref{3.7} guarantees the existence of a \emph{bounded} solution $\psi$ (since $s = \frac23 < 1$). As $q_1(0) = 0 = q_2(1)$, the solution $\psi$ also vanishes on the boundary of $\mathsf{X}$: $\psi(0) = 0 = \psi(1)$.

In order to obtain a \emph{continuous} solution, equations \eqref {joinup} must be satisfied. In our current setting, as all the functions involved are continuous on $\mathsf{X}$, we obtain
\[
\lim_{x\to 0-} q_2(x) + s_2(x)\psi(x) = q_1(1) + s_1(1)\psi(1),
\]
which is identical to \eqref{3.7}.

The solution/fixed point $\psi$ is therefore continuous. The graph of $\psi$ is depicted in Figure \ref{fig2}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=7cm, height= 4cm]{ex2.pdf}
\caption{A continuous solution/fixed point $\psi$.}\label{fig2}
\end{center}
\end{figure}

In \cite{SB}, a similar setting is considered.
\begin{remark}
For many applications, in particular the setting where $\mathsf{E} := \mathbb{R}$, the partition $(\mathsf{X}_i : i\in \mathbb{N}_n)$ is induced by a given finite set of data points $\{(x_i, y_i)\in \mathbb{R}\times\mathsf{F} : x_0 < x_1 < \cdots < x_n\}$ with $\mathsf{X}_i := [x_{i-1},x_i]$. The injective contractions $l_i$ are then affine mappings $[x_0,x_n]\to [x_{i-1},x_i]$. For $\mathsf{F} := \mathbb{R}$, this type of continuous fractal interpolation was first introduced in \cite{B2} and the fixed point $\psi$ was termed a \emph{fractal interpolation function} as it was also required that $y_j = \psi(x_j)$, $j\in\{0,1,\ldots, n\}$.
\end{remark}
\section{Local Fractal Interpolation}
In this section, we introduce a generalization of global fractal interpolation. Previously, we considered a fixed nonempty bounded subset $\mathsf{X}$ of a normed space $\mathsf{E}$ generating a partition of $\mathsf{X}$ via a finite family of injective contractions. Now, we replace the single subset $\mathsf{X}$ by a finite family of subsets $\mathsf{X}_i$ of $\mathsf{X}$.

More precisely, let $\{\mathsf{X}_i : i \in \mathbb{N}_n\}$ be a family of nonempty subsets of a fixed subset $\mathsf{X}$ of a normed space $\mathsf{E}$. Suppose $\{l_i\}_{i = 1}^{n}$ is a collections of injective mappings from $\mathsf{X}_i\to \mathsf{X}$ generating a partition of $\mathsf{X}$ in the sense that
\begin{align}
&\forall\;i, j\in \mathbb{N}_n, i\neq j: l_i(\mathsf{X}_i)\cap l_j(\mathsf{X}_j) = \emptyset;\label{cl1}\\
&\mathsf{X} = \bigcup_{i=1}^n l_i(\mathsf{X}_i).\label{cl2}
\end{align}
Note that the $l_i$ need not be contractive mappings here.
\begin{remark}
One can actually have $m<n$ subsets of $\mathsf{X}$ and still be able to define $n$ injective mappings satisfying \eqref{cl1} and \eqref{cl2}. In this case some of the injections $l_i$ share the same domain $X_j$, i.e., $\mathsf{X}_j$ is repeated a certain number of times $\mathsf{X}_j = \mathsf{X}_{j+1} = \cdots = \mathsf{X}_{j+r}$. This situation occurs in Section \ref{sect4.3}.
\end{remark}
Local fractal interpolation looks for local solutions $\psi: \mathsf{X} = \bigcup\limits_{i\in \mathbb{N}_n} l_i (\mathsf{X}_i)\to \mathsf{F}$ of functional equations or for fixed points of RB operators of the form
\begin{equation}\label{localpsi}
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad\text{$x\in\mathsf{X}_i$,\;\;$i\in \mathbb{N}_n$},
\end{equation}
respectively,
\begin{equation}\label{Tlocal}
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f_i\circ l_i^{-1})(x), \quad\text{$x\in l_i(\mathsf{X}_i)$,\;\; $i\in \mathbb{N}_n$},
\end{equation}
where $f_i := f\vert_{\mathsf{X}_i}$, on appropriate function spaces.
\subsection{Bounded Local Solutions.}
The extension of the results presented in Section \ref{GFI} carry over to the setting of local fractal interpolation. However, special care must be taken when considering the domains of the functions involved.

Besides the function space $\mathcal{B} (\mathsf{X},\mathsf{F})$ already introduced in the previous section, we also need the local version of this space, namely, $\mathcal{B} (\mathsf{X}_i,\mathsf{F})$, $i\in \mathbb{N}_n$. To this end and in view of the result below, we now assume that
\begin{enumerate}
\item $s_i\in \mathcal{B}(\mathsf{X}_i,\mathbb{R})$ and
\item $q_i\in \mathcal{B}(\mathsf{X}_i,\mathsf{F})$.
\end{enumerate}
Then the RB operator $T$ defined in \eqref{Tlocal} maps $\mathcal{B} (\mathsf{X},\mathsf{F})$ into itself. Hence, we arrive at the local version of Theorem \ref{sol}.
\begin{theorem}
The system of functional equations \eqref{localpsi} has a unique bounded solution $\psi: \mathsf{X}\to\mathsf{F}$, respectively, the RB operator \eqref{Tlocal} has a unique bounded fixed point $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item conditions \eqref{cl1} and \eqref{cl2} are satisfied, and
\item $s:= \max\limits_{i\in \mathbb{N}_n} \sup\limits_{x\in \mathsf{X}_i} |s_i(x)| < 1$.
\end{enumerate}
\end{theorem}

\begin{proof}
The proof follows the same arguments as those given in the proof of Theorem \ref{sol}; replace $\mathsf{X}_i$ and $\mathsf{X}$ there by $l_i(\mathsf{X}_i)$ and $\mathsf{X}_i$ here. We leave the details to the reader who may also consult \cite{SB}.
\end{proof}

As in the global setting of fractal interpolation, a few remarks are in order.
\begin{remarks}\hfill
\begin{enumerate}
\item The solution/fixed point $\psi$ is referred to as a \emph{bounded local fractal function}.
\item The function $\psi$ depends on $n$, the family of subsets $\{\mathsf{X}_i : i\in\mathbb{N}_n\}$, the partition $(l_i(\mathsf{X}) : i\in\mathbb{N}_n)$ induced by the injective mappings $l_i$, and the now locally defined functions $s_i$ and $q_i$.
\end{enumerate}
\end{remarks}

As observed above, conditions \eqref{cl1} and \eqref{cl2} cannot be relaxed without adding some compatibility conditions to guarantee that the RB operator $T$ has the form given by Eqn.~ \eqref{Tlocal}. Should Eqn. ~\eqref{cl1} not be satisfied, one would have to impose in our current setting the following, now local, compatibility conditions:
\begin{align}
\forall i,j\in &\mathbb{N}_n\,\forall x_1\in\mathsf{X}_i\,\forall x_2\in\mathsf{X}_j:\nonumber\\
& \; l_i (x_1) = l_j (x_2)\;\;\Longrightarrow\;\; q_i(x_1) + s_i (x_1) \psi (x_1) = q_j(x_2) + s_ j(x_2) \psi (x_2).
\end{align}
We again refer to \cite{SB,SB1} for more details regarding this issue.
\subsection{$L^p$ Local Solutions}
As an application of local fractal functions for a particular setting, we again consider local fractal functions in $L^p$ spaces for $p\in[0,\infty]$.

For this purpose, we choose the following set up. Let $\mathsf{E} := \mathbb{R} := \mathsf{Y}$ with the canonical Euclidean norm and let $\mathsf{X}:= [0,1]$. Suppose that we are given a partition of $\mathsf{X}$ of the form $\Delta := (0 =: x_0 < x_1 < \cdots < x_{n-1} < x_n := 1)$, for some integer $n>1$. Furthermore, suppose that $\{\mathsf{X}_i  :  i \in \mathbb{N}_n\}$ is a family of half-open intervals of $[0,1]$.

We define affine mappings $l_i:\mathsf{X}_i$ onto $[x_{i-1}, x_i)$, $i = 1, \ldots, n-1$, and from $\mathsf{X}_n^+ := \mathsf{X}_n\cup l_n^{-1}(1-)$ onto $[x_{n-1},x_n]$, where $l_n$ maps $\mathsf{X}_n$ onto $[x_{n-1}, x_n)$.

We have the following result for RB-operators defined on the Lebesgue spaces $L^p[0,1]$, $1\leq p \leq \infty$. (See, also, \cite{bhm}.)

\begin{theorem}\label{th4.2}
Assume that $q_i\in L^p (\mathsf{X}_i, [0,1])$ and $s_i\in L^\infty (\mathsf{X}_i, \mathbb{R})$, $i \in\mathbb{N}_n$.
The system of functional equations
\[
\psi (l_i (x)) = q_i (x) + s_i (x) \psi (x), \quad x\in [0,1],\; i\in \mathbb{N}_n,
\]
has a unique solution $\psi\in L^p[0,1]$, $1\leq p \leq \infty$, respectively, the RB operator
\begin{equation}\label{Top}
T f (x) = (q_i\circ l_i^{-1})(x) + (s_i\circ l_i^{-1})(x)\cdot (f_i\circ l_i^{-1})(x),\quad{x\in\mathsf{X}_i},\;i\in \mathbb{N}_n,
\end{equation}
has a unique fixed point $\psi\in L^p(\Omega, \mathbb{R}^k)$, where $f_i = f\vert_{\mathbb{X}_i}$, provided that
\begin{equation}\label{condition}
\begin{cases}
\left(\displaystyle{\sum_{i=1}^n}\, a_i \,\|s_i\|_{\infty, \mathsf{X}_i}^p\right)^{1/p} < 1, & p\in[1,\infty);\\
\max\limits_{i\in\mathbb{N}_n}\|s_i\|_{\infty,\mathsf{X}_i} < 1, & p = \infty,
\end{cases}
\end{equation}
where $a_i$ denotes the Lipschitz constant of $(l_i^{-1})'$. Here, we wrote $\|s_i\|_{\infty,\mathsf{X}_i}$ for $\sup\limits_{x\in \mathsf{X}_i} |s_i(x)|$.
\end{theorem}

\begin{proof}
Note that under the hypotheses on the functions $q_i$ and $s_i$ as well as the mappings $l_i$, $T f$ is well-defined and an element of $L^p[0,1]$. It remains to be shown that under condition \eqref{condition}, $T$ is contractive on $L^p[0,1]
$.

To this end, let $g,h \in L^p [0,1]$ and let $p\in [0,\infty)$. Then
\begin{align*}
\|T g - T h\|^p_{p} & = \int\limits_{[0,1]} |T g (x) - T h (x)|^p dx\\
& = \int\limits_{[0,1]} \left|\sum_{i=1}^{n} (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\,\chi_{l_i(\mathbb{X}_i)}(x)\right|^p\, dx\\
& = \sum_{i=1}^{n}\,\int\limits_{[x_{i-1},x_i]}\left| (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\right|^p\,dx\\
& = \sum_{i=1}^{n}\,a_i\,\int\limits_{\mathbb{X}_i} \left| s_i (x) [g_i(x)- h_i(x)]\right|^p\,dx\\
&  \leq \sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\,\int\limits_{\mathbb{X}_i} \left| g_i(x) - h_i(x)\right|^p\,dx = \sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\,\|f_i - g_i\|^p_{p,\mathbb{X}_i}\\
& = \sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\,\|g_i - h_i\|^p_{p} \leq \left(\sum_{i=1}^{n}\,a_i\,\|s_i\|^p_{\infty, \mathbb{X}_i}\right) \|g - h\|^p_{p}.
\end{align*}
Now let $p= \infty$. Then,
\begin{align*}
\|T g - T h\|_{\infty} & = \left\|\sum_{i=1}^{n} (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\,\chi_{l_i(\mathbb{X}_i)}(x)\right\|_\infty\\
& \leq \max_{i\in\mathbb{N}_n}\,\left\| (s_i\circ l_i^{-1})(x) [(g_i\circ l_i^{-1})(x) - (h_i\circ l_i^{-1})(x)]\right\|_{\infty,\mathbb{X}_i}\\
& \leq \max_{i\in\mathbb{N}_n}\|s_i\|_{\infty,\mathbb{X}_i} \left\|g_i - h_i]\right\|_{\infty,\mathbb{X}_i} = \max_{i\in\mathbb{N}_n}\|s_i\|_{\infty,\mathbb{X}_i} \left\|g_i - h_i]\right\|_{\infty}\\
&  \leq \left(\max_{{i\in\mathbb{N}_n}}\,\|s_i\|_{\infty,\mathbb{X}_i}\right) \left\|g - h]\right\|_{\infty}
\end{align*}
These calculations prove the claims.
\end{proof}
\begin{remark}
Although conditions \eqref{condition} resembles those presented in Theorem \ref{th3.2}, they are -- because of the local nature of the estimates -- more subtle than the former ones.
\end{remark}
\subsection{Continuous Local Solutions}\label{sect4.3}
As in the global case, we can consider continuous local solutions to the given set of functional equations, respectively, an RB operator. Here, we present the analog to Theorem \ref{th3.3} but refer to the literature for the proof. A good reference in the functional equation setting is again \cite{SB}. In the RB operator setting, we refer to \cite{m2}.

\begin{theorem}\label{th4.3}
The system of functional equations \eqref{localpsi} has a unique continuous solution $\psi: \mathsf{X}\to\mathsf{F}$ provided that
\begin{enumerate}
\item $\mathsf{X} = \bigcup\limits_{i=1}^n l_i(\mathsf{X}_i)$,
\item the functions $l_i$, $q_i$, and $s_i$ are continuous,
\item and $\forall i,j\in \mathbb{N}_n$, $i\neq j$, $\forall x_1\in\mathsf{X}$, $\forall x_2\in \mathsf{X}_i$:
\begin{equation}\label{joinup}
\lim_{\substack{x\to x_1\\ x\in\mathsf{X}_j}} f_j(x) = f_i(x_2)\;\;\Longrightarrow\;\; \lim_{\substack{x\to x_1\\ x\in\mathsf{X}_j}} q_j(x) + s_j(x) \psi(x) = q_i(x_2) + s_i(x_2) \psi(x_2).
\end{equation}
\end{enumerate}
\end{theorem}
As an application and an example for a continuous local fractal interpolation, we present the following set up which plays an important role in fractal-based numerical analysis as discussed in \cite{bhm}.

Suppose that $\mathsf{E} := \mathbb{R} =: \mathsf{F}$ and $\mathsf{X}:= [0,1]$. For an even integer $n\in\mathbb{N}$, define subsets
\begin{equation}\label{Xs}
\mathsf{X}_{2j-1} := \mathsf{X}_{2j} := \left[\tfrac{2j-2}{n},\tfrac{2j}{n}\right], \quad j \in \{1,\ldots, \tfrac{n}{2}\}
\end{equation}
and affine mappings $l_i: \mathsf{X}_i\to [0,1]$ by
\begin{equation}\label{ls}
l_{2j-1} (x) := \tfrac{x}{2} + \tfrac{j-1}{n}\quad\text{and}\quad l_{2j} (x) := \tfrac{x}{2} + \tfrac{j}{n}, \quad x\in \mathsf{X}_{2j-1} = \mathsf{X}_{2j}.
\end{equation}
Note that
\[
l_i (\mathsf{X}_i) = \left[\tfrac{i-1}{n},\tfrac{i}{n}\right], \quad i\in\mathbb{N}_n.
\]
Further, let $x_i := l_i(\mathsf{X}_i)\cap l_{i+1}(\mathsf{X}_{i+1}) = \{\frac{i}{n}\}$, $i\in\mathbb{N}_{n-1}$, $x_0:= 0$, and $x_n := 1$. In the terminology of \cite{SB}, the elements of $\{x_i : i\in \mathbb{N}_{n-1}\}$ are called \emph{contact points}.

We denote the distinct endpoints of the partitioning intervals $\{l_i(\mathbb{X}_{i})\}$ by $\{x_0 < x_1 <\ldots < x_{N}\}$ where $x_0 = 0$ and $x_N = 1$ and refer to them as {knots}.

Furthermore, we assume that we are given interpolation values at the endpoints of the intervals $\mathbb{X}_{2j-1} = \mathbb{X}_{2j}$:
\begin{equation}\label{Delta}
\Delta := \left\{(x_{2j}, y_j) :  j = 0, 1, \ldots, n/2\right\}.
\end{equation}
Let
$$
\mathcal{C}_{\Delta}(\mathsf{X}) := \{f\in \mathcal{C}(\mathsf{X}) : f(x_{2j}) = y_j, \,\forall\, j = 0, 1, \ldots, n/2\}.
$$
Here, $\mathcal{C}(\mathsf{X}) := \mathcal{C}(\mathsf{X},\mathbb{R})$ denotes the Banach space of all continuous functions $\mathsf{X}\to\mathbb{R}$ endowed with the supremum norm $\n{\cdot}$. Note that $\mathcal{C}_{\Delta}(\mathsf{X})$ is a closed, hence complete, \emph{metric} subspace of $\mathcal{C}(\mathsf{X})$ to which we can apply the Banach fixed point theorem. To this end, consider an RB operator $T$ of the form \eqref{Top} acting on $\mathcal{C}_{\Delta}(\mathsf{X})$.

In order for $T$ to map $\mathcal{C}_{\Delta}(\mathsf{X})$ into itself we require that
\begin{equation}\label{qsandss}
q_i, s_i \in \mathcal{C}(\mathsf{X}_i) := \mathcal{C}(\mathsf{X}_i,\mathbb{R}) := \{f: \mathsf{X}_i\to \mathbb{R} :  \text{$f$ continuous}\}
\end{equation}
and that
\begin{equation}\label{intcon}
y_{j-1} = (T f) (x_{2(j-1)}) \quad \wedge \quad y_{j} = (T f) (x_{2j}), \quad j = 1, \ldots, n/2,
\end{equation}
where $x_{2j} := \frac{2j}{n}$.

We remark that the preimages of the knots $x_{2(j-1)}$ and $x_{2j}$ are the endpoints of $\mathsf{X}_{2j-1} = \mathsf{X}_{2j}$. Substitution of $T$, as given in \eqref{Top}, into \eqref{intcon} and simplification results in
\begin{equation}\label{intcons}
\begin{split}
q_{2j-1} (x_{2(j-1)}) + \left(s_{2j-1}(x_{2(j-1)}) - 1\right) y_{j-1} & = 0,\\
q_{2j} (x_{2j}) + \left(s_{2j}(x_{2j}) - 1\right) y_{j} & = 0.
\end{split}
\end{equation}

To ensure global continuity of $T f$ on $\mathsf{X}=[0,1]$, we also have to impose the following join-up conditions at the oddly indexed knots. (Note that these oddly indexed knots are the images of the midpoints of the intervals $\mathsf{X}_{2j-1} = \mathsf{X}_{2j}$.)
\begin{equation}\label{prejoin}
(T f) (x_{2j-1}-)  = (T f) (x_{2j-1}+) , \quad j = 1, \ldots, n/2.
\end{equation}
These join-up conditions imply that
\begin{equation}\label{4.16}
q_{2j} (x_{2(j-1)}) + s_{2j} (x_{2(j-1)}) y_{j-1} = q_{2j-1} (x_{2j}) + s_{2j-1} (x_{2j})  y_j.
\end{equation}
In the case that all functions $q_i$ and $s_i$ are constant, \eqref{4.16} reduces to the condition given in \cite[Example 2]{bhm}.

We summarize these results in the next theorem.
\begin{theorem}\label{thm8}
Let $\mathsf{X}:= [0,1]$ and let $n\in 2\mathbb{N}$. Suppose that subsets of $\mathsf{X}$ are given by \eqref{Xs} and the associated mappings $l_i$ by \eqref{ls}. Further suppose that the functions $q_i$ and $s_i$ satisfy \eqref{qsandss} and that the join-up conditions \eqref{intcon} and \eqref{prejoin} hold. Then, an RB operator $T$ of the form \eqref{Top} maps $\mathcal{C}_\Delta (\mathsf{X})$ into itself and is well-defined.

If, in addition,
\begin{equation}\label{condfors}
\max\left\{\|s_i\|_{\infty,\mathsf{X}_i} : i\in\mathbb{N}_n\right\} < 1,
\end{equation}
then $T$ is a contraction on $\mathcal{C}_\Delta (\mathsf{X})$ and thus possesses a unique continuous fixed point $\psi:[0,1]\to \mathbb{R}$ satisfying $\psi\in \mathcal{C}_\Delta (\mathsf{X})$.
\end{theorem}

This unique fixed point $\psi$ is called a \emph{continuous local fractal function}.

\begin{proof}
It remains to show that under the condition \eqref{condfors}, $T$ is contractive on $\mathcal{C}_\Delta (\mathsf{X})$. This, however, follows immediately from the case $p=\infty$ in the proof of Theorem \ref{th4.2}.
\end{proof}

For RB operators of the form \eqref{Top} mapping function spaces such as H\"older, Sobolev, Besov or Triebel Lizorkin into themselves, we refer the interested reader to \cite{m2,m3}.
\section{Non-stationary Fractal Interpolation}\label{sect5}
In this section, we extend the notion of global fractal interpolation to a non-stationary setting. In other words, we no longer assume that we keep the functions $q_i$ and $s_i$ the same at each level of iteration in the construction of the fixed point $\psi$.
\subsection{The non-stationary setting}
To this end, consider a doubly-indexed family of injective contractions $\{l_{i_k,k} : i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}$ from $\mathsf{X}\to \mathsf{X}$ where $\mathsf{X}$ is nonempty bounded subset of a normed space $\mathsf{E}$ generating a partition of $\mathsf{X}$ for each $k\in \mathbb{N}$ in the sense of \eqref{cl1} and \eqref{cl2}.

Suppose that $\mathsf{F}$ is a Banach space, $\{q_{i_k,k}: i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\} \subset \mathcal{B}(\mathsf{X},\mathsf{F})$, and $\{s_{i_k,k}: i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}\subset \mathcal{B}(\mathsf{X},\mathbb{R})$ are such that
\[
s := \sup\limits_{k\in \mathbb{N}} \max\limits_{i_k\in \mathbb{N}_k} \|s_{i_k,k}\|_\infty < 1.
\]
For each $k\in \mathbb{N}$, we define an RB operator $T_k: \mathcal{B}(\mathsf{X},\mathsf{F})\to \mathcal{B}(\mathsf{X},\mathsf{F})$ by
\begin{align}
(T_k f)(l_{i_k,k}(x)) &:= q_{i_k,k}(x) + s_{i_k,k}(x) \cdot f(x), \quad\forall x\in\mathsf{X}.\label{nonstat1}
\end{align}
It is not difficult to verify that each $T_k$ is a contraction on $\mathcal{B}(\mathsf{X},\mathsf{F})$ with Lipschitz constant
\begin{equation}\label{Lip}
\Lip (T_k) = \max\limits_{i_k\in \mathbb{N}_k} \|s_{i_k,k}\|_\infty \leq s < 1.
\end{equation}
In order to continue, we require the following definition and result from \cite{LDV} adapted to our current setting.

\begin{definition}{\cite[Definition 3.6]{LDV}}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of transformations $T_k:\mathcal{B}(\mathsf{X},\mathsf{F})\to \mathcal{B}(\mathsf{X},\mathsf{F})$. A subset $\mathscr{I}$ of $\mathcal{B}(\mathsf{X},\mathsf{F})$ is called an invariant set of the sequence $\{T_k\}_{k\in \mathbb{N}}$ if
\[
\forall\,k\in \mathbb{N}\;\forall\,x\in \mathscr{I}: T_k (x)\in \mathscr{I}.
\]
\end{definition}
A criterion for obtaining an invariant domain for a sequence $\{T_k\}_{k\in \mathbb{N}}$ of transformations on $\mathcal{B}(\mathsf{X},\mathsf{F})$ is also given in \cite{LDV}.

\begin{proposition}{ \cite[Lemma 3.7]{LDV}}\label{prop2.1}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of transformations on $\mathcal{B}(\mathsf{X},\mathsf{F})$. Suppose there exists a $g\in \mathcal{B}(\mathsf{X},\mathsf{F})$ such that for all $f\in\mathcal{B}(\mathsf{X},\mathsf{F})$
\[
\n{T_k(x) - g} \leq \mu\,\n{f - g} + M,
\]
for some $\mu\in [0,1)$ and $M> 0$. Then the ball $B_r(g)$ of radius $r = M/(1-\mu)$ centered at $g$ is an invariant set for $\{T_k\}_{k\in \mathbb{N}}$.
\end{proposition}
\begin{proof}
The proof, although in a more general setting, is found in \cite{LDV}.
\end{proof}

\begin{proposition}\label{prop5.2}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of RB operators of the form \eqref{nonstat1} on $(\mathcal{B}(\mathsf{X},\mathsf{F}),\n{\cdot})$. Suppose that the elements of $\{q_{i_k,k}: i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}$ satisfy
\begin{equation}\label{condq}
\sup\limits_{k\in \mathbb{N}} \max\limits_{i_k\in \mathbb{N}_k} \n{q_{i_k,k}} \leq M,
\end{equation}
for some $M > 0$. Then the ball $B_r(0)$ of radius $r=M/(1-s)$ centered at $0\in \mathcal{B}(\mathsf{X},\mathsf{F})$ is an invariant set for $\{T_k\}_{k\in \mathbb{N}}$.
\end{proposition}
\begin{proof}
Let $x\in \mathsf{X}$. Then there exists an $i_k\in \mathbb{N}_{n_k}$ with $x\in l_{i_k,k}(\mathsf{X})$. Thus, for any $f\in \mathcal{B}(\mathsf{X},\mathsf{F})$,
\begin{align*}
\n{T_k f (x)}_\mathsf{F} &\leq \n{s_{i_k,k}\circ l_{i_k,k}^{-1}(x) \cdot f\circ l_{i_k,k}^{-1}(x)}_\mathsf{F} + \n{T_k(0)}_\mathsf{F}
\end{align*}
By \eqref{condq},  $T_k(0)$ is uniformly bounded on $\mathcal{B}(\mathsf{X},\mathsf{F})$ by $M > 0$. Hence,
\[
\n{s_{i_k,k}\circ l_{i_k,k}^{-1}(x) \cdot f\circ l_{i_k,k}^{-1}(x)}_\mathsf{F} \leq s \n{f\circ l_{i_k,k}^{-1}(x)}_\mathsf{F},
\]
which gives, after taking the supremum over $x\in\mathsf{X}$, $\n{T_k f} \leq s\,\n{f} + M$. Proposition \eqref{prop2.1} now yields the statement.
\end{proof}

To arrive at the main result of this section, we require yet another definition. (See also \cite[Section 4]{LDV}.)
\begin{definition}
Suppose that $f_0\in\mathcal{B}(\mathsf{X},\mathsf{F})$ and that $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of RB operators on $\mathcal{B}(\mathsf{X},\mathsf{F})$. The sequences
\begin{equation}
\Phi_k (f_0) := T_k\circ T_{k-1} \circ \cdots \circ T_1 (f_0)
\end{equation}
and
\begin{equation}
\Psi_k (f_0) := T_1\circ T_2 \circ \cdots \circ T_k (f_0)
\end{equation}
are called the forward, respectively, backward trajectory of $f_0$.
\end{definition}

In connection to the above definition, we need the following result from \cite[Corollary 11]{LDV} adapted to our setting.

\begin{theorem}\label{th5.1}
Suppose that $\{T_k\}_{k\in \mathbb{N}}$ is a sequence of RB operators of the form \eqref{nonstat1} on $\mathcal{B}(\mathsf{X},\mathsf{F})$. Further suppose that
\begin{enumerate}
\item there exists a nonempty closed invariant set $\mathscr{I}\subseteq\mathcal{B}(\mathsf{X},\mathsf{F})$ for $\{T_k\}_{k\in \mathbb{N}}$;
\item and
\begin{equation}\label{eq5.6}
\sum_{k=1}^\infty\prod_{j=1}^k \Lip (T_j) < \infty.
\end{equation}
\end{enumerate}
Then the backward trajectories $\Psi_k(f_0)$ converge for any initial $f_0\in\mathscr{I}$ to a unique function $\psi\in\mathscr{I}$.
\end{theorem}

With these preliminaries, we obtain the main result of this section.

\begin{theorem}\label{thm4.1}
The backwards trajectories $\{\Psi_k\}_{k\in \mathbb{N}}$ converge for any initial $f_0\in \mathscr{I}$ to a unique function $\psi \in \mathscr{I}$, where $\mathscr{I}$ is the closed ball in $\mathcal{B}(\mathsf{X},\mathsf{F})$ of radius $M/(1-s)$ centered at $0$.
\end{theorem}
\begin{proof}
By Proposition \ref{prop5.2} and Theorem \ref{th5.1} it remains to be shown that \eqref{eq5.6} is satisfied. This, however, follows directly from \eqref{Lip}:
\[
\prod_{j=1}^k \Lip (T_j) \leq s^k\quad\text{and}\quad\sum_{k=1}^\infty s^k = \frac{s}{1-s} < \infty.\qedhere
\]
\end{proof}

A fixed point $\psi$ generated by a sequence $\{T_k\}$ of non-stationary RB operators will be called a \emph{non-stationary fractal function (of class $\mathcal{B}(X,Y)$).}

Employing a non-stationary sequence of RB operators $\{T_k\}$ allows the construction of more general fractal functions exhibiting different local behavior at different scales. This is illustrated by the following example which is taken from \cite{m4}.

\begin{example}
Let $\mathsf{X}:= [0,1]\subset\mathbb{R}$ and $\mathsf{F}:= \mathbb{R}$. Consider the two RB operators
\[
T_1 f (x) := \begin{cases} 2x + \frac12 f(2x), & x\in [0,\frac12),\\
2 - 2x + \frac12 f(2x-1), & x\in [\frac12,1],
\end{cases}
\]
and
\[
T_2 f (x) := \begin{cases} 2x + \frac14 f(2x), & x\in [0,\frac12),\\
2 - 2x + \frac14 f(2x-1), & x\in [\frac12,1].
\end{cases}
\]
For both operators, $l_i (x) := \frac12 (x + i-1)$, $i=1,2$.

The sequence $T_1^k f \to \tau$, where $\tau$ denotes the Takagi function \cite{Tak} and $T_2^k\to q$, where $q(x)= 4x(1-x)$.

Consider the alternating sequence $\{T_k\}_{k\in \mathbb{N}}$ of RB operators given by
\[
T_k := \begin{cases} T_1, & 10(j-1) < k \leq 10j-5,\\
T_2, & 10j-5 < k \leq 10j,
\end{cases}\quad j\in \mathbb{N}.
\]
Two images of this hybrid attractor of the backward trajectory $\Psi_k$ starting with $f_0 \equiv 0$ are shown in Figure \ref{fig11}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=3cm]{frac1.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=3cm]{frac2.pdf}
\end{center}
\caption{The hybrid $\tau-q$ attractor. It is smooth at one scale but fractal at another.}\label{fig11}
\end{figure}
\end{example}

\subsection{Non-stationary Fractal Interpolation}

In this subsection, we consider the case $\mathsf{X} := [0,1]$ and $\mathsf{F} := \mathbb{R}$. In the following, we use a sequence $\{T_k\}$ of RB operators of a particular form to obtain a continuous fixed point $\psi$.

To this end, consider an RB operator defined by \eqref{eq3.17}, choose two arbitrary functions $f,b\in \mathcal{B}(\mathsf{X}) := \mathcal{B}(\mathsf{X}, \mathbb{R})$ and set
\begin{equation}\label{qs}
q_i := f\circ l_i - s_i\cdot b.
\end{equation}
The (stationary) RB operator associated with this particular choice is then given by
\begin{equation}
(Tg)(x) = f + s_i(l_i^{-1}(x))\cdot(g-b)(l_i (x)), \quad x\in l_i(\mathsf{X})
\end{equation}
\begin{remark}
The fixed point of the above defined RB operator may be thought of as the ``fractalization" of a given function $f$.
\end{remark}

To work in a non-stationary setting, let $k\in \mathbb{N}$ and let $\{l_{i_k,k} : i_k\in \mathbb{N}_{n_k}, \, k\in \mathbb{N}\}$ be family of injections from $[0,1]\to [0,1]$ generating a partition of $[0,1]$ in the sense of \eqref{cl1} and \eqref{cl2}. We assume w.l.o.g. that $l_{1,k}(0) = 0$ and $l_{n_k,k}(1) = 1$ and define
\begin{align*}
x_{i_k-1,k} := l_{i_k,k}(0), \quad x_{i_k,k} := l_{i_k,k}(1), \quad i_k\in \mathbb{N}_{n_k}
\end{align*}
where $x_{0,k} :=0$ and $x_{n_k,k} := 1$. By relabelling -- if necessary -- we may further assume that $0 = x_{0,k}  < \cdots < x_{i_k-1,k} < x_{i_k,k} < \cdots x_{n_k,k} = 1$.

Let $f\in \mathcal{C}(\mathsf{X})$ be arbitrary. Define a metric subspace of $\mathcal{C}(\mathsf{X})$ by
\[
\mathcal{C}_*(\mathsf{X}) := \{g\in \mathcal{C}(\mathsf{X}) : g(0) = f(0) \wedge g(1) = f(1)\}
\]
and note that $\mathcal{C}_*(\mathsf{X})$ becomes a complete linear metric space when endowed with the metric $d$ that is induced by the sup-norm from $\mathcal{C}(\mathsf{X})$.

Furthermore, let $b\in \mathcal{C}_*[0,1]$ be the unique affine function whose graph connects the points $(0, f(0))$ and $(1, f(1))$:
\begin{equation}\label{b}
b(x) = (f(1) - f(0)) x + f(0).
\end{equation}
Let $\{\mathcal{P}_k\}_{k\in \mathbb{N}}$ be a family of sets of points in $\mathsf{X}\times \mathsf{F}$ where
\[
\mathcal{P}_k :=\{(x_{j_k}, f(x_{j,k})\in \mathsf{X}\times \mathsf{F} : j = 0,1,\ldots, n\}.
\]
For $k\in \mathbb{N}$, define an RB operator $T_k: \mathcal{C}_*[0,1]\to \mathcal{C}_*[0,1]$ by
\begin{equation}\label{eq5.1}
T_k g = f + \sum_{i_k=1}^{n_k} s_{i_k,k}\circ l_{i_k,k}^{-1} \cdot (g-b)\circ l_{i_k,k}^{-1}\, \chi_{l_{i_k,k} [0,1]},
\end{equation}
where $\{s_{i_k,k}\}_{i_k=1}^{n_k}\subset \mathcal{C}[0,1]$ such that
\[
\sup_{k\in \mathbb{N}} \max_{i_k\in \mathbb{N}_{i_k}} \|s_{i_k,k}\|_\infty < 1.
\]
Note that
\begin{align*}
T_k g (x_{i_k,k}-) &= f(x_{i_k,k}-) + s_{i_k,k}\circ l_{i_k,k}^{-1}(x_{i_k,k}-) \cdot (g-b)\circ l_{i_k,k}^{-1}(x_{i_k,k}-) \\
&= f(x_{i_k,k}) + s_{i_k,k} (1)\cdot (f - b)(1) = f(x_{i_k,k})
\end{align*}
and
\begin{align*}
T_k g (x_{i_k,k}+) &= f(x_{i_k,k}+) + s_{i_k+1,k}\circ l_{i_k+1,k}^{-1}(x_{i_k,k}+) \cdot (g-b)\circ l_{i_k+1,k}^{-1}(x_{i_k,k}+) \\
&= f(x_{i_k,k}) + s_{i_k+1,k} (0)\cdot (f - b)(0) = f(x_{i_k,k}).
\end{align*}
implying that $T_k g$ is continuous at the points $x_{i_k,k}\in [0,1]$:
\[
T_k g (x_{i_k,k}-) = T_k g (x_{i_k,k}+), \quad\forall\,i_k\in\{1, \ldots, n-1\}.
\]
Hence,  $T_k g\in \mathcal{C}_*[0,1]$ and $T_k g$ interpolates $\mathcal{P}_k$ in the sense that
\[
T_k g (x_{i_k,k}) = f(x_{i_k,k}), \quad \forall\,i_k\in \mathbb{N}_{n_k}.
\]

\begin{proposition}
A nonempty closed invariant set for $\{T_k\}_{k\in \mathbb{N}}$ is given by the closed ball in $\mathcal{C}_*(\mathsf{X})$,
\begin{equation}\label{I}
\mathscr{I} = \left\{g\in \mathcal{C}_*[0,1] : \|g\| \leq \frac{\|f\| + s \|b\|}{1-s}\right\},
\end{equation}
where $s$ is defined by \eqref{Lip}.
\end{proposition}

\begin{proof}
From the form \eqref{qs} of the functions $q_{i_k,k}$, we obtain from \eqref{condq} the estimate $\|q_{i_k,k}\| \leq \|f\| + s \|b\|$, which by Proposition \ref{prop2.1} yields the result.
\end{proof}

Together with Theorem \ref{thm4.1}, the above arguments prove the next theorem.

\begin{theorem}
Let $\{T_k\}_{k\in \mathbb{N}}$ be a sequence of RB operators of the form \eqref{eq5.1} each of whose elements acts on the complete metric space $(\mathcal{C}_*(\mathsf{X}), d)$ where $f\in \mathcal{C}_*(\mathsf{X})$ is arbitrary and $b$ is given by \eqref{b}. Furthermore, let the family of functions $\{s_{i_k,k}\} \subset \mathcal{C}(\mathsf{X})$ satisfy \eqref{Lip}. Then, for any $f_0\in \mathscr{I}$, the backward trajectories $\Psi_k (f_0)$ converge to a function $\psi \in \mathscr{I}$ which interpolates $\mathcal{P}_k$.
\end{theorem}

We refer to the fixed point $\psi\in \mathcal{C}_*(\mathsf{X})$ as a \emph{continuous non-stationary fractal interpolation function}.

\begin{remark}
As $f_0$ one may choose $f$ or $b$.
\end{remark}

To illustrate the above results, we present the following example from \cite{m4}.

\begin{example}
Consider the two RB operators $T_i: C[0,1]\to C[0,1]$, $i=1,2$, defined by
\[
(T_1 f)(x) = \begin{cases} -\frac{1}{2}\,f(4x), & x\in[0,\frac{1}{4}),\\
 -\frac{1}{2} + \frac{1}{2}\,f(4x-1), & x\in[\frac{1}{4},\frac{1}{2}),\\
 \frac{1}{2}\,f(4x-2), & x\in[\frac{1}{2},\frac{3}{4}),\\
 \frac{1}{2} + \frac{1}{2}\,f(4x-3), & x\in[\frac{3}{4},1],\end{cases}
\]
and
\[
(T_2 f)(x) := \begin{cases} \frac34 f(2x), & x \in [0,\frac12),\\
\frac34 + \frac14 f(2x-1), &  x \in [\frac12,1].
\end{cases}
\]
The RB operators $T_1$ and $T_2$ generate \emph{Kiesswetter's fractal function} \cite{Kiess}, respectively, a \emph{Casino function} \cite{DS}.

Consider again the alternating sequence $\{T_i\}_{i\in \mathbb{N}}$ of RB operators given by
\[
T_k := \begin{cases} T_1, & 10(j-1) < k \leq 10j-5,\\
T_2, & 10j-5 < k \leq 10j,
\end{cases}\quad j\in \mathbb{N}.
\]
Two images of the hybrid attractor of the backward trajectory $\Psi_k$ starting with the function $f_0 (x) = x$, $x\in\mathsf{X}$, are shown below in Figure \ref{fig22}.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=4cm]{frac4.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{frac3.pdf}
\end{center}
\caption{The hybrid Kiesswetter-Casino attractor.}\label{fig22}
\end{figure}

\end{example}

\section{Quaternionic Fractal Interpolation}
In this section, we extend fractal interpolation to a quaternionic setting. As quaternions from a non-commutative division algebra, the non-commutativity generates more intricate fractal patterns.

First, we give a short introduction to quaternions and present some of those properties that are relevant for the remainder of this section. The interested reader is referred to the literature on the subject, a short and subjective list of which is \cite{Bour,GHS,GS,Krav,MGS}.

\subsection{A Brief Introduction to Quaternions}
Let $\{e_1, e_2, e_3\}$ be the canonical basis of the Euclidean vector space $\mathbb{R}^3$. We call $\{e_1, e_2, e_3\}$ imaginary units and require that the following multiplication rules hold:
\begin{gather}
e_1^2=e_2^2=e_3^2=-1,\label{6.1}\\
e_1e_2=e_3=-e_2e_1,\quad e_2e_3=e_1=-e_3e_2,\quad e_3e_1=e_2=-e_1e_3.\label{6.2}
\end{gather}
\begin{remark}
Note that \eqref{6.2} is equivalent to $e_1e_2e_3 = -1$.
\end{remark}
A \emph{real quaternion $q$} is then an expression of the form
\[
q = a+\sum_{i=1}^3v_ie_i, \quad a, v_1, v_2, v_3\in\mathbb{R}.
\]
The addition and multiplication of two quaternions $q_1 = a+\sum_{i=1}^3v_ie_i$ and $q_2 = b+\sum_{i=1}^3w_ie_i$ is defined by
\begin{align*}
q_1 + q_2 &:= (a+b) + \sum_{i=1}^3 (v_i+w_i) e_i\\
q_1 q_2&:= (ab - v_1w_1 - v_2w_2 - v_3w_3)\\
&\;\;+(a {w_1}+b {v_1}+{v_2} {w_3}-{v_3} {w_2})e_1\\
&\;\;+(a {w_2}+b {v_2}-{v_1} {w_3}+{v_3} {w_1})e_2\\
&\;\;+(a {w_3}+b {v_3}+{v-1}{w_2}-{v_2} {w_1})e_3.
\end{align*}
Each quaternion $q=a+\sum\limits_{i=1}^3v_ie_i$ may be decomposed as $q=\Sc (q)+\Ve (q)$ where $\Sc (q)=a$ is the {\it scalar part} of $q$ and $\Ve (q)=v=\sum\limits_{i=1}^3v_ie_i$ is the {\it vector part} of $q$. $v=\Ve (q)$ is also called a \emph{quaternionic vector}.

The {\it conjugate} $\overline{q}$ of the real quaternion $q=a+v$ is the quaternion $\overline{q}=a-v$. Note that $q\overline{q}=\overline{q}q=a^2+|v|^2=a^2+\sum\limits_{i=1}^3v_i^2$. Therefore, we can define a norm on $\mathbb{H}$ by setting
\[
\abs{q} := \sqrt{q\overline{q}}.
\]
The inverse of a quaternion $q$ is given by
\[
q^{-1} = \frac{\overline{q}}{\abs{q}^2}.
\]

It is straight-forward to establish that the collection of all real quaternions
\[
\mathbb{H} := \mathbb{H}_\mathbb{R} := \left\{a+\sum_{i=1}^3v_ie_i : a, v_1, v_2, v_3\in{\mathbb R}\right\},
\]
is a four-dimensional associative normed division algebra over $\mathbb{R}$. Due to the multiplication rules \eqref{6.1} and \eqref{6.2}, $\mathbb{H}$ is not commutative. We also note that $\mathbb{H}$, as already indicated above, is a four-dimensional vector space over $\mathbb{R}$ with basis $\{e_0,e_1,e_2,e_3\}$, where $e_0 := 1$.

\begin{remark}
We note that $\mathbb{H}$ and $\mathbb{R}^4$ are identical as \emph{point sets} but differ in their algebraic structure. The vector space $\mathbb{R}^4$ is not an algebra whereas $\mathbb{H}$ is one, albeit non-commutative.
\end{remark}

\begin{remark}
As we can define a norm and therefore a metric on $\mathbb{H}$, $\mathbb{H}$ becomes a topological space where open sets are defined via the metric. All other topological concepts such as limits, convergence, compactness etc. then follow. $\mathbb{H}$ thus becomes a topological space and, moreover, also a complete metric space.
\end{remark}

We also remark that if $v=\sum\limits_{j=1}^3v_je_j$ and $w=\sum\limits_{j=1}^3w_je_j$ are quaternionic vectors, then
\begin{equation}
vw=-\langle v,w\rangle +v\wedge w\label{vec mult},
\end{equation}
where $\langle v,w\rangle :=\sum\limits_{j=1}^3v_jw_j$ is the scalar product of $v$ and $w$ and
$$v\wedge w :=(v_2w_3-v_3w_2)e_1+(v_3w_1-v_1w_3)e_2+(v_1w_2-v_2w_1)e_3$$
is the vector (cross) product of $v$ and $w$.

For our purposes, we need to introduce the analog of a Banach space in the quaternionic setting. To this end, we begin with the following definitions. (See, also \cite{Bour}.)
\begin{definition}
A real vector space $\mathsf{V}$ is called a left quaternionic vector space if it is a left $\mathbb{H}$-module, i.e., if there exists a mapping $\mathbb{H}\times\mathsf{V}\to \mathsf{V}$, $(q,v)\mapsto q v$ which satisfies
\begin{enumerate}
\item $\forall v\in \mathsf{V}\,\forall q_1,q_2\in \mathbb{H}:$ $(q_1+q_2)v = q_1 v + q_2 v$.
\item $\forall v_1,v_2\in \mathsf{V}\,\forall q\in \mathbb{H}:$ $q (v_1+v_2) = q v_1 + q v_2$.
\item $\forall v\in \mathsf{V}\,\forall q_1,q_2\in \mathbb{H}:$ $q_1(q_2 v) = (q_1 q_2) v$.
\end{enumerate}
\end{definition}

\begin{remark}
In analogous fashion, one defines a \emph{right} quaternionic vector space as a right $\mathbb{H}$-module where the mapping is now $\mathsf{V}\times\mathbb{H}\to V$, $(v,q)\mapsto v q$.
\end{remark}

A \emph{two-sided quaternionic vector space} $\mathsf{V}$ is a left and right quaternionic vector space such that $\lambda v = v \lambda$, for all $\lambda\in\mathbb{R}$ and for all $v\in\mathsf{V}$. An example of a two-sided quaternionic vector space is given by $\mathbb{H}$ itself.

\begin{remark}
A quaternionic vector space becomes a real vector space when its scalar multiplication is restricted to $\mathbb{R}$.
\end{remark}

One can start with any real vector space $\mathsf{V}_\mathbb{R}$ and construct a two-sided quaternionic vector space by setting
\[
\mathsf{V}_\mathbb{H} := \left\{\sum_{i=0}^3 v_i\otimes e_i : v_i\in \mathsf{V}_\mathbb{R}\right\},
\]
where $\otimes$ denotes the algebraic tensor product. (See, for instance, \cite{Ng}.)

On the other hand, given any two-sided quaternionic vector space and defining
\[
\mathsf{V}_\mathbb{R} := \{v\in \mathsf{V} : \lambda v = v \lambda,\;\forall\lambda\in\mathbb{H}\},
\]
then $\mathsf{V}_\mathbb{R}$ is a real vector space called the \emph{real part of $\mathsf{V}$}.

The proof of the next result can be found in, i.e., \cite{Ng}.

\begin{proposition}
Let $\mathsf{V}$ be a two-sided quaternionic vector space and let $\mathsf{V}_\mathbb{R}$ denote its real part. Then $V\cong \mathsf{V}_\mathbb{R}\otimes\mathbb{H}$.
\end{proposition}

Next, we introduce quaternionic normed spaces.

\begin{definition}
Let $\mathsf{V}$ be a left quaternionic vector space. A function $\n{\cdot}: \mathsf{V}\to \mathbb{R}_0^+$ is called a norm on $V$ if
\begin{enumerate}
\item $\n{v} = 0$ iff $v=0$.
\item $\n{q v} = \abs{q} \n{v}$, for all $v\in \mathsf{V}$ and $q\in \mathbb{H}$.
\item $\n{v+w} \leq \n{v} + \n{w}$, for all $v,w\in \mathsf{V}$.
\end{enumerate}
A left quaternionic vector space endowed with a norm will be called a left quaternionic normed space.
\end{definition}

A left quaternionic normed space $\mathsf{E}$ is called \emph{complete} if it is a complete metric space with respect to the metric $d(x,y) = \n{x-y}$ induced by the norm $\n{\cdot}_\mathsf{E}$. In this case, we refer to $\mathsf{E}$ as \emph{left quaternionic Banach space}.

\begin{remark}
A left (or right) quaternionic Banach space becomes a real Banach space if the left (right) scalar multiplication is restricted to $\mathbb{R}$. (Cf., \cite[Section 2.3]{CGK}.)
\end{remark}

\begin{example}
The space $\mathbb{H}^k$ consisting of $k$-tuples of quaternions is both a left and a right quaternionic vector space. We represent elements $\xi\in\mathbb{H}^k$ as column vectors and define the quaternionic conjugate ${}^*$ of $\xi$ by
\[
\begin{pmatrix} \xi_1\\ \vdots\\ \xi_k
\end{pmatrix}^* := \begin{pmatrix} \overline{\xi_1} & \cdots & \overline{\xi_k}
\end{pmatrix},
\]
where each $\xi_j\in\mathbb{H}$. When endowed with the norm
\begin{equation}\label{dist}
\n{\xi}_k := \sqrt{\xi^*\xi} = \sqrt{\sum\limits_{l=1}^k \abs{\xi_j}^2},
\end{equation}
$\mathbb{H}^k$ becomes a two-sided quaternionic Banach space as $\lambda v = v \lambda$, for all $\lambda\in\mathbb{R}$ and for all $v\in\mathbb{H}^k$.

Note that under the norm $\n{\cdot}_k$, $\mathbb{H}^k$ becomes a topological space and also a complete metric space.
\end{example}

The final concept we need to introduce is that of left linear mapping between left quaternionic vector spaces.

\begin{definition}
Let $\mathsf{V}_1$ and $\mathsf{V}_2$ be left quaternionic vector spaces. A mapping $f:\mathsf{V}_1\to\mathsf{V}_2$ is called left linear if
\[
f(q v + w) = q f(v) + f(w), \quad\forall v,w\in \mathsf{V}, \forall q\in \mathbb{H}.
\]
A left linear mapping is called bounded if
\[
\n{f} := \sup_{x,y \in \mathsf{V}_1, x \neq y} \frac{\n{f(x)-f(y)}_{\mathsf{V}_2}}{\n{x-y}_{\mathsf{V}_1}} < \infty.
\]
\end{definition}
\subsection{Quaternionic Fractal Interpolation}
In this section, we introduce the novel concept of quaternionic fractal interpolation. For illustrative purposes, we do not choose the most general set up but restrict ourselves to the case where the left quaternionic Banach spaces are $\mathsf{E} := \mathbb{H}^k =:\mathsf{F}$, $k\in \mathbb{N}$.

For our purposes, we need the following function space. Let $\mathsf{X}\subset\mathbb{H}^k$ be compact (as defined via the norm $\n{\cdot}_k$) and let
\[
\mathcal{B}(\mathsf{X},\mathbb{H}^k) := \left\{f:\mathsf{X}\to\mathbb{H}^k : \text{$f$ is bounded}\right\}.
\]
A function $f:\mathsf{X}\to\mathbb{H}^k$ is called \emph{bounded} if there exists a real number $M>0$ such that $\n{f}_k \leq M$. If we define for $x\in\mathsf{X}$ and $\lambda\in\mathbb{H}$
\[
(f+g) (x) := f(x) + g(x)\quad\text{and}\quad (\lambda\cdot f)(x) := \lambda\cdot f(x)
\]
then $\mathcal{B}(\mathsf{X},\mathbb{H}^k)$ becomes a left quaternionic vector space. Setting for each $f\in \mathcal{B}(\mathsf{X},\mathbb{H}^k)$
\[
\n{f} := \sup_{x\in\mathsf{X}} \n{f(x)}_k,
\]
then $\mathcal{B}(\mathsf{X},\mathbb{H}^k)$ becomes a left Banach space. (See, for instance, \cite[Chapter IV.E.1.]{J} for the case $k=1$. The extension to $k > 1$ is straight-forward.)

For the sake of simplicity and the purpose of understanding the underlying issues, we concentrate on the special case that
\[
\mathsf{X} := \left\{q\in \mathbb{H} : \max\limits_{i=0,1,2,3} \abs{q_i} = 1\right\} \cong [-1,1]^4
\]
and that $k:=1$. The interested reader is encouraged to consider the general set-up.

We take the four-dimensional cube $\mathsf{X}$ and divide it into $n:=2^4$ congruent four-dimensional subcubes $\mathsf{X}_i$ each similar to $\mathsf{X}$ and such that $\{\mathsf{X}_i\}_{i=1}^{n}$ forms a partition of $\mathsf{X}$ in the sense of \eqref{c1} and \eqref{c2}. We leave it to the diligent reader to derive closed expressions for the $n$ injections $l_i$.

On the left Banach space $\mathcal{B}(\mathsf{X},\mathbb{H})$, we consider the RB operator $T: \mathcal{B}(\mathsf{X},\mathbb{H})\to\mathcal{B}(\mathsf{X},\mathbb{H})$ given by
\begin{equation}\label{6.5}
Tf(l_i(x)) := q_i(x) + s_i(x) f(x), \quad x\in \mathsf{X}, \;\;i \in \mathbb{N}_n,
\end{equation}
where $q_i, s_i:\mathsf{X}\to\mathbb{H}$ are bounded functions. Clearly, $T f\in \mathcal{B}(\mathsf{X},\mathbb{H})$ under these assumptions. Note that we left-multiply $f$ by $s_i$.

\begin{remark}
In the case of real Banach spaces, the RB operator is affine, i.e., $T - T(0)$ is a linear operator. In the quaternionic setting this is no longer true: $T - T(0)$ is not a left linear operator. It is only if $s_i (x) \in \mathbb{R}$.
\end{remark}

To simplify notation, we introduce the following abbreviation for the $m$-fold composition of functions from an IFS. Let $\mathcal{F}:=\{f_1, \ldots, f_n\}$. We write
\[
f_{i_m i_{m-1} \cdots i_1} := f_{i_m}\circ f_{i_{m-1}}\circ f_{i_1},
\]
where each $i_j\in \mathbb{N}_n$.

For each $m\in \mathbb{N}$, the $m$-fold application of $T$, $T^m f := T(T^{m-1} f)$, can be written as
\[
T^m f (l_{i_m i_{m-1} \cdots i_1}(x)) = \sum_{k=1}^{m} \prod_{j=1}^{k-1} s_{i_j} (x) q_{k} (x)+ \prod_{k=1}^m s_{i_k}(x) f(x),
\]
where the factors in the products $\prod$ are left-multiplied and where we set the empty product equal to 1. Notice the reverse order of the indices $i_m i_{m-1} \cdots i_1$ on the left- and right-hand side.

Continuing as in the previous section to ensure the existence of a fixed point $\psi$ for the RB operator \eqref{6.5} or a bounded solution of the associated system of functional equations, we obtain in summary the next result.
\begin{theorem}
For the above setting, the RB operator $T$ defined in \eqref{6.5} has a unique fixed point $\psi\in\mathcal{B}(\mathsf{X},\mathbb{H})$, i.e.,
\begin{equation}
T\psi = \psi\quad\Longleftrightarrow\quad\psi (l_i(x)) = q_i(x) + s_i(x) \psi(x), \quad x\in \mathsf{X}, \;\; i \in \mathbb{N}_n,
\end{equation}
provided that
\[
\max\limits_{i\in \mathbb{N}_n} \sup\limits_{x\in \mathsf{X}_i} |s_i(x)| < 1.
\]
\end{theorem}
The fixed point $\psi$ is then called a \emph{bounded quaternionic fractal function}.

\begin{remark}
A similar result exists of course for right Banach spaces and RB operators where the functions $s_i$ are right-multiplied onto $f$. Note that the fixed point $\psi$ depends on this right or left multiplication.
\end{remark}

\begin{example}
This example shows the versatility of quaternionic fractal interpolation. We choose as $\mathsf{X} := \{q\in\mathbb{H} : \Sc q \in [0,1) \wedge \Ve q = 0\}$. Note that $\mathsf{X} \cong [0,1] \subset\mathbb{R}$. Define injections $l_i:\mathsf{X}\to\mathsf{X}$ as follows:
\[
l_1 (x) : = \tfrac12 x\quad \text{and} \quad l_2 (x):= \tfrac12(x+1).
\]
Clearly, $\{l_1(\mathsf{X}), l_2(\mathsf{X})\}$ is a partition of $\mathsf{X}$ satisfying \eqref{c1} and \eqref{c2}.

Moreover, let $q_1 := e_0 + 2e_1-e_3+3e_4$ and $q_2 := -e_0 - 2e_1 + 2e_3+e_4$ be two quaternions. We set
\[
q_1 (x) := (1-q_1) x \quad \text{and} \quad q_2(x) := q_2 x^2.
\]

Then, $q_1,q_2\in \mathcal{B}(\mathsf{X},\mathbb{H})$. We define an RB operator $T$ by
\[
Tf (x) := \begin{cases} 2(1-q_1) x + s_1 f(2x), & x\in [0,\frac12),\\
q_2(2x-1)^2 + s_2 f(2x-1), & x\in [\frac12,1),
\end{cases}
\]
where $s_1 := \frac{1}{10}e_0 + \frac12 e_1 - \frac15 e_2 - \frac{1}{10} e_3$ and $s_1 := -\frac{1}{5}e_0 + \frac15 e_1 - \frac35 e_2 + \frac{1}{10} e_3$. Note that $\abs{s_1} = \frac{1}{10}\sqrt{31}$ and $\abs{s_2} = \frac{1}{10}\sqrt{45}$. Thus, $\max\{s_1,s_2\} < 1$.

Therefore, the RB operator $T$ is a contraction on $\mathcal{B}(\mathsf{X},\mathbb{H})$ and has a unique fixed point $\psi\in\mathcal{B}(\mathsf{X},\mathbb{H})$. The following Figure \ref{fig5} shows the projection of the graph of $\psi$ onto the $(e_0,e_i)$-plane for $i=0,1,2,3$.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=4cm]{qf1.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{qf2.pdf}\\
\includegraphics[width=5cm, height=4cm]{qf3.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{qf4.pdf}
\end{center}
\caption{The projections of $\psi$ onto the $(e_0,e_i)$-planes.}\label{fig5}
\end{figure}
As $\psi$ can be written as $\psi = \sum\limits_{i=0}^3 \psi_i e_i$, we display in Figure \ref{fig6} the parametric plots $(\psi_0, \psi_1, \psi_2)$ and $(\psi_0, \psi_2, \psi_4)$.
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5cm, height=4cm]{qf123.pdf}\hspace*{1cm}\includegraphics[width=5cm, height=4cm]{qf124.pdf}\end{center}
\caption{Some parametric plots of the components of $\psi$.}\label{fig6}
\end{figure}
\end{example}
\section{Summary}
In this chapter, we considered several aspects of fractal interpolation and identified their commonalities. Fractal interpolation is based either on solving a system of functional equations \cite{SB,SB1} or on obtaining the fixed point of an operator \cite{massopust,massopust1}. Both approaches yield the same results and are interchangeable.

In particular, we presented
\begin{itemize}
\item global fractal interpolation. This was the original approach undertaken in \cite{B2} which is based on a geometric setting. The analytic setting based on an RB operator commenced in \cite{bedford,dubuc1,dubuc2} and produced numerous generalizations and extensions. The main idea is to use a fixed bounded subset $\mathsf{X}$ of a normed space $\mathsf{E}$ to produce a partition of $\mathsf{E}$ over which copies of functions defined on the global set $\mathsf{X}$ reside. Each of the partitioning sets is then iteratively partitioned into subsets thus producing the fractal nature of the limiting object.
\item local fractal interpolation. Here, we allow subsets of a bounded subset $\mathsf{X}$ of a normed space $\mathsf{E}$ to produce a partition of $\mathsf{X}$. These subsets can be repeated in the construction of the partition and need no longer be defined by contractive injections. Here, functions defined on the local subsets reside over partitions of $\mathsf{X}$. This type of local interpolation was successfully employed in fractal compression and image analysis. (Cf., for instance, \cite{BH}.)
\item non-stationary fractal interpolation. In this approach, we do no longer require that we keep the same quantities at each level of iteration but allow them to vary. This set-up is very similar to non-stationary subdivision and was investigated in this respect in \cite{DLM,LDV,m4}. This new methodology shows great potential for future investigations into the subject.
\item quaternionic fractal interpolation. In this novel approach, we kept the set-theoretic structure but allowed for a change in the algebraic structure of the underlying sets. By requiring that we are using a division algebra such as the quaternions or more generally hypercomplex numbers, we obtain even more flexibility in the construction of fractal interpolants. The investigation into this type of setting commenced in \cite{m5} and promises to be a fruitful field of investigation in years to come.
\end{itemize}


\section{INTRODUCTION}
Robot manipulation systems (e.g. bin-picking, depalletization, loading/unloading) are rapidly gaining attention in manufacturing sites and logistic industries
due to the high potential upon revolutionizing productivity and labor efficiency.
The robot vision of such systems often perceives the workplace, localizes the target object and detects its pose, where such information is utilized to decide the robot operation policy.
Hence, when objects become occluded, the robot may overlook or estimate wrong postures. If the robot fails to grasp, not only it will fail to complete the task but may also damage the object by dropping or having collisions.
Therefore, occlusions are one of the biggest issues in robot vision systems.
\begin{figure}[!]
\centering
  \includegraphics[width=0.45\textwidth]{figure/fig1.png}
  \captionsetup{font=footnotesize}
  \caption{Our proposed concept of a robot vision system with tilt-type mirror reflection is shown. While the conventional direct sensing suffers from the occlusion caused by the robot arm, by concatenating with the reflection sensing results, we perceive all objects in the workplace with a single sensor.}
  \label{fig-concept}
\end{figure}

Robot vision systems to date have attempted to reduce occlusions as much as possible.
As a rule-of-a-thumb, occlusions in the scene can be minimized by mounting depth sensors on the zenith of the robot workplace, which many commercial robot systems follow \cite{mujin}\cite{righthand}\cite{toshiba}.
However, even with careful arrangements, an occlusion is always created directly below the robot arm, hiding objects underneath.
In our paper we define sensor blind spots as occlusions, which may be created by the robot itself or by another working robot.
One of the future goals in robot automation is a configuration in which multiple robots cooperate to perform one or more tasks simultaneously, which can significantly boost the work speed \cite{dual1}\cite{dual2}.
The impact of occlusions is further critical in such cases and is essential to construct a robot vision system enabling a further reduction of occlusions.

The most effective way to reduce the occlusion is to simply increase the number of sensors and sense from multiple angles.
However, since high-precision depth sensors (active projection cameras \cite{ensenso}, LiDARs \cite{velodyne}\cite{yoshioka2018}) are far more expensive than 2D cameras, the number of depth sensors should be minimized.
Also, since the optimal sensor placement will depend on the configuration of the cooperative robots, it is necessary to change the sensor placement every time the robot work content changes which lack flexibility.

In this paper, we propose a robot vision system which diminishes occlusions with tilt-type mirror reflection sensing. By combining both the conventional direct sensing result from the zenith sensor and the non-line-of-sight (NLoS) sensing result with mirror reflections, we acquire sensing results with significantly reduced occlusions. The two sensing results are obtained by rapidly configuring the sensor angle with a tilt-unit. Our system detects the occlusions within the workplace and the mirror reflection scan angle is adaptively set to effectively diminish the occlusions. The only required additional hardware is the tilt-unit and a commercially available mirror: compared to adding a depth sensor, the system cost can be greatly reduced.

Our original contributions are summarized bellow:
\begin{itemize}
\item To the best of our knowledge, this is the first robotic vision system that can adaptively remove occlusions utilizing mirror reflection and sensor tilting.
\item We propose a signal processing pipeline which achieves the concatenation of direct and mirror reflection sensing data with tilt-type sensing.
\item We propose a mirror displacement calibration that can be performed fully automatically without human interaction or oracles.
\item To evaluate the effectiveness of reflection sensing, we create a custom dataset which envisions multi-arm package picking. Compared to sensing only with a single zenith sensor, we confirm that reflection sensing improves the detection accuracy by 24\% and the accuracy is almost a par with the results obtained by using two sensors.
\end{itemize}

\section{Related Researches}
\textbf{Robot vision systems tackling occlusions:}
There are two major approaches upon reducing occlusions in robot vision systems, where one approach is to install multiple sensors and sense from multiple angles \cite{multicam}\cite{multicam2}\cite{multicam3}. Although this method do not prolong the sensing time, the sensor costs increase significantly.
Another approach is by attaching the sensor to the end-effector and controlling the arm to sense the work space from multiple angles, occlusions can be greatly reduced \cite{multiangle1}\cite{multiangle2}\cite{Holz2015}. However, e.g. ref.\cite{multiangle1} senses from 15 predefined viewpoints per scene; performing multiple-sensing may become time-consuming and reduce work efficiency. To summarize, conventional approaches results in a tradeoff between work efficiency and sensor cost.

Our proposed tilt-type system has the following advantages:
1) Compared to end-effector configurations, which moves the robot arm to several positions to capture multiple-sensing results, tilting the sensor itself can be enabled with significantly smaller movements which improve work efficiency. 2) Since the only required additional components are the tilt-unit and an ordinary mirror, the hardware cost increase can be marginal.

\textbf{Non-line-of-sight(NLoS) imaging:}
Our work belongs to the family of NLoS imaging, which uses reflection or diffusion to sense occluded objects\cite{nlosreview}. Since first demonstrated by ref.\cite{kirmani2009looking}, which proposed the concept of recovering occluded objects from time-of-flight data, NLoS have been applied to various medium such as laser light, acoustic beam\cite{acoustic}, and mm-wave radars\cite{scheiner2020seeing}. While it has been known that NLoS imaging can be realized via mirror reflections, to our best knowledge, we are the first to apply NLoS to adaptively reduce occlusions of the robot workspace.

\textbf{Processing depth sensor data under existence of mirrors:}
Mirrors or glass corrupt the sensing results of depth sensors, because the light unintentionally reflects on such surfaces but still returns to the sensor.
A "virtual" point cloud is generated as if an object exists behind the mirror and cause a negative impact on SLAM and mapping tasks. There have been researches to compensate for such effects \cite{Yang2008}\cite{Yang2011}\cite{Kim2016}\cite{Yun2018}\cite{Kashammer2015}, where a typical approach is to utilize the mirror symmetry and match the virtual point clouds to the real point clouds.

In our problem setup, the mirror is installed at a known location and thus, detection of the reflection data is not required. Also, while the recovery of the virtual reflection data is necessary in our system, our signal processing is original because previous works do not utilize sensor tilting.

\textbf{Improving robot vision system performance with mirror reflections:}
Our work is similar to ref.\cite{kulpate2005eye}\cite{Marchand2017} which realize visual servoing through mirror reflections and ref.\cite{okumura20151} \cite{Noguchi2012} which uses a mirror to realize a imaging system with a virtually wide FoV. However, while prior works aim to expand the FoV with fixed camera and mirror, our goal is to \emph{adaptively} reduce occlusions in robot vision systems by tilting the sensor and is clearly different.

\section{Tilt-type Robot Vision System with Mirror Reflections}
\subsection{System Concept}
Installing a depth sensor at the zenith of the robot workplace is a common configuration and provides fine sensing quality.
However, occlusions cannot be completely eliminated with a single sensor setup and the remaining occlusions may cause a major issue.
For example, if multiple robots work cooperatively in a single workplace, occlusions will continuously appear and may cause oversights. Fig.\ref{fig-concept} illustrates the occlusion caused by a robot arm: from the zenith sensor, the tomato below the arm is completely occluded and such object is impossible to detect.

We propose a robot vision system which diminishes occlusions with tilt-type mirror reflection sensing, and Fig.\ref{fig-concept} shows our main concept.
We define the sensing strategy where the workplace is directly sensed from the depth sensor as "direct sensing", and "reflection sensing" as to the NLoS sensing via the tilted depth sensor and the reflected light ray from the mirror.
By concatenating these two sensing results, we realize a robot vision system with highly diminished occlusions, since the "reflection sensing" provide results as if an extra sensor is installed in the mirror position.
Actual sensing results are shown in the bottom of Fig.\ref{fig-concept}, where the two target objects in the workplace (tomato and green pepper) are clearly visible.
Our system attaches a commercial tilt-unit to the sensor, enabling the "direct" and "reflection sensing" to be freely configured at high speeds by tilting.

\begin{figure}[!]
\centering
  \includegraphics[width=0.5\textwidth]{figure/fig2.png}
  \captionsetup{font=footnotesize}
  \caption{The signal processing pipeline for the tilt-based mirror reflection system is shown.
  First, the direct sensing is performed and from its results, the occlusion area is detected. Utilizing the detection results, we set the tilt angle so that the occlusion can be efficiently reduced.}
 \label{fig-signal}
 \end{figure}

 \begin{figure}[!]
\centering
  \includegraphics[width=0.5\textwidth]{figure/detect_occlusions.png}
  \captionsetup{font=footnotesize}
  \caption{The occlusion detection pipeline is shown. The occlusion area created by the robot arm is detected based on regions exceeding the height threshold.}
 \label{fig-detect}
 \end{figure}

\subsection{Signal Processing Pipeline}
In actual robot automation systems, locations of occlusions vary greatly with the work content ( e.g. operation scheme, placements) and "reflection sensing" with fixed tilt angles cannot remove occlusions sufficiently.
To handle such problems, we propose a signal processing pipeline which adaptively detect and diminish occlusions: the occlusion area is detected from the "direct sensing" results and the sensor tilt angle is adaptively set to minimize occlusions during "reflection sensing".

Fig.\ref{fig-signal} shows our overall signal processing pipeline.
Firstly, "direct sensing" is performed to sense the workplace, where the occlusions are  detected from its results. To mitigate the occlusions, the tilt-unit is rotated to perform "reflection sensing".
After converting the "reflection sensing" data to world coordinates, the mirror virtual images are converted to the real image to enable the concatenation of the direct and reflection data.
Finally, object detection is carried out using the concatenated data with diminished occlusions.
The following sections describe the details of each block.

\textbf{Detecting occlusions:}
We utilize a simple yet effective rule-based approach to detect the occluded area from the direct sensing data as shown in Fig.\ref{fig-detect}.
A height threshold is defined as the maximum height of an object allowed in the work place with some margins, and any object surpassing the threshold is classified as an occlusion area.
Then, $N$ largest occlusion areas are determined as the occlusions and others are filtered out to suppress sensor noise effects, where $N$ is the known number of robots in the workplace.
Finally, the centroid of the occlusion is derived to set the optimum sensor tilt angle.

In our system, we assume a robot manipulation task where the robot grasps the target object from above and then move to an another location.
In such tasks, the robot arm should approach from a position sufficiently higher than the object itself to avoid collisions. The arm elbow and wrists, which mainly generates the occlusions, will also take a sufficiently higher position than the surrounding objects and thus, occlusions can be sufficiently detected with height thresholds.
While image recognition techniques (e.g. DNNs) can be utilized to detect occlusions, its execution time is overwhelmingly long.
Since the latency of occlusion detections directly prolongs the sensing time, it may lead to degrading the robot working efficiency.

\begin{figure}[!]
\centering
  \includegraphics[width=0.35\textwidth]{figure/scanangle.png}
  \captionsetup{font=footnotesize}
  \caption{The derivation of $\theta_{tilt}$ with mirror reflections are shown. By converting the target coordinates to virtual coordinates $X_{tv}$, the optimal $\theta_{tilt}$ can be calculated by connecting the sensor and $X_{tv}$ with a straight line.}
  \label{fig-scan}
\end{figure}

\textbf{Setting the optimal tilt angle $\theta_{tilt}$:}
Fig.\ref{fig-scan} shows the method to set the tilt angle so that occlusions can be removed effectively.
In the case of conventional direct sensing, the target object is most efficiently captured when the tilt angle is set to form a straight line between the sensor and the target coordinate.
For simplicity, let us consider a two-dimensional space, where the sensor position is [$x_s$, $y_s$] and the target coordinate is [$x_t$, $y_t$]:
\begin{eqnarray}
    \centering
    \theta = \arctan(\frac{y_t - y_s}{x_t-x_s})
    \label{scanangle}
\end{eqnarray}
On the other hand, with mirror reflections three elements are required to set $\theta_{tilt}$: the sensor coordinate, the transformation matrix $H_m$ which converts the virtual reflection data to real data, and the target coordinate.
Let the target coordinate be $X_t$, the "virtual" mirror-imaged target coordinate $X_{tv}$. Such conversion can be realized with the inverse matrix $H_m^{-1}$:
\begin{eqnarray}
    \centering
    X_{tv} = H_m^{-1}X_t
\end{eqnarray}
Thus, as shown in Fig.\ref{fig-scan}, the optimal $\theta_{tilt}$ is set when $X_{tv}$ and the sensor is directly connected which can be calculated in the same manner as Eq.(\ref{scanangle}).

\begin{figure}[!]
\centering
  \includegraphics[width=0.5\textwidth]{figure/tilt-sensing.png}
  \captionsetup{font=footnotesize}
  \caption{(a) Definition of the world coordinate in our system. (b) The state of TiltedSensor coordinate with $\theta_{tilt}$ tilt.}
  \label{fig-tilt}
\end{figure}

\textbf{Converting tilt sensing reflection data to world coordinates:}
In our system, two coordinate translations are required: (1) translating tilted sensing data to the world coordinates, (2) translating virtual mirror reflection data to real data. Firstly, the translation of the world coordinates are explained.
In order to simplify the calculation, we define the world coordinate as shown in Fig.\ref{fig-tilt}(a). Note that the tilt rotation is around the Y-axis with the tilt-unit at the center.
Fig.\ref{fig-tilt}(b) illustrates an example when the sensor is tilted by $\theta_{tilt}$.  Since the data obtained after tilting is in the TiltedSensor coordinate system, it must be converted to the world coordinate system for processing. Such transformation matrix $H_S$ can be formulated where the the rotation radius of the tilt-unit is $r$:
\[
  H_{S} = \left(
    \begin{array}{cccc}
      \cos\theta_{tilt} & 0 & -\sin\theta_{tilt} & -r\sin\theta_{tilt} \\
      0 & 1 & 0 & 0 \\
      \sin\theta_{tilt} & 0 & \cos\theta_{tilt} & -r(1-\cos\theta_{tilt}) \\
      0 & 0 & 0& 1
    \end{array}
  \right)
\]
Since $\theta_{tilt}$ is given from the system and $r$ is determined by the mechanical of the tilt-unit, both of the values are known beforehand.

In addition, as shown in Fig.\ref{fig-scan}, the data sensed through the mirror forms a "virtual image" and must be converted to a "real image" by mirror-image transformation. As formulated with the householder transformation\cite{householder1958unitary}, when the mirror plane is represented by
\begin{eqnarray}
    \centering
    ax + by + cz + d = 0
\end{eqnarray}
and the L2 norm of a, b and c is unity, the mirror-image transformation can be expressed as:
\[
  H_{m} = \left(
    \begin{array}{cccc}
      1-2a^2 & -2ab & -2ac -2ad & -2ad \\
      -2ab   & 1-2b^2 & -2bc & -2bd \\
      -2ac   & -2bc   & 1-2c^2 & -2cd \\
      0 & 0 & 0& 1
    \end{array}
  \right)
\]

\textbf{Sensing error analysis for reflection sensing:}
In reflection sensing, the working distance of the sensor equivalently increases and increases the sensor measurement error. When the working distance of $\theta_{tilt}=0$ is 1, the distance to the object for reflection sensing ($D_{reflect}$) is:
\begin{eqnarray}
    \centering
    D_{reflect} = \frac{1}{\cos\theta_{tilt}}
\end{eqnarray}
While measurement errors of active 3D sensors depend on its architecture (e.g. stereo vision, dToF, iToF) and its lighting conditions, we can expect that the measurement error increase linearly or square to the working distance \cite{ensenso,yoshioka2018,he2017depth}.

Since the light decays with mirror reflection, the mirror quality can also become an error source. When the mirror reflectance is $\alpha$, the light emitted from the sensor will be attenuated by $1/\alpha^2$ compared to the direct path. Note that typical commercial mirrors have a high reflectance of $>$90\%; the distance-dependent error is dominate in our system.

\section{Automatic Mirror Displacement Calibration}

\begin{figure}[!t]
  \begin{algorithm}[H]
  \caption{Psuedocode to obtain the calibration-optimal robot pose}
\begin{algorithmic}
 \For{$joint=Joint_1,\cdots,Joint_n$}
   \State $maxpoints \leftarrow 0$
   \For{$angle=Angle_1, \cdots, Angle_n$}
     \State $SetRobotArmPose(joint, angle)$
     \State Obtain Direct and Reflection Sensing results
     \State Calculate $N_{points}$
     \If{$N_{points} > maxpoints$}
         \State $BestAngle \leftarrow angle$
         \State $maxpoints \leftarrow N_{points}$
     \EndIf
  \EndFor
  \State $SetRobotArmPose(joint, BestAngle)$
\EndFor
\end{algorithmic}
\end{algorithm}
\end{figure}

As described in section III, our mirror reflection vision system rely on the transformation matrix $H_m$ to process mirror reflection data and to set $\theta_{tilt}$.
However, with disturbance such as factory vibrations, the position/angle of the mirror may shift, altering $H_m$ from the predefined value.
Even with a mirror angle shift of only few degrees, significant error is introduced to the sensing results; an automatic mirror displacement calibration is necessary.
Moreover, while human-guided calibration (e.g. covering the mirror by non-reflective cloth) are easy to conduct, it is not desirable since it will persistently interrupt the robot operation.

To achieve automatic calibrations, we propose a calibration method which uses the robot itself as the calibration target.
$H_m$ is derived by registering the "virtual" point clouds (obtained from reflection sensing) to "real" point clouds (from direct sensing). However, to achieve high registration accuracy, it is important to configure a robot pose which have a large number of data in common between the two sensing results.
Since it is challenging to perceive the actual amount of common data between the two sensing results, we follow the intuition: when the same object is captured, the number of common data should scale with the total number of the sensed data.
Algorithm 1 shows the flow upon deciding the calibration robot pose, where the robot's joints are configured to find the posture which maximizes the total number of robot arm point cloud data ($N_{points}$). Note that $N_{points}$ is the sum of direct ($N_{directpoints}$) and reflection sensing results ($N_{reflectpoints}$), and the robot arm point clouds are obtained by height threshold filtering explained in Fig.\ref{fig-detect}.
Since the ray trajectory is longer in reflection sensing, the point cloud is inherently sparse and we aid $N_{reflectpoints}$ by $\alpha$ to mitigate the sparsity effects.
\begin{eqnarray}
    \centering
    N_{points} = N_{directpoints} + \alpha * N_{reflectpoints}
\end{eqnarray}
During experiments, we set $\alpha = 2$.
Importantly, only robot arm joints which have significant impact to $N_{points}$ should be searched to speed up the calibration process. During experiments, we search the optimal pose using 2 largest joints (shoulder and elbow) and 10 angle settings to secure sufficient calibration accuracy.

\section{Experiments}
\begin{figure}[!]
\centering
  \includegraphics[width=0.4\textwidth]{figure/experiments_v2.png}
  \captionsetup{font=footnotesize}
  \caption{ (a) The experimental setup for the tilt-type mirror reflection robot vision system.
(b) The two sensor experimental setup used for comparison.
(c) Close up view of the sensor mounted on the tilt-unit.
}
  \label{fig-setup}
\end{figure}

Fig. \ref{fig-setup} shows our experimental setup and the setup is used for all experiments. The sensor is placed 2.1m above ground, where the FLIR PTR-E46 was used as the tilt-unit, one DENSO COBOTTA was used as a robot arm to create occlusions, and one commercially available mirror was placed 1.2m away from the sensor in front of the robot. We arbitrary move the robot location during experiments to confirm the robustness of our system.  Note that while we chose stereo camera based EnsensoN35 as the depth sensor due to its precision, we confirmed that LiDARs and iToF cameras can be used in our reflection system as well.

\subsection{Object detection under occlusions}
\begin{figure}[!]
\centering
  \includegraphics[width=0.5\textwidth]{figure/pipeline_dl.png}
  \captionsetup{font=footnotesize}
  \caption{The object detection pipeline is shown. The point clouds are first filtered by height to remove the robot arm and then converted to birds-eye-view images. Then, instance segmentation is performed to detect the cardboard.}
  \label{fig-pipe}
\end{figure}
\begin{figure*}[ht!]
\centering
  \includegraphics[width=0.95\textwidth]{figure/eyecatch_v3_light.png}
  \captionsetup{font=footnotesize}
  \caption{Qualitative results of the object detection. With only either direct or reflection sensing, occlusions still remain which cause detection errors (e.g. one box recognized as two boxes, wrong box size). The detection result with direct+reflection sensing has significantly lower detection errors and is comparable to the result using two sensors.}
  \label{fig-eyecatch}
\end{figure*}

\begin{table}[!]
\captionsetup{font=footnotesize}
\caption{The object detection accuracy and the target data coverage is evaluated on our custom dataset, where the robot continuously creates occlusions. The accuracy is significantly improved by concatenating mirror reflection sensing and the accuracy is a par with the results obtained with two sensors.}
\large
\scalebox{0.75}{
\begin{tabular}{|c|c|c|c|c|c|}
\hline
\multicolumn{2}{|c|}{\multirow{2}{*}{Dataset}}                                                                                & \multicolumn{4}{c|}{Sensing Strategy}                                                                                \\ \cline{3-6}
\multicolumn{2}{|c|}{}                                                                                & \begin{tabular}[c]{@{}c@{}}Direct \\ only\end{tabular} & \begin{tabular}[c]{@{}c@{}}Mirror \\ only\end{tabular} & \textbf{\begin{tabular}[c]{@{}c@{}}Direct\\ \&Mirror\end{tabular}} & \begin{tabular}[c]{@{}c@{}}Two \\ sensors\end{tabular} \\ \hline
\multirow{2}{*}{Easy}                                  & \begin{tabular}[c]{@{}c@{}}mAP\\ @50\end{tabular}                             & \begin{tabular}[c]{@{}c@{}}73.6\\ (-25.0)\end{tabular} & \begin{tabular}[c]{@{}c@{}}86.5\\ (-12.1)\end{tabular} & \textbf{98.6}                                                      & \begin{tabular}[c]{@{}c@{}}98.7\\ (+0.1)\end{tabular}  \\ \cline{2-6}
                                                       & \begin{tabular}[c]{@{}c@{}}mAP\\ @75\end{tabular}                             & \begin{tabular}[c]{@{}c@{}}55.6\\ (-39.8)\end{tabular} & \begin{tabular}[c]{@{}c@{}}83.5\\ (-11.9)\end{tabular} & \textbf{95.4}                                                      & \begin{tabular}[c]{@{}c@{}}96.4\\ (+1.0)\end{tabular}  \\ \hline
\multirow{2}{*}{Hard}                                  & \begin{tabular}[c]{@{}c@{}}mAP\\ @50\end{tabular}                             & \begin{tabular}[c]{@{}c@{}}66.9\\ (-24.7)\end{tabular} & \begin{tabular}[c]{@{}c@{}}85.1\\ (-6.3)\end{tabular}  & \textbf{91.4}                                                      & \begin{tabular}[c]{@{}c@{}}91.6\\ (+0.2)\end{tabular}  \\ \cline{2-6}
                                                       & \begin{tabular}[c]{@{}c@{}}mAP\\ @75\end{tabular}                             & \begin{tabular}[c]{@{}c@{}}39.2\\ (-51.4)\end{tabular} & \begin{tabular}[c]{@{}c@{}}74.3\\ (-16.3)\end{tabular} & \textbf{90.6}                                                      & \begin{tabular}[c]{@{}c@{}}90.6\\ (+0)\end{tabular}    \\ \hline
\multirow{2}{*}{Averaged}                              & \begin{tabular}[c]{@{}c@{}}mAP\\ @50\end{tabular}                             & \begin{tabular}[c]{@{}c@{}}70.2\\ (-24.7)\end{tabular} & \begin{tabular}[c]{@{}c@{}}85.8\\ (-9.2)\end{tabular}  & \textbf{95.0}                                                      & \begin{tabular}[c]{@{}c@{}}95.1\\ (+0.1)\end{tabular}  \\ \cline{2-6}
                                                       & \begin{tabular}[c]{@{}c@{}}mAP\\ @75\end{tabular}                             & \begin{tabular}[c]{@{}c@{}}47.4\\ (-45.6)\end{tabular} & \begin{tabular}[c]{@{}c@{}}78.9\\ (-14.1)\end{tabular} & \textbf{93.0}                                                      & \begin{tabular}[c]{@{}c@{}}93.5\\ (+0.5)\end{tabular}  \\ \hline
\multicolumn{2}{|c|}{\begin{tabular}[c]{@{}c@{}}Mean target data\\ coverage (as compared \\ to two sensors) {[}\%{]}\end{tabular}} & 70.1                                                   & 83.9                                                   & \textbf{96.9}                                                      & 100                                                    \\ \hline
\end{tabular}
}
  \label{table}

\end{table}

The main purpose of this experiment is to demonstrate that our proposed mirror reflection vision system 1) can diminish occlusions and 2) can achieve similar detection accuracy as using two sensors.
The setup using two sensors is shown in Fig. \ref{fig-setup} (b), where the angle of sensor 2 was carefully adjusted to effectively remove the occlusion.

Here, we evaluate the object detection accuracy of the system by modeling a cardboard box manipulation task with collaborative robots.
For evaluation purposes, we create a dataset where we intentionally produce occlusion by randomly moving the robot arm over the cardboard boxes.
Under such circumstances, we acquire data following four sensing strategies: 1) direct sensing only, 2) reflection sensing only, 3) direct + reflection sensing, and 4) two sensors. We collect 8 scenes with different cardboard arrangements with each scene containing 50 point cloud data.
The 4 scenes have only 1 to 3 cardboard boxes (named "easy scenes") and the other 4 scenes containing 4 to 6 boxes (named "hard scenes").
To investigate the effect of occlusions on the object detection accuracy, we construct an object detection pipeline for evaluation. Fig. \ref{fig-pipe} shows our object detection pipeline. Firstly, since the robot arm is not an object of interest, we filter the robot point clouds by utilizing the height threshold. Then we project the point clouds into a 2D birds-eye-view image, where the height is mapped as the image intensity.
Finally, we use resnet50 Mask-RCNN \cite{Mask} to conduct instance segmentation to detect the cardboard, which is fine-tuned with 800 cardboard images.

We summarize the qualitative object detection results in Fig.\ref{fig-eyecatch} and the quantitative results in Table \ref{table}.
The instance segmentation mean average precision (mAP)\cite{coco} are calculated, and results for IoU=50\% and IoU=75\% are reported respectively.
The effect of occlusions were very large with either direct or reflection sensing only and held poor accuracy.
Some failure cases were: the object was completely hidden by the robot arm and overlooked, or a single object was misinterpreted as 2 objects due to cavities caused by occlusions.
By the proposed direct + reflection sensing, the occlusion can be greatly reduced, boosting the detection accuracy. Importantly, our proposed method achieves a similar accuracy with results using two sensors, where the accuracy difference was only 1\%. We also report the mean target data coverage respect to the two sensors result as well (computed from point cloud geometry), and we confirm that the coverage ratio correlates well with the detection accuracy.

\subsection{Mirror displacement calibrations}
\begin{figure}[!]
\centering
  \includegraphics[width=0.4\textwidth]{figure/cal_points.png}
  \captionsetup{font=footnotesize}
  \caption{(a) Robot pose used for evaluation versus number of data. Here, Pose6 was obtained by our calibration method.
(b) Robot pose versus translational error. (c) Robot pose versus rotational error. }
  \label{fig-calpoints}
\end{figure}

\begin{figure}[!]
\centering
  \includegraphics[width=0.45\textwidth]{figure/cal_results.png}
  \captionsetup{font=footnotesize}
  \caption{We report the calibration accuracy when the mirror angle is swept to -5 to 5 degrees. By using the proposed optimal pose, the calibration converges on various conditions with high accuracy. On the right, we show point cloud data before and after the displacement calibration.}
  \label{fig-calresults}
\end{figure}
Finally, the proposed mirror displacement calibration method is evaluated by altering the mirror angles from the initial condition.
Here, we follow ref.\cite{geiger2012we} and evaluate the average calibration accuracy (translational error and rotational error) of ten runs. The initial transformation matrix $H_{minit}$ is set when the mirror angle is placed perpendicular to the ground, and our calibration is carried out utilizing $H_{minit}$ as the initial value. For registration, we use standard RANSAC+ICP registration pipeline powered by Open3D\cite{open3d}.

Fig.\ref{fig-calpoints} reports the robot pose's impact to calibration accuracy.
Pose 6 was obtained using our proposed method, by configuring the robot pose to maximize $N_{points}$. On the other hand, Pose 1-5 are random robot postures with different $N_{points}$, shown for comparison.
As shown in Fig.\ref{fig-calpoints}(a), Pose 6 contains most point cloud data, and $N_{points}$ are smaller in Pose 1-5.
Fig.\ref{fig-calpoints}(b)(c) shows the calibration accuracy: Pose 1-3 which have insufficient $N_{points}$ concluded with degraded calibration accuracy. On the other hand, as $N_{points}$ increase in Pose 4,5, we saw an improvement in accuracy. The best calibration accuracy was achieved with Pose 6, showing that our $N_{points}$-based criteria contributes to improving the calibration accuracy.

Fig.\ref{fig-calresults} reports the mirror displacement calibration accuracy with various conditions, where the mirror angle was altered from -5 to 5 degrees. Here, we compare the calibration accuracy between the optimal pose (Pose6) and Pose 4.
In Pose 4, the calibration convergence is unstable: for some mirror angle situations, the calibration cannot converge (mirror angle $<$-2).
On the other hand, by using the optimal pose, regardless of various mirror angle situations, high calibration accuracy was achieved.

\section{CONCLUSIONS}
We proposed the first robot vision system with tilt-type mirror reflection sensing.
By concatenating the direct sensing and NLoS mirror reflection sensing results, the occlusions are greatly diminished. Our signal processing pipeline detects the occlusion area and dynamically configures the sensor tilt angle, allowing the system to adaptively remove the workplace occlusion. Through experiments, we confirmed that our system can achieve the same detection accuracy as that of multiple sensors, regardless  of  the  single-sensor  implementation.

In future works, we plan to configure a cooperative robot environment with multiple robots and evaluate our proposed vision system under such scenes. Since occlusions created in such environments are more complex, multiple mirrors should be installed and techniques to select the optimal mirror upon sensing will be required.
Moreover, we plan to evaluate the impact to robot grasp accuracy with the use of mirror reflections.

\bibliographystyle{IEEEbib}