\section{Introduction}\label{sec:intro}
\input{intro-paths.tex}

\section{Complex hypergraphs}\label{sec:complex}
\input{main-thm.tex}
\input{condensation.tex}

\section{Tools}\label{sec:tools}
\input{tools.tex}

\section{Shrinking the hypergraph}\label{sec:shrink}
\input{shrink.tex}

\section{Path-graphs}\label{sec:path-graphs}
\input{conflict-2.tex}
\input{rainbow-trees.tex}

\section{Proof of a key lemma}\label{sec:core}
\input{core-2.tex}

\section{Proof of Theorem~\ref{thm:main}}\label{sec:paths}
\input{main-proof.tex}

\section*{Acknowledgements}
We would like to thank to Julien Bensmail, Ararat Harutyunyan and Tien-Nam Le
for helpful discussions.

The work was done while the first author was a postdoc at Laboratoire d’Informatique du Parall\'elisme,
\'Ecole Normale Sup\'erieure de Lyon,
69364 Lyon Cedex 07, France.

\section{Graph fractions and other tools}
In this section, we introduce several tools used in the rest of the paper. First, we focus on connectivity in hypergraphs and state generalizations of some classical results for graphs.

The following Lemma~\ref{cor:sparse-tree} was essentially proven in~\cite{bib:thomassen-paths}. Our version (which we prove in~\cite{bib:our-trees}) differs in some details.

\begin{lemma}\label{cor:sparse-tree}
For every $\varepsilon>0$ and integer $m$, there exists $L$ such that every $2^n$-edge-connected stub graph $G$ with minimum degree at least $L$, where $n=2+m+ \lceil\log (1/\varepsilon)\rceil$,
has $2^{m}$ edge-disjoint spanning trees $T_1,\ldots, T_{2^m}$ such that
\[\sum_{i\in [2^m]}\deg_{T_i} v \leq \varepsilon \deg_{G} v\]
 for every $v\in V(G)$.
\end{lemma}

A {\em bridge} in a graph or a hypergraph is an edge such that after removing it, the graph or the hypergraph becomes disconnected.

A hypergraph $H$ is {\em $k$-partition-connected} if for any partition of $V(H)$ into $p$ non-empty parts, the number of hyperedges
of $H$ containing vertices in more than one part is at least $k(p-1)$. The following easy observation links partition-connectivity to edge-connectivity.

\begin{observation}\label{obs:part-con}
A $3k$-edge-connected hypergraph $H$ with edges of size at most three is $k$-partition-connected.
\end{observation}

\begin{proof}
Let $\HH$ be a partition of $H$ into $p$ nonempty parts. Since $H$ is $3k$-edge-connected, every part of $\HH$ is connected to other parts by at least $3k$ edges. Each such edge contains vertices in at most three parts of $\HH$, therefore in total there are at least $3kp/3=kp$ edges containing vertices in more than one part of $\HH$.
\end{proof}

A {\em spanning hypertree} of a hypergraph $(V,E)$ is a hypergraph $(V,F)$, where $F\subseteq E$, $|F|=|V|-1$, $\bigcup_{e\in F}e=V$ and there are at most $|X|-1$ hyperedges $e\in F$ such that $e\subseteq X$ for every $X\subseteq V$.
The classical result of Nash-Williams and Tutte about the number of edge-disjoint spanning trees
 was generalized to hypergraphs in~\cite{bib:nwh} as follows.

\begin{theorem}[Nash-Williams for hypergraphs~\cite{bib:nwh}]\label{thm:nwh}
A hypergraph $H$ is $k$-partition-connected if and only if $H$ contains $k$-hyperedge-disjoint spanning hypertrees.

\end{theorem}

Thus, every $1$-partition-connected hypergraph contains a spanning hypertree. By a result of Lov\'asz~\cite{bib:lovasz}, it is possible to shrink a spanning hypertree in such a way that the resulting graph is a spanning tree.

Let $H$ be a hypergraph. We call an ordered partition $(A,B)$ of the vertex set of $H$ into two parts a {\em cut} in $H$. We denote by $E(A,B)$ the set of hyperedges of $H$ incident with vertices in both $A$ and $B$ and call $|E(A,B)|$ the {\em order of the cut} $(A,B)$.
We define the {\em value} of the cut $(A,B)$ as $val(A,B)=\sum_{e\in E(A,B)}|e\cap A|$. Note that this definition is asymmetric, i.e., $val(A,B)$ and $val(B,A)$ are not equal in general. Next, we prove an analogue of Lemma~\ref{lem:cut} for hypergraphs.

\begin{lemma}\label{lem:hycut}
Let $k$ be an integer and $H$ be a hypergraph with edges of size at most three. Then there is a cut $(A,B)$ in $H$ of value at most $3k$ (and thus, of order at most $3k$) such that $H[A]$ is $k$-edge-connected or has only one vertex.
\end{lemma}

\begin{proof}
If $H$ is $k$-edge-connected, we are done. Otherwise, we consider a cut $(A,B)$ with value at most $3k$ such that $A$ is inclusion-wise minimal. (Such a cut exists in $H$, because $val(A,B)\leq 2|E(A,B)|$ for every cut $(A,B)$.) Consider a cut $(A_1,A_2)$ in $H[A]$ and let $r$ be its order. We show that $r>k$. We have
\[ val(A_1,V\setminus A_1) + val(A_2,V\setminus A_2)=val(A,B)+\sum_{e\in E(A_1,A_2)}|e|.\]
Note that $val(A,B)\leq 3k$ and $\sum_{e\in E(A_1,A_2)}|e|\leq 3r$. Moreover, $val(A_1,V\setminus A_1), val(A_2,V\setminus A_2)>3k$ by minimality of $A$.  It follows that $\sum_{e\in E(A_1,A_2)}|e|>3k$ and thus, $r>k$. Since this is true for any cut in $H[A]$, $H[A]$ is $k$-edge-connected.
\end{proof}

Let $G$ be a hypergraph (or a multigraph). We denote the minimum number of hyperedges in a cut separating vertices $u$ and $v$ in $G$ by $\lambda_{G}(u,v)$.

{\em Splitting off} a pair of edges $su$, $sv$ in a multigraph $G$ is replacing these two edges by a new edge $uv$. We define splitting off a pair of hyperedges of size two in a hypergraph in the same way.

\begin{theorem}[Mader splitting off theorem~\cite{bib:mader}]\label{thm:mader}

Let $G=(V,E)$ be a connected multigraph and let $s$ be a vertex such that the degree of $s$ is not three and $s$ is not incident with bridges.
Then there are two edges $e,f$ incident with $s$ such that the multigraph $G_{e,f}$ obtained by splitting off $e$ and $f$ satisfies $\lambda_{G_{e,f}}(u,v)=\lambda_{G}(u,v)$ for any $u,v\in V\setminus \{s\}$.
\end{theorem}

\begin{corollary}\label{cor:mader}
Let $H$ be a connected hypergraph with at least two vertices and hyperedges of size at most three and let $s$ be a vertex that is not contained in any hyperedge of size three and not incident with bridges. Then there are two hyperedges $e,f$ of size two incident with $s$ such that the hypergraph $H_{e,f}$ obtained from $H$ by splitting off $e$ and $f$ satisfies $\lambda_{H_{e,f}}(u,v)=\lambda_{H}(u,v)$ for any $u,v\in V\setminus \{s\}$.
\end{corollary}

\begin{proof}
Let $G$ be a hypergraph with hyperedges of size at most three.
Let $\tau$ be the operation creating a multigraph $\tau(G)$ from $G$ by replacing each hyperedge $h$ of size $3$ containing vertices $u_1,u_2$ and $u_3$ by a vertex $h$ adjacent to $u_1$, $u_2$ and $u_3$ (with possible multiplicities if $u_i=u_j$ for some $i\not= j$).
Observe that $\lambda_{G}(u,v)=\lambda_{\tau(G)}(u,v)$ for every $u,v\in V(G)$. In particular, if $s$ is not incident with a bridge in $G$, it is not incident with a bridge in $\tau(G)$.

Therefore, we can apply Theorem~\ref{thm:mader} to $\tau(H)$ and the vertex $s$. Let $e,f$ be a pair of edges incident with $s$ in $\tau(H)$ as in Theorem~\ref{thm:mader} and let $\tau(H)_{e,f}$ be the multigraph obtained from $\tau(H)$ by splitting off $e$ and $f$. Then, $\lambda_{\tau(H)_{e,f}} (u,v)=\lambda_{\tau(H)}(u,v)$ for every $u,v\in V(\tau(H))\setminus \{s\}$.

Note that since $s$ is not contained in any hyperedge of size three in $H$, $e$ and $f$ are hyperedges of size two in $H$. Moreover, $\tau(H)_{e,f}$ is isomorphic to $\tau(H_{e,f})$. It follows that $\lambda_{H_{e,f}}(u,v)=\lambda_{H}(u,v)$ for any $u,v\in V(H)\setminus \{s\}$.
\end{proof}

Let $G$ be a complex graph without stubs. We say that a spanning subgraph $H$ of $G$ is an {\em $\varepsilon$-fraction of $G$} if it satisfies $|\deg_H v- \varepsilon\deg_G v|\le 2$ for every vertex $v$ of $G$.
The following lemma shows that there exists an $\varepsilon$-fraction of $G$ for any  $\varepsilon=1/2^c$, $c\in \NN$.

\begin{lemma} \label{lem:fraction}
Let $m_1, \ldots, m_k$ be such that $\sum_{i\in [k]} 2^{-m_i}\leq 1$ and let $G$ be a complex graph without stubs. Then, there exist edge-disjoint subgraphs $G_1,\ldots,G_k$ of $G$ such that $G_i$ is a $2^{-m_i}$-fraction of $G$ for every $i\in [k]$.
\end{lemma}

\begin{proof}

Let $m=\max_{i\in [k]} m_i$
and let $H$ be the underlying multigraph of $G$.

Let $\bar{H}$ be the multigraph obtained from $H$ by adding an auxiliary vertex $s$ adjacent to all vertices of odd degree in $H$. Since the number of such vertices is even, all vertices of $\bar{H}$ have even degree.
Consider an Eulerian tour in each connected component of $\bar{H}$ and color its edges by two colors, $0$ and $1$, in an alternating way, starting from $s$ if the connected component contains $s$.

Let $H_0$ be the subgraph of $H$ of color $0$ and $H_1$ the subgraph of $H$ of color $1$. We apply the procedure described above for $H_0$ and $H_1$ to obtain graphs $H_{00}$, $H_{01}$, $H_{10}$ and $H_{11}$ and repeat it recursively until we obtain a collection $\{H_Q\}_{Q\in \{0,1\}^m}$ of edge-disjoint spanning subgraphs of $H$.

 Let $v$ be a vertex of $H$ and let $d$ be its degree. Observe that $\deg_{H_0} v,\deg_{H_1} v \in [d/2-1,d/2+1]$ and this is true for every iteration of the procedure. Thus,
  \[\deg_{H_Q} v\leq \underbrace{(\ldots (}_{i \times} d/\underbrace{2 +1)/2 +1)/2 \ldots )/2+1}_{i \times}\leq d/2^i+\sum_{j=0}^{i-1} 1/2^{j}\leq d/2^i +2\]
for every $H_{Q}$ with $Q\in \{0,1\}^{i}$.
 Analogously, $\deg_{H_Q} v \geq d/2^i-2$ for every $H_{Q}$ with $Q\in \{0,1\}^{i}$.

We choose $Q_i\in \{0,1\}^{m_i}$ for every $i\in [k]$ so that the graphs $H_{Q_i}$ are pairwise edge-disjoint. We let $G_i$ be the subgraph of $G$ with underlying multigraph $H_{Q_i}$ for every $i\in [k]$.
Consequently, all $G_i$'s are edge-disjoint $2^{-m_i}$-fractions of $G$.

\end{proof}

We conclude the section by stating the Asymmetric Lov\'asz Local Lemma, which we extensively use in the following sections.

