Face Boundary Formulation for Harmonic Models: Face Image Resembling

Abstract

This paper is devoted to numerical algorithms based on harmonic transformations with two goals: (1) face boundary formulation by blending techniques based on the known characteristic nodes and (2) some challenging examples of face resembling. The formulation of the face boundary is imperative for face recognition, transformation, and combination. Mapping between the source and target face boundaries with constituent pixels is explored by two approaches: cubic spline interpolation and ordinary differential equation (ODE) using Hermite interpolation. The ODE approach is more flexible and suitable for handling different boundary conditions, such as the clamped and simple support conditions. The intrinsic relations between the cubic spline and ODE methods are explored for different face boundaries, and their combinations are developed. Face combination and resembling are performed by employing blending curves for generating the face boundary, and face images are converted by numerical methods for harmonic models, such as the finite difference method (FDM), the finite element method (FEM) and the finite volume method (FVM) for harmonic models, and the splitting–integrating method (SIM) for the resampling of constituent pixels. For the second goal, the age effects of facial appearance are explored to discover that different ages of face images can be produced by integrating the photos and images of the old and the young. Then, the following challenging task is targeted. Based on the photos and images of parents and their children, can we obtain an integrated image to resemble his/her current image as closely as possible? Amazing examples of face combination and resembling are reported in this paper to give a positive answer. Furthermore, an optimal combination of face images of parents and their children in the least-squares sense is introduced to greatly facilitate face resembling. Face combination and resembling may also be used for plastic surgery, finding missing children, and identifying criminals. The boundary and numerical techniques of face images in this paper can be used not only for pattern recognition but also for face morphing, morphing attack detection (MAD), and computer animation as Sora to greatly enhance further developments in AI.

1. Introduction

1.1. Review of Previous Work

The face is the most prominent feature of a human being. Therefore, studies regarding face fusion and morphing remain active. In particular, face representation has been studied by Gabor wavelets in Du and Ward [1], and face recognition and ear recognition have been also reported by Chang et al. [2]. For face recognition and face morphing, we also cite some recent reports: Aloraibi [3], Indrawal and Sharma [4], Patel and Lapsiwala [5], Scherhag et al. [6,7], Tuncer et al. [8], Guo and Zhang [9], You et al. [10], and Venkatesh et al. [11]. This paper is a continued study of Li, Chiang, and Suen [12] who studied face image transformations using three steps of numerical techniques.

• Step I:

  Find the face boundary.

• Step II:

  Carry out harmonic transformations by FEM, FDM, FVM, etc.

• Step III:

  Transmit facial grayness by splitting algorithms.

We have devoted great efforts to splitting algorithms in Step III (see [13,14]) and harmonic transformations by numerical methods in Step II (see [12]). This paper will focus on Step I, which is to find face boundaries and their conversions during harmonic transformations.

Numerical methods are the main tools to perform face image transformation. For image transformation, new numerical methods have been studied since Li et al. [14], and a number of papers have been published. Among them, the splitting–integrating method (SIM) and harmonic transformation are the most powerful and advanced. The boundary is most important in describing and recognizing the pattern, the individual person, and the transformed models. Harmonic transformations consist of a couple of Laplace’s equations in the face domains and face boundaries. In fact, once the interior and exterior boundaries of face images are known, face domains are defined. Moreover, for each person, the formulation of face boundaries in digital images may be discriminative, complicated, and laborious.

For the boundary of the harmonic model, we solicit the blending curves. Surfaces and curves identified by interpolation are well-known and described in many reports [15,16]. In this paper, new methods for blending curves will be explored for generating face boundaries based on ordinary differential equations (ODE) and their numerical solutions [17]. It is worth pointing out that since the publication of [14], the basic numerical methods for geometric image transformations have been studied and developed, and the face transformations in this paper and [12] are important applications and developments of our basic research. The image examples of face combinations and face resembling in this paper can also find applications in many areas, such as criminal identification, finding missing children, plastic surgery, marriage counseling by predicting possible offspring images, etc. Evidently, the numerical techniques of face boundary in this paper can be used not only for pattern recognition but also for face fusion, face resampling, imaging morphing, morphing attack detection (MAD), and computer animation as Sora. It is also worth pointing out that the numerical algorithms in this paper are beneficial not only to face transformations but also to image geometric transformations.

1.2. Outlines of Numerical Face Transformations

Consider the nonlinear transformation as the harmonic transformation in [14]:

(1)$\begin{array}{ccc}\hfill & & \mathrm{\Omega }\stackrel{T}{\to }S,\partial \mathrm{\Omega }\stackrel{T}{\to }\partial S,\hfill \end{array}$

where $T:(\xi ,\eta )\to (x,y)$, and the harmonic functions

(2)$\begin{array}{ccc}\hfill & & x=x(\xi ,\eta ),y=y(\xi ,\eta )\hfill \end{array}$

satisfy the Laplace equations with the following Dirichlet conditions:

(3)$\begin{array}{cc}\hfill & \Delta x={\displaystyle \frac{{\partial }^{2}x}{\partial {\xi }^2}}+{\displaystyle \frac{{\partial }^{2}x}{\partial {\eta }^2}}{=0,x|}_{\partial \mathrm{\Omega }}=\overline{x},\end{array}$

(4)$\begin{array}{cc}\hfill & \Delta y={\displaystyle \frac{{\partial }^{2}y}{\partial {\xi }^2}}+{\displaystyle \frac{{\partial }^{2}y}{\partial {\eta }^2}}{=0,y|}_{\partial \mathrm{\Omega }}=\overline{y},\end{array}$

and $\overline{x}$ and $\overline{y}$ are given. For Steps I–III, the detailed techniques are described in the following sections.

Step I: Generate blending curves for mapping source and target face boundaries.

Boundary images have been treated by segmentation techniques, where superquadric models are explored, and the simple boundary is studied in [18] by affine kernel transformations. Furthermore, the symmetry of the boundary was explored in [19], and the piecewise linear approximation was used for planar curves for pattern recognition. These studies are confined to some special and simple types of image boundaries. This paper will deal with a rather general boundary of face images.

Let us briefly address the ideas of the blending techniques used in this paper. Based on the scattered 2D points, formulating a smooth curve is one of the basic topics in numerical analysis. Interpolation techniques, such as the Lagrange, Hermite, and spline interpolations, can be used. However, in this paper, we mimic a thin and flexible beam to be blended to pass through these given points exactly. Such a blending beam can be described as an ODE of order four with suitable boundary conditions. We choose the piecewise cubic Hermite polynomials as the admissible functions for the fourth-order ODE. The advantages of this technique are as follows: (1) the global curvature of the blending curve is minimal, and (2) the blending boundary has continuous derivatives (i.e., continuous slopes), thus fitting better and more naturally to the face boundary. The blending techniques are used for shaped evolution with structural and topological changes in [20] to “glue” two or several cutting surfaces together as one uniform subject. In fact, traditional techniques of blending surfaces in [15,21] can also be employed. However, we may also solicit the numerical PDE of biharmonic equations to better “glue” the cutting surfaces. Details are omitted.

First, the following preliminary work is prepared:

• (a). Input a source, which is a 2D face image.

• (b). Choose a target face frame to be transformed (or resembled).

• (c). Locate the control and characteristic points along open or closed curves on the face boundary and the important facial features.

In (c), the characteristic nodes on the face are related not only to the geometrical properties of facial features, such as the eyes, eyebrows, nose, mouth, etc., but also to the muscle structure of the human face. The significant characteristic nodes are provided in [22]. Those nodes may be found using skeleton techniques and by locating the boundary nodes, the sharp nodes of curves, and the nodes of inflection. Since the number of pixels contained in an image boundary is large, it is laborious to manually establish pixel correspondence between two boundaries to be matched, as shown in Figure 1 and Figure 2. New methods are solicited to formulate and map the two boundary curves only by a few matching characteristic points given.

Suppose that the given $(n+1)$ matching pairs of characteristic points $A_i$ and $B_i$ satisfy

$\begin{array}{ccc}\hfill & & A_i\stackrel{T}{\to }{B}_i,i=0,1,\dots ,n,\hfill \end{array}$

where $A_i\in \partial \mathrm{\Omega }$, and $B_i\in \partial S$. In Section 2, we will describe the numerical methods used to generate $\partial \mathrm{\Omega }$, $\partial S$, and their mapping in (1). We will use two methods: (M1) the cubic spline in [21] and (M2) ordinary differential equations (ODE) using Hermite interpolation by following [17]. Method (M2) is more flexible and advanced because it can be applied in parametric forms and is well suited to complicated curves, such as the contours of the head, eyes, eyebrows, nose, and ears. The new techniques of ODE can handle complicated boundary conditions for given slopes and curvatures and allow cubic splines to be embedded in the blending curves.

In [12], the boundary curves of face images of a male (man) and a female (woman) are formulated, as shown in Figure 2, and the face transformation from the female to the male in shape is given in the center of Figure 3. Face boundary formulation and mapping are completed based on the techniques described in this paper.

Step II: Carry out the harmonic model.

Since the image pixels are distributed uniformly, we may simply choose pixel points as the different nodes $(i,j)=({\xi }_i,{\eta }_j)=(iH,jH)$, where H is the pixel size. Denote $x_{i,j}=x(i,j)$ and $y_{i,j}=y(i,j)$. Take the Laplace equation $\Delta x=0$ in (3) as an example. The simple interior difference equations are given by the following:

(5)$\begin{array}{ccc}\hfill & & x_{i+1,j}+x_{i-1,j}+x_{i,j+1}+x_{i,j-1}-4x_{i,j}=0,(i,j)\in \mathrm{\Omega },\hfill \end{array}$

(6)$\begin{array}{ccc}\hfill & & x(i,j)={\overline{x}}_{ij},(i,j)\in \partial \mathrm{\Omega },\hfill \end{array}$

where ${\overline{x}}_{ij}$ is known. Note that the number of (5) is nearly that of the total pixels of the original image used to solve the above difference equations. The successive over-relaxation iteration (SOR) in [23] can be used to solve (5) and (6):

$\begin{array}{ccc}\hfill & & x_{ij}^{(k+1)}=x_{ij}^{(k)}-{\displaystyle \frac{w_{opt}}{4}}\left\{4x_{ij}^{(k)}-(x_{i-1,j}^{(k+1)}+x_{i,j-1}^{(k+1)}+x_{i+1,j}^{(k)}+x_{i,j+1}^{(k)})\right\},\hfill \end{array}$

where $w_{opt}$ is the optimal relaxation parameter, and $x_{ij}^{(0)}$ are the given initial values. For the arbitrary face domain $\mathrm{\Omega }$, although the optimal parameter $w_{opt}$ in [23] is unknown, some new numerical techniques in Chen [24] are provided to seek the optimal $w_{opt}$ easily by trial computation.

