A bisector line field approach to interpolation of orientation fields

TL;DR

This paper introduces a geometric approach using bisector line fields for global orientation field reconstruction, enabling polynomial interpolation without doubling phase steps, demonstrated through fingerprint analysis examples.

Contribution

The paper presents a novel bisector line field method for orientation field interpolation that simplifies the process by avoiding doubling phase, with applications in fingerprint analysis.

Findings

Effective polynomial interpolation of orientation fields.

Bypasses the doubling phase step in orientation reconstruction.

Demonstrated success in fingerprint analysis examples.

Abstract

We propose an approach to the problem of global reconstruction of an orientation field. The method is based on a geometric model called "bisector line fields", which maps a pair of vector fields to an orientation field, effectively generalizing the notion of doubling phase vector fields. Endowed with a well chosen energy minimization problem, we provide a polynomial interpolation of a target orientation field while bypassing the doubling phase step. The procedure is then illustrated with examples from fingerprint analysis.

Tables

Table 1. Table 1: Distribution of final RMSD.

	RMSD mean	RMSD standard deviation
$\left|ℐ\right|=10$	$3.809\times {10}^{-1}$	$1.121\times {10}^{-1}$
$\left|ℐ\right|=20$	$2.637\times {10}^{-1}$	$8.150\times {10}^{-2}$
$\left|ℐ\right|=40$	$1.519\times {10}^{-1}$	$5.978\times {10}^{-2}$
$\left|ℐ\right|=80$	$7.873\times {10}^{-2}$	$2.346\times {10}^{-2}$

Equations

$$V_i=\left(\begin{array}{c}\cos 2{\theta }_i\\ \sin 2{\theta }_i\end{array}\right).$$

$$\mathrm{B}(\mathbf{X},\mathbf{Y}):\mathcal{D}\setminus ({\mathcal{Z}}_{\mathbf{X}}\cup {\mathcal{Z}}_{\mathbf{Y}})\to {\mathbb{P}}^1$$

$$\mathbf{X}(p)=\Vert \mathbf{X}(p)\Vert \left(\begin{array}{c}\cos {\theta }_{\mathbf{X}}(p)\\ \sin {\theta }_{\mathbf{X}}(p)\end{array}\right)$$

$$\mathbf{Y}(p)=\Vert \mathbf{Y}(p)\Vert \left(\begin{array}{c}\cos {\theta }_{\mathbf{Y}}(p)\\ \sin {\theta }_{\mathbf{Y}}(p)\end{array}\right).$$

$$\mathrm{B}(\mathbf{X},\mathbf{Y})(p)=\frac{1}{2}\left({\theta }_{\mathbf{X}}(p)+{\theta }_{\mathbf{Y}}(p)\right)(\mathrm{mod}\pi ).$$

$${\mathrm{ind}}_p(\mathrm{B}(\mathbf{X},\mathbf{Y}))=\frac{1}{2}({\mathrm{ind}}_p(\mathbf{X})+{\mathrm{ind}}_p(\mathbf{Y})).$$

$$d({\theta }_1,{\theta }_2)=\{\begin{array}{rlrl}& {\theta }_1-{\theta }_2& \text{ if }& \left|{\theta }_1-{\theta }_2\right|\le \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{array}$$

$$\mathrm{J}(\mathbf{M},\mathbf{L})=\sum _{i\in \mathcal{I}}d{\left(\mathbf{M}(x_i,y_i),\mathbf{L}(x_i,y_i)\right)}^2.$$

$$\nabla d^2({\theta }_1,{\theta }_2)=2d({\theta }_1,{\theta }_2)({\partial }_{{\theta }_1}-{\partial }_{{\theta }_2})$$

$$d_m(({\theta }_1,\dots ,{\theta }_m),({\theta }_1^{\prime },\dots ,{\theta }_m^{\prime }))={\left(\sum _{i=1}^{m}d({\theta }_i,{\theta }_i^{\prime }{)}^2\right)}^{1/2}$$

$$({\theta }_1,\dots ,{\theta }_m)=(\mathbf{M}(x_1,y_1),\dots ,\mathbf{M}(x_m,y_m))$$

$$({\theta }_1^{\prime },\dots ,{\theta }_m^{\prime })=(\mathbf{L}(x_1,y_1),\dots ,\mathbf{L}(x_m,y_m)),$$

$$\mathrm{J}(\mathbf{M},\mathbf{L}{)}^{1/2}=d_m(({\theta }_1,\dots ,{\theta }_m),({\theta }_1^{\prime },\dots ,{\theta }_m^{\prime })).$$

$$\mathbf{X}(x,y)=\sum _{k=0}^n\sum _{j=0}^k\left(\begin{array}{c}{\alpha }_{k,j}\\ {\beta }_{k,j}\end{array}\right)x^{k-j}{y}^j,$$

$$\sum _{k=0}^n\sum _{j=0}^k{\alpha }_{k,j}^2+{\beta }_{k,j}^2=1.$$

$$\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\}.$$