\begin{lemma}[Lov\'asz Local Lemma (asymmetric version)~\cite{bib:lll}]\label{LLL}
Let ${\mathcal {A}}=\{A_{1},\ldots, A_{n}\}$ be a finite set of events in a probability space $\Omega$. For each $A\in {\mathcal {A}}$, we let $\Gamma (A)$ denote a subset of ${\mathcal {A}}$ such that $A$ is independent from the collection of events ${\mathcal {A}}\setminus (\{A\}\cup \Gamma (A))$. If there exists an assignment of reals $x:{\mathcal {A}}\to (0,1)$ to the events such that  for every $A\in {\mathcal {A}}$:
\[\Pr(A)\leq x(A)\prod_{B\in \Gamma (A)}(1-x(B)),\]
then the probability of avoiding all events in ${\mathcal {A}}$ is positive, in particular

    \[\Pr \left(\bigwedge_{A\in {\mathcal {A}}}{\overline {A}}\right)\geq \prod \nolimits _{A\in {\mathcal {A}}}(1-x(A)).\]
\end{lemma}


\section{Introduction} \label{sec:intro}

The present article deals with the question of global reconstruction of orientation fields on the basis of a discrete dataset.
The aim is to present an alternative way of modeling orientation fields that allows to
use a natural energy.

{ As continuous mathematical objects, orientation fields adequately model texture patterns predominantly displaying orientation information. They provide a  unifying framework for various patterns observed in nature such as fingerprints \cite{Mal-Mai2009,Kas-Wit1987,Wan-Hu2007a,Gup-Tiw2019}, liquid crystals arrangements in their nematic phase \cite{Cha1992,Dec-Red2015,Deg-Pro1995,Li-Shi2019} or the pinwheel structure of the visual cortex V1 of mammals \cite{pet2017,Hub-Wie1962,Bos-Che2014,Bos-Dup2012,Cit-Sar2006,Pra-Gau2018}.}

The model we want to discuss ties in with classical techniques used in the field of fingerprint reconstruction and authentication, therefore, the problem of the estimation of fingerprint ridge topologies is used to illustrate this approach. Indeed, as it is emphasized in \cite{Wan-Hu2007a,Mal-Mai2009,Wan-Hu2011,Jos-Dey2009,Gup-Tiw2019,Kho-Kho2017}, the estimation of fingerprint ridge topologies can be a necessary step before the use of high-level classification algorithms. In the
 review article \cite{Bia-Xu2019}, the authors proposed a classification of estimation methods into three broad categories: gradient-based methods, mathematical model-based methods, and learning-based methods. The matter discussed in the present paper falls into the second one, in particular, in the sub-class of methods that don't require prior heuristic knowledge.

The classical procedure consists in first obtaining a coarse estimation of the orientation field in the form of a discrete dataset. This  step can be achieved with gradient-based methods. Then modelling choices are made and an optimisation algorithm is applied to fit a model to the dataset. Finally, the orientation field is reconstructed with the help of the fully identified mathematical model.

On the first hand, large discrete datasets can be constructed by means of a rough method such as the computation of the gradients of the fingerprint image gray intensity changes \cite{Baz-Ger2002,Kas-Wit1987}. However, such a strategy is prone to introduce significant noise into the data \cite{Wan-Hu2007a}. Note that enhancement techniques can be applied in order to improve the situation \cite{Bia-Xu2019,Bia-Din2017,Zhu-Yin2006}.
 On the other hand, smaller but more reliable datasets can be obtained by using techniques focusing on the detection of stable and highly distinctive fingerprints features, such as minutiae \cite{Li-Kot2012,Jos-Dey2009}. In this case, the difficulty lies in the fact that the global reconstruction of the orientation field has to be performed on the basis of scarcer information.

From the modeling point of view, following \cite{She-Mon1993,God1983,Bos-Sac2016}, an orientation field can be regarded as a mapping from $\mathbb{R}^2$ to the orientation set ---that is, the interval $[0,\pi]$ where $\pi$ is equivalent to $0$--- which complicates the parametrization task. Indeed, consider the case of vector fields, which are mappings from $\mathbb{R}^2$ to $\mathbb{R}^2$. For a sufficiently regular vector field, singularities are identified by a rather simple criteria that can be computationally handled ---{\it i.e.} points where the vector field vanishes. Meanwhile, singularities of orientation fields correspond to discontinuities which, for the more classic ones, translate into points where all  orientations accumulate.
As a consequence, in order to efficiently handle orientation fields, one needs to propose a model that re-introduces some continuity property. One of the most popular solutions to this problem, which we call  \emph{doubling phase step}, sends each dataset angle from the orientation set to $\mathbb R^2$ by doubling its value and taking its \emph{sine} and \emph{cosine}, {\it e.g.} \cite{Kas-Wit1987,Wan-Hu2007a,Li-Kot2012,Wan-Hu2011,Yau-Li2004,Gup-Gup2016,Ram-Bischof2010,Got-Mih2009,Bia-Luo2014,Jir-Hou2011}. This resolves the orientation set's cyclicity, however, the discontinuity issue remains.

Starting from this basis, several models have been proposed in the literature, in order to reconstruct the doubling phase vector field and consequently, the orientation field \cite{She-Mon1994,Zhou-Gu2004a,Zhou-Gu2004b,Dass2004,Viz-Ger1996,Cap-Mai2007,Ros-Sha2007,Feng-Jain2011,Bia-Xu2019}. For instance, in the seminal paper \cite{Wan-Hu2007a}, the authors propose a method based on 2D Fourier Expansions to interpolate the doubling phase field, which helps address the discontinuity of the target without prior information. However, one can remark the two following points: first, this approach does not fully resolves the discontinuity issues, and second, the targeted field is not  the true orientation field.

The present paper proposes to address the problem of global reconstruction of orientation fields by means of bisector line fields. This model, discussed in \cite{Bos-Sac2016} in the framework of differential geometry,  appears as a natural extension of the doubling phase step.

 The \emph{bisector line field} is an orientation field constructed from  two vector fields according to the following procedure. At each point in $\mathbb{R}^2$,  the two vector fields define two directions (in $[0,2\pi]$). The mean value of these two directions, taken modulo $\pi$, belongs to the orientation set and corresponds to the orientation of the line bisecting the angle between the two vector fields.  This concept displays many properties that makes it an adequate tool for the global reconstruction task in the sense that, for sufficiently regular generating vector fields, the behavior of the bisector line field is very tractable. For instance, in the practical cases considered in this paper, singularities of bisector line fields happen whenever one of the two generating vector fields vanishes. As a consequence, this model organically solves the discontinuity issue discussed earlier.

The second contribution of this article is to rely on
an energy functional that measures the error between two orientation fields directly in the orientation set. In fact, it is very similar to the \emph{Root Mean Square Deviation} for orientation fields, which is acknowledged as well suited in order to measure the distance between the original ({\it i.e.} the ground truth) and the reconstructed orientation field \cite{Tur-Mal2011,Bia-Xu2019}.
In other words, the doubling phase step is discarded, and the reconstruction task is performed by minimising a classic comparison index.
Moreover, when this energy is associated with the bisector line field model, the resulting optimization problem is particularly suited to the gradient method since the complexity of the gradient only depends on the underlying vector field structure. Therefore, although polynomial interpolation is performed in the present paper, the proposed approach might also be used with the techniques presented in \cite{Wan-Hu2007a,Ram-Bischof2010,Tas-Hel2010,Liu-Liu2014,Gup-Gup2016,Gup-Gup2015}.

The rest of the article is organised as follows. Section~\ref{sec:bisector} is dedicated to the definition of the bisector line field model. A general presentation of the proposed approach for orientation field reconstruction is also exposed. Important mathematical properties that justifies the approach are presented in Section~\ref{sec:foundations}. It can be skipped by readers not interested in these aspects. Section~\ref{sec:energy} further details the interpolation method with the definition of the energy function and the calculation of the gradient. Finally, results are displayed and discussed in Section~\ref{sec:discussion}.

\section{Strategy overview}\label{sec:bisector}

The goal of orientation field interpolation is to transform a set of discretely localized orientation information into a global model, as illustrated in Figure~\ref{fig:ReconExpl}. Starting from a fingerprint image, a local estimation of the target orientation field is obtained. Enhancement techniques that improve quality of the dataset, {\it e.g.} \cite{Bia-Xu2019,Bia-Din2017,Zhu-Yin2006, Gup-Gup2016,Zhou-Gu2004a}, can  classically be used
at this point.

Let us stress that the amount of available data can be rather small in some cases, for instance in minutiae-based extraction \cite{Ros-Sha2007} or in low quality fingerprint impression.

\begin{definition}Let  $\mathcal{D}$ denote a subset of $\mathbb R^2$.

\begin{itemize}
    \item  $ \S^{n}=\left\{ x \in \mathbb{R}^{n+1} \mid
\sum_{k=1}^{n+1} x_k^2 = 1   \right\}$ denotes the \emph{unit n-sphere}.
In particular $\S^1$ is interpreted as the interval $[0,2\pi]$ where $0$ is identified with $2\pi$.
    \item The \emph{orientation set}, denoted by $\P^1$ is the interval $[0,\pi]$ where $0$ is identified with $\pi$.
    \item  A \emph{vector field} is a map $\mathbf{X}:\mathcal{D}\to \mathbb R^2$. We denote its \emph{zero set} by $\mathcal{Z}_\mathbf{X} = \left\{ p \in \mathcal{D}:\mathbf{X}(p) = (0,0) \right \}$.
    %
    \item  An \emph{orientation field} is a map $\L:\mathcal{D}\to \P^1$.
\end{itemize}
\end{definition}

A classical strategy, \cite{Kas-Wit1987,Wan-Hu2007a,Li-Kot2012,Wan-Hu2011,Yau-Li2004,Gup-Gup2016,Ram-Bischof2010,Got-Mih2009,Bia-Luo2014,Jir-Hou2011},  corresponding to the bottom branch  in  Figure~\ref{fig:TheStrategy},  is to not treat the estimation of the fingerprint orientation field as such.
Let $\left ( \theta_i \right)$ denote a collection of elements of $\P^1$, then a collection of elements $\left ( V_i \right)$ of $\mathbb R^2$ is obtained by:
$$ V_i = \begin{pmatrix} \cos{  2 \theta_i} \\ \sin{ 2 \theta_i  } \end{pmatrix} .$$
Next, data are fitted. This requires choosing both a mathematical model for the global vector field and a fitting technique. Finally, the interpolated data is taken back into the orientation set.

As emphasized in Figure~\ref{fig:TheStrategy}, we propose a more straightforward strategy based on a model called the {\it bisector line field}.

\begin{figure}
\begin{center}
\subfloat[Fingerprint]{\includegraphics[width = 0.3 \linewidth]{Demo_1.pdf} \label{fig:ReconExpl_Fingerprint}} \quad
\subfloat[Estimated OF]{\includegraphics[width = 0.3 \linewidth]{Demo_3.pdf}} \quad
\subfloat[Reconstructed OF]{\includegraphics[width = 0.3 \linewidth]{Demo_2.pdf}}
\caption{Reconstruction example. In (b) the grey level scale is such that black pixels correspond to $0$ and white ones to $\pi$.}
\label{fig:ReconExpl}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.65\columnwidth]{TheStrategy2}
\caption{FlowGraph of the proposed approach}\label{fig:TheStrategy}
\end{center}
\end{figure}

Similarly to vector fields, orientation fields admit an intrinsic definition from the differential topology point of view, see \cite{Bos-Sac2016} for a thorough discussion. However, such a definition lacks efficiency if one hopes to work with orientation fields and apply techniques from differential calculus. The introduction of bisector line fields is meant to bridge this gap.
\begin{definition}
Let  $\mathcal{D}$ denote a subset of $\mathbb R^2$.
Let $\mathbf{X}$ and $\mathbf{Y}$ be two vector fields on $\mathcal{D}$.
The \emph{bisector line field of $(\mathbf{X},\mathbf{Y})$} is the orientation field
$$
\mathrm{B}(\mathbf{X},\mathbf{Y}):\mathcal{D}\setminus(\mathcal{Z}_\mathbf{X}\cup \mathcal{Z}_\mathbf{Y})\to \P^1
$$
 such that at any point $p\in\mathcal{D}$, $\mathrm{B}(\mathbf{X},\mathbf{Y})(p)$ is the orientation of the line bisecting the oriented pair of vectors $(\mathbf{X}(p),\mathbf{Y}(p))$.
    %
