Nasal Patches and Curves for Expression-robust 3D Face Recognition

TL;DR

This paper introduces a novel nasal region-based 3D face recognition method that uses keypoint detection, surface normal features, and genetic algorithms to achieve high accuracy across multiple datasets, outperforming many existing methods.

Contribution

The paper presents a new five-step nasal region-based recognition algorithm that is robust to expressions and does not require complex alignment or training.

Findings

Achieves high recognition accuracy on FRGC, Bosphorus, and BU-3DFE datasets.

Outperforms many state-of-the-art 3D face recognition algorithms.

Does not rely on sophisticated alignment or denoising steps.

Abstract

The potential of the nasal region for expression robust 3D face recognition is thoroughly investigated by a novel five-step algorithm. First, the nose tip location is coarsely detected and the face is segmented, aligned and the nasal region cropped. Then, a very accurate and consistent nasal landmarking algorithm detects seven keypoints on the nasal region. In the third step, a feature extraction algorithm based on the surface normals of Gabor-wavelet filtered depth maps is utilised and, then, a set of spherical patches and curves are localised over the nasal region to provide the feature descriptors. The last step applies a genetic algorithm-based feature selector to detect the most stable patches and curves over different facial expressions. The algorithm provides the highest reported nasal region-based recognition ranks on the FRGC, Bosphorus and BU-3DFE datasets. The results are comparable with, and in many cases better than, many state-of-the-art 3D face recognition algorithms, which use the whole facial domain. The proposed method does not rely on sophisticated alignment or denoising steps, is very robust when only one sample per subject is used in the gallery, and does not require a training step for the landmarking algorithm. https://github.com/mehryaragha/NoseBiometrics

Tables

Table 1. TABLE I: Landmarking consistency error in mm.

Dataset	Nasal root
(L1)	Left eye corner
(L2)	Left alar groove
(L3)	Subnasale
(L5)	Right alar groove
(L6)	Right eye corner
(L7)
Bosphorus	1.06 $\pm$ 0.58	1.76 $\pm$ 1.03	1.06 $\pm$ 0.62	1.11 $\pm$ 0.38	1.19 $\pm$ 0.60	2.12 $\pm$ 1.14
FRGC	2.04 $\pm$ 1.09	2.95 $\pm$ 1.61	1.29 $\pm$ 0.82	1.86 $\pm$ 0.85	1.22 $\pm$ 0.62	2.91 $\pm$ 1.53

Table 2. TABLE II: Landmarking precision accuracy from the ground truth over the Bosphorus dataset samples.

Algorithm	Threshold (mm)	L3	L6	L2	L7	Nose tip L4
Proposed
method	10
12
15
20	99.55%
99.62%
99.66%
99.69%	99.35%
99.62%
99.66%
99.66%	96.69%
97.73%
99.04%
99.59%	94.59%
96.56%
98.38%
99.62%	97.52%
99.04%
99.66%
99.79%
Creusot et al. [37]	10
12
15
20	97.96%
99.18%
99.82%
99.90%	98.43%
99.71%
99.86%
99.90%	98.82%
99.65%
99.93%
99.93%	98.50%
99.43%
99.75%
99.86%	95.47%
98.15%
98.97%
99.33%

Table 3. TABLE III: R 1 RR performance for varying the training size per subject, when all samples of the FRGC dataset are merged from the three seasons.

Feature	Number of gallery samples per subject (No. gallery samples/ No. probe samples)
descriptors	1
(482/4330)	2
(880/3848)	3
(1206/3408)	4
(1432/3006)	5
(1610/2648)	6
(1752/2326)	7
(1757/2034)
Spherical patches	96.2%	98.9%	99.4%	99.6%	99.6%	99.7%	99.75%
Curves	91.6%	96.8%	98.1%	98.8%	98.9%	99.3%	99.4%

Table 4. TABLE IV: Performance comparison on the FRGC dataset.

Algorithm	Modality	Rank-one
FRGC v2.0	EER ROC III	0.1% FAR
ROC III	Neutral vs.
Neutral	Neutral vs.
Non-neutral
Spherical patches
Curves	3D Nose	97.9%
94.1%	2.4%
4.9%	93.5%
80.0%	98.45% (R1RR)
95.8% (R1RR)	98.5% (R1RR)
97.5% (R1RR)
Smeets et al. [34]	3D Face	89.6%	3.8%	77.2%	-	-
Osaimi et al. [20]	3D Face	96.5%	-	94.05%	98.35% (0.1% FAR)	97.8% (0.1% FAR)
Spreeuwers et al. [21]	3D Face
3D Nose	99.0%
94.5%	-	94.6%
83.7%	-	-
Drira et al. (2013) [24]	3D Face	97.0%	-	97.1%	99.2% (R1RR)	96.8% (R1RR)
Alyüz et al. (2010) [6]	3D Face
3D Nose	97.5%
91.81	1.91%
-	85.6%
-	98.39% (R1RR)
-	96.40% (R1RR)
-
Wang et al. (2010) [23]	3D Face	98.39%	-	98.04%	99.2% (0.1% FAR)	97.7% (0.1% FAR)
Wang et al. (2008) [5]	3D Nose	95% (44mm)
78% (24mm)	-	-	-	-
Drira et al. (2009) [7]	3D Face/Nose
(125 gallery)
(125 probe)	88% (Face)
77.5% (Nose)	-	-	-	-
Chang et al. [1]	3D Nose	-	Neutral 12%
Non-neutral 23%	-	97.1% (R1RR)	86.1% (R1RR)
Emambakhsh et al. [4]	3D Nose	89.61%	Neutral 8%
Non-neutral 18%	-	90.87% (R1RR)	81.61% (R1RR)
Li et al. (2014) [15]	3D Face	96.3%	-	-	98.0% (R1RR)	94.2% (R1RR)
Queirolo et al. [22]	3D Face	99.6%	-	96.6%	99.5% (R1RR)	94.8% (R1RR)
Berretti et al. (2013) [17]	3D Face	95.6%	-	-	97.3% (R1RR)	92.8% (R1RR)
Berretti et al. (2010) [16]	3D Face	94.15%	-	-	$\approx$97.3% (R1RR)	$\approx$91.0% (R1RR)
Mohammadzade et al. (2013) [14]	3D Face	-	-	99.2%	-	-
Mian et al. (2008) [33]	3D Face	93.5%	-	Neutral 99.9%
Non-neutral 92.7%	99%	86.7%
Mian et al. (2007) [9]	2D+3D Face
2D+3D Nose	95.91%
$\approx$92.2%	-	99.3%
92.5%	99.2%
$\approx$94.9%	95.37%
$\approx$80.0%