$$\omega =\left({\left({\alpha }_{k,j},{\beta }_{k,j}\right)}_{0\le j\le k\le m_{\mathbf{X}}},{\left({\gamma }_{k,j},{\delta }_{k,j}\right)}_{0\le j\le k\le m_{\mathbf{Y}}}\right);$$

$$\min \left\{\mathrm{J}(\mathbf{M},{\mathbf{L}}_{\omega })\mid \omega \in {\mathbb{S}}^{n_{\mathbf{X}}}\times {\mathbb{S}}^{n_{\mathbf{Y}}}\right\}.$$

$${\mathbf{L}}_{\omega }(x,y)=\frac{1}{2}\left({\theta }_{{\mathbf{X}}_{\omega }}(x,y)+{\theta }_{{\mathbf{Y}}_{\omega }}(x,y)\right),$$

$${\nabla }_{\omega }{\mathbf{L}}_{\omega }=$$

$$+\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}),$$

$$\begin{gathered} {\nabla }_{\omega }\mathrm{J}(\mathbf{M},{\mathbf{L}}_{\omega })= \\ -2\sum _{i\in \mathcal{I}}d\left(\mathbf{M}(x_i,y_i),{\mathbf{L}}_{\omega }(x_i,y_i)\right){\nabla }_{\omega }{\mathbf{L}}_{\omega }(x_i,y_i). \end{gathered}$$

$$\omega +\rho {\nabla }_{\omega }\mathrm{J}(\mathbf{M},{\mathbf{L}}_{\omega })\notin {\mathbb{S}}^{n_{\mathbf{X}}}\times {\mathbb{S}}^{n_{\mathbf{Y}}}.$$

$$\begin{array}[]{cccc}p:&\mathbb{R}^{n_{\mathbf{X}}}\times\mathbb{R}^{n_{\mathbf{Y}}}&\longrightarrow&\mathbb{S}^{n_{\mathbf{X}}}\times\mathbb{S}^{n_{\mathbf{Y}}}\\ &(\omega_{1},\omega_{2})&\longmapsto&\left(\dfrac{\omega_{1}}{\|\omega_{1}\|},\dfrac{\omega_{2}}{\|\omega_{2}\|}\right).\end{array}$$

$${\omega }_{k+1}\leftarrow p\left({\omega }_k+\rho {\nabla }_{\omega }\mathrm{J}\left(\mathbf{M},{\mathbf{L}}_{{\omega }_k}\right)\right).$$

$$RMSD(\overline{\mathbf{M}},\overline{\mathbf{L}})=\sqrt{\frac{{\sum }_{j\in \mathcal{J}}d{\left(\mathbf{M}({\bar{x}}_j,{\bar{y}}_j),\mathbf{L}({\bar{x}}_j,{\bar{y}}_j)\right)}^2}{|\mathcal{J}|}}$$

$${\mathrm{ind}}_p\mathrm{B}(\mathbf{X},\mathbf{Y})=\frac{1}{2}\mathrm{sign}(\det {\mathrm{D}}_p\mathbf{X}).$$

Full text

1 Introduction

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 [33, 28, 49, 22], liquid crystals arrangements in their nematic phase [10, 14, 13, 30] or the pinwheel structure of the visual cortex V1 of mammals [34, 25, 6, 7, 11, 35].

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 [49, 33, 48, 27, 22, 29], the estimation of fingerprint ridge topologies can be a necessary step before the use of high-level classification algorithms. In the review article [5], 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 [2, 28]. However, such a strategy is prone to introduce significant noise into the data [49]. Note that enhancement techniques can be applied in order to improve the situation [5, 3, 52]. 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 [31, 27]. 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 [38, 18, 8], 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 [math]— 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 —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 doubling phase step, sends each dataset angle from the orientation set to $\mathbb{R}^{2}$ by doubling its value and taking its sine and cosine, e.g. [28, 49, 31, 48, 47, 21, 37, 19, 4, 26]. 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 [39, 50, 51, 12, 46, 9, 36, 16, 5]. For instance, in the seminal paper [49], 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 [8] in the framework of differential geometry, appears as a natural extension of the doubling phase step.

The 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 Root Mean Square Deviation for orientation fields, which is acknowledged as well suited in order to measure the distance between the original (i.e. the ground truth) and the reconstructed orientation field [45, 5]. 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 [49, 37, 43, 32, 21, 20].