\end{definition}

In practice, this definition implies the following.
Let $p\in \mathcal{D}\setminus(\mathcal{Z}_\mathbf{X}\cup \mathcal{Z}_\mathbf{Y})$, there exist unique $\theta_\mathbf{X}(p)$ and $\theta_\mathbf{Y}(p)$ in $\S^1$ such that
$$
\mathbf{X}(p)=\|\mathbf{X}(p)\|
\begin{pmatrix}
\cos\theta_\mathbf{X}(p)
\\
\sin\theta_\mathbf{X}(p)
\end{pmatrix}$$
and
$$
\mathbf{Y}(p)=\|\mathbf{Y}(p)\|
\begin{pmatrix}
\cos\theta_\mathbf{Y}(p)
\\
\sin\theta_\mathbf{Y}(p)
\end{pmatrix}.
$$
Then, as illustrated in  Figure~\ref{F:bissection},
$$
\mathrm{B}(\mathbf{X},\mathbf{Y})(p)=\frac{1}{2} \left(\theta_\mathbf{X}(p)+\theta_\mathbf{Y}(p)\right)\pmod{\pi}.
$$

\begin{figure}[htb]
\begin{center}
\begin{tikzpicture}[line join=round,>= latex,x=2.0cm,y=2.0cm]

\clip(1.4,1.8) rectangle (5.,3.8);

\draw [->,line width=1.pt,color=black] (2.8326419725368748,2.2164129218750697) -- (4.685595354668596,3.420621357556255);

\draw [->,line width=1.pt,color=black] (2.8326419725368748,2.2164129218750697) -- (1.498242498948219,3.359383032014635);

\draw [line width=1.pt,domain=1.4:5.,color=black] plot(\x,{(-2.6803154837669454--0.9978232810591935*\x)/0.0659446720839969});

\draw [shift={(2.8326419725368748,2.2164129218750697)},->,line width=1.pt,color=black] (0.5553920758325129:0.6) arc (0.5553920758325129:86.21890474934047:0.6);

\draw [->,line width=1.pt,color=black] (2.8326419725368748,2.2164129218750697) -- (3.832594991705581,2.2261061904357087);

\draw [shift={(2.8326419725368748,2.2164129218750697)},->,line width=1.pt,color=black] (0.5553920758325129:0.8) arc (0.5553920758325129:33.019275032959115:0.8);

\draw [shift={(2.8326419725368748,2.2164129218750697)},->,line width=1.pt,color=black] (0.5553920758325129:0.4) arc
(0.5553920758325129:139.4185344657218:0.4) ;

\draw[color=black] (4.7,3.2) node {$\mathbf{X}$};
\draw[color=black] (1.65,3.45) node {$\mathbf{Y}$};
\draw[color=black] (3.8,2.05) node {$\mathbf{Z}$};
\draw[color=black] (3.3,3.6) node {$\mathrm{B}(\mathbf{X},\mathbf{Y})$};
\draw[color=black] (3.25,2.9) node {$\frac{\theta_\mathbf{X}+\theta_\mathbf{Y}}{2}$};
\draw[color=black] (3.8,2.5) node {$\theta_\mathbf{X}$};
\draw[color=black] (2.65,2.75) node {$\theta_\mathbf{Y}$};