Table 5. TABLE V: Comparison of some of the previous works on the 105 subjects in the Bosphorus dataset: the numbers in parentheses show the No. gallery samples/ No. of probe samples.

Algorithm	Modality and size	R1RR
Spherical patches
Curves	3D Nose (105/2797)	95.35%
86.1%
Li et al. (2014) [15]	3D Face (105/2797)	95.4%
Dibeklioğlu [8]	3D Nose (47/1527)
(47/423) [rotated]	89.2%
62.6%
Li et al. (2011) [35]	3D Face (105/4561)	94.1%
Alyüz et al. (2008) [39]	3D Face (34/441)
(47/1508)	95.9%
95.3%

Table 6. TABLE VI: R 1 RR for different expression types from the Bosphorus dataset used as probes, with 105 neutral samples used as the gallery.

Algorithm	Facial expression (No. of probe samples)
	Happy (106)	Surprise (71)	Fear (70)	Sadness (66)	Anger (71)	Disgust (69)	Neutral (194)
Spherical patches
Curves	98.08%
85.85%	100%
92.96%	98.55%
87.14%	96.92%
92.31%	94.12%
84.06%	88.24%
69.12%	98.96%
96.88%
Li et al. [35]	95.28%	98.59%	92.86%	95.45%	88.73%	76.81%	100%

Table 7. TABLE VII: R 1 RR for different expression types from the BU-3DFE dataset used as probes, with 100 neutral samples used as the gallery.

Algorithm	Facial expression (No. of probe samples)
	Happy (400)	Surprise (400)	Fear (400)	Sadness (400)	Anger (400)	Disgust (400)
Spherical patches
Curves	88.5%
81.8%	91.0%
87.8%	89.8%
85.3%	92.3%
87.6%	90.1%
83.2%	81.8%
69.8%
Hajati et al. [36]	86.0%	84.0%	82.0%	85.0%	93.0%	79.0%

Equations