The rest of the article is organised as follows. Section 2 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 3. It can be skipped by readers not interested in these aspects. Section 4 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 5.

2 Strategy overview

The goal of orientation field interpolation is to transform a set of discretely localized orientation information into a global model, as illustrated in Figure 1. Starting from a fingerprint image, a local estimation of the target orientation field is obtained. Enhancement techniques that improve quality of the dataset, e.g. [5, 3, 52, 21, 50], 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 [36] or in low quality fingerprint impression.

Definition 1**.**

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

• •

$\mathbb{S}^{n}=\left\{x\in\mathbb{R}^{n+1}\mid\sum_{k=1}^{n+1}x_{k}^{2}=1\right\}$* denotes the unit n-sphere. In particular $\mathbb{S}^{1}$ is interpreted as the interval $[0,2\pi]$ where [math] is identified with $2\pi$.*

• •

The orientation set, denoted by $\mathbb{P}^{1}$ is the interval $[0,\pi]$ where [math] is identified with $\pi$.

• •

A vector field is a map $\mathbf{X}:\mathcal{D}\to\mathbb{R}^{2}$. We denote its zero set by $\mathcal{Z}_{\mathbf{X}}=\left\{p\in\mathcal{D}:\mathbf{X}(p)=(0,0)\right\}$.

• •

An orientation field is a map $\mathbf{L}:\mathcal{D}\to\mathbb{P}^{1}$.

A classical strategy, [28, 49, 31, 48, 47, 21, 37, 19, 4, 26], corresponding to the bottom branch in Figure 2, is to not treat the estimation of the fingerprint orientation field as such. Let $\left(\theta_{i}\right)$ denote a collection of elements of $\mathbb{P}^{1}$, then a collection of elements $\left(V_{i}\right)$ of $\mathbb{R}^{2}$ is obtained by:

[TABLE]

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 2, we propose a more straightforward strategy based on a model called the bisector line field.

Similarly to vector fields, orientation fields admit an intrinsic definition from the differential topology point of view, see [8] 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.

Definition 2**.**

Let $\mathcal{D}$ denote a subset of $\mathbb{R}^{2}$. Let $\mathbf{X}$ and $\mathbf{Y}$ be two vector fields on $\mathcal{D}$. The bisector line field of $(\mathbf{X},\mathbf{Y})$ is the orientation field

[TABLE]

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))$.

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 $\mathbb{S}^{1}$ such that

[TABLE]

and

[TABLE]

Then, as illustrated in Figure 3,

[TABLE]

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 [51], Taylor expansions [17, 47], orthogonal/Legendre/Chebychev polynomials [43, 32, 37, 4], Fourier expansions [49, 42, 41], and others [5]. For the sake of simplicity, we chose to implement a Taylor expansion model.

Remark 1**.**

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” [24, 8]. We chose the terminology that is more popular in the fingerprint analysis community.

3 Bisector line fields

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

3.1 Singularities of line fields

In our application case, singularities of orientation fields should be understood as unresolvable discontinuities of the map $\mathbf{L}:\mathbb{R}^{2}\to\mathbb{P}^{1}$. To this extent, we refer in the following to continuous orientation fields $\mathbf{L}:\mathcal{D}\to\mathbb{P}^{1}$, where $\mathcal{D}$ is the continuity domain of $\mathbf{L}$, and singularities of $\mathbf{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}$.

Definition 3**.**

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

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 $\mathbf{L}:\mathcal{D}\to\mathbb{P}^{1}$ be a continuous orientation field and let $p$ be an isolated singularity of $\mathbf{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}(\mathbf{L})$, is a half integer quantifying the winding of $\mathbf{L}$ around $p$ (for a precise construction of this object see, for instance, [8, Section 3.2]).

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

Proposition 1** ([8]).**

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 $\mathbf{L}:\mathcal{D}\to\mathbb{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 $\mathbf{L}=\mathrm{B}(\mathbf{X},\mathbf{Y})$.

Proposition 2**.**

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

[TABLE]

Remark 2**.**

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.

3.2 Generic properties of bisector line fields

Following Thom [44], it is understood, philosophically 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 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, e.g. [23, 1]).

When generic features of vector fields are considered, more can be said [8]. 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 1, we get the following property for generic smooth bisector line fields, as illustrated in Figure 4.

Proposition 3**.**

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 $\pm 1/2$.

Remark 3**.**

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 [8]. Another popular model relies on differential 2-forms on surfaces [40].

4 Interpolation as an energy minimization problem

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.

4.1 Discrete energy

The Root Mean Square Deviation for orientation fields is acknowledged as well suited to measure the distance between two orientation fields [45, 5]. As seen in Section 2, 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