\end{tikzpicture}
\caption{ Representation of the bisection operation on an oriented pair of vectors, where $\mathbf{Z}=(1,0)'$ marks the reference vector to measure angles.}
\label{F:bissection}
\end{center}
\end{figure}
Our strategy consists in building $\mathbf{X}$ and $\mathbf{Y}$ such that $B\left(\mathbf{X},\mathbf{Y}  \right)$ interpolates the data. The data fitting step then relies on the above formula. In this regard, the modelling of the vector fields $\mathbf{X}$ and $\mathbf{Y}$ is completely open. Numerous strategies have been identified for this purpose, be it  rational complex functions \cite{Zhou-Gu2004b}, Taylor expansions \cite{For-Str1993,Yau-Li2004}, orthogonal/Legendre/Chebychev polynomials \cite{Tas-Hel2010,Liu-Liu2014,Ram-Bischof2010,Bia-Luo2014}, Fourier expansions \cite{Wan-Hu2007a,Tas-Hel2009,Tao-Yan2010}, and others \cite{Bia-Xu2019}. For the sake of simplicity, we chose to implement a Taylor expansion model.

\begin{remark}
As described in the introduction, the concept of orientation field arises in different fields of applications. However, as it is a rather uncommon notion,  terminology can vary between fields of study. In the context of differential geometry, the term line field is sometimes preferred, following Hopf's terminology ``fields of line elements'' \cite{Hopf,Bos-Sac2016}. We chose the terminology that is more popular in the fingerprint analysis  community.
\end{remark}

\section{Bisector line fields}\label{sec:foundations}

{
In this section, we expose properties of  bisector line fields  that highlight the suitability of this structure for the  interpolation of orientation fields.
}
\subsection{Singularities of line fields}

In our application case, singularities of orientation fields should be understood as unresolvable discontinuities of the map $\L:\mathbb R^2\to \P^1$. To this extent, we refer in the following to continuous orientation fields $\L:\mathcal{D}\to \P^1$, where $\mathcal{D}$ is the continuity domain of $\L$, and singularities of $\L$ are located in $\bar{\mathcal{D}}\setminus \mathcal{D}$, the border of $\mathcal{D}$, where $\bar{\mathcal{D}}$ denotes the topological closure of $\mathcal{D}$.

\begin{definition}
A continuous orientation field $\L:\mathcal{D}\to \P^1$ is said to be \emph{singular} at a point $p\in \bar{\mathcal{D}}\setminus \mathcal{D}$ if it cannot be uniquely continuously extended at $p$.
\end{definition}

In the following, we focus on isolated singularities, which are sufficient to our framework.
 In particular, we can introduce the topological index, which is a useful tool for the study of such singularities of line fields.

Let $\L:\mathcal{D}\to \P^1$ be a continuous orientation field and let $p$ be an isolated singularity of $\L$ such that $p$ belongs to the interior of $\bar{\mathcal{D}}$ (for instance there exists an open subset $V$ of $\mathbb R^2$   such that $V=\{p\}\cup(V\cap\mathcal{D})$). Then the index of $V$ at $p$, denoted $\mathrm{ind}_p(\L)$, is a half integer quantifying the winding of $\L$ around $p$  (for a precise construction of this object see, for instance, \cite[Section 3.2]{Bos-Sac2016}).

In the context of the study of orientation fields singularities, the following properties of bisector line fields arise.

\begin{proposition}[\cite{Bos-Sac2016}]\label{P:continuity}
 If $\mathbf{X}$ and $\mathbf{Y}$ are continuous vector fields over $\mathbb R^2$ then $\mathrm{B}(\mathbf{X},\mathbf{Y})$ is continuous on $\mathcal{D}=\mathbb R^2\setminus(\mathcal{Z}_X\cup \mathcal{Z}_\mathbf{Y})$.

On the other hand, let $\L:\mathcal{D}\to \P^1$ be a $C^k$ orientation field, for some $k\in \mathbb N\cup\{\infty\}$. Then there exist two $C^k$ vector fields $\mathbf{X},\mathbf{Y}$ over $\mathbb R^2$ such that $\L=\mathrm{B}(\mathbf{X},\mathbf{Y})$.
\end{proposition}

\begin{proposition}
Let $(\mathbf{X}, \mathbf{Y} )$ be a pair of continuous vector fields over $\mathbb R^2$. Given an isolated point $p$ of $\mathcal{Z}_\mathbf{X} \cup \mathcal{Z}_\mathbf{Y}$, we have
$$
\mathrm{ind}_p(\mathrm{B}(\mathbf{X}, \mathbf{Y} )) =  \frac{1}{2}(\mathrm{ind}_p(\mathbf{X}) + \mathrm{ind}_p(\mathbf{Y} )).
$$
\end{proposition}

\begin{remark}
An immediate consequence of this model is that singularities, initially characterized by an analytical property, now coincide with the zero sets of regular functions.
In particular, isolated singularities of the bisector line field corresponds to  isolated singularities of either one of the two generating vector fields.
\end{remark}

\subsection{Generic properties of bisector line fields}

Following Thom \cite{thom1972stabilite}, it is understood, philosophi\-cally speaking, that ``typical'' behaviors of mathematical objects should be the only ones visible in nature. From the point of view of transversality theory, these typical properties are known as \emph{generic}, in the sense that a property is generic on a topological set if it is satisfied on a residual subset.
A residual subset is understood to be a large dense set
in the following sense:
it is a countable intersection of open and dense subsets.

The classical application case of this theory is geared towards regular maps, endowed with the $C^k$ Whitney topology, for some $k\in \mathbb N$ (see, {\it e.g.} \cite{Hir1976,Abr-Rob1967}).

When generic features of vector fields are considered, more can be said \cite{Bos-Sac2016}. Example of such features are isolation of singularities of vector fields or the fact that indices of such singularities must be  $\pm 1$.

Let $\mathbf{X}$ and $\mathbf{Y}$ be two smooth vector fields. Generically with respect to the $C^2$ Whitney topology, $\mathcal{Z}_\mathbf{X}$ and $\mathcal{Z}_\mathbf{Y}$ are discrete collections of points that do not accumulate.
Thus, singularities of $\mathbf{X}$ and $\mathbf{Y}$ are isolated, and moreover have $\pm 1$ index.
Furthermore, when considering $(\mathbf{X},\mathbf{Y})$ as a pair, we get  $\mathcal{Z}_\mathbf{X}\cap \mathcal{Z}_\mathbf{Y}=\emptyset$.

As a consequence of these facts and Proposition~\ref{P:continuity}, we get the following property for generic smooth bisector line fields, as illustrated in Figure~\ref{F:Bissection_example}.

\begin{proposition}
Let $(\mathbf{X},\mathbf{Y})$ be a generic pair of vector fields.
Singularities of $\mathrm{B}(\mathbf{X},\mathbf{Y})$ are isolated and have index $\pm1/2$.
\end{proposition}

\begin{remark}
The set of orientation fields is not endowed with a topology allowing study of genericity or stability. It is then necessary to provide an extrinsic structure and this was the motivation for the study of bisector line fields in \cite{Bos-Sac2016}. Another popular model relies on differential 2-forms on surfaces \cite{Sot-Gar2008}.
\end{remark}

{
\begin{figure}[hbtp]
\begin{center}
\subfloat[Vector field $\mathbf{X}$]{
\includegraphics[width=.3\linewidth,trim=1mm 1mm 1mm 1mm, clip=true]{ChampX.pdf}}
\quad
\subfloat[Vector field $\mathbf{Y}$]{
\includegraphics[width=.3\linewidth,trim=1mm 1mm 1mm 1mm, clip=true]{ChampY.pdf}}
\quad
\subfloat[Orientation field $\mathrm{B}(\mathbf{X},\mathbf{Y})$]{
\includegraphics[width=.3\linewidth,trim=1mm 1mm 1mm 1mm, clip=true]{ChampBXY.pdf}}
\end{center}
\caption{
Example of the bisection of two vector fields.
Singularities of positive index have been marked with a { $\oplus$} symbol, while singularities of negative index have been marked with a $\ominus$ symbol.
}
\label{F:Bissection_example}
\end{figure}

}

\section{Interpolation as an energy minimization problem}\label{sec:energy}

This section presents the elements needed the implementation of the orientation field reconstruction task. That is, the energy function, the polynomial model used for vector fields, and the gradient descent scheme.

\subsection{Discrete energy}\label{subsec:DiscEner}

The \emph{Root Mean Square Deviation} for orientation fields is acknowledged as well suited to measure the distance between two orientation fields \cite{Tur-Mal2011,Bia-Xu2019}. As seen in Section~\ref{sec:bisector}, the bisector line field model allows to work in the orientation set. Therefore, the root mean square deviation can now be directly used to perform the reconstruction task. This section details this quantity in the energy formalism, which is suitable for gradient descent schemes. Let us underline that, for a given number of data, minimizing the energy is equivalent to minimizing the root mean square deviation.

Let $d:[0,\pi]\times [0,\pi]\to [-\pi/2,\pi/2] $ be piecewise defined by
\begin{equation}
 d(\theta_1,\theta_2)=\left\{
\begin{aligned}
 &\theta_1-\theta_2 & \text{ if }& \left|\theta_1-\theta_2\right| \leq \frac{\pi }{2}, \\
&\theta_1-\theta_2-\pi  & \text{ if }&\theta_1-\theta_2>\frac{\pi }{2} ,
\\
&\theta_1-\theta_2+\pi  & \text{ if }&\theta_1-\theta_2<-\frac{\pi }{2} .
\end{aligned}
 \right.
 \label{equ:distance}
\end{equation}
 Let $\mathcal{I} = \left \{1,\hdots,m \right \} $ where $m$ is the number of data, and let $(x_i,y_i)_{i\in \mathcal{I}}$ be a  collection of points in $\mathbb R^2$.
For $\mathbf{M},\L:\mathcal{D}\subset \mathbb R^2\to \P^1$, such that $(x_i,y_i)_{i\in \mathcal{I}}\subset \mathcal{D}$,
we set  the least-squares energy functional to be
$$
\mathrm{J}(\mathbf{M},\L)=\sum_{i\in \mathcal{I}} d\left(\mathbf{M}(x_i,y_i),\L(x_i,y_i)\right)^2.
$$

The relevance of $d$ is further justified by the following observations.

\begin{remark}\label{R:torus}

\begin{enumerate}
\item
Consider the energy given by $d^2$ on the torus $\P^1\times \P^1$. The gradient of $d^2$ is
$$
\nabla d^2(\theta_1,\theta_2)= 2 d(\theta_1,\theta_2)(\partial_{\theta_1}-\partial_{\theta_2})
$$
As a consequence, one can check that the gradient flow of $d^2$ is actually parallel to the geodesic flow on the flat torus starting from the affine sets $\{\theta_1-\theta_2=\pm \pi/2\}$ to the diagonal $\{\theta_1=\theta_2\}$ as it is illustrated in
 Figure~\ref{F:torus}.
 \item
 %
 The two connected components of $\P^1\times \P^1\setminus\{\theta_1 =\theta_2 \pmod{\pi/2}\}$ are geodesically convex once $\P^1\times \P^1$ has been endowed with the flat torus Riemannian metric. Furthermore,
 the maps $d_{|C}$ and $d^2_{|C}$ are then geodesically convex  and strictly geodesically convex respectively.
\item Notice that $\mathrm{J}^{1/2}$
is a pseudo-metric (it is symmetric and satisfies the triangular inequality) on the space of $\P^1$-valued maps over a domain containing $(x_i,y_i)_{i\in \mathcal{I}}$.

Indeed $|d|$ is a distance over $\P^1$ and we {\it de facto} have the classical product metric over $\left(\P^1\right)^{m}$ given by
$$
d_m((\theta_1,\dots,\theta_m),(\theta_1',\dots,\theta_m'))=
\left(
\sum_{i=1}^m
d(\theta_i,\theta_i')^2
\right)^{1/2}
$$
Then with $\mathcal{I} = \left \{1,...,m\right\}$,
$$
(\theta_1,\dots,\theta_m)=(\mathbf{M}(x_1,y_1),\dots,\mathbf{M}(x_{m},y_{m}))
$$ and
$$
(\theta_1',\dots,\theta_m')=(\L(x_1,y_1),\dots,\L(x_{m},y_{m})),
$$
one has
$$
\mathrm{J}(\mathbf{M},\L)^{1/2}
=
d_m((\theta_1,\dots,\theta_m),(\theta_1',\dots,\theta_m')).
$$
\end{enumerate}
\end{remark}

\begin{figure}[hbtp]
\centering
\hfill
\subfloat[Flat representation]{
\includegraphics[width=.4\linewidth,trim=1mm 1mm 1mm 1mm, clip=true]{gradient_d2_1.pdf}}
\hfill
\subfloat[Embedded 3d representation]{
\includegraphics[width=.4\linewidth,trim=1mm 1mm 1mm 1mm, clip=true]{gradient_d2_2.pdf}}
\hfill{}
\caption{
Representation of the gradient flow of $d^2$ on the flat torus $\P^1\times \P^1$. The set $\{\theta_1-\theta_2=\pm \pi/2\}$ is denoted with a dashed line, the set  $\{\theta_1=\theta_2\}$ is denoted with a continuous line.}
\label{F:torus}
\end{figure}

\subsection{Polynomial bisector line fields} \label{subsec:PolyModel}

For a given target $\mathbf{M}$, we minimize $\mathrm{J}(\mathbf{M},\L)$ over the set of polynomial bisector line fields of a fixed maximal degree.

Notice that for a given vector field $\mathbf{X}$, $\lambda \mathbf{X}$ yields the same direction $\theta_\mathbf{X}$ for all $\lambda>0$. Hence we can compactify the set of polynomial vector fields under consideration.
To this extent, let us define for all $n\in \mathbb N$ the set $\mathfrak{X}^n$ of polynomial vector fields of degree $n$ such that $\mathbf{X}\in \mathfrak{X}^n$ if there exist $(\alpha_{k,j})_{0\leq j\leq k\leq n}$, $(\beta_{k,j})_{0\leq j\leq k\leq n}$
such that
$$
\mathbf{X}(x,y)=\sum_{k=0}^{n}\sum_{j=0}^k
\begin{pmatrix}
\alpha_{k,j}
\\
\beta_{k,j}
\end{pmatrix}
x^{k-j}y^j,
$$
and
$$
\sum_{k=0}^{n}\sum_{j=0}^k
\alpha_{k,j}^2
+
\beta_{k,j}^2
=1.
$$
In particular $\left(\alpha_{k,j},\beta_{k,j}\right)_{0\leq j\leq k}\in \S^{(n+1)(n+2)}$ and the space $\mathfrak{X}^n$ is compact for any given $n\in \mathbb N$.

The polynomial interpolation of the target $\mathbf{M}$ is performed by finding $\mathbf{X}$ and $\mathbf{Y}$ solving, for $m_\mathbf{X}, m_\mathbf{Y}\in \mathbb N$,
\begin{equation}\label{E:min_problem}
    \min
    \left\{
    \mathrm{J}(\mathbf{M},\mathrm{B}(\mathbf{X},\mathbf{Y}))
    \mid
    \mathbf{X}\in \mathfrak{X}^{m_\mathbf{X}}, \mathbf{Y}\in \mathfrak{X}^{m_\mathbf{Y}}
    \right\}.
\end{equation}

Let us denote
\begin{enumerate}
 \item $n_\mathbf{X} =  (m_\mathbf{X}+1)(m_\mathbf{X}+2)$ and
$ n_\mathbf{Y} =  (m_\mathbf{Y}+1)(m_\mathbf{Y}+2)$;
    \item
     $ \omega \in \S^{n_\mathbf{X}}\times \S^{n_\mathbf{Y}}$ such that
      $$\omega=\left(\left(\alpha_{k,j},\beta_{k,j}\right)_{0\leq j\leq k\leq m_\mathbf{X}},\left(\gamma_{k,j},\delta_{k,j}\right)_{0\leq j\leq k\leq m_\mathbf{Y}}\right);$$
    %
    \item $\mathbf{X}_\omega(x,y)= \displaystyle\sum_{k=0}^{m_\mathbf{X}}\displaystyle\sum_{j=0}^k
\begin{pmatrix}
\alpha_{k,j}
\\
\beta_{k,j}
\end{pmatrix}
x^{k-j}y^j$

 and

 $
\mathbf{Y}_\omega(x,y)= \displaystyle \sum_{k=0}^{m_\mathbf{Y}} \displaystyle\sum_{j=0}^k
\begin{pmatrix}
\gamma_{k,j}
\\
\delta_{k,j}
\end{pmatrix}
x^{k-j}y^j;
$
\item $\L_\omega=\mathrm{B}(\mathbf{X}_\omega,\mathbf{Y}_\omega)$.
\end{enumerate}

Problem~\eqref{E:min_problem} is then equivalent to solving
\begin{equation}
\label{E:min_problem_compact}
\min
\left\{
\mathrm{J}(\mathbf{M},\L_\omega)\mid \omega \in  \S^{n_\mathbf{X}}\times \S^{n_\mathbf{Y}}
\right\}.
\end{equation}

The existence of minimizers is guaranteed by the compacity of $\S^{n_\mathbf{X}}\times \S^{n_\mathbf{Y}}$.

\subsection{Projected gradient descent approach}  \label{subsec:GradientDesc}
The function $d$, given in Equation~\eqref{equ:distance}, and Problem~\eqref{E:min_problem_compact} have been designed to allow the use of gradient based optimisation algorithms. Indeed, since
$$
\L_{\omega}(x,y)=\frac{1}{2}\left(\theta_{\mathbf{X}_\omega}(x,y)+\theta_{\mathbf{Y}_\omega}(x,y)\right),
$$
we have on $\mathbb R^2\setminus \mathcal{Z}_{\mathbf{X}_\omega}\cup \mathcal{Z}_{\mathbf{Y}_\omega}$
\begin{align}
\nonumber
\nabla_\omega \L_\omega
= &
\frac{1}{2}
\left(\frac{\mathbf{X}_{\omega,1} \nabla_\omega \mathbf{X}_{\omega,2}-\mathbf{X}_{\omega,2} \nabla_\omega \mathbf{X}_{\omega,1}}{\mathbf{X}_{\omega,1}^2+\mathbf{X}_{\omega,2}^2} \right.
 \\ &
 %
 \hspace{30pt}
 %
  \left. +
\frac{\mathbf{Y}_{\omega,1} \nabla_\omega \mathbf{Y}_{\omega,2}-\mathbf{Y}_{\omega,2} \nabla_\omega \mathbf{Y}_{\omega,1}}{\mathbf{Y}_{\omega,1}^2+\mathbf{Y}_{\omega,2}^2}
\right),
\label{E:derivative_B}
\end{align}
where $(\cdot)_{\omega,i}$ denotes the $i^{th}$ component of $(\cdot)_{\omega}$.

As a consequence, we can immediately compute by the  chain rule
\begin{multline*}
    \nabla_\omega \mathrm{J}(\mathbf{M},\L_\omega)=\\-2\sum_{i\in \mathcal{I}} d\left(\mathbf{M}(x_i,y_i),\L_\omega(x_i,y_i)\right)
\nabla_\omega \L_\omega(x_i,y_i).
\end{multline*}
For $\rho>0$, it  should be noted  that in general:
$$ \omega+\rho \nabla_\omega \mathrm{J}(\mathbf{M},\L_\omega)\notin \S^{n_\mathbf{X}}\times \S^{n_\mathbf{Y}}.$$ However, we solve this issue with the  projection
$$
\begin{array}{cccc}
    p :& \mathbb R^{n_\mathbf{X}}\times \mathbb R^{n_\mathbf{Y}}&\longrightarrow&\S^{n_\mathbf{X}}\times \S^{n_\mathbf{Y}} \\
    &(\omega_1,\omega_2) &\longmapsto&
    \left(
    \dfrac{\omega_1}{\|\omega_1\|}
    ,
    \dfrac{\omega_2}{\|\omega_2\|}
    \right).
\end{array}
$$

Thus in the gradient descent scheme, we iterate the recursive transform
$$
\omega_{k+1}\leftarrow p
\left(
\omega_k+\rho \nabla_\omega \mathrm{J}\left(\mathbf{M},\L_{\omega_k}\right)
\right).
$$

\section{Examples and discussion}\label{sec:discussion}

\subsection{Experimental setting}

In order to illustrate the paper's approach,
simulations have been performed through an implementation of a constant step-size gradient descent method in Matlab.
In line with the main application of this theory ---fingerprint analysis---  we tested the method on elements of the  FVConGoing Initiative data set  \cite{Dor-Cap2009}.

More precisely, a target has been obtained with the classical elementary method of lifting the orientations of the finger ridges from the gradient of a grey-scale image of a fingerprint \cite{Kas-Wit1987,Mal-Mai2009}. In order to discard parts of the image that don't correspond to the fingerprint, we neglected orientations
where the gradient's norm was too low. As a consequence,
a few data corresponding to the actual fingerprint might be missing without significant effect on the result.

In the following, we present two reconstruction examples obtained with datasets of orientations of large size but low quality ---{\it i.e.} $245 \times145-$pixels images leading to approximately 30000 data-points. One has been performed on a classical loop fingerprint, the other on a classical whorl fingerprint. A third reconstruction example has been obtained with a scarcer dataset of high fidelity, similar to the minutiae based methods found in \cite{Feng-Jain2011,Ros-Sha2007}. This reconstruction, based on a smaller dataset made of $40$ elements,  also serves to illustrate situations when latent fingerprints are treated, or the image corrupted. The scope of the present paper is not to propose a method to extract datasets in these degraded situations. However, as it is explained in Section~\ref{sec:bisector} and schematized in Figure~\ref{fig:TheStrategy}, any method that efficiently extracts data can be used prior to the bisector line field based reconstruction.

Additional reconstruction experiments with datasets made of $10$, $20$, $40$ and $80$ orientations values have been performed. In each of these four cases, ten datasets where randomly generated. The results are compared to a reference set by computing the root mean square deviation \cite{Tur-Mal2011,Gup-Gup2016,Bia-Xu2019}, defined as follows.

Let $\mathcal{J}\subset \mathbb N$ and  $(\bar x_j,\bar y_j)_{j \in \mathcal{J}}$ be the collection of points ranging all the relevant pixel positions of the fingerprint image. Let  $\overline{\mathbf{M}}=\left(\mathbf{M}(\bar x_j,\bar y_j)\right)_{j \in \mathcal{J}}$ and $\overline{\L}=\left(\L(\bar x_j,\bar y_j)\right)_{j \in \mathcal{J}}$ be two matrices of orientations. Then, the root mean square deviation (RMSD) is:
 $$ RMSD (\overline{\mathbf{M}}, \overline{\L}) = \displaystyle \sqrt{  \frac{\sum_{j\in \mathcal{J}} d \left({\mathbf{M}}(\bar x_j,\bar y_j),{\L}(\bar x_j,\bar y_j)\right)^2}{\vert \mathcal{J} \vert   }  }  $$
The considered reference set is the orientation field that was reconstructed in the first round of experiments. This procedure is presented for the loop fingerprint only. The measured RMSD in these experiments serve as a basis for comparison with the third reconstruction example, where the data-points were picked to be meaningful.

\subsection{Results}\label{subsec:results}
The result of the interpolation experiments with large datasets are presented in Figures~\ref{fig:Whorl-FatSet}~and~\ref{fig:Loop-FatSet}.
The first image is the input of the algorithm, a grey-scale image of a fingerprint. The second is a representation of  the reconstructed orientation field, as a field of lines, that has been superimposed on the input.

The third and fourth images correspond to a representation of the phases of the two orientation fields, that is, grey-scale images where angles from $0$ to $\pi$ are mapped to a light intensity (near $0$, dark, near $\pi$, bright). Image (c) is then a representation of the  matrix of targeted orientations $\overline{\mathbf{M}}=\left(\mathbf{M}(\bar x_j,\bar y_j)\right)_{j \in \mathcal{J}}$, while image (d) is the matrix $\overline{\L}=\left(\L(\bar x_j,\bar y_j)\right)_{j \in \mathcal{J}}$ of the interpolated field.
Notice that the hard lines separating black and white correspond to the location of transitions from $0$ to $\pi$ in the orientation. Likewise, singularities are points where all grey levels accumulate (similar to the pinwheel singularities observed in the visual cortex \cite{pet2017}).

In the case of a scarce dataset, the input is a small collection of triples $(x,y,\theta)\in \mathbb R^2\times \P^1$, hence, the results are presented in  Figure~\ref{fig:Loop-SparseSet} in a slightly different way. Indeed, the representation of the matrix $\overline{\mathbf{M}}$ is no longer appropriate: in the third experiment, the data-set $(x_i,y_i)_{i\in \mathcal{I}}$ does not coincide with $(\bar x_j,\bar y_j)_{j \in \mathcal{J}}$. Therefore, we added a representation of the $40$ inputs as line elements on the plane, cf.~Figure~\ref{Fig:ScarceInput}.

\begin{figure*}
\centering
\subfloat[Fingerprint]{\frame{\includegraphics[width = 0.22\textwidth,trim=6.5cm 9.1cm 5.8cm 8.6cm,clip]{Whorl_1.pdf}}}
\quad
\subfloat[Output]{\frame{\includegraphics[width = 0.22\textwidth,trim=6.5cm 9.1cm 5.8cm 8.6cm,clip]{Whorl_2.pdf}}}
\quad
\subfloat[Matrix $\overline{\mathbf{M}}$]{\frame{\includegraphics[width = 0.22\textwidth,trim=6.5cm 9.1cm 5.8cm 8.6cm,clip]{Whorl_3.pdf}}}
\quad
\subfloat[Matrix $\overline{\L}$]{\frame{\includegraphics[width = 0.22\textwidth,trim=6.5cm 9.1cm 5.8cm 8.6cm,clip]{Whorl_4.pdf}}}
\quad
\caption{``Whorl'' Experiment ---see details in Section~\ref{subsec:results}}
\label{fig:Whorl-FatSet}
\end{figure*}

\begin{figure*}
\centering
\subfloat[Fingerprint]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{Loop_1}}}
\quad
\subfloat[Output]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{Loop_2}}}
\quad
\subfloat[Matrix $\overline{\mathbf{M}}$]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{Loop_3}}}
\quad
\subfloat[Matrix $\overline{\L}$]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{Loop_4}}}
\quad
\caption{``Loop'' Experiment  ---see details in Section~\ref{subsec:results}}
\label{fig:Loop-FatSet}
\end{figure*}