To reduce computation complexity, we may choose some pixels $(\overline{i},\overline{j})$ as the difference nodes and establish the non-uniform difference equations. We define $(i,j)=({\xi }_i,{\eta }_j)$, $x_{ij}=x({\xi }_i,{\eta }_j)$, where ${\xi }_i=\overline{i}H$ and ${\eta }_j=\overline{j}H$ are not uniform. The mesh spacings are denoted by $h_i={\xi }_{i+1}-{\xi }_i$ and $k_j={\eta }_{j+1}-{\eta }_j$. Hence, for ${\displaystyle \frac{{\partial }^{2}x}{\partial {\xi }^2}}+{\displaystyle \frac{{\partial }^{2}x}{\partial {\eta }^2}}=f$ in Poisson’s model, the interior difference equations are obtained by direct difference approximations (see Li [25]):

(7)$\begin{array}{ccc}\hfill & & {\displaystyle \frac{2}{h_i+h_{i-1}}}({\displaystyle \frac{x_{i+1,j}-x_{i,j}}{h_i}}-{\displaystyle \frac{x_{i,j}-x_{i-1,j}}{h_{i-1}}})+{\displaystyle \frac{2}{k_j+k_{j-1}}}({\displaystyle \frac{x_{i,j+1}-x_{i,j}}{k_j}}-{\displaystyle \frac{x_{i,j}-x_{i,j-1}}{k_{j-1}}})=f_{ij}.\hfill \end{array}$

We multiply two sides of (7) by the factor $-({\displaystyle \frac{h_i+h_{i-1}}{2}})({\displaystyle \frac{k_j+k_{j-1}}{2}})$ to yield the symmetric difference equations:

$\begin{array}{ccc}& & -({\displaystyle \frac{k_j+k_{j-1}}{2}})({\displaystyle \frac{x_{i+1,j}-x_{i,j}}{h_i}}-{\displaystyle \frac{x_{i,j}-x_{i-1,j}}{h_{i-1}}})-({\displaystyle \frac{h_i+h_{i-1}}{2}})({\displaystyle \frac{x_{i,j+1}-x_{i,j}}{k_j}}-{\displaystyle \frac{x_{i,j}-x_{i,j-1}}{k_{j-1}}})\hfill \\ & & =-({\displaystyle \frac{h_i+h_{i-1}}{2}})({\displaystyle \frac{k_j+k_{j-1}}{2}})f_{ij},(i,j)\in S.\hfill \end{array}$

FDM is simpler than FEM and FVM, and its programming is easy. As a result, FDM is beneficial to the images of medium size converted by harmonic transformations. For example, in Section 4, with near $500\times 500$ pixels, the computation of FDM by the SOR with the optimal parameter $w_{opt}$ can be computed very quickly by using a personal computer.

After the solutions $(x_{i,j},y_{i,j})$ have been obtained from FDM, the piecewise bilinear transformations are given by $x=\widehat{x}(\xi ,\eta )$ and $y=\widehat{y}(\xi ,\eta )$. The functions $(\widehat{x}(\xi ,\eta ),\widehat{y}(\xi ,\eta ))$ are approximated to $(x(\xi ,\eta ),y(\xi ,\eta ))$ of the harmonic transformation with $A_0^i\stackrel{T}{\to }{B}_0^i,i=0,1,2,3$. The functions $\widehat{x}$ and $\widehat{y}$ inside are obtained by the bilinear interpolations in Figure 4 or simply by the linear interpolations in Figure 5. Note that $A_0^i$ should be located at the pixel points in $\xi O\eta$, but it is not necessary for $B_0^i$ in $XOY$.

Step III: Use combinations of the SIM in Li [13] to obtain the distorted and restored images.

It is worth pointing out that the numerical algorithms via Steps I–III in this paper are beneficial not only to face transformations but also to image geometric transformations. We cite the following related reports: Castleman [26], Lakemond et al. [27], Rosa et al. [28], Luhmann et al. [29], Holden [30], Li et al. [31], Ma et al. [32], Gao et al. [33], Fang et al. [34], Chen et al. [35], Tian et al. [36], You et al. [37], and Püspöki [38].

This paper is organized as follows: in Section 2, the blending techniques and their applications are explored for Step I. In Section 3, the numerical methods for Steps II and III are addressed. In Section 4, some image examples of face transformations and resembling are provided. In Section 5, a few concluding remarks are addressed. In Appendix A, a proof of Theorem 1 is presented, and in Appendix B, the combinations of cubic splines and ODE approaches are discussed.

2. Curves for Face Boundaries

For simplicity, we consider the right boundary of a face in $XOY$ in Figure 6. We choose the nose’s central point as the origin of polar coordinates, where $\overrightarrow{z}$ is parallel to direction $\overrightarrow{x}$. The right-half boundary $\partial S$ of the half-face contour can be denoted by

$\begin{array}{ccc}\hfill & & r=r(\theta ),-{\displaystyle \frac{\pi }{2}}\le \theta \le {\displaystyle \frac{\pi }{2}}.\hfill \end{array}$

We assume the ending conditions of $\partial S$ are just horizontal,

(8)$\begin{array}{ccc}\hfill & & r^{\prime }(-{\displaystyle \frac{\pi }{2}})=r^{\prime }({\displaystyle \frac{\pi }{2}})=0.\hfill \end{array}$

Equation (8) makes the curve figuration easier to connect the right half-face contour, see Figure 1.

Suppose that $B_i$ is located in $(x_i,y_i)$ in $XOY$, where $x_i=k_{i}H$, and $y_i=j_{i}H$. The origin of the (Cartesian) coordinates is located at $({\overline{x}}_0,{\overline{y}}_0)$, with ${\overline{x}}_0=\overline{i}H$ and ${\overline{y}}_0=\overline{j}H$, where H is the mesh resolution. Hence, the polar coordinates of point $B_i=(r,\theta )$ are given by

$\begin{array}{ccc}\hfill & & r=r(B_i)=\sqrt{{(x_i-{\overline{x}}_0)}^2+{(y_i-{\overline{y}}_0)}^2}=\sqrt{{(k_i-\overline{i})}^2+{(j_i-\overline{j})}^2}H,\hfill \\ & & \theta =\theta (B_i)=arctan{\displaystyle \frac{y_i-\overline{y}}{x_i-\overline{x}}}=arctan{\displaystyle \frac{j_i-\overline{j}}{k_i-\overline{i}}}.\hfill \end{array}$

For simplicity, we only choose the beginning and end points located at the vertical line: $x(B_0)=x(B_n)={\overline{x}}_0$, see Figure 6. Two interpolation methods are provided below. The cubic splines are well-known, and the ODE using the Hermite interpolation for blending curves is new and developed from [17].

2.1. Cubic Spline Interpolation

Consider the piecewise cubic polynomials $s(x)\in {\mathbf{C}}^2[a,b]$, where ${\mathbf{C}}^k[a,b]$ denotes the set of functions having k-order continuous derivatives. Choose $a=x_0<x_1<\dots <x_n=b$ and $h_i=x_i-x_{i-1}$. Assume that $(x_i,y_i),i=0,1,\dots ,n$ are given and the derivatives $y_0^{\prime }$ and $y_n^{\prime }$ on two boundary nodes are known (e.g., $y_0^{\prime }=y_n^{\prime }=0$ in Figure 7). For the function $s(x)$, denote

$\begin{array}{ccc}& & s(x_i)=y_i,s^{\prime }(x_i)=m_i,s^{″}(x_i)=M_i,i=0,1,2,\dots ,n.\hfill \end{array}$

Then, the piecewise cubic polynomials $s(x)$ are given by (see Su and Liu [21] p. 13)

(9)$\begin{array}{ccc}\hfill & & s(x)={\displaystyle \frac{M_{i-1}}{6h_i}}{(x_i-x)}^3+{\displaystyle \frac{M_i}{6h_i}}{(x-x_{i-1})}^3+\left({\displaystyle \frac{y_{i-1}}{h_i}}-{\displaystyle \frac{h_i{M}_{i-1}}{6}}\right)(x_i-x)\hfill \\ & & +\left({\displaystyle \frac{y_i}{h_i}}-{\displaystyle \frac{h_i{M}_i}{6}}\right)(x-x_{i-1}),\mathrm{in}[x_{i-1},x_i].\hfill \end{array}$

For the given boundary derivatives $m_0$ and $m_n$, the values of $m_i,i=1,\dots ,n-1$ can be obtained from the linear algebraic system,

(10)$\begin{array}{ccc}\hfill & & {\lambda }_i{m}_{i-1}+2m_i+{\mu }_i{m}_{i+1}=C_i,i=1,2,\dots ,n-1,\hfill \end{array}$

where the constants are

$\begin{array}{ccc}\hfill & & {\lambda }_i={\displaystyle \frac{h_{i+1}}{h_i+h_{i+1}}},{\mu }_i={\displaystyle \frac{h_i}{h_i+h_{i+1}}},C_i=3\left({\lambda }_i{\displaystyle \frac{y_i-y_{i-1}}{h_i}}+{\mu }_i{\displaystyle \frac{y_{i+1}-y_i}{h_{i+1}}}\right).\hfill \end{array}$

Equations (10) can be easily solved due to the tri-diagonal coefficient matrix. Then, the values $M_i,i=1,\dots ,n-1$ of second-order derivatives can be found by

$\begin{array}{ccc}\hfill & & M_i=s^{″}(x_i^{-})={\displaystyle \frac{2m_{i-1}}{h_i}}+{\displaystyle \frac{4m_i}{h_i}}-6{\displaystyle \frac{y_i-y_{i-1}}{h_i^2}},\hfill \end{array}$

or

$\begin{array}{ccc}\hfill & & M_i=s^{″}(x_i^{+})=-{\displaystyle \frac{4m_i}{h_{i+1}}}-{\displaystyle \frac{2m_{i+1}}{h_{i+1}}}+6{\displaystyle \frac{y_{i+1}-y_i}{h_{i+1}^2}},\hfill \end{array}$

where $x_i^{-}$ and $x_i^{+}$ are the one-side limit of $x_i$ from the left and right, respectively. Once $m_i$ and $M_i$ are given, the entire function $s(x)$ in (9) is provided.

2.2. Ordinary Differential Equations Using the Hermite Interpolation

For general 2D curves, since suitable Cartesian coordinates (or polar coordinates) cannot be found, the parametric functions are chosen instead (see Figure 8).

$\begin{array}{ccc}\hfill & & x=x(S),y=y(S),0\le S\le 1.\hfill \end{array}$

In [21] (p. 69), the cumulative chord length can be chosen as the following variables:

$\begin{array}{ccc}\hfill & & t_i=\left|{\overline{B_{i}B}}_{i-1}\right|=\sqrt{{(x_i-x_{i-1})}^2+{(y_i-y_{i-1})}^2},1\le i\le n.\hfill \end{array}$