[TABLE]

Let $\mathcal{I}=\left\{1,\ldots,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},\mathbf{L}:\mathcal{D}\subset\mathbb{R}^{2}\to\mathbb{P}^{1}$, such that $(x_{i},y_{i})_{i\in\mathcal{I}}\subset\mathcal{D}$, we set the least-squares energy functional to be

[TABLE]

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

Remark 4**.**

Consider the energy given by $d^{2}$ on the torus $\mathbb{P}^{1}\times\mathbb{P}^{1}$. The gradient of $d^{2}$ is

[TABLE]

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 5. 2. 2.

The two connected components of $\mathbb{P}^{1}\times\mathbb{P}^{1}\setminus\{\theta_{1}=\theta_{2}\pmod{\pi/2}\}$ are geodesically convex once $\mathbb{P}^{1}\times\mathbb{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. 3. 3.

Notice that $\mathrm{J}^{1/2}$ is a pseudo-metric (it is symmetric and satisfies the triangular inequality) on the space of $\mathbb{P}^{1}$-valued maps over a domain containing $(x_{i},y_{i})_{i\in\mathcal{I}}$.

Indeed $|d|$ is a distance over $\mathbb{P}^{1}$ and we de facto have the classical product metric over $\left(\mathbb{P}^{1}\right)^{m}$ given by

[TABLE]

Then with $\mathcal{I}=\left\{1,...,m\right\}$,

[TABLE]

and

[TABLE]

one has

[TABLE]

4.2 Polynomial bisector line fields

For a given target $\mathbf{M}$, we minimize $\mathrm{J}(\mathbf{M},\mathbf{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

[TABLE]

and

[TABLE]

In particular $\left(\alpha_{k,j},\beta_{k,j}\right)_{0\leq j\leq k}\in\mathbb{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}$,

[TABLE]

Let us denote

$n_{\mathbf{X}}=(m_{\mathbf{X}}+1)(m_{\mathbf{X}}+2)$ and $n_{\mathbf{Y}}=(m_{\mathbf{Y}}+1)(m_{\mathbf{Y}}+2)$; 2. 2.

$\omega\in\mathbb{S}^{n_{\mathbf{X}}}\times\mathbb{S}^{n_{\mathbf{Y}}}$ such that

[TABLE] 3. 3.

$\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};$ 4. 4.

$\mathbf{L}_{\omega}=\mathrm{B}(\mathbf{X}_{\omega},\mathbf{Y}_{\omega})$.

Problem (2) is then equivalent to solving

[TABLE]

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

4.3 Projected gradient descent approach

The function $d$, given in Equation (1), and Problem (3) have been designed to allow the use of gradient based optimisation algorithms. Indeed, since

[TABLE]

we have on $\mathbb{R}^{2}\setminus\mathcal{Z}_{\mathbf{X}_{\omega}}\cup\mathcal{Z}_{\mathbf{Y}_{\omega}}$

[TABLE]

where $(\cdot)_{\omega,i}$ denotes the $i^{th}$ component of $(\cdot)_{\omega}$.

As a consequence, we can immediately compute by the chain rule

[TABLE]

For $\rho>0$, it should be noted that in general:

[TABLE]

However, we solve this issue with the projection

[TABLE]

Thus in the gradient descent scheme, we iterate the recursive transform

[TABLE]

5 Examples and discussion

5.1 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 [15].

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 [28, 33]. 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 —i.e. $245\times 145-$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 [16, 36]. 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 2 and schematized in Figure 2, 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 [45, 21, 5], 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{\mathbf{L}}=\left(\mathbf{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:

[TABLE]

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.

5.2 Results

The result of the interpolation experiments with large datasets are presented in Figures 6 and 7. 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 [math] to $\pi$ are mapped to a light intensity (near [math], 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{\mathbf{L}}=\left(\mathbf{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 [math] to $\pi$ in the orientation. Likewise, singularities are points where all grey levels accumulate (similar to the pinwheel singularities observed in the visual cortex [34]).

In the case of a scarce dataset, the input is a small collection of triples $(x,y,\theta)\in\mathbb{R}^{2}\times\mathbb{P}^{1}$, hence, the results are presented in Figure 8 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 8b.

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 9a-9b and Figures 9c-9d). 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 [49]). 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 7, in particular the highlighted region in Figures 9e-9f). 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 8, 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 1. Representative examples of the batches are shown in Figure 10. With only 10 data-points, relevant information is easily missed, and singularities are not reconstructed, see Figure 10e. The singularities appear with more datapoints, but their distribution influences the quality of the placement, see Figure 10f-10g-10h.

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 8). 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.

5.3 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

[TABLE]

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.

6 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.