\begin{figure*}
\centering
\subfloat[Fingerprint]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{SPARSE_Loop_1.pdf}}}
\quad
\subfloat[Input data]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{SPARSE_Loop_2.pdf}\label{Fig:ScarceInput}}}
\quad
\subfloat[Output]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{SPARSE_Loop_3.pdf}}}
\quad
\subfloat[Matrix $\overline{\L}$]{\frame{\includegraphics[width = 0.22\textwidth,trim=7.4cm 9.1cm 6.6cm 8.6cm,clip]{SPARSE_Loop_4.pdf}}}
\quad
\caption{``Loop'' Experiment with a dataset made of 40 selected elements  ---see details in Section~\ref{subsec:results}}
\label{fig:Loop-SparseSet}
\end{figure*}

Let us briefly comment on these results. As expected, singularities are elegantly fitted without prior knowledge of their positions { (see, for instance, the zoomed in details of Figures~\ref{sfig:whorl_detail1}-\ref{sfig:whorl_detail2} and Figures~\ref{sfig:loop_detail1}-\ref{sfig:loop_detail2})}. Interpolations of degree $(3,3)$ seem to be sufficient to obtain these results and higher order does not improve the final energy significantly. This is in line with other similar interpolation methods (see, for instance, the discussion on the matter in \cite{Wan-Hu2007a}).
It also appears that inflexions in the orientation fields seem to be hard to fit properly
{ (such as some areas of the Loop experiment, Figure~\ref{fig:Loop-FatSet}, in particular the highlighted region in Figures~\ref{sfig:inflexion_detail1}-\ref{sfig:inflexion_detail2}).}
This is a weakness of the method that seems to be a general difficulty  observed across the different interpolation methods we  encountered in the literature. However, this does not seem to be an obstacle for the later use of finer analysis and reconstruction methods.

In the case of a scarcer dataset, a moderate number of  meaningful points is sufficient to recover some global information on the structure of the orientation field, such as index and curvature of certain regions. However, precise placement of the singularities is dependent on the position of elements chosen as inputs.
Indeed, in the  example case shown in Figure~\ref{fig:Loop-SparseSet},
a majority of data points tightly fit the singular regions which allows a good reconstruction of the underlying orientation field.

In contrast, we ran a series of experiments with data-points positions randomly picked in the domain and corresponding orientations taken from the results of the first two experiments. As such, these are considered to be high-fidelity data in the sense that the first reconstruction filtered out the noise. Four series of experiments have been performed with $10$, $20$, $40$ and $80$ data-points. Each batch was made of $10$ datasets. The average RMSD and the standard deviation for each category are reported in Table~\ref{tab:SCARCE_RMSD}. { Representative examples of the batches are shown in Figure~\ref{fig:Loop-SparseSet-Random}.
With only 10 data-points, relevant information is easily missed, and singularities are not reconstructed, see Figure~\ref{fig:SCARCE10_C}.  The singularities appear with more datapoints, but their distribution influences the quality of the placement, see Figure~\ref{fig:SCARCE20_C}-\ref{fig:SCARCE40_C}-\ref{fig:SCARCE80_C}.
}

As one can expect, the mean RMSD decreases with an increasing number of data-points, and the standard deviation narrows down. However, the RMSD of the third experiment is 0.07278 (experiment corresponding to Figure~\ref{fig:Loop-SparseSet}).
This value, obtained with 40 hand-picked data-points, turns out to be more in line with the values observed with 80 points, highlighting the influence of the singularities on the global pattern.

\begin{figure}
    \centering
    \begin{minipage}{.3\linewidth}
    \subfloat[Whorl details\label{sfig:whorl_detail1}]{
        \frame{\includegraphics[width = \linewidth,trim=8cm 12.5cm 7cm 11.2cm,clip]{Whorl_1.pdf}}}
    \end{minipage}
    \qquad
    \begin{minipage}{.3\linewidth}
    \subfloat[Interpolation\label{sfig:whorl_detail2}]{
        \frame{\includegraphics[width = \linewidth,trim=8cm 12.5cm 7cm 11.2cm,clip]{Whorl_2.pdf}}}
    \end{minipage}

    \begin{minipage}{.3\linewidth}
    \subfloat[Loop details\label{sfig:loop_detail1}]{
        \frame{\includegraphics[width = \linewidth,trim=8.3cm 10.7cm 8.2cm 14.5cm,clip]{Loop_1.pdf}}}
    \end{minipage}
    \qquad
    \begin{minipage}{.3\linewidth}
    \subfloat[Interpolation\label{sfig:loop_detail2}]{
        \frame{\includegraphics[width = \linewidth,trim=8.3cm 10.7cm 8.2cm 14.5cm,clip]{Loop_2.pdf}}}
    \end{minipage}

    \begin{minipage}[t]{.3\linewidth}
    \subfloat[Region of inflexion\label{sfig:inflexion_detail1}]{
        \frame{\includegraphics[width = \linewidth,trim=8cm 12.cm 9.5cm 14.2cm,clip]{Loop_1.pdf}}}
    \end{minipage}
    \qquad
    \begin{minipage}[t]{.3\linewidth}
    \subfloat[Imprecise interpolation\label{sfig:inflexion_detail2}]{
        \frame{\includegraphics[width = \linewidth,trim=8cm 12cm 9.5cm 14.2cm,clip]{Loop_2.pdf}}}
    \end{minipage}

    \caption{Details of interpolation figures.}
    \label{fig:details}

\end{figure}