Define

$\begin{array}{ccc}\hfill & & S_i={\displaystyle \frac{{\sum }_{j=1}^i{t}_j}{{\sum }_{j=1}^n{t}_j}};\hfill \end{array}$

then, $0\le S_i\le 1$. We follow [17] to form the Hermite functions satisfying ODE. Consider $x(S)$ and $y(S)$ on $[0,1]$ with the nodes $S_i$, $0=S_0<S_1<\dots <S_n=1$. On $[S_i,S_{i+1}]$, we choose the Hermite interpolations

(11)$\begin{array}{ccc}\hfill & & x_h^i(S)=x_i{\varphi }_0({\displaystyle \frac{S-S_i}{h_{i+1}}})+x_{i+1}{\varphi }_1({\displaystyle \frac{S-S_i}{h_{i+1}}})+h_{i+1}\left[x_i^{\prime }{\psi }_0({\displaystyle \frac{S-S_i}{h_{i+1}}})+x_{i+1}^{\prime }{\psi }_1({\displaystyle \frac{S-S_i}{h_{i+1}}})\right],\hfill \end{array}$

where $h_{i+1}=S_{i+1}-S_i$, $x_i=x(S_i)$, and $x_i^{\prime }=x^{\prime }(S_i)$. The Hermite basis functions are

(12)$\begin{array}{ccc}\hfill & & {\varphi }_0(\overline{S})=2{\overline{S}}^3-3{\overline{S}}^2+1,{\varphi }_1(\overline{S})=-2{\overline{S}}^3+3{\overline{S}}^2,{\psi }_0(\overline{S})={\overline{S}}^3-2{\overline{S}}^2+\overline{S},{\psi }_1(\overline{S})={\overline{S}}^3-{\overline{S}}^2,\hfill \end{array}$

where $\overline{S}\in [0,1]$. Assume that

$\begin{array}{ccc}\hfill & & {\left.{\displaystyle \frac{dy}{dx}}\right|}_{S=0,1}={\displaystyle \frac{y^{\prime }(S)}{x^{\prime }(S)}}=0.\hfill \end{array}$

Then, we have the following boundary conditions:

(13)$\begin{array}{ccc}\hfill & & (1)x_0,x_n,y_0,\mathrm{and}{y}_n\mathrm{are}\mathrm{given},y_0^{\prime }=y^{\prime }(0)=0,y_n^{\prime }=y^{\prime }(1)=0.\hfill \end{array}$

(14)$\begin{array}{ccc}\hfill & & (2)x_0^{\prime }\mathrm{and}{x}_n^{\prime }\mathrm{arbitrary}.\hfill \end{array}$

Equations (13) and (14) enable us to separate the ODE.

The blending curves are modeled as follows. Consider the single equation of the fourth-order ODE

$\begin{array}{ccc}& & {\displaystyle \frac{d^2}{d{S}^2}}\left(p^x(S){\displaystyle \frac{d^2}{d{S}^2}}x(S)\right)=f^x(S),0\le S\le 1,\hfill \end{array}$

where $p^x(S)(>0)$ is the reflection coefficient of the material, and $f^x(S)$ denotes the outward force. The boundary conditions $x_0$ and $x_n$ are given, but $x_0^{\prime }$ and $x_n^{\prime }$ may be zero or arbitrary. Hence, the variational description is employed. Denote $V_h$ as the space of the piecewise Hermite functions (11) satisfying the given boundary conditions. Also, define the energy as follows:

$\begin{array}{ccc}\hfill & & I_x(x(S))={\displaystyle \frac{1}{2}}{\int }_0^1{p}^x(S){\left({\displaystyle \frac{d^{2}x}{d{S}^2}}\right)}^{2}dS-{\int }_0^1{f}^x(S)x(S)dS.\hfill \end{array}$

We seek the minimum solution $x_h$ by

(15)$\begin{array}{ccc}\hfill & & I_x(x_h(S))=\underset{\forall x(S)\in V_h}{min}{I}_x(x(S)).\hfill \end{array}$

For conditions (13) and (14), $x_i$ and $y_i$ are given, but $x_i^{\prime }(i=0,1,\dots ,n)$ and $y_i^{\prime }(i=1,\dots ,n-1)$ are unknown, for which the difference equations will be established. Take the difference equations for $x_i^{\prime }$, for example, since the equations for $y_i^{\prime }$ are similar. The equations for $x_i^{\prime }$ can be obtained from (15) (see [17]) and given as follows. We choose Simpson’s rule to approximate the integration and denote $\widetilde{I_x}(\approx I_x)$. The interior difference equations for $i=1,2,\dots ,n-1$ are given by the following:

(16)$\begin{array}{ccc}\hfill & & 0={\displaystyle \frac{\partial }{\partial x_i^{\prime }}}\widetilde{I_x}(\overrightarrow{X})={\overline{a}}_i^x{x}_{i-1}+{\overline{b}}_i^x{x}_i+{\overline{c}}_i^x{x}_{i+1}+{\overline{d}}_i^x{x}_{i-1}^{\prime }+{\overline{e}}_i^x{x}_i^{\prime }+{\overline{g}}_i^x{x}_{i+1}^{\prime }+{\overline{r}}_i^x,\hfill \end{array}$

where the constants are as follows:

$\begin{array}{ccc}& & {\overline{a}}_i^x={\displaystyle \frac{2}{h_i^2}}\left(2p_i^x+p_{i-1}^x\right),{\overline{b}}_i^x=-\left({\overline{a}}_i^x+{\overline{c}}_i^x\right),\hfill \\ \hfill & & {\overline{c}}_i^x=-{\displaystyle \frac{2}{h_{i+1}^2}}\left(p_{i+1}^x+2p_i^x\right),{\overline{d}}_i^x={\displaystyle \frac{4}{3}}{\displaystyle \frac{1}{h_i}}\left(p_i^x+p_{i-1}^x\right)-{\displaystyle \frac{2}{3}}{\displaystyle \frac{p_{i-\frac{1}{2}}^x}{h_i}},\hfill \\ \hfill & & {\overline{e}}_i^x={\displaystyle \frac{1}{3h_{i+1}}}\left(2p_{i+1}^x+8p_i^x\right)+{\displaystyle \frac{1}{3h_i}}\left(8p_i^x+2p_{i-1}^x\right)+{\displaystyle \frac{2}{3}}\left({\displaystyle \frac{p_{i+\frac{1}{2}}^x}{h_{i+1}}}+{\displaystyle \frac{p_{i-\frac{1}{2}}^x}{h_i}}\right),\hfill \\ \hfill & & {\overline{g}}_i^x={\displaystyle \frac{4}{3}}{\displaystyle \frac{1}{h_{i+1}}}\left(p_{i+1}^x+p_i^x\right)-{\displaystyle \frac{2}{3}}{\displaystyle \frac{p_{i+\frac{1}{2}}^x}{h_{i+1}}},{\overline{r}}_i^x=-{\displaystyle \frac{1}{12}}(f_{i+\frac{1}{2}}^x{h}_{i+1}^2-f_{i-\frac{1}{2}}^x{h}_i^2).\hfill \end{array}$

Next, the boundary difference equation for $x_0^{\prime }$ is given by

(17)$\begin{array}{ccc}\hfill & & 0={\displaystyle \frac{\partial }{\partial x_0^{\prime }}}\widetilde{I_x}(\overrightarrow{X})={\overline{b}}_0^x{x}_0+{\overline{c}}_0^x{x}_1+{\overline{e}}_0^x{x}_0^{\prime }+{\overline{g}}_0^x{x}_1^{\prime }+{\overline{r}}_0^x,\hfill \end{array}$

where the constants are as follows:

$\begin{array}{ccc}\hfill {\overline{b}}_0^x& =& {\displaystyle \frac{1}{h_1^2}}\left(4p_0^x+2p_1^x\right),{\overline{c}}_0^x=-{\overline{b}}_0^x,{\overline{e}}_0^x={\displaystyle \frac{2}{3h_1}}(4p_0^x+p_1^x+p_{\frac{1}{2}}^x),\hfill \\ \hfill {\overline{g}}_0^x& =& {\displaystyle \frac{2}{3h_1}}(2p_0^x+2p_1^x-p_{\frac{1}{2}}^x),{\overline{r}}_0^x=-{\displaystyle \frac{1}{12}}{f}_{\frac{1}{2}}^x{h}_1^2.\hfill \end{array}$

Last, the boundary difference equation for $x_n^{\prime }$ is given by

(18)$\begin{array}{ccc}\hfill & & 0={\displaystyle \frac{\partial }{\partial x_n^{\prime }}}\widetilde{I_x}(\overrightarrow{X})={\overline{a}}_n^x{x}_{n-1}+{\overline{b}}_n^x{x}_n+{\overline{d}}_n^x{x}_{n-1}^{\prime }+{\overline{e}}_n^x{x}_n^{\prime }+{\overline{r}}_n^x,\hfill \end{array}$

where the constants are as follows:

$\begin{array}{ccc}& & {\overline{a}}_n^x={\displaystyle \frac{1}{h_n^2}}\left(2p_{n-1}^x+4p_n^x\right),{\overline{b}}_n^x=-{\overline{a}}_n^x,{\overline{d}}_n^x={\displaystyle \frac{2}{3h_n}}(2p_{n-1}^x+2p_n^x-p_{n-\frac{1}{2}}^x),\hfill \\ \hfill & & {\overline{e}}_n^x={\displaystyle \frac{2}{3h_n}}(p_{n-1}^x+4p_n^x+p_{n-\frac{1}{2}}^x),{\overline{r}}_n^x={\displaystyle \frac{1}{12}}{h}_n^2{f}_{n-\frac{1}{2}}^x.\hfill \end{array}$

For the right half boundary, we may use the ODE techniques as in Figure 9 to seek more intermediate points. Once $x_i$ and $x_i^{\prime }$ are obtained, we find the intermediate points from the Hermite functions in (11). Note that the curves of the Hermite functions by ODE have minimal energy, which grants the curves good smooth properties. For the given slope and curvature at the curve boundary in Figure 10, the numerical techniques are provided by Chen [24].

We provide the following theorem, whose proofs are given in Appendix A.

Theorem 1.

Denote the energy

(19) $\begin{array}{ccc}\hfill & & E^{*}(y(S))=\frac{1}{2}{\int }_0^{1}p(S){\left({\displaystyle \frac{d^{2}y}{d{S}^2}}\right)}^{2}dS-{\int }_0^{1}f(S)y(S)dS.\hfill \end{array}$