$$\left\{\begin{array}[]{l}{\bf{P}_{r}}={\bf{P}}\times\left[\begin{array}[]{cc}\cos\alpha&\sin\alpha\\ -\sin\alpha&\cos\alpha\end{array}\right]=[{\bf{X_{f}^{\alpha}}},{\bf{Y_{f}^{\alpha}}}]\\ {\bf{MIN}}=V_{n,\alpha}({\bf{P}_{r}})\end{array}\right.,$$

$${\mathbf{SMIN}}_m=V_{1,\beta }({\mathbf{S}}_{{\gamma }_m}^m)$$

$${\theta }_z^k=\arccos \frac{\left|[L{4}_x^j-L{1}_x^i,L{4}_y^j-L{1}_y^i]\cdot {\hat{a}}_y\right|}{\left|[L{4}_x^j-L{1}_x^i,L{4}_y^j-L{1}_y^i]\right|},$$

$$E^k({L1}^i,{L4}^j)=\underset{y\in {\mathbf{S}}_y^k}{\max }\left(\sum _x\left|{\mathbf{Z}}_L^k-{\mathbf{Z}}_R^k\right|\right)$$

$${\theta }_z^{opt}=\underset{{\theta }_z^k}{\mathrm{argmin}}{E}^k({L1}^i,{L4}^j).$$

$$r=\left\{\begin{array}[]{cc}r_{0}\cos^{a_{1}}(\theta)&0\leq\theta<\pi\\ r_{0}\cos^{a_{2}}(\theta)&\pi\leq\theta<2\pi\end{array}\right.$$

$$r^{\prime}=\left\{\begin{array}[]{cc}r_{0}^{\prime}\cos^{a_{3}}(\theta^{\prime})&0\leq\theta^{\prime}<\pi\\ r_{0}^{\prime}\cos^{a_{4}}(\theta^{\prime})&\pi\leq\theta^{\prime}<2\pi\end{array}\right.,$$

$$\left\{\begin{array}[]{l}{{\bf{RMIN}}_{j}}=V_{3,\zeta}({{\bf{Q}}^{j}}_{\zeta})\\ {{\bf{LMIN}}_{j}}=V_{3,\eta}({{\bf{Q}}^{j}}_{\eta})\end{array}\right.$$

$$\left\{\begin{array}[]{l}{{\bf{RMIN}}_{k}^{\prime}}=V_{3,\zeta}({{\bf{Q}}^{\prime k}}_{\zeta})\\ {{\bf{LMIN}}_{k}^{\prime}}=V_{3,\eta}({{\bf{Q}}^{\prime k}}_{\eta})\end{array}\right.,$$

$$\left\{\begin{array}[]{l}{\bf{n}}_{s}=\nabla\left[{\bf{N_{x}}},{\bf{N_{y}}},{\bf NGm}_{s}\right]\\ s=1,2,\ldots,s_{m}\end{array}\right.$$

$$\left\{\begin{array}[]{l}{\bf{F}}=\left[{\bf{F}}_{s_{1}},{\bf{F}}_{s_{2}},\ldots,{\bf{F}}_{s_{n}}\right],\\ {\bf{F}}_{s_{k}}=\left[{\bf{Fx}}_{s_{k}},{\bf{Fy}}_{s_{k}},{\bf{Fz}}_{s_{k}}\right],\end{array}\right.$$

$$\left\{\begin{array}[]{l}{\bf{Fx}}_{s_{k}}=\left[{\bf{Hx}}_{1,s_{k}},{\bf{Hx}}_{2,s_{k}},\ldots,{\bf{Hx}}_{K,s_{k}}\right],\\ {\bf{Fy}}_{s_{k}}=\left[{\bf{Hy}}_{1,s_{k}},{\bf{Hy}}_{2,s_{k}},\ldots,{\bf{Hy}}_{K,s_{k}}\right],\\ {\bf{Fz}}_{s_{k}}=\left[{\bf{Hz}}_{1,s_{k}},{\bf{Hz}}_{2,s_{k}},\ldots,{\bf{Hz}}_{K,s_{k}}\right].\end{array}\right.$$

$$\left\{\begin{array}[]{l}{\bf{B}}_{s_{k}}=\left[{\bf{B}}_{1,s_{k}},{\bf{B}}_{2,s_{k}},\ldots,{\bf{B}}_{K,s_{k}}\right],\\ {\bf{B}}_{i,s_{k}}=\left\{\begin{array}[]{l}\left[0,0,\ldots,0\right]\text{ if }{\bf{Bn}}(i)=0\\ \left[1,1,\ldots,1\right]\text{ if }{\bf{Bn}}(i)=1\end{array}\right..\end{array}\right.$$

$$\mathbf{B}=\left[\underset{\text{length}=s_n\times 3\times K\times h_l}{\underbrace{\stackrel{\text{for all normals in scale 1}}{\overbrace{[{\mathbf{B}}_{s_1},{\mathbf{B}}_{s_1},{\mathbf{B}}_{s_1}]}},\dots ,\stackrel{\text{for all normals in scale n}}{\overbrace{[{\mathbf{B}}_{s_n},{\mathbf{B}}_{s_n},{\mathbf{B}}_{s_n}]}}}}\right].$$

$${\mathbf{D}}_{\mathbf{g},\mathbf{p}}=-\left(\frac{{\mathbf{X}}_{\mathbf{g}}}{\sqrt{|{\mathbf{X}}_{\mathbf{g}}{\Sigma }^{-1}{{\mathbf{X}}_{\mathbf{g}}}^{⊺}|}}\right){\Sigma }^{-1}{\left(\frac{{\mathbf{X}}_{\mathbf{p}}}{\sqrt{|{\mathbf{X}}_{\mathbf{p}}{\Sigma }^{-1}{{\mathbf{X}}_{\mathbf{p}}}^{⊺}|}}\right)}^{⊺}.$$

$${\mathbf{Bn}}_{opt}=\underset{\mathbf{Bn}}{\mathrm{argmax}}\{R_1\}$$

Full text

https://github.com/mehryaragha/NoseBiometrics

I Introduction

While much previous research on expression invariant 3D face recognition has focused on modelling expressions and detecting expression insensitive facial parts, there have been relatively few studies evaluating the potential of the nasal region for addressing this issue. Despite this, the nose has a number of salient features that make it suitable for expression robust recognition. It can be easily detected, due to its discriminant curvature and convexity [1], is difficult to hide without attracting suspicion [2, 3], is relatively stable over various facial expressions ([1, 4, 5, 6, 7, 8, 9]) and is rarely affected by unintentional occlusions caused by hair and scarves. Although it has been reported that the 2D image of the nose has too few discriminant features to be used as a reliable region for human identification [10], its 3D surface has much undiscovered potential. This paper further investigates the 3D nasal region for human identity authentication and verification purposes and presents a novel algorithm that provides very high discriminant strength, comparable with recent 3D face recognition algorithms, which use the whole facial domain.

The proposed approach is based on a very consistent and accurate landmarking algorithm, which overcomes the issue of robust segmentation of the nasal region. The algorithm first finds an approximate location of the nose tip and then finely tunes its location, while accurately determining the position of the nasal root and detecting the symmetry plane of the face. Next, the locations of three sets of landmarks are found: subnasale, eye corners and nasal alar groove. These landmarks are utilised on feature maps created by applying multi-resolution Gabor wavelets to the surface normals of the depth map. Two types of feature descriptors are used: spherical patches and nasal curves. Feature selection is then performed using a heuristic genetic algorithm (GA) and, finally, the expression-robust feature descriptors are applied to the well-known and widely used 3D Face Recognition Grand Challenge (FRGC) [11], Bosphorus [12] and Binghamton University 3D Facial Expression (BU-3DFE) [13] datasets.

Results show the algorithm’s high potential to recognise nasal regions, and hence faces, over different expressions, with very few gallery samples per subject. The highest rank-one recognition rates (R1RR) achieved are: 1) a R1RR of 97.9% and equal error rate (EER) of 2.4% for FRGC v2.0 and receiver operator characteristic (ROC) III experiments, respectively; 2) a R1RR of 98.45% and 98.5% for FRGC’s neutral vs. neutral and neutral vs. non-neutral samples, respectively; 3) a R1RR of 96.2% when one gallery sample per subject is used for the FRGC dataset (482 gallery samples (subjects) vs. 4330 probe samples); 4) a R1RR of 95.35% for the Bosphorus dataset when 2797 scans of 105 subjects are used as probes and the set of 105 neutral scans (one per subject) is used as the galley.

The remainder of the paper is organized as follows. After the literature review provided in section II, the alignment and nasal region cropping steps, followed by the nasal region landmarking, are detailed in section III. The feature extraction algorithm is described in section IV and section V explains the feature descriptors used. The feature selection algorithm is detailed in section VI and experimental results, including a thorough comparison with previous work, is provided in section VII. Finally, conclusions are given in section VIII.

I-A Scientific contribution and comparison with previous work

The major contribution of this paper is a novel surface normal-based recognition algorithm that provides a thorough evaluation of the recognition potential of the 3D nasal region. The results achieved are not only better than previous 3D nose recognition algorithms but also higher than many recognition algorithms that employ the whole face. The algorithm employs a novel, training-free, highly consistent and accurate landmarking algorithm for the nasal region and a robust feature space, based on the response of Gabor wavelets to surface normal vectors, is also introduced. To localise the expression robust regions on the nose a heuristic GA feature selection is applied to two different geometrical feature descriptors. Because of the smoothing effects of the Gabor wavelets, there is no need for sophisticated denoising algorithms. Indeed, only simple median filtering is required for the surface normals, even with noisy datasets such as the FRGC Spring 2003 folder. An additional advantage of the proposed approach is that a fast Principal Component Analysis (PCA)-based self-dependent method can be employed for facial pose correction. This eliminates the need for sophisticated pose correction algorithms or reference faces for fine tuning the alignment.

The proposed approach significantly extends our previous work [4] in which the nasal landmarking and recognition was performed on the depth map. This paper increases the number of landmarks and their detection accuracy and presents new feature extraction and selection algorithms. The work is inspired by recent algorithms on utilising facial normal vectors in 3D [14] and regional normal vectors [15]. To compare the new algorithms with previous approaches which used similar methodologies, the application of normals, computed over the nasal surface, is used for identification as well as the verification scenario. By using multi-resolution Gabor wavelets the ability of the algorithm to handle more noisy samples is enhanced, providing higher R1RR than the approach of Li et al. [15], which excluded the noisy FRGC Spring 2003 samples. This work also extends the application of facial curves, introduced as feature descriptors by Berretti et al. ([16] and [17]), to nasal spherical patches, producing a R1RR increase of $>2\%$, and showing a higher class separability for the spherical patches than for curves for 3D face recognition.

II Recent literature review

Robustness against the deformations caused by facial expressions has been a popular research topic in 3D face recognition. The face is a non-rigid object and therefore 3D matching techniques for rigid objects, such as the iterative closest point (ICP) algorithm [18], can become trapped in local minima and fail to provide accurate matching scores.

An empirical approach to deal with the variations caused by expressions is to capture a range of facial expressions for each subject and store them in the gallery [19]. Then, the facial biometric features of each test subject can be compared with all the stored expressions and a decision made on the identity of the subject. This method has numerous disadvantages: capturing a range of facial expressions for each subject is not always straightforward and requires a high storage capacity per subject. In addition, facial expressions will not necessarily remain constant and may differ between the test and gallery captures [19].

One approach to overcome this problem is to use computer graphics algorithms to artificially create different expressions for each facial capture. In [20], expressions are learned using PCA eigenvectors and then used to re-generate the expressions on the probe samples. Although this approach does not require multiple samples per subject in the gallery, it is still vulnerable to the number of training samples used to model the facial expressions. Also, a universal definition of facial expression for all subjects still remains to be found [19] and the need to classify the expression types prior to face recognition increases the computational complexity.

Another approach is to employ region-based methods, in which the least variant parts of the face over different expressions are detected using facial segmentation [1, 21, 6, 22] or extracted using their expression invariant capabilities [14, 23, 5, 15]. Spreeuwers proposes a multiple regional approach based on a PCA-Linear Discriminant Analysis (LDA) feature extraction method [21]. In regional recognition, scores from a combination of different masks on the nose, cheek, forehead, chin and mouth are fused to finalise the decision making.

Alyüz et al. use a regional registration algorithm in conjunction with LDA classifiers, giving an expression robust 3D face recognition approach [6]. They also demonstrate that the nasal region has a high discriminatory power. A focus on integrating multiple regions is provided by Queirolo et al. [22] in which four regions (the upper face image, the whole face and two nasal regions) are segmented and stored for the gallery sessions before matching is performed using a novel matching criterion, called the surface interpenetration measure, and simulated annealing.

Using facial curves is another popular approach to 3D face recognition that can be categorised as a subset of regional algorithms. Drira et al. use the intersections of planes with the facial surface to define a set of radial curves which pass through the nose tip, and then perform a quality assessment in order to handle missing data and occlusions [24]. Another curve-based algorithm is proposed by Berretti et al. [17]. First, keypoints are detected on the facial surface and then the least variant curves on the face are selected using a statistical model and matched with those in the gallery. As an extension to curves, isogeodesic stripes centralised on the nose tip are used in an expression invariant 3D face recognition method that employs a novel descriptor, termed the 3D weighted walkthroughs, to quantify the differences between corresponding stripes [16]. In another curve-based approach, Drira et al. find geodesic curves on the nasal region for a subset of the FRGC dataset [7].

To overcome the sensitivity of holistic face recognition algorithms to expression variations, Mian et al. propose a landmark-based method, in conjunction with a localised feature descriptor that incorporates the 2D texture and 3D point clouds. In an alternative approach, Wang et al. apply shape difference boosting to the Bosphorus dataset to learn the expressions and identify those facial regions which remain constant over different expressions [23]. Instead of using depth or the point coordinates for 3D registration, Mohammadzade et al. use the surface normals of the points in conjunction with a Fisher’s discriminant paradigm [14]. This approach selects the normals which maximise the concentration of within-class scatter while simultaneously maximising the between-class distribution. Recently, Li et al. proposed local normals histograms, captured from multiple rectangular regions on the face, to set up an expression-robust feature space and use a novel sparse classifier to perform the matching [15].

Despite the robustness of these algorithms against facial expressions, they often rely on accurate and consistent facial segmentation, which is not a straightforward task in 3D. To address this issue, some researches have focused on the nasal region, which shows high consistency over different expressions. For example, in one of the first investigations on 3D nose recognition, Chang et al. initially segment the face into different non-overlapping regions, using the curvature information [1]. Then, three overlapping nasal regions are detected and stored in the gallery. The same regions are segmented in the probe images and matched using the ICP algorithm. Wang et al. propose the use of local shape difference boosting for 3D face recognition and also apply the boosting algorithm to different nasal regions [5]. The regions are cropped using the intersection of spheres of radius $r$, centred on the nose tip, with the face surface. When the value of $r$ was increased, the recognition ranks reached a maximum and then plateaued. A combination of the nasal region, forehead and eyes are used for a 2D/3D face recognition by Mian et al. [9]. A modified ICP algorithm is used for matching, in conjunction with a pattern rejector based on spherical face representation (SFR) and shift-invariant feature transform (SIFT) [25], producing high recognition ranks on the FRGC dataset, in particular for the neutral probes. Dibeklioğlu et al. used the Dijkstra algorithm to segment the nose and evaluated the performance using a subset of the Bosphorus dataset [8].

III Preprocessing and nasal region landmarking

The algorithm explained in [4] is used to crop the face. Next, median filtering with a $2.5\times 2.5$ mm2 mask size is applied twice on the cropped face. The image is then resampled to a uniform grid with $0.5$ mm/pixel horizontal and vertical resolutions using Delaunay triangulation and aligned using the iterative PCA algorithm [9]. The aligned face is then intersected with three cylinders to crop the nasal region, according to [4]. The depth map of the cropped nasal region is again median filtered with a $2.5\times 2.5$ mm2 mask to further smooth its surface and decrease the spike noise effects. The block diagram in Fig. 1-a shows how the landmarks in Fig. 1-b are detected.

III-A Local minima detector, nose tip re-localisation, nasal root and subnasale detection

First, an initial position of the nasal root ($\bf{L1^{0}}$) is detected by [4]. Then, the location of the nose tip ($\bf{L4^{0}}$), found in section III, is more finely tuned. Various planes, passing through $\bf{L4^{0}}$ with normals $\cos(\theta_{i})\hat{a}_{x}+\sin(\theta_{i})\hat{a}_{y}$ are intersected with the nose surface, where $\theta_{i}$ is the angle of the $i^{th}$ plane with the $y$-axis, and $\hat{a}_{x}=[1,0]$ and $\hat{a}_{y}=[0,1]$ are the unit vectors along the $x$ and $y$ axes, respectively. This process results in several curves on the nasal region, shown in Fig. 3-a.

The proposed landmarking algorithm relies on a minima detector, which finds a set of minima on rotated versions of the curves and then maps them to the original curve. The rotation is required because some of the original curves are strictly decreasing functions that do not have an actual minimum. Assuming ${\bf{P}}=[{\bf{X_{f}}},{\bf{Y_{f}}}]$ is a $K\times 2$ matrix representing points of a curve, instead of directly differentiating ${\bf{P}}$ to find the minima on the curves, as proposed by Segundo et al. [26], the curves are first rotated in the $z$-axis (roll direction) by an angle $\alpha$ around a given point on the curve. This operation is given by,

[TABLE]

where ${\bf{P}_{r}}$ is the rotated version of ${\bf{P}}$ and the function $V_{n,\alpha}(.)$ finds the location of the $n$ smallest local minima on ${\bf{P}_{r}}$ and then remaps them to the original curve $\bf P$ using the rotation angle $\alpha$. The output ${\bf{MIN}}$ is an $n\times 2$ matrix, containing the locations of the $n$ local minima. $V_{n,\alpha}(.)$ computes the first order differentiation (first order difference in discrete space) of $\bf P_{r}$, which is then given to the signum function to detect its sign changes. This finds the locations of all the local minima in ${\bf{P}_{r}}$, which are then sorted based on their value in ascending order and the $n$ lowest are selected and rotated back to the original curve using $\alpha$. Fig. 2 shows an example of this procedure for $n=1$ (the global extremum).

The value allocated to $\alpha$ should be small enough in order to preserve the single-valued functionality of ${\bf{X_{f}^{\alpha}}}\rightarrow{\bf{Y_{f}^{\alpha}}}$, i.e. each projection of any point on ${\bf{P}_{r}}$ to the horizontal axis should correspond to only one point on the vertical axis. Based on the type of landmark to be extracted, the value of $\alpha$ is chosen using trial and error. The minima detector operator of (1) is applied to each curve in Fig. 3-a as follows,

[TABLE]

in which, ${{\bf{S}}^{m}_{\gamma_{m}}}$ represents the $m^{th}$ curve, which is rotated by $\gamma_{m}$ around $\bf{L4^{0}}$. ${{\bf{SMIN}}_{m}}$ represents the location of the global minimum for each curve, which itself, constitutes a curve whose global maximum gives $\bf{L1^{0}}$, the initial location of the nasal root. Figure 3 shows the set of curves ${{\bf{S}}^{m}_{\gamma}}$ (in blue), their minima (in red) and the maximum of the minima (green).

The nasal root and tip locations ($\bf{L1^{0}}$ and $\bf{L4^{0}}$) may be slightly inaccurate due to the depth variations caused by the noise and facial expressions. In order to improve the accuracy of their locations, for the points situated on a 5 $\times$ 5 mm2 area (shown in Fig. 3-b) around the nasal root and saddle, the following angular deviation is calculated,

[TABLE]

in which, $[L1^{i}_{x},L1^{i}_{y}]$ and $[L4^{j}_{x},L4^{j}_{y}]$ are the $[x,y]$ projections of the two pairs of points ${\bf{L1}}^{i}$ and ${\bf{L4}}^{j}$ ($i=1,\ldots,I$, $j=1,\ldots,J$ and $k=1,\ldots,J\times I$), which are selected from the overall $I$ and $J$ points on the region of interest (RoI) from the nasal root and tip regions, respectively, see Fig. 3-b. $\theta_{z}^{k}$ is used to rotate the nose region in the roll direction and around ${\bf{L4}}^{j}$. Then the image is divided into the left and right halves. Assuming the rotated nasal region is translated so that the nose tip ${\bf{L4}}^{j}$ is at the origin, for the $y$-axis indices within the strip ${\bf{S}}^{k}_{y}$ shown in Fig. 3-c (computed using $\theta_{z}^{k}$), the objective function $E^{k}$ is calculated by

[TABLE]

in which ${\bf{Z}}^{k}_{L}$ and ${\bf{Z}}^{k}_{R}$ are the depth maps of the flipped and cropped left and right sides of ${\bf{L4}}^{j}$ and ${\bf{L1}}^{i}$, respectively. The two points ${\bf{L1}}^{opt}$ and ${\bf{L4}}^{opt}$ that minimise $E^{k}$ have the most similar values of ${\bf{Z}}^{opt}_{L}$ and ${\bf{Z}}^{opt}_{R}$ and their projections onto the $x$ axis (${L1^{opt}_{x}}$ and ${L4^{opt}_{x}}$) correspond to the $x$ values of the accurate nose tip and root locations such that,

[TABLE]

This is an example of a ”min-max” optimisation, which finds the best worst case for the optimum [27]. A plane passing through ${\bf{L1}}^{opt}$ and ${\bf{L4}}^{opt}$ is then intersected with the nose surface, with normal vector $[\cos\theta_{z}^{opt},\sin\theta_{z}^{opt}]$, see Fig. 4-b. The locations of the maximum and minimum of the resulting curve are the positions of ${L1^{opt}_{y}}$ and ${L4^{opt}_{y}}$. This procedure is illustrated in Fig. 4. ${\bf{L4}}=[{L4^{opt}_{x}},{L4^{opt}_{y}},{L4^{opt}_{z}}]$ and ${\bf{L1}}=[{L1^{opt}_{x}},{L1^{opt}_{y}},{L1^{opt}_{z}}]$ give the final locations of the nose tip and nasal root, respectively.

The points on the same curve, which are located below the nasal tip $\bf{L4}$ (shown in Fig. 4-c) are then rotated around $\bf{L4}$ by an angle $\phi$. The location of the lowest minimum of the resulting curve ($\bf{S}_{\phi}$) provides the subnasale $\bf{L5}$ after applying (1), i.e. ${\bf{L5}}=V_{1,\phi}({\bf{S}}_{\phi})$. Finally, $\theta_{z}^{opt}$ is used to update the nose region and correct the pose by applying a roll directional rotation around ${\bf{L4}}$.

III-B Nose alar groove and eye corners localisation

The location of the nose tip ($\bf{L4}$) is moved to the origin and an RoI defined to detect the nasal alar grooves (Fig. 5-a) by,

[TABLE]

where $r_{0}$, $a_{1}$ and $a_{2}$ are scalar constants determining the length and directivity of the lobes in the RoI. These are chosen to be able to crop the nasal alar region, while avoiding redundant parts (in subsequent experiments, $r_{0}$ = 30 mm, $a_{1}=4$ and $a_{2}=0.75$). $r$ and $\theta$ are the distance from the nose tip and angular rotation from the horizontal axis passing the nose tip location, respectively. Similarly, the RoI used to detect the eye corners ($\bf{L2}$ and $\bf{L7}$) is depicted in Fig. 5-c and is found using,

[TABLE]

in which $r_{0}^{\prime}=45$ mm, $a_{3}=4$ and $a_{4}=0.75$. The polar coordinate system is characterised by $r^{\prime}$ and $\theta^{\prime}$, which are the distance from the nasal root ($\bf L1$) and the angular rotation from the horizontal axis passing $\bf L1$, respectively.

Planes parallel with the $xz$-plane are then intersected with each row of the RoIs. For the $j^{th}$ and $k^{th}$ intersections over the nasal alar groove and eye corner RoIs, curves ${\bf Q}^{j}$ and ${\bf Q^{\prime}}^{k}$ are found, respectively, and (1) is used to find three minima for the nasal alar groove,

[TABLE]

and the eye corners,

[TABLE]

where $\{{{\bf{RMIN}}_{j}},{{\bf{LMIN}}_{j}}\}$ and $\{{{\bf{RMIN}}_{k}^{\prime}},{{\bf{LMIN}}_{k}^{\prime}}\}$ are matrices with three rows, in which each row has the location of the lowest minimum found from the $j^{th}$ and $k^{th}$ rows of the RoIs, for the right and left sides of $\bf{L4}$ and $\bf{L1}$, respectively. Then, for each row, ${{\bf{RMIN}}_{j}}$ and ${{\bf{LMIN}}_{j}}$ are compared and the pairs with the most similar Euclidean distances to the nose tip $\bf{L4}$ selected. Using a similar approach, the distance of ${{\bf{RMIN}}_{k}^{\prime}}$ and ${{\bf{LMIN}}_{k}^{\prime}}$ to $\bf{L1}$ for each row is computed and those pairs with the most similar distances kept. Figures 5-b and -c illustrate these processes.

The points found as candidates for the nasal alar groove and eye corners might contain some outliers. This is because of the imaging noise and deformations on the face due to the facial expressions. To remove the outliers, an iterative approach is used. First, the 3D Euclidean distances between the points on each consecutive row are computed. Then the standard deviation ($\sigma$) of the resulting vector is used to reject the points whose $\sigma$ is higher than a given threshold (in mm). This process continues until the number of inliers remains unchanged. Compared to the outlier removal method of [4], which uses $K$-means clustering as a criterion to localise the outliers, this approach is deterministic and, unlike $K$-means, is not vulnerable to empty clusters. The outlier removal algorithm results in the green points labelled as the inliers in Fig. 5-b and -d. The left and right pairs, which have the closest value of $y$ to that of the nose tip are selected as $\bf{L3}$ and $\bf{L6}$. Also, the points amongst the inliers in Fig. 5-d, with the smallest depth values, are detected and the pair with the most similar distance to $\bf{L4}$ are selected as the eye corners ($\bf{L2}$ and $\bf{L7}$). The eye corners and nasal alar groove landmarks are the red points in Fig. 5-b and -d, respectively.

IV Feature extraction

The proposed feature space is based on surface normals. For an aligned depth map of the nasal region, represented by its point clouds as ${\bf{N}}=[{\bf{N_{x}}},{\bf{N_{y}}},{\bf{N_{z}}}]$ the normals are $\bf{n}=[{\bf{n_{x}}},{\bf{n_{y}}},{\bf{n_{z}}}]=\nabla N$, where ${\bf{n_{x}}}\circ{\bf{n_{x}}}+{\bf{n_{y}}}\circ{\bf{n_{y}}}+{\bf{n_{z}}}\circ{\bf{n_{z}}}={\bf{1}}$ ($\circ$ and $\bf 1$ represent the Hadamard product operator and matrix of ones, respectively). In order to reduce the sensitivity of the normal vectors to noise and enable the extraction of multi-resolution directional region-based information from the nasal region, instead of calculating the normal vectors directly from the nose surface, they are derived from the Gabor wavelet [28] filtered depth map. The algorithm proposed by Manjunath et al. is used to minimise the wavelets overlap and redundancy in the filtered images [29].

The discrete Fourier transform of the resampled Gabor wavelet ${\bf{G}}_{s,o}$ for the $s^{th}$ scale and $o^{th}$ orientation level ($s=\{1,2,\ldots,s_{m}\}$ and $o=\{1,2,\ldots,o_{m}\}$) is computed and its zero frequency component is set to zero. The Hadamard product of the resulting ${\bf{G}}^{f}_{s,o}$ and the Fourier transform of ${\bf{N_{z}}}$ is then calculated and the absolute value of its inverse Fourier transform is computed for each scale and orientation, i.e. ${\bf N^{f}_{z}}_{s,o}=\left|\mathcal{F}^{-1}\left\{\mathcal{F}\{{\bf N_{z}}\}\circ{\bf{G}}^{f}_{s,o}\right\}\right|$. The maximum of all the corresponding elements of the filtered images is computed over all orientations for each scale $s$: $\left\{{\bf NGm}_{s}|\forall i,j,o:{\bf NGm}_{s}(i,j)\geq{\bf N^{f}_{z}}_{s,o}(i,j)\right\}$. In other words, ${\bf NGm}_{s}=\underset{o}{\text{max}}\left({\bf N^{f}_{z}}_{s,o}\right)$, where $\underset{o}{\text{max}}\left(.\right)$ computes the maximum of the corresponding elements along orientations $o$. Finally, the normal vectors of the resulting per scale maximal map ${\bf NGm}_{s}$ is calculated using the aligned nose coordinate maps ${\bf{N_{x}}}$ and ${\bf{N_{y}}}$,

[TABLE]

where ${\bf{n}}_{s}=[{\bf{n_{x}}}_{s},{\bf{n_{y}}}_{s},{\bf{n_{z}}}_{s}]$ is a block matrix containing the normal vectors for the $s^{th}$ scale level.

V Localised feature descriptors using spherical patches and curves

The feature descriptors are used to define a part of the nasal region, containing a set of normal vectors from the Gabor wavelets filters. Histograms of the resulting feature vectors for the $\bf X$, $\bf Y$ and $\bf Z$ maps are concatenated to create the feature space. This procedure is illustrated in Fig. 6 for $s_{m}=3$ and $o_{m}=4$. The feature descriptors are used to reduce the dimensionality of the feature space, decrease the redundancy and enable the use of probabilistic feature selection to lower the sensitivity to facial expressions while maintaining the most discriminative parts.

The basic landmarks previously identified, see Fig. 1-b, are used to create the new keypoints shown in Fig. 7-a. These new landmarks are easily obtained by dividing the horizontal and vertical lines that connect the landmarks. A sphere centralised on each point is then intersected with the nasal surface and its inner parts are cropped. Then, the histogram of the normals of Gabor-wavelet filtered depth images are computed, based on the procedure explained in section IV. The intersection process is depicted in Fig. 7-b. A set of spheres of identical radii (in this case 7 mm) are intersected with the nose surface. These spherical feature descriptors provide the capability to evaluate the potential of overlapping spherical regions on the nasal surface, when used as feature vectors.

Alternatively, using different pairs of landmarks, a set of orthogonal planes to the nasal region can be found. Intersecting the planes with the nose surface results in a set of curves on the nasal region. For example, the normal vector of a plane passing through two nasal landmarks ${\bf{A_{1}}}$ and ${\bf{A_{2}}}$, and orthogonal to the $xy$ plane can be defined by $\hat{p}_{{\bf{A_{1}}}{\bf{A_{2}}}}=\hat{a}_{z}\times\frac{[{\bf{A_{1}}}-{\bf{A_{2}}}]}{\sqrt{[{\bf{A_{1}}}-{\bf{A_{2}}}]}}$, where $\hat{a}_{z}=[0,0,1]$ is the unit vector along the $z$-axis. When ${\bf{A_{1}}}$ and ${\bf{A_{2}}}$ are selected from the set of landmarks shown in Fig. 7-c, they can be used to create the set of curves shown in Fig. 7-d, which provide the feature descriptors. For each curve, the concatenated histograms of the $x$, $y$ and $z$ components of the normal vectors from the Gabor wavelet filters outputs are computed, giving the feature vector.

VI Feature selection using GA

The feature selection step selects those subsets of feature vectors extracted from the curves and spherical patches that are more robust against facial expressions. For a given feature descriptor and $n$ different Gabor wavelets scales $\{s_{1},s_{2},\ldots,s_{n}\}$, the feature vector is computed by,

[TABLE]

where ${\bf{Fx}}_{s_{k}}$, ${\bf{Fy}}_{s_{k}}$ and ${\bf{Fz}}_{s_{k}}$ are the features of the $s_{k}^{th}$ scale, for the $x$, $y$ and $z$ surface normal components, respectively. For $K$ feature descriptors, each feature set of the normal maps is represented by the concatenation of $K$ different histograms, of length $h_{l}$ from the feature descriptors, giving

[TABLE]

In (12), ${\bf{Hx}}_{i,s_{k}}$, ${\bf{Hy}}_{i,s_{k}}$ and ${\bf{Hz}}_{i,s_{k}}$ are the normalised histograms computed using the $i^{th}$ feature descriptor ($i=1,\ldots,K$) for the $s_{k}^{th}$ scale ($k=1,\ldots,n$) on the normal map ${\bf{n}}_{s_{k}}$, which is computed using (10).

Here the aim is to find a binary vector to be used as a switch to select the most robust and remove the vulnerable feature descriptors to facial expressions. Using a $1\times K$ binary vector ${\bf{Bn}}$, the vector ${\bf{B}}_{s_{k}}$, whose length is equal to the length of ${\bf{Fx}}_{s_{k}}$, ${\bf{Fy}}_{s_{k}}$ or ${\bf{Fz}}_{s_{k}}$ can be computed for the $s_{k}^{th}$ scale by,

[TABLE]

The elements of ${\bf{B}}_{i,s_{k}}$ ($i=1,\ldots,K$) are set to zero or one, depending on the value of the $i^{th}$ element of ${\bf{{Bn}}}$. Finally, ${\bf{B}}_{s_{k}}$ is concatenated over all scales to create a binary vector $\bf B$, whose length is equal to the feature space dimensionality,

[TABLE]

The value of each element of $\bf B$ can be altered using the nucleus binary vector ${\bf{Bn}}$. A curve or patch is selected or omitted based on the value of ${\bf{Bn}}$ elements. If ${\bf{Bn}}(i)=1$, then the $i^{th}$ curve or patch is selected, otherwise it is omitted. By grouping the neutral samples for the gallery and the non-neutral samples for the test phase, and varying ${\bf{Bn}}$, the most expression robust curves and patches can be selected. As shown in (13) and (14) when ${\bf{B}}_{s_{1}}={\bf{B}}_{s_{2}}=\ldots={\bf{B}}_{s_{n}}$, all curves and patches are selected or removed, simultaneously, for all scales. The resulting low dimensional samples are matched with those in the gallery using the Mahalanobis cosine distance,

[TABLE]

The Kernel Fisher’s analysis (KFA) algorithm with polynomial kernel is applied to the feature space to project the features to a lower dimensional space using a supervised approach. If $S_{g}$ and $S_{p}$ are the number of gallery and probe samples, respectively, and $d_{p}$ is the dimension of the projected subspace (in all subsequent experiments, $d_{p}=S_{g}$), ${\bf{X_{g}}}$ and ${\bf{X_{p}}}$ will be matrices of dimensions $S_{g}\times d$ and $S_{p}\times d$, respectively. $\Sigma$ is the $d_{p}\times d_{p}$ covariance matrix computed over ${\bf{X_{g}}}$, and ${\bf{D_{g,p}}}$ is a $S_{g}\times S_{p}$ distance matrix containing the matching errors. To maximise the probability of assigning the test samples ${\bf{X_{t}}}$ to their corresponding classes (subjects), when compared with the gallery samples ${\bf{X_{g}}}$, ${\bf{Bn}}$ can be varied and its optimum found by,

[TABLE]

in which $R_{1}$ is the average probability that the label corresponding to the smallest matching error, found by (15), is the same as the label of the probe sample. In other words, $R_{1}$ is the rank one recognition rate which is maximised as ${\bf{Bn}}$ is changed. The excellent capability of GA in high dimensional binary parameter space [22, 30] make it well suited for this non-convex optimisation problem. The GA used in this work is a modified Non-dominated Sorting Genetic Algorithm-II (NSGA-II) [31] which, in comparison with NSGA [32], is an elitism-based approach, relies on an improved sorting algorithm, has lower computational complexity and does not require sharing parameter assignment [31]. The modified NSGA-II incorporates elitism over the individuals that increase the diversity of the population in addition to those with better fitness output. The parameter assignments for the GA are explained in section VII-C.