\begin{table}[tb]
    \centering
    \begin{tabular}{| c|c | c|}
    \hline
           &  RMSD mean  & RMSD standard deviation  \\
           \hline
           \hline
        $\left \vert \mathcal{I}  \right \vert = 10 $ &  $3.809  \times 10^{-1}$ & $1.121 \times 10^{-1}$  \\
 $\left \vert \mathcal{I}  \right \vert = 20 $ &  $2.637 \times 10^{-1} $ & $8.150 \times 10^{-2}$ \\
  $\left \vert \mathcal{I}  \right \vert = 40 $ & $1.519 \times 10^{-1}$ &  $5.978 \times 10^{-2}$ \\
   $\left \vert \mathcal{I}  \right \vert = 80 $ & $ 7.873\times 10^{-2}$  & $2.346 \times 10^{-2}$
   \\ \hline
    \end{tabular}
    \caption{Distribution of final RMSD.}
    \label{tab:SCARCE_RMSD}
\end{table}

\begin{figure*}
\begin{center}

\subfloat[10 data points]{\frame{\includegraphics[width = 0.22\textwidth, trim=7.4cm 9.1cm 6.5cm 8.5cm,clip]{SCARCE10_B_Wsquare.pdf}}}
\quad
\subfloat[20 data points]{\frame{\includegraphics[width = 0.22\textwidth, trim=7.4cm 9.1cm 6.5cm 8.5cm,clip]{SCARCE20_B_Wsquare.pdf}}}
\quad
\subfloat[40 data points]{\frame{\includegraphics[width = 0.22\textwidth, trim=7.4cm 9.1cm 6.5cm 8.5cm,clip]{SCARCE40_B_Wsquare.pdf}}}
\quad
\subfloat[80 data points]{\frame{\includegraphics[width = 0.22\textwidth, trim=7.4cm 9.1cm 6.5cm 8.5cm,clip]{SCARCE80_B_Wsquare.pdf}}}

\subfloat[RMSD = $3.366\times 10^{-1}$\label{fig:SCARCE10_C} ]{\frame{\includegraphics[width = 0.22\linewidth, trim=8.cm 11cm 8.cm 13cm, clip]{SCARCE10_C.pdf}}}
\quad
\subfloat[RMSD = $2.893\times 10^{-1}$\label{fig:SCARCE20_C}  ]{\frame{\includegraphics[width = 0.22\textwidth, trim=8cm 11cm 8cm 13cm, clip]{SCARCE20_C.pdf}}}
\quad
\subfloat[RMSD = $1.301 \times 10^{-1}$ \label{fig:SCARCE40_C}]{\frame{\includegraphics[width = 0.22\textwidth, trim=8cm 11cm 8cm 13cm, clip]{SCARCE40_C.pdf}}}
\quad
\subfloat[RMSD = $0.744\times 10^{-1}$ \label{fig:SCARCE80_C}]{\frame{\includegraphics[width = 0.22\textwidth, trim=8cm 11cm 8cm 13cm, clip]{SCARCE80_C.pdf}}}

\caption{``Loop” Experiment with typical randomly generated small datasets. The first line represents the input datasets. The corresponding outputs are shown on the second line.
In the outputs, we focus on the delimited regions in order to display the improvements of the interpolation.}
\label{fig:Loop-SparseSet-Random}
\end{center}
\end{figure*}

\subsection{Discussion}

As illustrated in the simulations,
the algorithm we designed is satisfactory.
Through a rudimentary gradient based method, it is possible to recover a smooth interpolation of the targeted orientation field while discarding some of the shortcomings of prior methods. For instance the continuity issue is solved with the introduction of the bisector model and a proper energy functional.
Indeed, the doubling phase step is based on extracting a pair of discontinuous real valued functions which are, in a second step, interpolated with smooth functions.
Conversely, we were able to provide a procedure that performs a smooth interpolation with the true orientation field as direct target.

Furthermore, let us mention some of the applications of the bisector line field interpolation in the context of the study of orientation fields. This generalization of the doubling phase with smooth functions allows the study of the line field with methods geared towards smooth functions.
Regarding the study of singularities, it is well known that interpolation can be used to locate and describe singularities of the orientation field. In the case of bisector line fields, this observation still holds true and we can illustrate it with an example. Singularities of $\mathrm{B}(\mathbf{X},\mathbf{Y})$ correspond to points of $\mathbb R^2$ such that either $\mathbf{X}$ or $\mathbf{Y}$ vanish. For instance, if $p\in \mathbb R^2$ is such that $\mathbf{X}(p)=0$ and $\mathbf{Y}(p)\neq 0$, the bisector line field $q\mapsto \mathrm{B}(\mathrm{D}_p\mathbf{X}\cdot q,\mathbf{Y}(p))$ acts as a linearization of $\mathrm{B}(\mathbf{X},\mathbf{Y})$
 at $p$ as soon as $\mathrm{D}_p\mathbf{X}$ is invertible (which can be assumed as this is the generic case). As a consequence, for instance,  one has  in this case
$$
\mathrm{ind}_p \mathrm{B}(\mathbf{X},\mathbf{Y})=\frac{1}{2}\mathrm{sign}(\det \mathrm{D}_p \mathbf{X}).
$$

The energy we proposed is the  natural choice when considering this problem and the gradient descent method is one possible direction to optimize it. Furthermore, its definition is adaptable and more can be done when information on the target is known.  For instance, one can introduce  weights depending on the quality of the first lift or prior knowledge of the singularity locations.

Finally, solving this optimization problem on the set of polynomial vector fields of a fixed degree is not a requirement of the method.
What is actually necessary to solve the problem with this methodology is the introduction of a family of smooth functions that serves as a basis for the space of regular vector fields on a bounded domain of $\mathbb R^2$. Hence this entire method can be adapted to the classical case of trigonometric polynomials.

\section{Conclusion}
In this paper we proposed a solution to the problem of interpolation of orientation fields with smooth functions.
To this end, we introduced a methodology based on the bisector line field model associated with a well suited energy functional. On the one hand, the bisector model has the double virtue of generalizing known techniques from the field of fingerprint analysis while resolving the continuity issues from the classical approach. On the second hand, the energy is coherent with unique aspects of this problem on the space of orientation fields and facilitates to use of gradient descent methods. Finally, the procedure has been applied to perform polynomial interpolation of orientation fields in the framework of fingerprint analysis.

\bibliographystyle{plain}