Let $y_h(S)$ be the Hermite functions given by ${\displaystyle E^{*}(y_h(S))=\underset{\forall y(S)\in V_h}{min}{E}^{*}(y(S))}$ with $p(S)\equiv 1$ and $f\equiv 0$. Then, the minimal solution of $E^{*}(y)$ gives the exact cubic spline functions satisfying $y\in {\mathbf{C}}^2[0,1]$ and (13).

Furthermore, the new combinations of cubic splines and ODE approaches are explored in Appendix B.

Figure 9

Adding intermediate points.

[Figure omitted. See PDF]

Figure 10

A face contour (left) and nose contour (right).

[Figure omitted. See PDF]

2.3. Intermediate Pixel–Pixel Correspondence Between Two Face Boundaries

To map two face boundaries (1), the intermediate points between $A_i$ and $A_{i+1}$ and $B_i$ and $B_{i+1}$ can also be used (see Figure 11):

$\begin{array}{ccc}& & \xi =\xi (t),\eta =\eta (t),t_i<t_{i+1},A_i\to A_{i+1},\hfill \\ & & x=x(S),y=y(S),S_i<S_{i+1},B_i\to B_{i+1}.\hfill \end{array}$

Consider the function $S=S(t)$, $0<t_i<t_{i+1}<1,0<S_i<S_{i+1}<1$. The simplest function of $S(t)$ is the piecewise linear interpolation:

$\begin{array}{ccc}\hfill & & S=S_i+{\displaystyle \frac{S_{i+1}-S_i}{t_{i+1}-t_i}}(t-t_i).\hfill \end{array}$

However, the spline cubic interpolation and the ODE using the Hermite interpolation provide better intermediate pixels.

3. Splitting–Integrating Method for the Image Grayness of Harmonic Transformations

The splitting–integrating method (SIM) originated in [14] and was then developed in both algorithms and error analysis in [13]. The image grayness ${\mathrm{\Phi }}_{ij}$ and $B_{IJ}$ can be represented by the mean of continuous (or piecewise continuous) intensity functions $\varphi (\xi ,\eta )$ and $b(x,y)$:

(20)$\begin{array}{ccc}\hfill & & {\mathrm{\Phi }}_{ij}={\displaystyle \frac{1}{H^2}}\int {\int }_{{\square }_{ij}}\varphi (\xi ,\eta )d\xi d\eta ,B_{IJ}={\displaystyle \frac{1}{H^2}}\int {\int }_{{\square }_{IJ}}b(x,y)dxdy,\hfill \end{array}$

where $b(x,y)=\varphi (\xi (x,y),\eta (x,y))$, and

(21)$\begin{array}{ccc}\hfill & & {\square }_{ij}=\left\{(\xi ,\eta )|(i-\frac{1}{2})H\le \xi <(i+\frac{1}{2})H,(j-\frac{1}{2})H\le \eta <(j+\frac{1}{2})H\right\},\hfill \end{array}$

(22)$\begin{array}{ccc}& & {\square }_{IJ}=\left\{(x,y)|(I-\frac{1}{2})H\le x<(I+\frac{1}{2})H,(J-\frac{1}{2})H\le y<(J+\frac{1}{2})H\right\}.\hfill \end{array}$

Consider an image with 256 gray levels. Numerical methods are used to evaluate image grayness under harmonic transformation. In [12], the splitting–shooting method (SSM) and the splitting–integrating method (SIM) are used for harmonic transformation T and its inversion transformation $T^{-1}$, respectively. In this paper, we choose the SIM for both T and $T^{-1}$. The SIM is well suited to harmonic transformation since the sequential convergence of $O(N^{-2})$ is higher than $O(N^{-1.5})$ in the SSM (see [13]).

3.1. Splitting–Integrating Method (SIM) for $T^{-1}$

First, consider the inverse transformation $T^{-1}:(x,y)⟶(\xi ,\eta )$. From the known distorted image pixels $\left\{{\widehat{Z}}_{IJ}\right\}$, the approximate functions

$\begin{array}{ccc}& & b(x,y)\approx b_{\mu }(x,y)={\varphi }_{\mu }(\xi (x,y),\eta (x,y)),\hfill \end{array}$

where $\mu$ is the order of grayness interpolations, and $\mu =0$ and $\mu =1$ denote the constant and the linear interpolation, respectively (see [14]). Below, the composite centroid rule will be used to evaluate the integration values in (20). Let ${\square }_{ij}$ in (21) be split into $N\times N$ uniform squares ${\square }_{ij,k\ell }$, i.e., ${\square }_{ij}={\bigcup }_{k,\ell }{\square }_{ij,k\ell }$, where

$\begin{array}{c}\hfill {\square }_{ij,k\ell }=\left\{\begin{array}{cc}(\xi ,\eta ),\hfill & (i-{\displaystyle \frac{1}{2}})H+(k-1)h\le \xi <(i-{\displaystyle \frac{1}{2}})H+kh,\hfill \\ & (j-{\displaystyle \frac{1}{2}})H+(\ell -1)h\le \eta <(j-{\displaystyle \frac{1}{2}})H+\ell h\hfill \end{array}\right\},\end{array}$

where h is the boundary length of ${\square }_{ij,k\ell }$ given by $h={\displaystyle \frac{H}{N}}$. For small subpixels ${\square }_{ij,k\ell }$, the coordinates of the center of gravity are given by

$\begin{array}{ccc}\hfill & & {\xi }_{\dot{G}}={\xi }_{{\dot{G}}_{ij,k\ell }}=(i-\frac{1}{2})H+(k-\frac{1}{2})h,{\eta }_{\dot{G}}={\eta }_{{\dot{G}}_{ij,k\ell }}=(j-\frac{1}{2})H+(\ell -\frac{1}{2})h.\hfill \end{array}$

Based on the composite centroid rule, we have

(23)$\begin{array}{ccc}\hfill & & \int {\int }_{{\square }_{ij,k\ell }}\varphi (\xi ,\eta )d\xi d\eta =\int {\int }_{{\square }_{ij,k\ell }}b(x(\xi ,\eta ),y(\xi ,\eta ))d\xi d\eta \hfill \\ & & \approx h^{2}b\left(x({\xi }_{{\dot{G}}_{ij,k\ell }},{\eta }_{{\dot{G}}_{ij,k\ell }}),y({\xi }_{{\dot{G}}_{ij,k\ell }},{\eta }_{{\dot{G}}_{ij,k\ell }})\right).\hfill \end{array}$

Hence, the normalized image grayness of ${\widehat{W}}_{ij}$ at pixel $(i,j)$ is obtained by

$\begin{array}{ccc}\hfill & & {\mathrm{\Phi }}_{ij}\approx {\displaystyle \frac{h^2}{H^2}}\sum _{k,\ell =1}^N{b}_{\mu }\left(x({\xi }_{{\dot{G}}_{ij,k\ell }},{\eta }_{{\dot{G}}_{ij,k\ell }}),y({\xi }_{{\dot{G}}_{ij,k\ell }},{\eta }_{{\dot{G}}_{ij,k\ell }})\right).\hfill \end{array}$

Note that the computational algorithms (23) do not involve any nonlinear solutions, and the sequential errors as $\mu =1$ are proven to be $O(N^{-2})$ in [13].

3.2. Harmonic Face Image Transformation T

Let ${\square }_{IJ}$ in (22) be split into $N\times N$ uniform squares, i.e., ${\square }_{IJ}={\bigcup }_{k,\ell }{\square }_{IJ,k\ell }$. The pixel grayness via T is also evaluated by the composite centroid rule,

(24)$\begin{array}{ccc}\hfill B_{IJ}& =& {\displaystyle \frac{1}{H^2}}\sum _{k,\ell =1}^N\int {\int }_{{\square }_{IJ,k\ell }}b(x,y)dxdy\hfill \\ \hfill & \approx & {\displaystyle \frac{1}{H^2}}\sum _{k,\ell =1}^N\int {\int }_{{\square }_{IJ,k\ell }}{\varphi }_{\mu }(\xi (x,y),\eta (x,y))dxdy\hfill \\ \hfill & \approx & ({\displaystyle \frac{h}{H}}{)}^2\sum _{k,\ell =1}^N{\varphi }_{\mu }(\xi (x_{\dot{g}},y_{\dot{g}}),\eta (x_{\dot{g}},y_{\dot{g}})):=B_{IJ}^{(N)},\hfill \end{array}$

where $(x_{\dot{g}},y_{\dot{g}})$ is the gravity center of ${\square }_{IJ,k\ell }$. In (24), we do need the values ${\xi }_{\dot{g}}=\xi (x_{\dot{g}},y_{\dot{g}})$ and ${\eta }_{\dot{g}}=\eta (x_{\dot{g}},y_{\dot{g}})$, where $x_{\dot{g}}=x({\xi }_{\dot{g}},{\eta }_{\dot{g}})$ and $y_{\dot{g}}=y({\xi }_{\dot{g}},{\eta }_{\dot{g}})$. In general, to find $({\xi }_{\dot{g}},{\eta }_{\dot{g}})$, we need to solve the two nonlinear equations in (2) for the SIM. However, since the approximate harmonic functions are piecewise linear functions, as in Figure 5, the nonlinear solutions are bypassed. Hence, for the harmonic transformation T as $\mu =1$, we may retain the convergence rate $O(N^{-2})$ by the SIM in [13], higher than $O(N^{-1.5})$ by the SSM in [12].

Under the harmonic models, the approximate values ${\widehat{x}}_{\dot{g}}$ and ${\widehat{y}}_{\dot{g}}$, are obtained from the FDM in $\widehat{x}(\xi ,\eta )\approx x(\xi ,\eta )$ and $\widehat{y}(\xi ,\eta )\approx y(\xi ,\eta )$, see Section 1.2. Suppose that the transformation T is approximated by the piecewise bilinear transformation $\widehat{T}$. Then, when the values $({\widehat{x}}_{ij},{\widehat{y}}_{ij})$($\leftarrow ({\xi }_{ij},{\eta }_{ij}$)) are known, the values of $({\xi }_{\dot{g}},{\eta }_{\dot{g}})$ can be easily found. In fact, the FDM can be regarded as a special kind of finite element method using piecewise linear and bilinear interpolant functions and using special rules of integration approximations (see [25]). Therefore, the FDM solutions $({\widehat{x}}_{ij},{\widehat{y}}_{ij})$ may formulate the piecewise linear transformation $\widehat{T}$:

$\begin{array}{ccc}& & (\xi ,\eta )\to (\widehat{x}(\xi ,\eta ),\widehat{y}(\xi ,\eta )),({\xi }_{ij},{\eta }_{ij})\stackrel{\widehat{T}}{⟶}({\widehat{x}}_{ij},{\widehat{y}}_{ij}).\hfill \end{array}$

Based on the piecewise bilinear interpolation functions, $(\widehat{\xi },\widehat{\eta })$ are formulated from the known $({\widehat{\xi }}_{IJ},{\widehat{\eta }}_{IJ})$ by the following:

$\begin{array}{ccc}& & \xi (x,y)\approx \widehat{\xi }(x,y)={\widehat{\xi }}_{I+1,J+1}{\displaystyle \frac{(x-IH)(y-JH)}{H^2}}+{\widehat{\xi }}_{I,J+1}{\displaystyle \frac{((I+1)H-x)(y-JH)}{H^2}}\hfill \\ & & +{\widehat{\xi }}_{I+1,J}{\displaystyle \frac{(x-IH)(y-(J+1)H)}{H^2}}+{\widehat{\xi }}_{I,J}{\displaystyle \frac{((I+1)H-x)((J+1)H-y)}{H^2}}.\hfill \end{array}$

Then, the redesigned SIM of (24) is given by

(25)$\begin{array}{ccc}& & B_{IJ}\approx {\widehat{B}}_{IJ}^{(N)}={\displaystyle \frac{h^2}{H^2}}\sum _{k,\ell }{\widehat{\varphi }}_{\mu }(\widehat{\xi }(x_{\dot{g}},y_{\dot{g}}),\widehat{\eta }(x_{\dot{g}},y_{\dot{g}}))={\displaystyle \frac{h^2}{H^2}}\sum _{k,\ell }{\widehat{b}}_{\mu }(x_{\dot{g}},y_{\dot{g}}),\hfill \end{array}$

where $\dot{g}={\dot{g}}_{IJ,k\ell }$ and ${\widehat{b}}_{\mu }(x,y)={\widehat{\varphi }}_{\mu }(\widehat{\xi }(x,y),\widehat{\eta }(x,y))$. Note that no nonlinear solutions are needed in (25) either. Here, only one question remains: How do we find $({\widehat{\xi }}_{IJ},{\widehat{\eta }}_{IJ})$ from the FDM solutions $({\widehat{x}}_{ij},{\widehat{y}}_{ij})$? Now, we propose a new technique consisting of four mini steps.

• Step (a):. Compute all ${\widehat{x}}_{ij}$ and ${\widehat{y}}_{ij}$ by the FDM in Step II of Section 2.1,

  $\begin{array}{ccc}& & {\widehat{x}}_{ij}=\widehat{x}({\xi }_i,{\eta }_j)=\widehat{x}(iH,jH),{\widehat{y}}_{ij}=\widehat{y}({\xi }_i,{\eta }_j)=\widehat{y}(iH,jH).\hfill \end{array}$

• Step (b):. Find potential candidates $(I,J)\in {\widehat{\mathrm{\Omega }}}_{IJ}({\widehat{\mathrm{\Omega }}}_{ij}\stackrel{\widehat{T}}{⟵}{\square }_{ij})$, which may be determined by $I_{min}\le I\le I_{max}$ and $J_{min}\le J\le J_{max}$, where

  $\begin{array}{ccc}& & I_{max}=max\left\{\lfloor {\displaystyle \frac{{\widehat{x}}_{ij}}{H}}\rfloor ,\lfloor {\displaystyle \frac{{\widehat{x}}_{i+1,j}}{H}}\rfloor ,\lfloor {\displaystyle \frac{{\widehat{x}}_{i,j+1}}{H}}\rfloor ,\lfloor {\displaystyle \frac{{\widehat{x}}_{i+1,j+1}}{H}}\rfloor \right\},\hfill \\ & & J_{max}=max\left\{\lfloor {\displaystyle \frac{{\widehat{y}}_{ij}}{H}}\rfloor ,\lfloor {\displaystyle \frac{{\widehat{y}}_{i+1,j}}{H}}\rfloor ,\lfloor {\displaystyle \frac{{\widehat{y}}_{i,j+1}}{H}}\rfloor ,\lfloor {\displaystyle \frac{{\widehat{y}}_{i+1,j+1}}{H}}\rfloor \right\}.\hfill \end{array}$

  Here, $\lfloor x\rfloor$ is the largest integer $\le x$. The definitions of $I_{min}$ and $J_{min}$ are similar.

• Step (c):. Split ${\square }_{ij}$ into two triangles: ${▵}_{ij,1}$ and ${▵}_{ij,2}$ in $XOY$ in Figure 12, where ${\widehat{▵}}_{ij,\ell }$ are also triangles, and find all possible $(I,J)$ in $XOY$ such that

  (26)$\begin{array}{ccc}\hfill & & (IH,JH)\in {\widehat{▵}}_{ij,\ell }.\hfill \end{array}$

  The ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are obtained from the following equations:

  $\begin{array}{ccc}\hfill & & \left\{\begin{array}{c}{x}_I={\lambda }_1{x}_a+{\lambda }_2{x}_b+{\lambda }_3{x}_c,\hfill \\ y_J={\lambda }_1{y}_a+{\lambda }_2{y}_b+{\lambda }_3{y}_c,\hfill \\ 1={\lambda }_1+{\lambda }_2+{\lambda }_3,\hfill \end{array}\right.\hfill \end{array}$

  where $a,b,c$ are the vertices of triangles ${\widehat{▵}}_{ij,\ell },\ell =1,2$, and $(x_a,y_a)$ are the coordinates of the vertices in Figure 12. The sufficient and necessary conditions for (26) are $0\le {\lambda }_i\le 1,i=1,2,3$.

• Step (d):. Obtain the approximate values of $({\xi }_{IJ},{\eta }_{IJ})$ related to $(IH,JH)$ in (26) by

  $\begin{array}{ccc}& & {\widehat{\xi }}_{IJ}={\lambda }_1{\xi }_{\overline{a}}+{\lambda }_2{\xi }_{\overline{b}}+{\lambda }_3{\xi }_{\overline{c}},{\widehat{\eta }}_{IJ}={\lambda }_1{\eta }_{\overline{a}}+{\lambda }_2{\eta }_{\overline{b}}+{\lambda }_3{\eta }_{\overline{c}},\hfill \end{array}$

  where $\overline{a}$, $\overline{b}$, and $\overline{c}$ are the vertices of ${▵}_{ij,\ell }$, as shown in Figure 12.

Note that Steps (a)–(d) above are easy to perform. The computation complexity is also only $O(M^2)$, where $M^2$ is the total pixel number. More importantly, the computations for faces resembling under T and $T^{-1}$ do not involve any nonlinear solutions.

4. Image Experiments of Age Effects and Face Resembling

The second goal of this paper is to develop and extend the face transformations in [12]. In Figure 3, we converted the female image using the face boundary of the male. By harmonic transformations, the numerical algorithms are given in Steps I–III in Section 1.1, where the finite volume method (FVM) with Delaunay triangulation was used. An image of a teenage girl was produced and is shown in the center of Figure 3. She looks pretty, but she is virtual. No such girl exists in the world.

Encouraged by Figure 3, we propose a new challenging task. Based on the photos of the parents at 50 years old and the child at five years old at the top of Figure 13, we will seek the integrated image to resemble the current image of her at 20 years old at the bottom of Figure 13 as closely as possible. This face resembling is more challenging than the face transformation in [12]. Furthermore, age effects can also be found in Palumbo et al. [39] and Taskiran et al. [40]. Now, we report two image experiments.

• Experiment I: Age Effects of Face Appearance. First, we choose the images of the parents only. The images of the father and mother are converted to the same frame of the young girl, and the combined images with different ratios are provided in Figure 14, where “father:mother = 0.1:0.9” denotes the proportions of grayness values given by $10\%$ of the father’s image and $90\%$ of the mother’s image. Disappointedly, the combined images look like a 50-year-old, not like a 20-year-old. Second, we choose images of the mother and child; the combined images are provided in Figure 15. Surprisingly, the images with different ages between teenagers and adults at 50 years old can be observed clearly.

• Experiment II: Face Resembling. Based on Experiment I, we add the child’s image and choose three bases of images at the top of Figure 13 to resemble the young girl’s image. Since the young girl looks more like her mother, we may assume that the proportion (i.e., contribution) of her father is less, say about $10\%$ and $20\%$. Then, we may select the resembling images among many possible combinations, which are provided in Figure 16 with their grayness proportion values. However, this selection requires a great deal of computational work. Also, the resemblance selection purely by our eyes may not be trustworthy and reasonable. New numerical techniques are required to choose optimal images automatically.

When three face images of father, mother, and child have been reshaped to the girl’s face frame (see the bottom of Figure 13), we denote their pixel grayness as $u_{ij}$, $v_{ij}$, and $w_{ij}$ at $(i,j)$, respectively. Also, the pixel grayness of the young girl is denoted by ${\varphi }_{ij}$. The sets of $\left\{u_{ij}\right\}$, $\left\{v_{ij}\right\}$, and $\left\{w_{ij}\right\}$ are regarded as three bases, and a linear combination of them is given by $\left\{z_{ij}\right\}=\alpha \left\{u_{ij}\right\}+\beta \left\{v_{ij}\right\}+\gamma \left\{w_{ij}\right\}$, where $\alpha$, $\beta$, and $\gamma$ are positive proportions to be sought. The optimal proportions $\alpha$, $\beta$, and $\gamma$ should be chosen to achieve $\left\{{\varphi }_{ij}\right\}\approx \left\{z_{ij}\right\}$ as best as possible. To this end, the least squares method (LSM) may be employed. The LSM is a popular technique for geometric transformations.

Define the errors

$\begin{array}{ccc}& & E=E(\alpha ,\beta ,\gamma )=\sum _{ij}{\{{\varphi }_{ij}-z_{ij}\}}^2=\sum _{ij}{\{{\varphi }_{ij}-(\alpha u_{ij}+\beta v_{ij}+\gamma w_{ij})\}}^2.\hfill \end{array}$

The optimal $\alpha ,\beta ,\gamma$, satisfy the extreme conditions: ${\displaystyle \frac{\partial E}{\partial \alpha }}=0$, ${\displaystyle \frac{\partial E}{\partial \beta }}=0$, and ${\displaystyle \frac{\partial E}{\partial \gamma }}=0$. This gives a system of three linear algebraic equations:

(27)$\begin{array}{ccc}\hfill & & a_{11}\alpha +a_{12}\beta +a_{13}\gamma =b_1,a_{21}\alpha +a_{22}\beta +a_{23}\gamma =b_2,a_{31}\alpha +a_{32}\beta +a_{33}\gamma =b_3,\hfill \end{array}$

where the coefficients $a_{ij}=a_{ji}$, and the explicit constants are given by

$\begin{array}{ccc}& & a_{11}=\sum _{ij}{u}_{ij}^2,a_{22}=\sum _{ij}{v}_{ij}^2,a_{33}=\sum _{ij}{w}_{ij}^2,\hfill \\ & & a_{12}=\sum _{ij}{u}_{ij}{v}_{ij},a_{13}=\sum _{ij}{u}_{ij}{w}_{ij},a_{23}=\sum _{ij}{v}_{ij}{w}_{ij},\hfill \\ & & b_1=\sum _{ij}{u}_{ij}{\varphi }_{ij},b_2=\sum _{ij}{v}_{ij}{\varphi }_{ij},b_3=\sum _{ij}{w}_{ij}{\varphi }_{ij}.\hfill \end{array}$

Since the total number of $(i,j)$ is much larger than three, the linear Equation (27) are linearly independent. Hence, the values of $\alpha$, $\beta$, and $\gamma$ can be easily obtained and given by $\alpha =0.052$, $\beta =0.619$, and $\gamma =0.301$, as shown in Figure 17, with different hairstyles. Amazingly, two artificial integrated images with such optimal proportions reach a good resemblance to the girl’s original image in the center of Figure 18.

To close this section, we introduce two remarks.

Remark 1.

Face resembling in this paper is critical to finding missing children and identifying criminals via the police. Note that the face images in Figure 13 are standard photos required in a passport: shot from the front with no smile and a closed mouth. Suppose a target and potential patterns are given, where the patterns are of standard face images, but the face target suffers from geometric transformations and illumination effects since it may be taken from a cell phone by chance. The geometric transformations include basic geometric transformations (such as translation, rotation, and scaling) and prospective transformation. Our main efforts are paid to seek these geometric transformations and illumination effects. Then, the restored target can be found by the numerical algorithms in this paper, and the identification of the restored target with one of the standard patterns may be accomplished using the least squares method (LSM). Details are reported in another paper.

Remark 2.

In this remark, let us address the new techniques of face resembling proposed in this paper. The algorithms in Steps I–III in Section 1 are, indeed, of numerical interpolation, numerical integration, numerical ordinary differential equations (ODE), and numerical partial differential equations (PDE). These numerical algorithms were developed from our previous study [13,14,25] for geometric image transformations. From Figure 3, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18, the virtual face images were created by merging two (or more) face images; they belong to face morphing. Since numerical algorithms are used, numerical face resembling (or morphing) may be called. The numerical algorithms in this paper and [12] are different from those in the face recognition method in [4,8,9,10], where convolutional neural networks (CNN), textual transformations, etc., are used. Our key numerical algorithms are also different from those used in face morphing [3,5,6,11], where deep neural networks, feature comparisons, cross-dissolve techniques, etc., are used. The numerical interpolation in this paper and the Delaunay triangulation in [12] are found in [6]. Furthermore, the rapping techniques in [6] and the facial landmark techniques in [11] are also related to geometrical transformations. The advanced numerical algorithms in this paper and [12] may be used to improve the quality of face morphing. On the other hand, their effective techniques may also be employed for better identification, as in Figure 13 and Figure 16. When face morphing is used by an attacker in malice, face verification is vulnerable. Face-morphing attacks have appeared since 2014. Our numerical techniques can also be used for morphing attack detection (MAD). To describe new numerical algorithms, some computational formulas are essential. To display the efficiency of the algorithms, not only are some face images with grayness errors provided, but also the error analysis of algorithms as in [13] is needed. Nowadays, artificial intelligence (AI) is the hottest topic in research and applications since ChatGPT and Sora appeared. AI includes three basic elements: data, models, and algorithms. In face resembling (or morphing), geometrical and harmonic transformations are the models. We do not deal with massive data but seek good numerical algorithms. Hence, the CPU time can be saved greatly. In summary, the new and advanced numerical algorithms of this paper and [12] can also be applied to Sora to greatly enhance further developments of AI.

5. Concluding Remarks

To conclude, let us address the novelties in this paper.

1. The face transformation achieved by harmonic models is shown in three steps in Section 1; this paper focuses on Step I to generate the face boundary, which is needed in the Dirichlet boundary condition for Laplace’s equations in harmonic models.

2. We provide two methods for blending boundary curves: the cubic spline and the ODE using the Hermite interpolation. The latter is new and developed from our previous study [17]. For the interior and exterior boundary of the face contour, the curves are piecewise smooth only with second-order derivatives. Under some constraints of characteristic points, the remaining displacements and derivatives of the curves can be determined uniquely by minimizing the energy (15).

3. The obtained curves have small curvatures, and they are smooth (or piecewise smooth). Moreover, the popular cubic spline is proved in Appendix A to be a special case of the ODE using the Hermite interpolation. A combined technique in Appendix B is then explored to couple the cubic spline and the ODE using the Hermite interpolation.

4. By the SIM, for the general nonlinear functions $x(\xi ,\eta )$ and $y(\xi ,\eta )$ in (2), their solutions are needed. However, the approximate harmonic functions are piecewise linear so that the nonlinear solutions can be bypassed. Hence, the algorithms of the SIM with $\mu =1$ for harmonic transformations are simple, while the convergence rate $O(N^{-2})$ may remain higher than $O(N^{-1.5})$ by the SSM for T used.

5. Two image experiments are provided in Section 4. First, different ages of facial appearance can be produced by the images of adults and children. Then, a young boy’s/girl’s image can be resembled from his/her parents and those from their childhood.

6. Inspired by two successful experiments, many applications of face transformations may follow: modifying facial appearance, plastic surgery, generating possible offspring images for marriage counseling, and identifying a person for different purposes, such as finding missing children and confirming criminals.

7. The boundary techniques and the face combinations used in this paper can be applied not only to image processing and pattern recognition but also to face fusion and morphing, morphing attack detection (MAD), and computer animation, such as Sora, to greatly enhance further developments in AI. See Remarks 1 and 2.

Author Contributions

Methodology, Y.W.; Investigation, C.Y.S.; Writing—original draft, Z.-C.L.; Writing—review & editing, H.-T.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding authors.

Acknowledgments

We thank C.C. Chen for the image experiments in Section 4, and we are also grateful to our reviewers for their valuable comments and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Face transformation from a source face image to a target one.

Figure 2. The corresponding exterior and boundary curves of face images of the male (man) and female (woman), where the left shows the blending curves with characteristic nodes, and the right shows the blending curves (with the permission from Z.C. Li).

Figure 3. The face transformation from female to male, where the center one is the fusion image (with the permission from Z.C. Li).

Figure 4. Piecewise bilinear transformations.

Figure 5. Piecewise bilinear transformations on a square.

Figure 6. Right half of face contour.

Figure 7. Cubic spline interpolation.

Figure 8. General 2D curves.

Figure 11. Intermediate points.

Figure 12. The inverse transformation at [Forumla omitted. See PDF.] in [Forumla omitted. See PDF.].

Figure 13. The pictures of a young girl (original girl), her father (Father), her mother (Mother), and her childhood (Child) (with the permission from C.C. Chen).

Figure 14. Age effects of parents’ images (with the permission from C.C. Chen).

Figure 15. Age effects of the mother and child (with the permission from C.C. Chen).

Figure 16. The resembling images chosen by us (with the permission from C.C. Chen).

Figure 17. Images of father, mother, and child at the top, reshaped to the girl’s image at the bottom (with the permission from C.C. Chen).

Figure 18. The optimal resembling images obtained by the least squares method (with the permission from C.C. Chen).

Appendix A. Proofs for Theorem 1

In this first appendix, we will provide the direct but heuristic proofs of Theorem 1, although other proofs can be found based on Vandergraft [41] (Theorem 4.2, p. 136). We will derive the explicit equations for the Hermite function just as those of the cubic spline functions, by which we will explore the combined algorithms of the two methods given in Appendix B. First, let us derive the cubic splines from the ODE using the Hermite interpolation. For simplicity, we use the simple function $y=f(x)$ to denote the curve since the parametric functions $(x(S),y(S))$ can be proved similarly.

Let $\overline{S}={\displaystyle \frac{x-x_{i-1}}{h_i}}$. From [21] (p. 12), we obtain the ODE Hermite interpolation from (11) and (12), $$\begin{array}{ccc}& & y(x)=m_{i-1}{\displaystyle \frac{{(x_i-x)}^2(x-x_{i-1})}{h_i^2}}-m_i{\displaystyle \frac{{(x-x_{i-1})}^2(x_i-x)}{h_i^2}}\hfill \\ \hfill & & +y_{i-1}{\displaystyle \frac{{(x_i-x)}^2[2(x-x_{i-1})+h_i]}{h_i^3}}+y_i{\displaystyle \frac{{(x-x_{i-1})}^2[2(x_i-x)+h_i]}{h_i^3}}.\hfill \end{array}$$Then, we have the first- and second-order derivatives, $$\begin{array}{ccc}& & y^{\prime }(x)=m_{i-1}{\displaystyle \frac{(x_i-x)(2x_{i-1}+x_i-3x)}{h_i^2}}\hfill \\ \hfill & & -m_i{\displaystyle \frac{(x-x_{i-1})(2x_i+x_{i-1}-3x)}{h_i^2}}+6{\displaystyle \frac{y_i-y_{i-1}}{h_i^3}}(x_i-x)(x-x_{i-1}),\hfill \end{array}$$ (A1)$$\begin{array}{ccc}\hfill & & y^{″}(x)=2m_{i-1}{\displaystyle \frac{3x-2x_i-x_{i-1}}{h_i^2}}+2m_i{\displaystyle \frac{3x-x_i-2x_{i-1}}{h_i^2}}+6{\displaystyle \frac{y_i-y_{i-1}}{h_i^3}}(x_i+x_{i-1}-2x).\hfill \end{array}$$The two-sided limits of the second-order derivatives at $x_i$ are given by (A2)$$\begin{array}{ccc}\hfill & & y^{″}(x_i^{-})={\displaystyle \frac{2m_{i-1}}{h_i}}+{\displaystyle \frac{4m_i}{h_i}}-6{\displaystyle \frac{y_i-y_{i-1}}{h_i^2}},\hfill \end{array}$$ (A3)$$\begin{array}{ccc}\hfill & & y^{″}(x_i^{+})=-{\displaystyle \frac{4m_i}{h_{i+1}}}-{\displaystyle \frac{2m_{i+1}}{h_{i+1}}}+6{\displaystyle \frac{y_{i+1}-y_i}{h_{i+1}^2}}.\hfill \end{array}$$Moreover, the continuity of the second-order derivatives at an internal knot, i.e., $y^{″}(x_i^{-})=y^{″}(x_i^{+})$, yields (A4)$$\begin{array}{ccc}\hfill & & {\lambda }_i{m}_{i-1}+2m_i+{\mu }_i{m}_{i+1}=c_i,i=1,2,\dots ,n-1,\hfill \end{array}$$ where (A5)$$\begin{array}{ccc}\hfill & & {\lambda }_i={\displaystyle \frac{h_{i+1}}{h_i+h_{i+1}}},{\mu }_i={\displaystyle \frac{h_i}{h_i+h_{i+1}}},c_i={\displaystyle \frac{3}{h_i+h_{i+1}}}\left(h_i{\displaystyle \frac{y_{i+1}-y_i}{h_{i+1}}}+h_{i+1}{\displaystyle \frac{y_i-y_{i-1}}{h_i}}\right).\hfill \end{array}$$Equations (A4) and (A5) are the exact cubic interpolation (10). This implies that by adding the continuity of second-order derivatives at the nodes, the ODE Hermite solutions lead to the cubic spline interpolant.

On the other hand, we may also explicitly derive the equations for the minimum energy $E^{*}(y(S))$ in (19) to provide it exactly (A4). Since $p(S)\equiv 1$ and $f\equiv 0$, we consider the energy $E(y)={\sum }_{i=1}^n{E}_i(y)$, where $$\begin{array}{ccc}& & E_i(y)={\displaystyle \frac{1}{2}}{\int }_{x_{i-1}}^{x_i}{(y^{″}(x))}^{2}dx.\hfill \end{array}$$The derivatives $y^{″}(x)$ of order two in (A1) are simplified as $$\begin{array}{ccc}& & y^{″}(x)=g_0(x)m_{i-1}+g_1(x)m_i+g_2(x)(y_i-y_{i-1}),\hfill \end{array}$$ with the known functions $$\begin{array}{ccc}& & g_0(x)={\displaystyle \frac{2(3x-2x_i-x_{i-1})}{h_i^2}},g_1(x)={\displaystyle \frac{2(3x-x_i-2x_{i-1})}{h_i^2}},g_2(x)={\displaystyle \frac{6}{h_i^3}}(x_i+x_{i-1}-2x).\hfill \end{array}$$Therefore, we have (A6)$$\begin{array}{ccc}\hfill & & E_i(y)={\displaystyle \frac{1}{2}}(m_{i-1},m_i)\left[\begin{array}{cc}{a}_{11}\hfill & a_{12}\hfill \\ a_{21}\hfill & a_{22}\hfill \end{array}\right]\left(\begin{array}{c}{m}_{i-1}\hfill \\ m_i\hfill \end{array}\right)+(m_{i-1},m_i)\left(\begin{array}{c}{b}_1\hfill \\ b_2\hfill \end{array}\right)+c_0,\hfill \end{array}$$ where $c_0$ is a constant, and (A7)$$\begin{array}{ccc}\hfill & & a_{11}={\int }_{x_{i-1}}^{x_i}{(g_0(x))}^{2}dx,a_{22}={\int }_{x_{i-1}}^{x_i}{(g_1(x))}^{2}dx,a_{12}=a_{21}={\int }_{x_{i-1}}^{x_i}{g}_0(x)g_1(x)dx,\hfill \\ & & b_1=(y_i-y_{i-1}){\int }_{x_{i-1}}^{x_i}{g}_0(x)g_2(x)dx,b_2=(y_i-y_{i-1}){\int }_{x_{i-1}}^{x_i}{g}_1(x)g_2(x)dx.\hfill \end{array}$$

Next, we evaluate the integrals (A7) by calculus. By using the transformation $$\begin{array}{ccc}& & t={\displaystyle \frac{x-x_{i-1}}{h_i}},h_i=x_i-x_{i-1},dx=h_{i}dt,\hfill \end{array}$$ the explicit entry $a_{11}$ is obtained by $$\begin{array}{ccc}& & a_{11}={\int }_{x_{i-1}}^{x_i}{\left[{\displaystyle \frac{2(3x-2x_i-x_{i-1})}{h_i^2}}\right]}^{2}dx={\displaystyle \frac{4}{h_i^4}}{\int }_0^1{\left[h_i(3t-2)\right]}^2{h}_{i}dt={\displaystyle \frac{4}{h_i}}.\hfill \end{array}$$Also, we obtain other explicit entries similarly, $$\begin{array}{ccc}& & a_{22}={\displaystyle \frac{4}{h_i^4}}{\int }_0^1{(h_i(3t-1))}^2{h}_{i}dt={\displaystyle \frac{4}{h_i}},a_{12}={\displaystyle \frac{4}{h_i}}{\int }_0^1(3t-1)(3t-2)dt={\displaystyle \frac{2}{h_i}},\hfill \\ & & b_1={\displaystyle \frac{12}{h_i^2}}(y_i-y_{i-1}){\int }_0^1(3t-2)(1-2t)dt=-{\displaystyle \frac{6}{h_i^2}}(y_i-y_{i-1}),\hfill \\ & & b_2={\displaystyle \frac{12}{h_i^2}}(y_i-y_{i-1}){\int }_0^1(3t-1)(1-2t)dt=-{\displaystyle \frac{6}{h_i^2}}(y_i-y_{i-1}).\hfill \end{array}$$Hence, the matrix and vector in (A6) are expressed explicitly, (A8)$$\begin{array}{ccc}\hfill & & \left[\begin{array}{cc}{a}_{11}\hfill & a_{12}\hfill \\ a_{21}\hfill & a_{22}\hfill \end{array}\right]={\displaystyle \frac{2}{h_i}}\left[\begin{array}{cc}2\hfill & 1\hfill \\ 1\hfill & 2\hfill \end{array}\right],\left[\begin{array}{c}{b}_1\hfill \\ b_2\hfill \end{array}\right]=-{\displaystyle \frac{6(y_i-y_{i-1})}{h_i^2}}\left[\begin{array}{c}1\hfill \\ 1\hfill \end{array}\right].\hfill \end{array}$$From the derivatives $$\begin{array}{ccc}& & {\displaystyle \frac{\partial }{\partial x_i}}(E(y))={\displaystyle \frac{\partial }{\partial x_i}}\sum _{i=1}^n{E}_i(y)=0,\hfill \end{array}$$ we have the algebraic equations (A9)$$\begin{array}{ccc}\hfill & & \mathbf{A}\mathbf{x}=\mathbf{b}.\hfill \end{array}$$From (A6) and (A8), the interior equations are obtained, $$\begin{array}{ccc}& & {\displaystyle \frac{2}{h_i}}[m_{i-1}+2m_i]-{\displaystyle \frac{6}{h_i^2}}(y_i-y_{i-1})+{\displaystyle \frac{2}{h_{i+1}}}(2m_i+m_{i+1})-{\displaystyle \frac{6}{h_{i+1}^2}}(y_{i+1}-y_i)=0.\hfill \end{array}$$This is (A4) and completes the proof of Theorem 1. □

Appendix B. Combinations of Cubic Splines and ODE Approaches

In the second Appendix, new combinations of cubic splines and ODE solutions are explored, as in Figure A1. For simplicity, we only consider $y=y(x)$ since the case of parametric functions can be extended easily. Define an energy $$\begin{array}{ccc}& & E(y)={\displaystyle \frac{1}{2}}{\int }_a^{b}p(x){(y^{″}(x))}^{2}dx-{\int }_a^{b}f(x)y(x)dx,\hfill \end{array}$$ where $p(x)=1$ in $[a,c]$, $f(x)=0$ in $[a,c]$, with $a<c<b$. From Theorem 1 above, the minimum of $E^{*}(y)$ is exactly the cubic spline, where $$\begin{array}{ccc}\hfill & & E^{*}(y)={\displaystyle \frac{1}{2}}{\int }_a^c{(y^{\prime \prime }(x))}^{2}dx.\hfill \end{array}$$We conclude that the cubic spline is just a special case when $p(x)\equiv 1$ and $f(x)\equiv 0$. Now we may carry out a combination as in Figure A1, where point $c\in (a,b)$ is the interior point such that $y^{″}(x)\in {\mathbf{C}}^2(a,c)\wedge y^{\prime }(x)\in {\mathbf{C}}^1[a,b]$.

Figure A1. A combination of cubic splines and ODE solutions.

Now, we give the algorithms for a combination in $[a,b]$ in Figure A1 to couple the cubic spline interpolant in $(a,c)$ and the ODE approaches with the Hermite interpolant in $(c,b)$. Assume $(x_i,y_i),i=0,1,\dots ,n$, $x_k=x_c$ are known, but $y_i^{\prime }$ are determined in the following three cases:

• Case A:

  The equations $m_i(i=0,1,\dots ,k-1)$ are obtained from the cubic splines in Section 2.1.

• Case B:

  The equations $y_i^{\prime }=m_i(i=k+1,\dots ,n)$ are obtained from (16) in Section 2.2.

• Case C:

  For the equation $m_k$ at the intersection point c, there occurs the continuity of second-order derivatives. Comparing the right boundary equations in (A3) with the left boundary Equation (17), we obtain $$\begin{array}{ccc}\hfill & & {\displaystyle \frac{2}{h_k}}{m}_{k-1}+({\displaystyle \frac{4}{h_k}}+{\overline{e}}_0)m_k+{\overline{g}}_0{m}_{k+1}+{\overline{r}}_0-{\displaystyle \frac{6(y_k-y_{k-1})}{h_k}}=0,\hfill \end{array}$$ where $$\begin{array}{ccc}& & {\overline{e}}_0={\displaystyle \frac{2}{3h_{k+1}}}(4p_k^{+}+p_{k+1}^{+}+p_{k+{\displaystyle \frac{1}{2}}}),\hfill \\ \hfill & & {\overline{g}}_0={\displaystyle \frac{2}{3h_{k+1}}}(2p_k+2p_{k+1}-p_{k+{\displaystyle \frac{1}{2}}}),{\overline{r}}_0=-{\displaystyle \frac{1}{12}}{f}_{k+{\displaystyle \frac{1}{2}}}{h}_{k+1}^2.\hfill \end{array}$$

Note that by using this combination, the coefficient matrix $\mathbf{A}$ in (A9) is still sparse, positive definite, and symmetric. Therefore, the combined solutions can be easily obtained. It is worth pointing out that the equations in Case A are obtained from (16) in cases where $p_k^{-}=1,i\le k$ and $f_k^{-}=0,i\le k$.

To close this appendix, let us consider the simple support boundary condition. Suppose that the homogeneous simple support boundary condition is given by $y_0^{″}=\alpha \ne 0$, and that $y_0$, $y_k$ and $y_k^{\prime }$ are known. We may modify the energy by $$\begin{array}{ccc}& & E_1(y)={\displaystyle \frac{1}{2}}{\int }_a^c{(y^{″}(x))}^{2}dx+{\displaystyle \frac{\alpha y^{\prime }(a)}{2}}\hfill \end{array}$$ or (A10)$$\begin{array}{ccc}\hfill & & {\int }_a^c{y}^{″}(x)w^{″}(x)dx+\alpha ·w^{\prime }(a)=0,\hfill \end{array}$$ where $w(x)$ are the admissible functions satisfying $w(a)=w(c)=0$ and $w^{\prime }(c)=0$. From integration by parts, we obtain from (A10) $$\begin{array}{ccc}& & {\int }_a^c{y}^{(4)}wdx-(y^{″}(a)-\alpha )w^{\prime }(a)=0.\hfill \end{array}$$Since $w^{\prime }$ is arbitrary at $x=a$, we conclude $y^{″}(a)=\alpha$. From (A3), we have the left boundary equation $$\begin{array}{ccc}& & {\displaystyle \frac{2}{h_1}}[m_1+2m_0]-{\displaystyle \frac{6}{h_1^2}}(y_1-y_0)+\alpha =0,\hfill \end{array}$$ i.e., $$\begin{array}{ccc}& & 2m_0+m_1=3\left({\displaystyle \frac{y_1-y_0}{h_1}}\right)-{\displaystyle \frac{h_1{y}_0^{\prime \prime }}{2}},\hfill \end{array}$$ which is exactly the same as that proposed by Su and Liu [21] (p. 13).

References

1. Du, S.; Ward, R.K. Improved face representation by nonuniform multilevel selection of Gabor convolution features. IEEE Trans. Syst. Man and Cybern. Part B Cybern.; 2009; 39, pp. 1408-1419.

2. Chang, K.; Bowyer, K.W.; Sarkar, S.; Victor, B. Comparison and combination of ear and face images in appearance-based biometrics. IEEE Trans. Patt. Anal. Machine Intell.; 2003; 25, pp. 1160-1165. [DOI: https://dx.doi.org/10.1109/TPAMI.2003.1227990]

3. Aloraibi, A.Q. Image morphing techniques: A review. Technium; 2023; 9, pp. 41-53. [DOI: https://dx.doi.org/10.47577/technium.v9i.8699]

4. Indrawal, D.; Sharma, A. Multi-module convolutional neural network based optimal face recognition with minibatch optimization. Int. J. Image Graph. Signal Process.; 2022; 3, pp. 32-46. [DOI: https://dx.doi.org/10.5815/ijigsp.2022.03.04]

5. Patel, A.; Lapsiwala, P. Image morphing algorithm: A survey. Int. J. Comput. Appl.; 2015; 5, pp. 156-160.

6. Scherhag, U.; Rathgeb, C.; Merkle, J.; Breithaupt, R.; Bush, C. Face recognition systems under morphing attacks: A Survey. IEEE Access; 2019; 7, pp. 2169-3536. [DOI: https://dx.doi.org/10.1109/ACCESS.2019.2899367]

7. Scherhag, U.; Rathgeb, C.; Merkle, J.; Busch, C. Deep face representations for differential morphing attack detection. IEEE Trans. Inf. Forensics Secur.; 2020; 15, pp. 3625-3639. [DOI: https://dx.doi.org/10.1109/TIFS.2020.2994750]

8. Tuncer, T.; Dogan, S.; Subasi, A. Automated facial expression recognition using novel textural transformation. J. Ambient Intell. Human. Comput.; 2023; 14, pp. 9439-9449. [DOI: https://dx.doi.org/10.1007/s12652-023-04612-x]

9. Guo, G.; Zhang, N. A survey on deep learning based face recognition. Comput. Vis. Image Underst.; 2019; 189, 102805. [DOI: https://dx.doi.org/10.1016/j.cviu.2019.102805]

10. You, M.; Han, X.; Xu, Y.; Li, L. Systematic evaluation of deep face recognition methods. Neurocomputing; 2020; 388, pp. 144-156. [DOI: https://dx.doi.org/10.1016/j.neucom.2020.01.023]

11. Venkatesh, S.; Ramachandra, R.; Raja, K.; Busch, C. Face morphing attack generation and detection: A comprehensive survey. IEEE Trans. Technol. Soc.; 2021; 2, pp. 128-145. [DOI: https://dx.doi.org/10.1109/TTS.2021.3066254]

12. Li, Z.C.; Chiang, J.Y.; Suen, C.Y. Face transformation with harmonic models by finite volume method with Delaunay triangulation. IEEE Trans. Syst. Man and Cybern. Part B Cybern.; 2010; 40, pp. 1543-1554.

13. Li, Z.C. Analysis of discrete techniques for image transformations. Numer. Algor.; 1996; 13, pp. 225-263. [DOI: https://dx.doi.org/10.1007/BF02207696]

14. Li, Z.C.; Bui, T.D.; Tang, Y.Y.; Suen, C.Y. Computer Transformation of Digital Images and Patterns; World Scientific: Singapore, London, UK, 1989.

15. Farin, G. Curve and Surfaces for Computer Aided Geometric Design, A Practical Guide; 2nd ed. Academic Press: Cambridge, MA, USA, 1990.

16. Foley, J.D.; van Dam, A.; Feiner, S.K.; Hughes, J.K. Computer Graphics, Principles and Practice; 2nd ed. Addison-Wesley: Reading, MA, USA, 1990.

17. Li, Z.C.; Huang, H.T. Blending curves for landing problems by numerical differential equations, II. Numerical methods. Comput. Math. Appl.; 2000; 39, pp. 165-187. [DOI: https://dx.doi.org/10.1016/S0898-1221(00)00052-3]

18. Leichter, I.; Lindenbaum, M.; Rivlin, E. Tacking by afffine kernel transformations using color and boundary cues. IEEE Trans. Patt. Anal. Machine Intell.; 2009; 31, pp. 164-171. [DOI: https://dx.doi.org/10.1109/TPAMI.2008.194] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19029554]

19. Stahl, J.S.; Wang, S. Globally optimal grouping of symmetric closed boundaries by combining boundary and region information. IEEE Trans. Patt. Anal. Machine Intell.; 2008; 30, pp. 395-411. [DOI: https://dx.doi.org/10.1109/TPAMI.2007.1186]

20. DeCarlo, D.; Metaxas, D. Shape evolution with structural and topological changes using blending. IEEE Trans. Patt. Anal. Machine Intell.; 2003; 20, pp. 1186-1205. [DOI: https://dx.doi.org/10.1109/34.730554]

21. Su, B.Q.; Liu, D.Y. Computational Geometry Curve and Surface Modeling; Academic Press: Boston, MA, USA, 1989.

22. Ostermann, J. Face Animation in MPEG-4. MPEG-4 Facial Animation: The Standard, Implementation and Applications; Pandzic, I.S.; Forchheimer, R. John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2002; 21.

23. Hageman, L.A.; Young, D.M. Applied Iteration Method; Academic Press: Boston, MA, USA, 1981.

24. Chen, C.C. Face Transformation by Harmonic Model, Generating the Face Boundary. Master’s Thesis; Department of Applied Mathematics, National Sun Yat-sen University: Kaohsiung, Taiwan, 2002.

25. Li, Z.C. Combined Methods for Elliptic Equations with Singularities, Interface and Infinities; Kluwer Academic Publishers: Boston, MA, USA, Amsterdam, The Nerzerland, 1998.

26. Castleman, K.R. Geometric transformations. Microscope Image Processing; 2nd ed. Merchant, F.A.; Castleman, K.R. Academic Press/Elsevier: London, UK, 2023; pp. 47-54.

27. Lakemond, N.; Holmberg, G.; Pettersson, A. Digital transformation in complex systems. IEEE Trans. Eng. Manag.; 2024; 71, pp. 192-204. [DOI: https://dx.doi.org/10.1109/TEM.2021.3118203]

28. de la Rosa, F.L.; Gomez-Sirvent, J.L.; Sanchez-Reolid, R.; Morales, R.; Fernandez-Caballero, A. Geometric transformation-based data augmentation on defect classification of segmented images of semiconductor materials using a ResNet50 convolutional neural network. Expert Syst. Appl.; 2022; 206, 117731.

29. Luhmann, T.; Robson, S.; Kyle, S.; Boehm, J. Close-Range Photogrammetry and 3D Imaging; De Gruyter: Berlin, Germany, Boston, MA, USA, 2020; [DOI: https://dx.doi.org/10.1515/9783110607253]

30. Holden, M. A review of geometric transformations for nonrigid body registration. IEEE Trans. Med. Imaging; 2008; 27, pp. 111-128. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18270067]

31. Li, J.; Zhou, S.K.; Chellappa, R. Appearance modeling using a geometric transform. IEEE Trans. Image Process.; 2009; 18, pp. 889-902.

32. Ma, K.; Duanmu, Z.; Wang, Z. Geometric transformation invariance image quality assessment using convolutional neural networks. Proceedings of the 2018 IEEE International Conference on Acoustics; Calgary, AB, Canada, 15–20 April 2018; pp. 6732-6736.

33. Gao, X.; Lu, W.; Tao, D.; Li, X. Image quality assessment based on multiscale geometric analysis. IEEE Trans. Image Process.; 2009; 18, pp. 1409-1423. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19447715]

34. Fang, L.; Shi, Z.; Liu, Y.; Li, C.; Pang, M.; Zhao, E. A general geometric transformation model for linescan image registration. EURASIP J. Adv. Signal Process.; 2023; 78, [DOI: https://dx.doi.org/10.1186/s13634-023-01041-y]

35. Chen, J.; Wang, L.; Li, X.; Fang, Y. Arbitrary Continuous Geometric Transformation Networks for Image Registration. Proceedings of the Advances in Neural Information Processing Systems 32 (NeurIPS 2019); Vancouver, BC, Canada, 8–14 December 2019.

36. Tian, C.; Zheng, M.; Zuo, W.; Zhang, B.; Zhang, Y.; Zhang, D. Multi-stage image denoising with the wavelet transform. Patt. Recogn.; 2023; 134, 109050. [DOI: https://dx.doi.org/10.1016/j.patcog.2022.109050]

37. You, N.; Han, L.; Zhu, D.; Song, W. Research on image denoising in edge detection based on wavelet transform. Appl. Sci.; 2023; 13, 1837. [DOI: https://dx.doi.org/10.3390/app13031837]

38. Püspöki, Z. Local Geometric Transformations in Image Analysis; EPFL: Lausanne, Switzerland, 2016; [DOI: https://dx.doi.org/10.5075/epfl-thesis-7014]

39. Palumbo, R.; Adams, R.B., Jr.; Hess, U.; Kleck, R.E.; Zebrowitz, L. Age and gender differences in facial attractiveness, but not emotion resemblance, contribute to age and gender stereotypes. Front. Psychol.; 2017; 8, pp. 1-18. [DOI: https://dx.doi.org/10.3389/fpsyg.2017.01704] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29033881]

40. Taskiran, M.; Kahraman, N.; Erdem, C.E. Face recognition: Past, present and future (a review). Digit. Signal Process.; 2020; 106, pp. 1-28. [DOI: https://dx.doi.org/10.1016/j.dsp.2020.102809]

41. Vandergraft, J.S. Introduction to Numerical Computations; Academic Press: Cambridge, MA, USA, 1983.