\section{Introduction and motivations: presentation of the main theorem}\label{IntroductionSection}
The biological models presented by Keller and Segel  in their landmarking papers \cite{K-S-1970,Keller-1971-MC}, and describing chemotaxis phenomena, have  been lately inspiring many interesting investigations in the fields of both theoretical and applied mathematics.  In this sense, aim of the present research is focusing on a precise variant (which, as technically detailed in \cite{Luca2003Alzheimer}, fits with real applications concerning aggregation phenomena of microglia observed in Alzheimer's disease) of the aforementioned models  and enhance its underlying  mathematical theory.

To be precise, this paper deals with the analysis of the mathematical problem idealizing the motion of a certain cell density $u(x, t)$ at the position $x$ and at the time $t$, initially distributed according to the law of $u_0(x):=u(x,0)$, and moving in an insulated domain under a repulsion effect, from a certain chemical signal concentration $w(x, t)$, and an attraction impact, from another one $v(x, t)$, which is ``slightly'' weaker (as specified later) than the first. Mathematically, we will face this system
\begin{equation}\label{problem}
\begin{cases}
u_{ t}=\Delta u -\chi \nabla \cdot (u\nabla v)+\xi \nabla \cdot (u\nabla w) & \textrm{in }\Omega, t>0, \\
0=\Delta v+\alpha u^\rho-\beta v & \textrm{in } \Omega,  t>0,\\
0=\Delta w+\gamma u-\delta w & \textrm{in } \Omega,  t>0,\\
 u_{\nu}= v_{\nu}= w_{\nu}=0 & \textrm{on } \partial \Omega, t>0, \\
u_0(x):=u(x,0)\geq 0 & x \in  \Omega,
\end{cases}
\end{equation}
defined in $\Omega$, a bounded and smooth domain of $\R^n$ with $n\geq 2$, and where $\alpha,\beta,$ $\gamma,\delta,\chi,\xi >0,$ $0<\rho<1$,  whilst ${(\cdot)}_{\nu}$ indicates the outward normal derivative on $\partial \Omega$ and $u_0(x)$ is a nonnegative and sufficiently regular function over $\bar{\Omega}$.
 Physically, the parameters $\chi$ and $\xi$ measure the influences of the attraction and repulsion, whereas the second and third equations idealize that chemoattractant and chemorepellent, $v$ and $w$, are sublinearly and linearly released by cells, moreover decaying with rates $\beta$ and $\delta$, respectively.

Similarly to what happens in the abundantly studied parabolic-elliptic Keller--Segel system, obtained by \eqref{problem} letting $\xi=0$, $\rho=1$ and eliminating the third unknown $w$, also model \eqref{problem} itself is likely to manifest  the so called \textit{chemotactic collapse}, the mechanism according to which the movement of the cells may eventually degenerate into aggregation processes giving rise to $\delta$-formations.

As far as this research is concerned, our accurate bibliographic research did not show any result about prototypes as that in \eqref{problem} presenting sublinear production. Subsequently, since we establish here that such weakening in the rate growth of the chemical signal associated to the cells' gathering suffices to prevent any impulsive concentration of the same cells independently by the coalescence effects coming from other factors of the model, we believe that this work provides a more complete picture about attraction-repulsion chemotaxis problems.

Conversely, confining now our attention to the parabolic-elliptic-elliptic system  \eqref{problem} with $\rho=1$ (we just mention that other variants with nonlinear diffusion and/or chemo-sensitivities or logistic sources are available in the literature),  it is known that for $n\geq 2$ and for  $\xi\gamma-\chi\alpha>0$ (repulsion dominates attraction)  it only admits globally bounded classical solutions, whilst in the bi-dimensional setting and for $\xi\gamma-\chi\alpha<0$  unbounded solutions can be detected (see \cite{EspejoSuzukiAttractionRepulsion,TaoWangAttractionRepulsion,YuGuoZhengAttractionRepulsion}).
Additionally, in \cite{CriticalMassAttrRepulGuoEtAl} the authors establish that for $n=2$ and $\chi\alpha-\xi\gamma>0$,  the value $\frac{4\pi}{\chi\alpha-\xi\gamma}$ is the critical mass deciding whether all solutions are global or, on the contrary, certain ones may blow-up; if this last scenario occurs, lower bounds of the blow-up time are estimated in \cite{Viglialoro2019RepulsionAttraction}.

Now, by virtue of all of the considerations, we rigorously formulate  our main result.
\begin{theorem}\label{Maintheorem}
Let $\Omega$ be a bounded and smooth domain of $\R^n$, $n\geq 2$. Then, for any $\alpha,\beta,\gamma,\delta,\chi,\xi>0$, $0<\rho<1 $, and any nonnegative and nontrivial initial data $0 \leq u_0(x)\in C^0(\bar{\Omega})$, problem \eqref{problem} admits a unique solution $(u,v,w)$ of nonnegative and bounded functions in the class
\begin{equation*}
C^0([0,\infty);C^0(\Omega))\cap C^{2,1}(\bar{\Omega}\times (0,\infty))\times C^{2,0}(\bar{\Omega}\times (0,\infty))\times C^{2,0}(\bar{\Omega}\times (0,\infty)).
\end{equation*}
\end{theorem}
\section{From local to globally bounded solutions}
One of the first steps in dealing with solutions of \eqref{problem} is showing that they do exist, at least locally.
\begin{lemma}\label{LocalExistenceLemma}
Let $\Omega$ be a bounded and smooth domain of $\R^n$, $n\geq 2$. Then, for any $\alpha,\beta,\gamma,\delta\chi,\xi,\rho>0$, and any nonnegative and nontrivial initial data $0 \leq u_0(x)\in C^0(\bar{\Omega})$, problem \eqref{problem} admits a unique solution $(u,v,w)$ of nonnegative functions in the class
\begin{equation*}
C^0([0,T_{max});C^0(\Omega))\cap C^{2,1}(\bar{\Omega}\times (0,T_{max}))\times C^{2,0}(\bar{\Omega}\times (0,T_{max}))\times C^{2,0}(\bar{\Omega}\times (0,T_{max})).
\end{equation*}
Here $T_{max}\in(0,\infty]$, denoting the maximal existence time, is such that (dichotomy criterion) or $T_{max}=\infty$  (global-in-time classical solution) or if $T_{max}<\infty$ (local-in-time classical solution) then necessarily
\begin{equation}\label{extensibility_criterion_Eq}
\limsup_{t\nearrow T_{max}}\lVert u (\cdot,t)\rVert_{L^\infty(\Omega)}=\infty.
\end{equation}
Moreover,
\begin{equation}\label{Bound_of_u}
\int_\Omega u (\cdot ,t)  = m:=\int_\Omega u_0>0\quad \textrm{for all}\quad  t\in (0,T_{max}).
\end{equation}
\begin{proof}
The first statements can be shown by straightforward
 adaptations of well-established methods involving an appropriate fixed point framework and standard parabolic and elliptic regularity theory (see, for instance,
\cite[Lemma 2.1]{FujieWinklerYokotaSignalSensi}), as well as related comparison principles. On the other hand, relation \eqref{Bound_of_u}  directly comes by integrating over $\Omega$  the equation for $u$ in \eqref{problem}.
\end{proof}
\end{lemma}
Once the solvability (at least in the local sense) for problem \eqref{problem} is ensured, the bridge establishing the globability and boundedness  is achieved throughout some precise $L^p$-bound for these solutions. To be precise we have this
\begin{lemma}\label{FromLocalToGLobalBoundedLemma}
Under the assumptions of Lemma \ref{LocalExistenceLemma}, let  $(u,v,w)$  be the classical solution to problem \eqref{problem}. If for some $\frac{n}{2}<p<n$ the $u$-component belongs  to $L^\infty((0,T_{max});L^p(\Omega))$, then  $T_{max}=\infty$ and $u$ is uniformly bounded in $\Omega \times (0,\infty)$.
\begin{proof}
Well-known elliptic regularity theory in conjunction with  Sobolev embedding theorems infer, through the second equation of \eqref{problem} and  $u\in L^\infty((0,T_{max});L^p(\Omega))$, that $v\in L^\infty((0,T_{max});W^{2,p}(\Omega))$, $\nabla v\in  L^\infty((0,T_{max});W^{1,p}(\Omega))$, and finally $ v\in  L^\infty((0,T_{max});C^{[2-(n/p)]}(\bar{\Omega}))$ and $\nabla v\in  L^\infty((0,T_{max});L^{q}(\Omega))$ for all $n<q<p^*:=\frac{np}{n-p}$. In particular, since the same reasoning is valid for $w$, by posing $\tilde{v}=\chi  u-\xi  w$ we have that for some positive constant $C_q$
\begin{equation}\label{Bound_v_1-q}
\lVert  \tilde{v} (\cdot, t)\lVert_{L^{q}(\Omega)}+\lVert  \nabla \tilde{v} (\cdot, t)\lVert_{L^{q}(\Omega)}  \leq C_q\quad \textrm{for all}\quad t\in(0,T_{max}).
\end{equation}
Additionally, for any $(x,t)\in \Omega \times (0,T_{max})$, the first equation of  \eqref{problem} reads $u_t=\Delta u-\nabla \cdot (u \nabla \tilde{v})$ so that for $t_0:=\max\{0,t-1\}$ the representation formula yields
\begin{equation*
\begin{split}
u (\cdot,t) &\leq e^{(t-t_0)\Delta}u(\cdot,t_0)-\int_{t_0}^t e^{(t-s)\Delta}\nabla \cdot (u (\cdot,s) \nabla \tilde{v} (\cdot,s))ds=: u_{1}(\cdot,t)+u_{2}(\cdot,t).
\end{split}
\end{equation*}
In these circumstances, the rest of the proof follows that done in \cite[Lemma 4.1]{ViglialoroWoolleyApplicableanalysis}; precisely,  in order to control the $L^\infty(\Omega)$-norm of $u$ on $(0,T_{max})$, first one can control a suitable  norm of the cross diffusion term $u \nabla \tilde{v}$ by replacing relation (24) therein with bound \eqref{Bound_v_1-q}, then applications of known smoothing estimates for the Neumann heat semigroup entail such uniform bound. Finally, the conclusion is achieved by relying on the dichotomy criterion  \eqref{extensibility_criterion_Eq}.
\end{proof}
\end{lemma}
\begin{remark}
From the above lemma, the open question given in \cite[Remark 1]{Viglialoro2019RepulsionAttraction} has a response: indeed, in the context of \cite[Theorem 3.1]{Viglialoro2019RepulsionAttraction}, if the $u$-component of the solution $(u,v,w)$ to the bi-dimensional version of problem \eqref{problem} with $\rho=1$ becomes unbounded at some finite time $t^*$ (in the sense of the $L^\infty(\Omega)$-norm) it also blows-up in the $L^p(\Omega)$-norm for any $p>1$, since otherwise from Lemma \ref{FromLocalToGLobalBoundedLemma} with $n=2$ it would be globally bounded. In particular, that theorem continues valid also without the extra assumption that $E(t):=\int_\Omega u^p \nearrow \infty$ as $t\nearrow t^*$, therein required.
\end{remark}
\section{Some properties of classical solutions: proof of the main theorem}
\subsection{Deriving a proper absorptive differential inequality}
With the crucial implication of Lemma \ref{FromLocalToGLobalBoundedLemma} in our hands, in this section we aim at bounding on $(0,T_{max})$ the functional $E(t):=\int_\Omega u^p$, for $p>1$, by means of a time independent constant. This will be obtained by deriving a proper absorptive differential inequality for $E$, exactly with the aid of this sequel of lemmas.

\begin{lemma}\label{EllipticEhrlingSystemLemma}
Under the assumptions of Lemma \ref{LocalExistenceLemma}, let $(u,v,w)$ be the classical solution to problem \eqref{problem} and $m:=\int_\Omega u_0$. Then for any $p>1$ and $\sigma>0$ there exists $\tilde{c}=\tilde{c}(p,\sigma) >0$ such that the $w$-component
satisfies, for any $\hat{c}>0,$
\begin{equation}\label{EhrlingTypeInequality}
\hat{c}\int_\Omega w^{p+1} \leq \sigma \int_\Omega u^{p+1}+\tilde{c}m^{p+1}\quad \textrm{for all } t \in (0,T_{max})
.
\end{equation}
Additionally, for $0<\theta=\frac{\frac{p}{2}-\frac{1}{2}}{\frac{p}{2}+\frac{1}{n}-\frac{1}{2}}<1$ and some constant $C_{GN}>0$, the $u$-component  fulfills
\begin{equation}\label{InequalityTipoG-N}
\int_\Omega u^p \leq \frac{4(p-1)}{p}\int_\Omega \lvert \nabla u^\frac{p}{2}\rvert^2 +c^* \quad \textrm{for all } t \in (0,T_{max}),
\end{equation}
being $c^*=2m^pC_{GN}^2[(1-\theta)m^p(\frac{2(p-1)}{p \theta C_{GN}^2})^\frac{\theta}{\theta-1}+1].$
\begin{proof}
As to the first conclusion, for the sake of completeness and clarity we retrace in detail what presented in \cite[Lemma 2.2]{WinklerHowFar} and \cite[Lemma 2.2]{LankeitChemoPrevent}. A direct integration over $\Omega$ of the third equation in \eqref{problem} produces, for any $p>1$, and using \eqref{Bound_of_u}
\begin{equation}\label{EqualityComingConservationMass}
\Big(\int_\Omega w \Big)^{p+1}  =\Big(\frac{\gamma}{\delta}\Big)^{p+1}\Big( \int_\Omega u \Big)^{p+1} =\Big(\frac{\gamma}{\delta}\Big)^{p+1}m^{p+1}\quad \textrm{ on } (0,T_{max}),
\end{equation}
whilst testing procedures and Young's inequality on the same equation yield
\begin{equation*}
\begin{split}
p\int_\Omega w^{p-1}\lvert \nabla w \rvert^2 &+\delta \int_\Omega w^{p+1} =\gamma \int_\Omega w^p u\leq \frac{4p}{(p+1)^2}\int_\Omega w^{p+1}  \\ & +\frac{\gamma^{p+1}}{4^p}(p+1)^{p-1}\int_\Omega u^{p+1}\quad \textrm{ for all } (0,T_{max}).
\end{split}
\end{equation*}
This, through the identity $ \lvert \nabla w^\frac{p+1}{2} \rvert^2=\frac{(p+1)^2}{4}w^{p-1} \lvert \nabla w \rvert^2$, reads for all $\eta \in (0,\frac{1}{2})$
\begin{equation}\label{InequalityLikeEhrling_1}
\eta \int_\Omega \lvert \nabla w^\frac{p+1}{2} \rvert^2 \leq \eta\int_\Omega w^{p+1}  +\eta \frac{\gamma^{p+1}}{4^{p+1}p}(p+1)^{p+1}\int_\Omega u^{p+1} \quad \textrm{ on } (0,T_{max}).
\end{equation}
On the other hand, for the same $\eta \in (0,\frac{1}{2})$, by virtue of the inclusions
\[W^{1,2}(\Omega)\hookrightarrow \hookrightarrow L^2(\Omega)\hookrightarrow L^\frac{2}{p+1}(\Omega),\]
Ehrling's Lemma (see \cite[Lemma 1.1]{ShowalterBOOK}) yields a constant $c_E(\eta)>0$ such that
\begin{equation*}
\|V\|_{L^2(\Omega)}^2\leq \eta \|V\|_{W^{1,2}(\Omega)}^2+c_E(\eta)\|V\|_{L^\frac{2}{p+1}(\Omega)}^2\quad \textrm{for all } V \in W^{1,2}(\Omega);
\end{equation*}
subsequently, posing in this last relation $V=w^\frac{p+1}{2}$, and making use of \eqref{EqualityComingConservationMass} and \eqref{InequalityLikeEhrling_1},  as well as of the conservation of mass property \eqref{Bound_of_u}, we obtain for all $t\in (0,T_{max})$
\begin{equation*
(1-2\eta)\int_\Omega w^{p+1} \leq \eta \frac{(\gamma(p+1))^{p+1}}{4^{p+1}p}\int_\Omega u^{p+1} +\Big(\frac{\gamma}{\delta}\Big)^{p+1}c_E(\eta) m^{p+1}.
\end{equation*}
Finally, for any $\hat{c}>0$ we introduce the function $\sigma: (0,\frac{1}{2})\rightarrow (0,\infty)$ defined as $\sigma(\eta)=\frac{\eta}{1-2\eta}\frac{(\gamma(p+1))^{p+1}\hat{c}}{4^{p+1}p}$, and estimate \eqref{EhrlingTypeInequality} follows with the choice $$\tilde{c}=\tilde{c}(p,\sigma):=\big(\frac{\gamma}{\delta}\big)^{p+1}\Big(\hat{c}+\frac{2\sigma}{\frac{(\gamma(p+1))^{p+1}\hat{c}}{4^{p+1}p}}\Big)c_E\Big(\frac{\sigma}{2\sigma+\frac{(\gamma(p+1))^{p+1}\hat{c}}{4^{p+1}p}}\Big).$$

In turn, the proof of \eqref{InequalityTipoG-N} comes from an application of a general case of the Gagliardo--Nirenberg inequality: in particular, for any $p>1$, we can use \cite[(22) of Lemma 4]{ViglialoroApplMathOpt2019} with $f=u^\frac{p}{2}$, $\mathfrak{p}=\mathfrak{q}=2$ and $\mathfrak{r}=\frac{2}{p}$ so to explicitly have
\begin{equation*
\begin{split}
\int_\Omega u^p&= \|  u^\frac{p}{2} \|_{L^{2}(\Omega)}^2 \leq C_{GN}^2 ( \| \nabla u^\frac{p}{2} \|_{L^{2}(\Omega)}^{\theta} \| u^\frac{p}{2} \|_{L^{\frac{2}{p}}(\Omega)}^{1 - \theta}+  \| u^\frac{p}{2} \|_{L^{\frac{2}{p}}(\Omega)})^2 \\ &
= 2C_{GN}^2m^{p(1-\theta)}\Big(\int_\Omega |\nabla u^ \frac{p}{2}|^2\Big)^\theta+2C_{GN}^2m^{p}\quad \textrm{ on } (0,T_{max});
\end{split}
\end{equation*}
hence, we conclude invoking the Young inequality with exponents $\theta$ and $(1-\theta).$ (We remark that the proof of inequality \eqref{InequalityTipoG-N} does not rely on the fact that $(u,v,w)$ solves \eqref{problem}. Such an estimate, indeed, holds for any general function belonging to $W^{1,2}(\Omega)\cap L^\frac{2}{p}(\Omega)$; nevertheless,  to facilitate the reading we preferred to present it in this way.)
\end{proof}
\end{lemma}
\begin{lemma}\label{MainInequalityU^2Lemma}
Under the assumptions of Lemma \ref{LocalExistenceLemma}, but for $0<\rho<1$, let $(u,v,w)$ be the classical solution to problem \eqref{problem} and $m:=\int_\Omega u_0$. Then for any $p>1$  and $\tilde{c}=\tilde{c}(p,\frac{\gamma\xi(p-1)}{3})$ taken from Lemma \ref{EllipticEhrlingSystemLemma}, the $u$-component satisfies
\begin{equation}\label{MainInequalityU^2}
\frac{d}{d t}\int_\Omega u^p \leq -\frac{4(p-1)}{p}\int_\Omega \lvert \nabla u^\frac{p}{2} \rvert^2+\bar{c} \quad \textrm{ for all} \quad  t \in (0,T_{max}),
\end{equation}
being $\bar{c}=c_1+\tilde{c}m^{p+1}$, with $c_1=\frac{\alpha\chi(p-1)(1-\rho)}{p+1}(\frac{(p+1)\gamma\xi}{(p+\rho)3\alpha \chi})^\frac{p+\rho}{\rho-1}|\Omega|$.
 \begin{proof}
For any $p>1$ by using problem \eqref{problem} and the divergence theorem (this in particular twice in both cross-diffusion terms), we have for all $t\in (0,T_{max})$
\begin{equation}\label{DerivativeE-FirstStep}
\begin{split}
\frac{d}{dt}\int_\Omega u^p&=p\int_\Omega u^{p-1} u_t  =p\int_\Omega u^{p-1} [\Delta u -\chi \nabla \cdot (u\nabla v)+\xi \nabla \cdot (u\nabla w) ]\\&
=-p(p-1)\int_\Omega u^{p-2} |\nabla u|^2-\chi\beta (p-1)\int_\Omega u^p v \\& \quad+ \alpha\chi (p-1)\int_\Omega u^{p+\rho} + \xi\delta (p-1)\int_\Omega u^pw-\gamma\xi (p-1)\int_\Omega u^{p+1}.
\end{split}
\end{equation}
On the other hand, if we neglect the nonpositive term $-\chi\beta (p-1)\int_\Omega u^p v $, use the Young inequality and  \eqref{EhrlingTypeInequality} with  $\tilde{c}$  as in our hypotheses (corresponding to the choice $\sigma=\frac{\gamma\xi(p-1)}{3}$),  $c_1$ as established in this statement  and $\hat{c}=\frac{\xi\delta(p-1)}{p+1}(\frac{(p+1)\gamma}{3p\delta})^{-p}$, we have that for all $t\in(0,T_{max})$ these relations are complied:
\begin{equation*}
\xi\delta (p-1)\int_\Omega u^pw \leq \frac{\gamma\xi(p-1)}{3}\int_\Omega u^{p+1}+\hat{c}\int_\Omega w^{p+1} \leq \frac{2\gamma\xi(p-1)}{3}\int_\Omega u^{p+1}+\tilde{c}m^{p+1},
\end{equation*}
and
\begin{equation*}
 \alpha\chi (p-1)\int_\Omega u^{p+\rho}\leq \frac{\gamma\xi(p-1)}{3}\int_\Omega u^{p+1}+c_1
 \quad \textrm{ on } (0,T_{max}).
\end{equation*}
Then, by virtue of the pointwise identity $u^{p-2}|\nabla u|^2=\frac{4}{p^2}|\nabla u^\frac{p}{2}|^2$ and the previous two inequalities, estimate
\eqref{DerivativeE-FirstStep} actually reads as claim \eqref{MainInequalityU^2}, with $\bar{c}=c_1+\tilde{c}m^{p+1}.$
\end{proof}
\end{lemma}
\begin{lemma}\label{MainfinalLemma}
Under the assumptions of Lemma \ref{LocalExistenceLemma},  but for $0<\rho<1$, let $(u,v,w)$ be the classical solution to problem \eqref{problem}. Then for all $p>1$ there exists $C>0$ such that
\begin{equation*
\int_\Omega u^p \leq C \quad \textrm{ for all} \quad  t \in (0,T_{max}).
\end{equation*}
 \begin{proof}
Collecting \eqref{InequalityTipoG-N}  and \eqref{MainInequalityU^2}, for $E(t)=\int_\Omega u^p$ and $c_{*}=c^{*}+\bar{c}$, we get the absorptive differential inequality $E'(t)\leq c_{*}- E(t)$ on  $(0,T_{max})$ which, complemented with the natural initial condition $E(0)=\int_\Omega u_0^p$, manifestly leads to $E(t)\leq \max\{E(0),c_{*}\}=:C$, for all $t\in(0,T_{max}).$
 \end{proof}
 \end{lemma}
\subsection{Proof of Theorem \ref{Maintheorem}}\label{MainTheoremProofSection}
 For any $n\geq 2$, let $(u,v,w)$ be the classical solution  to \eqref{problem} provided by Lemma \ref{LocalExistenceLemma}. By choosing $\frac{n}{2}<p<n$, we have  from Lemma \ref{MainfinalLemma} that $u\in L^\infty((0,T_{max});L^p(\Omega))$, so that in turn an application of  Lemma \ref{FromLocalToGLobalBoundedLemma}
immediately concludes the proof.
 \qed


\section{Introduction} \label{sec:intro}
\input{sources/introduction}

\section{Related Work} \label{sec:related}
\input{sources/related}

\section{Problem Formulation} \label{sec:problem}
\input{sources/problem}

\section{Learning Algorithms and Regret Bounds}\label{sec:case1}

\textbf{A Learning Algorithm for Case 1:}
\input{sources/case1}

\textbf{A Learning Algorithm for Case 2:}
\input{sources/case2}

\textbf{A Learning Algorithm for Case 3:}
\input{sources/case3}

\section{Experiments} \label{sec:experiments}
\input{sources/experiments}

\section{Conclusion and Future Work} \label{sec:conc}
We proposed a new multiobjective MAB problem, called the Lex-MAB, where the learner aims to minimize its regret with respect to lexicographic optimal arms. We studied the Lex-MAB with and without prior information on the expected rewards. We proved that the regret is uniformly bounded when the learner knows either the lexicographic optimal expected rewards or near-lexicographic optimal expected rewards. \reva{We also made a connection between knowing near-lexicographic optimal expected rewards and satisficing, and proved that uniformly bounded regret can be achieved for the MAB with satisficing objectives.} Interestingly, we showed that the regret incurred in an objective also depends on the prior information and observations from the other objectives. We also showed that the learner can still achieve $\tilde{O}(T^{2/3})$ gap-free regret without prior information. An important open question is whether $O(\sqrt{T})$ gap-free regret can be achieved in the Lex-MAB. This will answer the question of whether learning in the Lex-MAB without prior information is fundamentally more difficult than learning in the classical stochastic MAB without prior information. The case where there is prior information only for a subset of objectives is also worth investigating in the future. Another important future research direction is to extend the current work in the multiobjective MAB to minimize the regret with respect to a target arm other than a lexicogrpahic optimal arm by exploiting the prior information about the expected rewards of the target arm.

\newpage

\bibliographystyle{IEEEtran}

\subsection{MAB with Prior Information}

\citet{lai1985asymptotically} shows that in the classical stochastic MAB problem, for any uniformly good policy the regret grows at least logarithmically over time. As opposed to this, \citet{lai1984optimal} proves for the two-armed stochastic bandit that when the learner has prior information on the maximum expected reward $\mu^*$ and the minimum nonzero suboptimality gap $\Delta$, there exists policies that can achieve uniformly bounded regret. This idea is further investigated in \citet{bubeck2013bounded}, which shows that bounded regret of order $1/\Delta$ is achieved for the case with finitely many arms when the learner knows $\mu^*$ and a positive lower bound on $\Delta$.
\rev{\citet{garivier} studies the case where only $\mu^*$ is known and proposes an algorithm with bounded regret of order $\log(1/\Delta) (1/\Delta)$, and also proves a lower bound of order $1/\Delta$. This paper also provides a generic tool to prove gap-dependent lower bounds on the regret.} Similarly, \citet{bubeck2013prior} considers Thompson sampling and shows that its regret is uniformly bounded when $\mu^*$ and a positive lower bound on $\Delta$ are known. On the other hand, \citet{vakili2013achieving} considers a weaker prior information model where the learner knows a near-optimal expected reward $\eta$, which can be computed using $\mu^*$ and a positive lower bound on $\Delta$. The proposed algorithm obtains $\sum_a \Delta_a/\delta^3$ regret, where $\delta=\mu^*-\eta<\Delta$ and $\Delta_a$ is the suboptimality gap of arm $a$.
\citet{mersereau2009structured} and \citet{lattimore2014bounded} consider as prior information the knowledge of parametrized expected reward functions for each arm. In these works, the only unknown is the true parameter, which can be estimated by using reward observations from all of the arms.

Different from the works mentioned above, in this paper we consider a multiobjective MAB problem with lexicographically ordered objectives. We design algorithms that exploit the prior information in all objectives simultaneously to rule out arms that are not lexicographic optimal. \rev{Our regret bounds match the one in \cite{garivier} and improve the one in \cite{vakili2013achieving} for the case with a single objective.}

Numerous works have investigated regret minimization in multiobjective variants of the MAB problem. For instance, \citet{drugan2013designing} defines for each suboptimal arm its distance to the Pareto front as the Pareto suboptimality gap and the regret as the sum of the Pareto suboptimality gaps of the arms chosen by the learner. It proposes a learning algorithm that achieves $O(\log T)$ gap dependent Pareto regret. \citet{turgay2018multi} considers the multiobjective contextual MAB problem with similarity information, and extends the contextual zooming algorithm in \citet{slivkins2014contextual} to minimize the Pareto regret while making fair selections among the estimated Pareto optimal arms. The proposed algorithm is shown to achieve $\tilde{O}(T^{(1+d_p)/(2+d_p)})$ Pareto regret where $d_p$ is the Pareto zooming dimension. In addition, \citet{tekin2018multi} considers a biobjective contextual MAB problem with lexicographically ordered objectives. Unlike this work, we study the general case with $D$ lexicographically ordered objectives and focus on the importance of prior information in learning.

Our problem is also related to the tresholding MAB studied in \citet{locatelli2016optimal} and satisficing MAB studied in \citet{reverdy2017satisficing}. Specifically, if the prior information is given with respect to a target arm which is not necessarily a lexicographic optimal arm, and the regret in each objective is defined with respect to the expected reward of the target arm in that objective, then our algorithms for the cases with prior information can easily be adapted to minimize the regret with respect to the target arm. \reva{Specifically, we adapt our second case to the satisfaction-in-mean-reward problem introduced in \cite{reverdy2017satisficing} and prove that our algorithm achieves bounded regret for this problem as well.}\footnote{\reva{See Section \ref{sec:satisficing} for our results on the multiarmed bandit with satisficing objectives.}}