# Intelligent Reflecting Surface-Aided Maneuvering Target Sensing: True
Velocity Estimation

Lei Xie, Xianghao Yu, and Shenghui Song
Dept. of ECE, The Hong Kong University of Science and Technology, Hong Kong
Email: $\\{$eelxie, eexyu<EMAIL_ADDRESS>This work was supported by the
Shenzhen Science and Technology Innovation Committee under Grant
SGDX20210823103201006 and the HKUST-BDR Joint Research Institute under Grant
OKT22EG04.

###### Abstract

Maneuvering target sensing will be an important service of future vehicular
networks, where precise velocity estimation is one of the core tasks. To this
end, the recently proposed integrated sensing and communications (ISAC)
provides a promising platform for achieving accurate velocity estimation.
However, with one mono-static ISAC base station (BS), only the radial
projection of the true velocity can be estimated, which causes serious
estimation error. In this paper, we investigate the estimation of the true
velocity of a maneuvering target with the assistance of an intelligent
reflecting surface (IRS). We propose an efficient velocity estimation
algorithm by exploiting the two perspectives from the BS and IRS to the
target. We propose a two-stage scheme where the true velocity can be recovered
based on the Doppler frequency of the BS-target link and BS-IRS-target link.
Experimental results validate that the true velocity can be precisely
recovered and demonstrate the advantage of adding the IRS.

###### Index Terms:

Intelligent reflecting surface, integrated sensing and communications, Doppler
effect, true velocity estimation.

##  I Introduction

Vehicular to everything (V2X) communications is becoming a key enabling
technology for many thrilling applications, e.g., smart city and intelligent
transportation systems. However, V2X poses stringent requirements on both
sensing accuracy and communication quality. To this end, the recently proposed
integrated sensing and communication (ISAC) provides a promising platform that
amalgamates the abilities of sensing and communication, especially for 5G and
beyond systems [1, 2, 3]. However, high-accuracy beam tracking is necessary to
compensate for the path loss imposed by the high-frequency signals and
guarantee the quality-of-service (QoS) for V2X, particularly in high-mobility
vehicular scenarios [4].

To achieve accurate beam-tracking, precise velocity estimation is required.
Doppler effect, which refers to the change in frequency of an electromagnetic
wave caused by the relative motion between the observer and the wave source,
is typically utilized for velocity estimation. Existing works on velocity
estimation mainly adopted matched filters (MFs) to estimate the motion
parameters [5, 6, 4]. Micro-Doppler frequency estimation was considered in [5]
for orthogonal frequency division multiplexing (OFDM) radar systems. In [6],
the authors proposed to mitigate the inter-carrier interference caused by the
Doppler effect to improve the performance of Doppler estimation by
oversampling. A predictive beamforming design was investigated for an ISAC
system in [4]. In general, the performance of the MF-based method is limited
by the inherent grid issue. To address this problem, a joint range-velocity
estimation was proposed in [7] based on the multiple signal classification
(MUSIC) algorithm in OFDM radar systems, which achieves a higher resolution
than the MF-based methods. However, due to the nature of the Doppler effect,
the conventional mono-static ISAC BS can only measure the radial projection of
the true velocity, and thus the true velocity cannot be uniquely recovered.
One solution is to obtain another perspective toward the maneuvering target so
that the true velocity can be uniquely determined.

In recent years, intelligent reflecting surfaces (IRSs) have been introduced
to enhance the performance of both communication and sensing [8, 9, 10]. Some
research efforts have been made to IRS-aided systems to tackle user equipment
(UE) mobility [11, 12, 13]. The authors in [11] considered a continuous model
for IRS-aided satellite communication, based on which the IRS phase shifters
are optimized to simultaneously maximize the received power and minimize the
delay and Doppler spread. In [12], a two-stage protocol for channel estimation
was proposed for an IRS-aided high-mobility communication system, which
deploys an IRS at a high-speed vehicle to enhance the QoS for passengers by
mitigating the Doppler effect. In [13], the authors considered the
localization problem in an IRS-aided single-input single-output (SISO) system
by taking into account the UE mobility and spatial wideband effects. However,
the potential of IRSs for true velocity estimation has not been fully
understood.

In this paper, by leveraging the additional perspective provided by the IRS,
we study the true velocity estimation in an IRS-aided ISAC system. By
exploiting the geometric information, the true velocity can be recovered based
on the Doppler frequency of the BS-target link and BS-IRS-target link111For
simplicity, we will refer to the BS-target link and BS-IRS-target link as the
direct link and IRS link, respectively.. To estimate the Doppler frequency of
both links, we propose a two-stage scheme. In the first stage, a coarse
estimation of radial projection of the true velocity is obtained by turning
off the IRS. In the second stage, a gridless estimation for the Doppler
frequency of the direct link and IRS link is performed by exploiting the
method of direction estimation (MODE) [14]. Compared with existing gridless
estimation methods such as root-multiple signal classification (root-MUSIC)
[15] and estimation of signal parameters via rotation invariance techniques
(ESPRIT) [16], the MODE approach is shown to provide more accurate parameter
estimation [14].

##  II Preliminary

Consider an ISAC system, as illustrated in Fig. 1, where the velocity vector
of the point-like maneuvering target, denoted by $\mathbf{v}$222The velocity
vector is defined by
$\mathbf{v}=\left[|\mathbf{v}|\cos\theta_{v},|\mathbf{v}|\sin\theta_{v}\right]^{\mathrm{T}}$,
where $\theta_{v}$ represents the direction of velocity as illustrated in Fig.
2., is estimated by one BS with the help of an IRS. The BS and IRS are
equipped with a uniform linear array (ULA) with $N$ and $M$ antennas,
respectively.

Figure 1: Illustration of the considered IRS-aided ISAC system.

###  II-A Signal Model

Assume that the directions of the target with respect to the BS and IRS,
denoted by ${\theta}_{tb}$ and ${\theta}_{it}$ in Fig. 2(b), are well
estimated in advance. To estimate the true velocity, the BS transmits a
probing beam

$\begin{split}\mathbf{f}&=\mathbf{a}\left(\theta_{tb}\right)\triangleq\left[1,e^{j\frac{2\pi
d}{\lambda}\cos\theta_{tb}},\cdots,e^{j\frac{2\pi
d(N-1)}{\lambda}\cos\theta_{tb}}\right]^{\mathrm{T}},\end{split}$ (1)

where $d$ and $\lambda$ represent the inner spacing of the ULA and the
wavelength, respectively. Here, $\mathbf{a}(\cdot)$ represents the steering
vector of the antenna array of the BS towards the target. Note that due to the
use of the narrow beam, the effect of the non-light-of-sight (NLoS) path is
considered negligible.

The echo of the $k$-th pilot symbol $s_{k}$ received by the BS is given by

$\begin{split}&\mathbf{y}_{r,k}=\underbrace{\alpha_{d}e^{j2\pi\mu_{d}(k-1)T_{s}}\mathbf{a}(\theta_{tb})\mathbf{a}^{\mathrm{H}}(\theta_{tb})\mathbf{f}s_{k}}_{\text{direct
link}}\\\
&+\underbrace{\alpha_{r}e^{j2\pi\mu_{r}(k-1)T_{s}}\mathbf{G}\mathbf{\Psi}\mathbf{b}(\theta_{it})\mathbf{a}^{\mathrm{H}}(\theta_{tb})\mathbf{f}s_{k}}_{\text{IRS
link}}+\mathbf{p}_{r,k},\end{split}$ (2)

where

$\begin{split}\mathbf{b}\left({\theta}_{it}\right)\triangleq\left[1,e^{j\frac{2\pi
d}{\lambda}\cos{\theta}_{it}},\cdots,e^{j\frac{2\pi
d(M-1)}{\lambda}\cos{\theta}_{it}}\right]^{\mathrm{T}}\end{split}$ (3)

represents the response vector of the IRS, $\mathbf{G}$ denotes the channel
between BS and IRS, which is modeled as a Rician channel comprising a light-
of-sight (LoS) path and a number of NLoS paths, $\mathbf{\Psi}$ denotes the
phase shifter which is designed to align toward the target333The optimized
design of phase shifter will be left as future work., $\alpha_{d}$ and
$\alpha_{r}$ denote the complex channel gain, accounting for the path loss and
target reflectivity444The fluctuations of target reflectivity are typically
modeled using the standard Swerling classes. In this paper, we adopt the
Swerling I model where the reflectivity is assumed to be constant within a
symbol [17].. $\mu_{d}$ and $\mu_{r}$ represent the Doppler frequencies
corresponding to the direct link and IRS link, respectively, $T_{s}$ denotes
the symbol period, and $\\{\mathbf{p}_{r,k}\\}_{k=1}^{N_{r}}$ denote the noise
vectors whose elements are drawn independently from a complex Gaussian
distribution with zero mean and covariance matrix
$\frac{\sigma_{r}^{2}}{N}\mathbf{I}$ [13].

###  II-B Geometrical Insights of Doppler Effect

(a)  (b)

Figure 2: Illustration of Doppler effect. (a) Direct link; (b) IRS link.

Since the Doppler effect depends on the mobility of the target, we propose to
estimate the true velocity $\mathbf{v}$ by exploiting the Doppler frequencies
$\\{\mu_{d},\mu_{r}\\}$.In the following, we first discuss the geometrical
insights of the Doppler effect to reveal the relation between true velocity
$\mathbf{v}$ and the Doppler frequencies $\\{\mu_{d},\mu_{r}\\}$.

The Doppler frequency of the relative motion between the BS and the
maneuvering target is given by

$f_{tb}=\frac{v\cos\left(\theta_{v}-\theta_{tb}\right)}{\lambda},$ (4)

where $v$ denotes the amplitude of the true velocity, and $\theta_{v}$
represents the angle between the velocity direction and the $x$-axis, which is
marked in red in Fig. 2. Note that $v\cos\left(\theta_{v}-\theta_{tb}\right)$
is the projection of velocity on the line across the BS and target, which is
marked in dark green in Fig. 2(a). Then, the Doppler frequency corresponding
to the direct link is given by

$\mu_{d}=2f_{tb}=\frac{2v\cos\left(\theta_{v}-\theta_{tb}\right)}{\lambda}.$
(5)

The deployment of the IRS creates another path, which allows the BS to observe
the velocity from an additional perspective. The Doppler frequency of the
relative motion between the IRS and the maneuvering target is given by

$f_{it}=\frac{v\cos\left(\theta_{it}-\theta_{v}\right)}{\lambda}.$ (6)

Then, the Doppler frequency corresponding to the IRS link is given by

$\begin{split}\mu_{r}&=f_{tb}+f_{it}=\frac{2v}{\lambda}\cos\theta_{1}\cos\theta_{2},\end{split}$
(7)

where $\theta_{1}={\frac{\theta_{tb}+\theta_{it}}{2}-\theta_{v}}$ and
$\theta_{2}=\frac{\theta_{it}-\theta_{tb}}{2}$ which are marked in dark green
and purple in Fig. 2(b), respectively.

From (7), $\mu_{r}$ comes from two projections:

  1. 1.

$\cos\theta_{1}$: This term projects the velocity to the angular bisector of
the angle $\angle_{\text{BS-Target-IRS}}$. Note that the angular bisector is
perpendicular to the tangent line of the ellipse with foci at the BS and IRS.

  2. 2.

$\cos\theta_{2}$: This term projects the projected speed $v\cos\theta_{1}$ to
the line connecting the BS and the target (or the line connecting the IRS and
the target). Since the direction of the projected speed $v\cos\theta_{1}$ is
fixed, there will be no ambiguity on the velocity at the projection
corresponding to $\cos\theta_{2}$.

Remark: The conventional mono-static radar can only access $\mu_{d}$, which
corresponds to the radial component of the velocity. As a result, the tangent
projection will cause ambiguity, e.g., different velocities may yield the same
observation at the BS. One example is illustrated in Fig. 2(a) where two
velocity vectors, i.e., $\mathbf{v}_{amb}$ (marked in gray) and $\mathbf{v}$
(marked in red), give the same projection. Similarly, for the IRS link, the
velocity ambiguity occurs at the tangent line of the ellipse, as illustrated
in Fig. 2(b). Thus, either of the two links, by itself, will cause velocity
ambiguity. Nevertheless, by jointly considering the two links, the true
velocity vector can be uniquely determined. Note that we assume that the
distance between the BS and the IRS should be comparable to the distance
between the target and the BS/IRS.

The following proposition provides the true velocity estimation in the IRS-
aided system by jointly leveraging $\mu_{d}$ and $\mu_{r}$.

###### Proposition 1

The direction and absolute value of the true velocity can be determined
uniquely as

$\begin{split}{\theta}_{v}=\arctan\left[\left(1-\frac{\mu_{d}}{\mu_{r}}\right)\cot\left(\frac{\theta_{it}-\theta_{tb}}{2}\right)\right]+\frac{\theta_{tb}+\theta_{it}}{2},\end{split}$
(8)
$\begin{split}{v}=\frac{\lambda}{2}\cdot\frac{\mu_{d}}{\cos\left(\theta_{v}-\theta_{tb}\right)}.\end{split}$
(9)

_Proof_ : See Appendix. $\blacksquare$

##  III True Velocity Estimation

According to _Proposition 1_ , to estimate the true velocity, we need to
estimate $\mu_{d}$ and $\mu_{r}$. In this section, we propose a two-stage
scheme to estimate the true velocity of the maneuvering target. Specifically,
in the first stage, we obtain a coarse estimation of $\mu_{d}$. In the second
stage, $\mu_{d}$ and $\mu_{r}$ are estimated accurately by utilizing the MODE.
The coarse estimation serves two purposes. First, it can be utilized to
initialize the iterative algorithm to speed up the iteration in the second
stage. Second, the coarse estimation can help match the estimated Doppler
frequencies to different links.

###  III-A Stage I: Coarse Estimation of $\mu_{d}$

We first turn off the IRS to obtain a coarse estimation of $\mu_{d}$ only
based on the direct link. This can be implemented through a receiving beam
with zero gain on the direction of IRS. In this case, the system degenerates
to a conventional mono-static one. Assume the BS transmits $N_{d}$ pilot
symbols in Stage I. Given the pilot symbol $s_{k}=1$, the received echo of the
$k$-th symbol at the BS is given by

$\begin{split}\mathbf{y}_{d,k}=&\alpha_{d}e^{j2\pi\mu_{d}(k-1)T_{s}}\mathbf{a}(\theta_{tb})\mathbf{a}^{\mathrm{H}}(\theta_{tb})\mathbf{f}+\mathbf{p}_{d,k},\end{split}$
(10)

where $\\{\mathbf{p}_{d,k}\\}_{k=1}^{N_{d}}$ denotes the noise vector whose
elements are drawn independently from a complex Gaussian distribution with
zero mean and covariance matrix $\frac{\sigma_{d}^{2}}{N}\mathbf{I}$

Define $\mathbf{w}_{d}=\mathbf{a}(\theta_{tb})$ as the combiner corresponding
to the received beam which points to the estimated target position. Then the
output for the $k$-th symbol is given by

$\begin{split}z_{d,k}=\mathbf{w}_{d}^{\mathrm{H}}\mathbf{y}_{d,k}=\beta_{d}e^{j2\pi\mu_{d}(k-1)T_{s}}&+n_{d,k},\end{split}$
(11)

where
$\beta_{d}=\alpha_{d}\mathbf{w}_{d}^{\mathrm{H}}\mathbf{a}(\theta_{tb})\mathbf{a}^{\mathrm{H}}(\theta_{tb})\mathbf{f}$
and $n_{d,k}=\mathbf{w}_{d}^{\mathrm{H}}\mathbf{p}_{d,k}$ denotes the combined
noise following the Gaussian distribution, i.e.,
$n_{d,k}\sim\mathcal{CN}(0,\sigma_{d}^{2})$. Stacking $z_{d,k}$ into a vector
yields

$\begin{split}\mathbf{z}_{d}&=\left[z_{d,1},z_{d,2},\cdots,z_{d,N_{d}}\right]^{\mathrm{T}}\in\mathbb{C}^{N_{d}\times
1}.\end{split}$ (12)

Define the steering vector in the Doppler domain as

$\begin{split}\mathbf{a}_{f}^{(N_{d})}(\mu)&\triangleq\left[1,e^{-j2\pi\mu
T_{s}},\cdots,e^{-j2\pi\mu(N_{d}-1)T_{s}}\right]^{\mathrm{T}}.\end{split}$
(13)

The coarse estimation of $\mu_{d}$ is obtained as

$\begin{split}\hat{\mu}_{d,\text{coarse}}=\arg\max_{\mu}|\mathbf{z}_{d}^{\mathrm{T}}\mathbf{a}_{f}^{(N_{d})}(\mu)|^{2}.\end{split}$
(14)

###  III-B Stage II: Refined Gridless Estimation of $\mu_{d}$ and $\mu_{r}$

At this stage, we turn on the IRS to obtain a refined estimation of $\mu_{d}$
and $\mu_{f}$. The BS transmits $N_{r}$ pilot symbols at this stage. Define
$\mathbf{w}_{r}=\mathbf{a}(\theta_{tb})+\mathbf{a}({\theta}_{bi})$ as the
combiner illuminated at the estimated target position and IRS, simultaneously,
where $\theta_{bi}$ denotes the direction of the IRS with respect to the BS.
Then, the output for the $k$-th symbol is given by

$\begin{split}z_{r,k}=\mathbf{w}_{r}^{\mathrm{H}}\mathbf{y}_{r,k}=\beta_{r}e^{j2\pi\mu_{r}(k-1)T_{s}}&+n_{r,k},\end{split}$
(15)

where
$\beta_{r}=\alpha_{r}\mathbf{w}_{r}^{\mathrm{H}}\mathbf{G}\mathbf{\Psi}\mathbf{b}(\theta_{it})\mathbf{a}^{\mathrm{H}}(\theta_{tb})\mathbf{f}$
and $n_{r,k}=\mathbf{w}_{r}^{\mathrm{H}}\mathbf{p}_{r,k}$ denotes the combined
noise, which follows the Gaussian distributions, i.e.,
$n_{r,k}\sim\mathcal{CN}(0,\sigma_{r}^{2})$.

Before proceeding, we first construct a set of $P\times 1$ signal vectors by
rearranging $\\{z_{r,k}\\}_{k=1}^{N_{r}}$. By doing that, the Doppler
frequency estimation problem can be solved by the MODE. Given a positive
integer $P$, the signal vector is defined as

$\begin{split}\mathbf{z}_{r,k}=\left[z_{r,k},z_{r,k-1},\cdots,z_{r,k-P+1}\right]^{\textrm{T}}&\in\mathbb{C}^{P\times
1}.\end{split}$ (16)

From (15), (16) can be rewritten as

$\mathbf{z}_{r,k}=\mathbf{A}\mathbf{x}_{r,k}+\mathbf{n}_{r,k},$ (17)

where

$\mathbf{A}=\left[\mathbf{a}_{f}^{(P)}(\mu_{d}),\mathbf{a}_{f}^{(P)}(\mu_{r})\right],$
(18)
$\begin{split}\mathbf{x}_{r,k}=[\beta_{d}e^{j2\pi\mu_{d}(k-1)T_{s}},\beta_{r}e^{j2\pi\mu_{r}(k-1)T_{s}}]^{\textrm{T}},\end{split}$
(19)
$\mathbf{n}_{r,k}=\left[n_{r,k},n_{r,k-1},\cdots,n_{r,k-P+1}\right]^{\textrm{T}}.$
(20)

The covariance matrix of $\mathbf{z}_{r,k}$ is

$\begin{split}\mathbf{R}_{r}=\mathbb{E}\left\\{\mathbf{z}_{r,k}\mathbf{z}_{r,k}^{\textrm{H}}\right\\}=\mathbf{A}\mathbf{\Phi}\mathbf{A}^{\textrm{H}}+\sigma_{r}^{2}\mathbf{I},\end{split}$
(21)

where $\mathbb{E}{\\{\cdot\\}}$ denotes the expectation over $k$ and
$\mathbf{\Phi}=\mathbb{E}{\\{\mathbf{x}_{r,k}\mathbf{x}_{r,k}^{\mathrm{H}}\\}}$.
The eigenvalue decomposition of $\mathbf{R}_{r}$ gives

$\begin{split}\mathbf{R}_{r}=\sum_{i=1}^{P}\lambda_{i}\mathbf{g}_{i}\mathbf{g}_{i}^{\mathrm{H}},\end{split}$
(22)

where
$\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}=\cdots=\lambda_{P}=\sigma_{r}^{2}$
are the eigenvalues,
$\mathbf{G}_{s}=\left[\mathbf{g}_{1},\mathbf{g}_{2}\right]\in\mathbb{C}^{P\times
2}$ and
$\mathbf{G}_{N}=\left[\mathbf{g}_{3},\mathbf{g}_{4},\cdots,\mathbf{g}_{P}\right]\in\mathbb{C}^{P\times(P-2)}$
correspond to the signal and noise subspace, respectively. Note that the
linear span of $\mathbf{G}_{s}$ is the same as that of $\mathbf{A}$, i.e.,

$\text{span}(\mathbf{G}_{s})=\text{span}(\mathbf{A}).$ (23)

Note that the maximum likelihood estimate (MLE) of
$\bm{\mu}=\\{\mu_{d},\mu_{r}\\}$ depends on $\mathbf{R}_{r}$, which is unknown
in practice. To address this, we first obtain the estimated eigenvectors
$\\{\widehat{\mathbf{g}}_{i}\\}_{i=1}^{P}$ and the estimated eigenvalues
$\\{\widehat{\lambda}_{i}\\}_{i=1}^{P}$ by performing the eigenvalue
decomposition on the sample covariance matrix of
$\\{\mathbf{z}_{r,k}\\}_{k=P}^{N_{r}}$, i.e.,

$\begin{split}\widehat{\mathbf{R}}_{r,\text{SCM}}\triangleq\frac{1}{N_{r}-P+1}\sum_{k=P}^{N_{r}}\mathbf{z}_{r,k}\mathbf{z}_{r,k}^{\textrm{H}}.\end{split}$
(24)

Then, an approximated maximum likelihood estimate (MLE) of
$\bm{\mu}=\\{\mu_{d},\mu_{r}\\}$ is given by [14]

$\begin{split}{\bm{\widehat{\mu}}=\arg\min_{\bm{\mu}}\mathrm{tr}\left\\{\left(\mathbf{I}-\Phi_{\mathbf{A}}(\bm{\mu})\right)\mathbf{\widehat{G}}_{s}\widehat{\mathbf{\Gamma}}\mathbf{\widehat{G}}_{s}^{\rm{H}}\right\\},}\end{split}$
(25)

where $\mathrm{tr}(\cdot)$ denotes the trace,
$\Phi_{\mathbf{A}}(\bm{\mu})=\mathbf{A}\left({\mathbf{A}}^{\mathrm{H}}\mathbf{A}\right)^{-1}{\mathbf{A}}^{\mathrm{H}}$
represents the orthogonal projector onto the range space of $\mathbf{A}$, and
$\mathbf{\widehat{G}}_{s}=[\widehat{\mathbf{g}}_{1},\widehat{\mathbf{g}}_{2}]$,
while

$\begin{split}\widehat{\mathbf{\Gamma}}\triangleq{\rm{diag}}\left\\{\widehat{\lambda}_{1}^{-1}\left(\widehat{\lambda}_{1}-\widehat{\sigma}_{r}^{2}\right)^{2},\widehat{\lambda}_{2}^{-1}\left(\widehat{\lambda}_{2}-\widehat{\sigma}_{r}^{2}\right)^{2}\right\\},\end{split}$
(26)

and $\widehat{\sigma}_{r}^{2}=\sum_{i=3}^{P}\widehat{\lambda}_{i}/(P-2)$.

However, the problem in (25) does not have an analytical solution and is thus
hard to solve directly. To solve (25), we then consider introducing some
tractable intermediate variables, based on which $\mu_{d}$ and $\mu_{r}$ can
be estimated. Construct a polynomial

$\begin{split}\left(1-e^{j2\pi\mu_{d}T_{s}}\omega\right)\left(1-e^{j2\pi\mu_{r}T_{s}}\omega\right)=1+c_{1}\omega+c_{2}\omega^{2},\end{split}$
(27)

which has $2$ roots in the complex plane. It can be validated that
$c_{1}=-\left(e^{j2\pi\mu_{d}T_{s}}+e^{j2\pi\mu_{r}T_{s}}\right)$ and
$c_{2}=e^{j2\pi(\mu_{d}+\mu_{r})T_{s}}$. Clearly,
$\omega_{d}=e^{-j2\pi\mu_{d}T_{s}}$ and $\omega_{r}=e^{-j2\pi\mu_{r}T_{s}}$
are the roots of this polynomial. Then, the strategy of solving (25) is to
obtain the intermediate variables $c_{1}$ and $c_{2}$, based on which
$\mu_{d}$ and $\mu_{r}$ can be estimated, instead of estimating $\mu_{d}$ and
$\mu_{r}$ directly.

Define a matrix

$\begin{split}\mathbf{C}=\left[\begin{matrix}1&c_{1}&c_{2}&0&\cdots&0\\\
0&1&c_{1}&c_{2}&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&0\\\
0&\cdots&0&1&c_{1}&c_{2}\end{matrix}\right]\in\mathbb{C}^{(P-2)\times
P}.\end{split}$ (28)

From the definition of $\mathbf{A}$ in (18), we have
$\mathbf{C}\mathbf{A}=\mathbf{0}$, which indicates

$\begin{split}\mathbf{I}-\Phi_{\mathbf{A}}(\bm{\mu})=\Phi_{\mathbf{C^{\mathrm{H}}}}=\mathbf{C}^{\mathrm{H}}\left(\mathbf{C}\mathbf{C}^{\mathrm{H}}\right)^{-1}\mathbf{C}.\end{split}$
(29)

Substituting (29) into (25) produces

$\begin{split}&\mathrm{tr}\left\\{\mathbf{C}^{\mathrm{H}}\left(\mathbf{C}\mathbf{C}^{\mathrm{H}}\right)^{-1}\mathbf{C}\mathbf{\widehat{G}}_{s}\widehat{\mathbf{\Gamma}}\mathbf{\widehat{G}}_{s}^{\rm{H}}\right\\}\\\
=&\mathrm{tr}\left\\{\left(\mathbf{C}\mathbf{C}^{\mathrm{H}}\right)^{-1}\mathbf{C}\mathbf{\widehat{G}}_{s}\widehat{\mathbf{\Gamma}}\mathbf{\widehat{G}}_{s}^{\rm{H}}\mathbf{C}^{\mathrm{H}}\right\\}\\\
=&\mathrm{vec}\left\\{\mathbf{C}\mathbf{\widehat{G}}_{s}\right\\}^{\mathrm{H}}\mathbf{W}\mathrm{vec}\left\\{\mathbf{C}\mathbf{\widehat{G}}_{s}\right\\}.\end{split}$
(30)

where $\mathrm{vec}(\cdot)$ denotes the vectorization operator and

$\begin{split}\mathbf{W}=\widehat{\mathbf{\Gamma}}\otimes\left(\mathbf{C}\mathbf{C}^{\mathrm{H}}\right)^{-1}.\end{split}$
(31)

Recall from (27) that
$\mathbf{c}\triangleq\left[c_{1},c_{2}\right]^{\mathrm{T}}$ is related to
$\bm{\mu}$. Then, the MLE of $\mathbf{c}$ is given by

$\begin{split}\mathbf{\widehat{c}}=\arg\min_{\mathbf{c}}\mathrm{vec}\left\\{\mathbf{C}\mathbf{\widehat{G}}_{s}\right\\}^{\mathrm{H}}\mathbf{W}\mathrm{vec}\left\\{\mathbf{C}\mathbf{\widehat{G}}_{s}\right\\}.\end{split}$
(32)

Define matrix

$\begin{split}\mathbf{\widehat{\Psi}}_{j}\triangleq\left[\begin{matrix}\widehat{g}_{2,j}&\widehat{g}_{3,j}\\\
\widehat{g}_{3,j}&\widehat{g}_{4,j}\\\ \vdots&\vdots\\\
\widehat{g}_{P-1,j}&\widehat{g}_{P,j}\\\
\end{matrix}\right]\in\mathbb{C}^{(P-2)\times 2},\quad j=1,2,\end{split}$ (33)

and vector

$\begin{split}\mathbf{\widehat{q}}_{j}\triangleq-\left[\widehat{g}_{1,j},\widehat{g}_{2,j},\cdots,\widehat{g}_{P-2,j}\right]^{\mathrm{T}},\quad
j=1,2,\end{split}$ (34)

where $\widehat{g}_{i,j}$ denotes the $(i,j)$-th entry of
$\widehat{\mathbf{G}}_{s}$. We have

$\begin{split}\mathrm{vec}\left\\{\mathbf{C}\mathbf{\widehat{G}}_{s}\right\\}&=\mathbf{\widehat{\Psi}}\mathbf{c}-\mathbf{\widehat{q}},\end{split}$
(35)

where

$\mathbf{\widehat{\Psi}}=\left[\mathbf{\widehat{\Psi}}_{1}^{\mathrm{T}},\mathbf{\widehat{\Psi}}_{2}^{\mathrm{T}}\right]^{\mathrm{T}},$
(36)
$\mathbf{\widehat{q}}=\left[\mathbf{\widehat{q}}_{1}^{\mathrm{T}},\mathbf{\widehat{q}}_{2}^{\mathrm{T}}\right]^{\mathrm{T}}.$
(37)

Substituting (35) into (32) leads to

$\begin{split}\mathbf{\widehat{c}}&=\arg\min_{\mathbf{c}}\left(\mathbf{\widehat{\Psi}}\mathbf{c}-\mathbf{\widehat{q}}\right)^{\mathrm{H}}\mathbf{W}\left(\mathbf{\widehat{\Psi}}\mathbf{c}-\mathbf{\widehat{q}}\right).\end{split}$
(38)

Note that (38) is a weighted least squares problem whose solution is given by

$\begin{split}\mathbf{\widehat{c}}_{\dagger}=\left(\mathbf{\widehat{\Psi}}^{\mathrm{H}}\mathbf{W}\mathbf{\widehat{\Psi}}\right)^{-1}\mathbf{\widehat{\Psi}}^{\mathrm{H}}\mathbf{W}\mathbf{\widehat{q}}.\end{split}$
(39)

However, $\bf{W}$ also depends on $\mathbf{C}$ through (31) and (28). In order
to tackle this challenge, we propose an iterative algorithm to find the
solution, which is summarized in Algorithm 1.

Algorithm 1 MODE-based velocity estimation algorithm

Initialize $\widehat{\mathbf{\Gamma}}$,$\mathbf{\widehat{\Psi}}$, and
$\mathbf{\widehat{q}}$ via (26), (36), and (37).
Initialize
$\mathbf{\widehat{c}}_{(0)}=[-\left(e^{j2\pi\hat{\mu}_{d,\text{coarse}}T_{s}}+1\right),e^{j2\pi\hat{\mu}_{d,\text{coarse}}T_{s}}]^{\mathrm{T}}$,
$t=0$.
Repeat:

  1. 1.

Update $\widehat{\mathbf{C}}_{(t+1)}$ using $\mathbf{\widehat{c}}_{(t)}$, via
(28)

  2. 2.

Update $\mathbf{\widehat{c}}_{(t+1)}$ using $\widehat{\mathbf{C}}_{(t+1)}$,
via (39)

  3. 3.

$t\leftarrow t+1$.

Until some termination condition is met.

Once $\mathbf{\widehat{c}}$ is obtained, we can solve the equation

$\begin{split}1+\widehat{c}_{1}\omega+\widehat{c}_{2}\omega^{2}=0\end{split}$
(40)

to obtain two roots

$\omega_{1}=\frac{-\widehat{c}_{1}-\sqrt{\widehat{c}_{1}^{2}-4\widehat{c}_{2}}}{2\widehat{c}_{2}},\omega_{2}=\frac{-\widehat{c}_{1}+\sqrt{\widehat{c}_{1}^{2}-4\widehat{c}_{2}}}{2\widehat{c}_{2}}.$
(41)

The corresponding $\mu_{1}$ and $\mu_{2}$ are obtained as

$\hat{\mu}_{1}=\frac{1}{2\pi
T_{s}}\mathcal{P}(\omega_{1}),\hat{\mu}_{2}=\frac{1}{2\pi
T_{s}}\mathcal{P}(\omega_{2}),$ (42)

where $\mathcal{P}(\cdot)$ denotes the phase of a complex number.

We can obtain the refined estimates of $\mu_{d}$ and $\mu_{r}$, i.e.,

$\begin{split}&\hat{\mu}_{d}=\mathop{\arg\min}_{\mu\in\\{\hat{\mu}_{1},\hat{\mu}_{2}\\}}|\hat{\mu}_{d,\text{coarse}}-\mu|,\\\
&\hat{\mu}_{r}=\mathop{\arg\max}_{\mu\in\\{\hat{\mu}_{1},\hat{\mu}_{2}\\}}|\hat{\mu}_{d,\text{coarse}}-\mu|.\end{split}$
(43)

By exploiting _Proposition 1_ , the direction and absolute value of true
velocity can be estimated.

##  IV Simulation Results

In this section, we show the performance of the proposed velocity estimator.
The number of antennas at the BS is set as $N=16$ and the number of elements
at the IRS is set as $M=32$. The carrier frequency is $f_{c}=3$ GHz and the
symbol period is $T_{s}=0.5$ ms. The number of symbols at stages I and II are
$N_{d}=16$ and $N_{r}=16$, respectively. Assume that the BS and IRS are
located at $(0,0)$ and $(20,0)$, respectively, where all coordinates are in
meters. The direction of the target with respect to the BS and IRS are
$\theta_{tb}=30^{\circ}$ and $\theta_{it}=120^{\circ}$, respectively.
$\mathbf{G}$ is modeled as a Rician channel comprising a LOS path and a number
of NLOS paths. The number of paths for $\mathbf{G}$ is set as $N_{p}=3$, where
the AoA and AoD parameters are uniformly generated from $[-\pi/2,\pi/2]$. The
Rician factor is set as $13.2$ dB. The signal-to-noise ratio (SNR) of the
direct link is defined as
$\text{SNR}=\frac{\mathbb{E}(\alpha_{d}^{2})}{\sigma_{r}^{2}}$.

Figure 3: Convergence of the proposed MODE-based method.  Figure 4: NMSE
performance versus SNR.

We first show the convergence performance of the proposed MODE-based method,
which involves an iterative process. For that purpose, we choose the average
discrepancy of $\widehat{\mathbf{c}}_{(t)}$ between two adjacent iterations as
the metric, i.e.,
$\mathcal{D}_{(t)}=||\widehat{\mathbf{c}}_{(t+1)}-\widehat{\mathbf{c}}_{(t)}||.$
For each abscissa, $10,000$ Monte-Carlo trials are carried out. From Fig. 3,
we can see that the iteration converges in about ten rounds.

We compare the proposed MODE-based method with two gridless alternatives,
i.e., root-MUSIC [15] and root-ESPRIT [16]. To quantify the estimation
accuracy, we evaluate the normalized mean squared error (NMSE), i.e.,
$\text{NMSE}=\sqrt{\frac{1}{N_{mc}}{\sum_{i=1}^{N_{mc}}\frac{||\mathbf{v}-\widehat{\mathbf{v}}_{i}||^{2}}{||\mathbf{v}||^{2}}}},$
where $N_{mc}$ denotes the number of Monte-Carlo trials and
$\widehat{\mathbf{v}}_{i}$ denotes the estimation result in the $i$-th Monte-
Carlo trial.

The NMSEs of different methods under different SNRs are shown in Fig. 4. The
velocity of target is $40$ m/s with $\theta_{v}=60^{\circ}$ and
$N_{mc}=10,000$. From Fig. 4, we observe that the true velocity can be
recovered with a considerable accuracy. As SNR increases, the NMSEs of all
these methods decrease. Furthermore, the proposed MODE can achieve higher
accuracy than MUSIC and ESPRIT. This is because the MODE can achieve the
minimum estimation variance among the class of methods based on the eigenspace
of the sample covariance matrix including MUSIC and ESPRIT [14].

Figure 5: NMSE performance versus velocity $|\mathbf{v}|$.

Next, we show the NMSE performance under different target velocities in Fig.
5. We set $\theta_{v}=60^{\circ}$, $N_{mc}=10,000$, and $\text{SNR}=0$ dB.
From Fig. 5, we can see that the NMSE performance is significantly improved
with the deployment of the IRS. This is because only a radial projection of
the true velocity can be estimated for a single BS without the IRS, such that
the tangent projection of the true velocity dominates the estimation error.
Meanwhile, with the IRS, the NMSE performance improves as the target velocity
increases. The reason for this phenomenon requires a further theoretical study
on the Cramér-Rao bound, which will be left to our future work. Moreover, the
gap between the curves with and without the IRS becomes larger as $v$
increases, which indicates that the benefit of the deployment of IRSs will be
more pronounced in the high-mobility case.

##  V Conclusion

This paper considered the true velocity estimation for a maneuvering target
with the assistance of an IRS. Thanks to the extra link provided by the IRS,
the BS can observe the true velocity from two different perspectives. In
particular, we proposed a two-stage algorithm to estimate the Doppler
frequency of the direct link and IRS link, based on which the true velocity
recovered. Compared with existing approaches, the proposed approach achieved
much better estimation performance. The work in this paper represents a small
step to fully exploit the power of IRS in ISAC applications.

## Proof of Proposition 1

(9) can be directly obtained from (5), thus we only prove (8) in the
following. Recalling the definition of $\theta_{1}$ and $\theta_{2}$, we have

$\theta_{v}-\theta_{tb}=\theta_{2}-\theta_{1}.$ (44)

Then, (5) can be reformulated as

$\mu_{d}=\frac{2v\cos\left(\theta_{2}-\theta_{1}\right)}{\lambda}.$ (45)

According to (7), we have

$\begin{split}\frac{\mu_{d}}{\mu_{r}}&=\frac{\cos\left(\theta_{2}-\theta_{1}\right)}{\cos\theta_{1}\cos\theta_{2}}=1+\tan\theta_{1}\tan\theta_{2}.\end{split}$
(46)

By rearranging (46), we have

$\theta_{v}-\frac{\theta_{tb}+\theta_{it}}{2}=-\theta_{1}=\arctan\left[\left(1-\frac{\mu_{d}}{\mu_{r}}\right)\cot\theta_{2}\right],$
(47)

which completes the proof of (8). $\blacksquare$

## References

  * [1] F. Liu, C. Masouros, A. P. Petropulu, H. Griffiths, and L. Hanzo, “Joint radar and communication design: Applications, state-of-the-art, and the road ahead,” _IEEE Trans. Commun._ , vol. 68, no. 6, pp. 3834–3862, Feb. 2020.
  * [2] L. Xie, P. Wang, S. Song, and K. B. Letaief, “Perceptive mobile network with distributed target monitoring terminals: Leaking communication energy for sensing,” _IEEE Trans. Wireless Commun._ , Jun. 2022.
  * [3] L. Xie, S. Song, Y. C. Eldar, and K. B. Letaief, “Collaborative sensing in perceptive mobile networks: Opportunities and challenges,” _arXiv preprint arXiv:2205.15805_ , 2022.
  * [4] F. Liu, W. Yuan, C. Masouros, and J. Yuan, “Radar-assisted predictive beamforming for vehicular links: Communication served by sensing,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 11, pp. 7704–7719, Aug. 2020.
  * [5] S. Sen, “Adaptive OFDM radar waveform design for improved micro-Doppler estimation,” _IEEE Sens. J._ , vol. 14, no. 10, pp. 3548–3556, Jul. 2014.
  * [6] J. Lim, S.-R. Kim, and D.-J. Shin, “Two-step Doppler estimation based on intercarrier interference mitigation for OFDM radar,” _IEEE Antennas Wirel. Propag. Lett._ , vol. 14, pp. 1726–1729, Apr. 2015.
  * [7] S. H. Dokhanchi, B. S. Mysore, K. V. Mishra, and B. Ottersten, “A mmwave automotive joint radar-communications system,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 55, no. 3, pp. 1241–1260, Feb. 2019.
  * [8] X. Yu, D. Xu, D. W. K. Ng, and R. Schober, “IRS-assisted green communication systems: Provable convergence and robust optimization,” _IEEE Trans. Commun._ , vol. 69, no. 9, pp. 6313–6329, Jun. 2021.
  * [9] X. Yu, V. Jamali, D. Xu, D. W. K. Ng, and R. Schober, “Smart and reconfigurable wireless communications: From IRS modeling to algorithm design,” _IEEE Wireless Commun._ , vol. 28, no. 6, pp. 118–125, Dec. 2021.
  * [10] J. He, H. Wymeersch, T. Sanguanpuak, O. Silvén, and M. Juntti, “Adaptive beamforming design for mmwave RIS-aided joint localization and communication,” in _Proc. IEEE Wireless Commun. Netw. Conf. Workshops (WCNCW)_. IEEE, Jun. 2020, pp. 1–6.
  * [11] B. Matthiesen, E. Björnson, E. De Carvalho, and P. Popovski, “Intelligent reflecting surface operation under predictable receiver mobility: A continuous time propagation model,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 2, pp. 216–220, Sep. 2020.
  * [12] Z. Huang, B. Zheng, and R. Zhang, “Transforming fading channel from fast to slow: Intelligent refracting surface aided high-mobility communication,” _IEEE Trans. Wireless Commun._ , Dec. 2021.
  * [13] K. Keykhosravi, M. F. Keskin, G. Seco-Granados, P. Popovski, and H. Wymeersch, “RIS-enabled SISO localization under user mobility and spatial-wideband effects,” _IEEE J. Sel. Topics Signal Process._ , May 2022.
  * [14] P. Stoica and K. Sharman, “Maximum likelihood methods for direction-of-arrival estimation,” _IEEE Trans. Acoustics, Speech, Signal Process._ , vol. 38, no. 7, pp. 1132–1143, Jul. 1990.
  * [15] R. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE Trans. Antennas Propag._ , vol. 34, no. 3, pp. 276–280, Mar. 1986.
  * [16] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” _IEEE Trans. Acoustics, Speech, Signal Process._ , vol. 37, no. 7, pp. 984–995, Jul. 1989.
  * [17] L. Xie, Z. He, J. Tong, and W. Zhang, “A recursive angle-Doppler channel selection method for reduced-dimension space-time adaptive processing,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 56, no. 5, pp. 3985–4000, Apr. 2020.


# Light Weight Color Image Warping with Inter-Channel Information

###### Abstract

Image warping is a necessary step in many multimedia applications such as
texture mapping, image-based rendering, panorama stitching, image resizing and
optical flow computation etc. Traditionally, color image warping interpolation
is performed in each color channel independently. In this paper, we show that
the warping quality can be significantly enhanced by exploiting the cross-
channel correlation. We design a warping scheme that integrates intra-channel
interpolation with cross-channel variation at very low computational cost,
which is required for interactive multimedia applications on mobile devices.
The effectiveness and efficiency of our method are validated by extensive
experiments.

Index Terms—  Image warping, inter-channel correlation, Laplacian filtering,
image enhancement

##  1 Introduction

Image warping is fundamental for a variety of multimedia applications such as
image resizing (e.g, [1]), texture mapping (e.g., [2]), image-based rendering
(e.g., [3]), panorama stitching (e.g., [4]), stereo reconstruction (e.g.,
[5]), optical flow computation (e.g., [6]), to name but a few. Image warping
can be briefly described as transforming a source image $I_{1}$ into a target
image $I_{2}$ under a geometric mapping $\mathbf{H}$, which can be linear,
affine, perspective, non-parametric etc.[7]. Particularly, assuming $I_{2}$
has $M$ rows and $N$ columns, each pixel $(i,j)$ in the target image lattice
$[1\quad M]\times[1\quad N]$ has a correspondence image point
$[x,y]=\mathbf{H}^{-1}([i,j])$ in the source image, from which the intensity
value $I_{1}[x,y]$ is assigned to $I_{2}[i,j]$. The central problem of image
warping is to estimate $I_{1}[x,y]$ from pixels in the vicinity of $[x,y]$,
which are generally not on the integer grid, and hence the true intensity
$I_{1}[x,y]$ is not immediately available.

The most popular estimation techniques for image warping in real practice are
nearest neighbour, bilinear and bicubic (spline) interpolation schemes [8],
which do not adapt well to irregular image regions. Therefore, academic
research effort has been mainly dedicated to edge guided interpolation, such
as anisotropic filtering (e.g., [9], [10]) and variational regularization
(e.g., [11]). Particularly on the sub-topic of Signal Image Super-resolution,
i.e., in the special case that $\mathbf{H}$ degenerates to a scaling function,
mathematical models such as compressive sensing (e.g., [12]), dictionary
learning (e.g., [13]) and Convolutional Neural Networks (CNN, e.g., [1]) have
also been intensively investigated.

Although numerous image warping schemes have been proposed, they are mostly
designed for monochrome images. This means that on warping color images, these
methods should be applied to each channel independently (e.g., [14]) or to the
luminance channel only (e.g.[15]), without exploiting the cross-channel
correlation. The only exception is end-to-end CNNs that take RGB images as
input. In this situation, the color correlation is captured during the
training process, but the learned correlation is CNN specific rather than
being general. In contrast, cross-channel correlation is thoroughly studied
for Image Demosaicking, which recovers the full RGB images from subsampled
color channels, with only one primary color component available at each pixel
[16]. Although each channel can be recovered by merely using the intra-channel
interpolation, it is commonly known that inter-channel correlation information
can significantly improve the demosaicking accuracy.

In this paper, we integrate intra-channel image warping with cross-channel
correlation, and show that our approach significantly improve the warping
accuracy at very low cost. The proposed integration scheme is inspired by the
Malvar-He-Cutler High Quality Linear Interpolation (HQLI) demosaicking method
[17], which is probably the fastest demosaicking algorithm in the literature
(except the Nearest Neighbour and Bilinear interpolation), but its trade-off
between accuracy and speed is superior to many sophisticated algorithms [18].
We first examine the potential of cross-channel correlation for image warping,
by applying HQLI to postprocess color images upsampled with intra-channel
interpolation. We then generalize this combination to formulate general color
image warping. The performance of our method is assessed in the scenario of
image super-resolution, as it is the most active sub-area of image warping.
Despite its very low cost, the proposed algorithm achieves comparable accuracy
to state of the art algorithms at about $40$ times faster speed. Compared to
CNN-based methods, our model has the flexibility of resizing the image to
arbitrary size. These features are desirable for real time interactive
applications on multimedia devices that have limited computation power (e.g.,
mobile smart phones).

##  2 Inter-Channel Correlation Gain

Image warping is generally implemented by backward mapping. Given $\mathbf{H}$
the geometric mapping from the source image $I_{1}$ to the target image
$I_{2}$, and $\mathbf{H}^{-1}$ be its inverse mapping, $\mathbf{H}^{-1}$ maps
the coordinate frame of $I_{2}$ back to the coordinate frame of $I_{1}$. Image
$I_{2}$ is generated in such a way that, for each pixel $(i,j)$ of $I_{2}$,
the intensity value $I_{2}[i,j]$ is copied from $I_{1}[x,y]$, where

$\displaystyle\left[\begin{array}[]{c}x\\\
y\end{array}\right]=\mathbf{H}^{-1}\left[\begin{array}[]{c}i\\\
j\end{array}\right].$ (5)

As $x$ and $y$ are generally not integers, $I_{1}[x,y]$ has to be estimated
from the surrounding pixels. Let $m=\left\lfloor x\right\rfloor$ and
$n=\left\lfloor y\right\rfloor$ be the floor integer parts of $x$ and $y$, and
define $s=x-m$, $t=y-n$. For example, the bilinear warping scheme estimates
$I_{2}[i,j]$ by

$\displaystyle\begin{split}\hat{I}_{1}[x,y]&=(1-s)(1-t){I}_{1}[m,n]+st{I}_{1}[m+1,n+1]\\\
&+s(1-t){I}_{1}[m+1,n]+(1-s)t{I}_{1}[m,n+1].\end{split}$ (6)

Today, multimedia images are usually chromatic rather than being grayscale.
Therefore, $I_{1}$ and $I_{2}$ generally have RGB three channels. Typically,
bilinear warping for color images is to replace $I_{1}$ in Eq.6 with its color
channels $R_{1}$, $G_{1}$ and $B_{1}$, to estimate $R_{2}[i,j]$, $G_{2}[i,j]$
and $B_{2}[i,j]$, without cross-channel interaction.

The importance of cross-channel correlation is well investigated for the
problem of Image Demosaicking, entailed in current consumer-level digital
cameras, which commonly use Color Filter Arrays (CFA) to sample one of the RGB
components for each pixel. Fig.1 shows an example of the Bayer pattern CFA,
adopted for most digital cameras.

(a) RGGB

(b) BGGR

(c) GRBG

(d) GBRG

Fig. 1: An illustration of the Bayer pattern for CFA.

Define $\mathcal{R}$, $\mathcal{G}$ and $\mathcal{B}$ to be the sets that
collect the pixels where the original red, green and blue values are captured
respectively. Malvar-He-Cutler propose estimating the missing color values by
a weighted combination of the local intra-channel average and the Laplacian of
another channel [17]. For example, at a pixel $[m,n]\in\mathcal{R}$, the green
value $G[m,n]$ is estimated by

$G[m,n]\approx\bar{G}+\alpha\Delta R[m,n],$ (7)

where constant $\alpha$ is learned from training images, $\bar{G}$ stands for
the average of $G[m,n-1]$, $G[m,n+1]$, $G[m-1,n]$ and $G[m+1,n]$, and

$\displaystyle\begin{split}\Delta
R[m,n]=R[m,n]-\frac{1}{4}&\left(R[m+2,n]+R[m-2,n]\right.\\\
&\left.+R[m,n-2]+R[m,n+2]\right).\\\ \end{split}$ (8)

The other missing color values are recovered in the same fashion. This
demosaicking method, namely HQLI, adds second-order details (the Laplacian) of
one channel to the average of another channel, thereby it performs image
enhancement.

To examine whether the inter-channel information supplements intra-channel
recovery in image warping. We conduct a preliminary experiment using image
super-resolution as an example application. We upsample test images by
bilinear and bicubic interpolation, and then refine the upsampled image values
by sequential HQLI demosaicking. In particular, we first apply the GRBG CFA
pattern to filter the upsampled color image, and estimate the filtered-out
color values by the remaining color values. Subsequently, we apply the RGGB
CFA to the demosaicked full RGB image, and perform the second-round
demosaicking. We then change the CFA pattern to BGGR, and repeat HQLI
demosaicking for the final round. In this way, each color component of a pixel
is refined by the intensity values in the local neighbourhood of the other two
channels. Table1 compares the reconstruction accuracy, measured by Peak Signal
to Noise Ratio (PSNR), before and after demosaicking (i.e., the cross-channel
information) on benchmark datasets Set5 [19], Set14 [20], BSD100 [21] and
Urban100 [22]. It can be seen that the inter-channel correlation introduced by
demosaicking notably improves the intra-channel interpolation, on all test
datasets. Numerically, the improvement is upto $1.57$ dB and $1.09$ dB for
bilinear and bicubic intra-channel interpolation respectively.

Table 1: The average PSNR before and after refining bilinear and bicubic upsampling by sequential HQLI demosaicking, tested on benchmark datasets. Dataset | Bilinear |  | Bilinear
---
+HQLI
Bicubic |  | Bicubic
---
+HQLI
Set5 | 28.86 | 30.43 | 28.64 | 29.73
Set14 | 26.55 | 27.50 | 26.23 | 26.85
BSD100 | 26.44 | 26.89 | 26.02 | 26.16
Urban100 | 23.57 | 24.12 | 23.17 | 23.45

Note that to simulate the real world applications, here the low resolution
images are obtained by Nearest Neighbour downsampling of the high resolution
test images. This is different from the common SISR algorithm evaluation
routine, which utilizes Bicubic downsampling. The Nearest Neighbour sampling
is known to suffer severe aliasing and blurring, whereas Bicubic downsampling
implicitly performs low pass filtering before sampling, thereby largely avoids
such artifacts [16]. However, using Bicubic downsampling also implicitly
assumes that the high resolution images are known, which is not the case for
real applications. Therefore we only apply bicubic downsampling in the
comparison experiments for fair conditions. Due to the coarse downsampling,
bicubic warping has lower accuracy than bilinear warping, which is seemingly
counter intuitive.

##  3 Cross-Channel Color Image Warping

###  3.1 General Formulation

Based on the analysis of the gain achieved by the inter-channel information,
we formulate general color image warping by

$\displaystyle G_{1}[x,y]$ $\displaystyle\approx$
$\displaystyle\tilde{G}_{1}[x,y]+\omega_{g,r}\Delta
R_{1}[x,y]+\omega_{g,b}\Delta B_{1}[x,y]$ (9a) $\displaystyle R_{1}[x,y]$
$\displaystyle\approx$ $\displaystyle\tilde{R}_{1}[x,y]+\omega_{r,g}\Delta
G_{1}[x,y]+\omega_{r,b}\Delta B_{1}[x,y]$ (9b) $\displaystyle B_{1}[x,y]$
$\displaystyle\approx$ $\displaystyle\tilde{B}_{1}[x,y]+\omega_{b,g}\Delta
G_{1}[x,y]+\omega_{b,r}\Delta R_{1}[x,y],$ (9c)

where $\tilde{G}_{1}[x,y]$, $\tilde{R}_{1}[x,y]$ and $\tilde{B}_{1}[x,y]$ can
be any intra-channel estimation of $G_{1}[x,y]$, $R_{1}[x,y]$ and
$B_{1}[x,y]$; also, $\Delta G_{1}[x,y]$, $\Delta R_{1}[x,y]$ and $\Delta
B_{1}[x,y]$ can be estimated by interpolating ${\Delta
G_{1}[k,l]}_{[k,l]\in\Omega(m,n)}$, ${\Delta
R_{1}[k,l]}_{[k,l]\in\Omega(m,n)}$ and ${\Delta
B_{1}[k,l]}_{[k,l]\in\Omega(m,n)}$ respectively, where $\Omega[m,n]$ stands
for the local neighbourhood of pixel $[m,n]$. The involved interpolation
schemes can be either linear or non-linear, isotropic or anisotropic,
depending on the particular requirements of the applications at hand. The
weights $\omega_{(\cdot,\cdot)}$are to be learned from training images, as
described in Section3.2.

###  3.2 Parameter Training

We illustrate the weight learning process using the Single Image Super
Resolution (SISR) application as an example, for which the geometric mapping
$\mathbf{H}$ and its inverse mapping $\mathbf{H}^{-1}$ are

$\displaystyle\mathbf{H}=\left[\begin{array}[]{cc}S&0\\\
0&S\end{array}\right]\quad\mathbf{H}^{-1}=\left[\begin{array}[]{cc}\frac{1}{S}&0\\\
0&\frac{1}{S}\end{array}\right],$ (14)

where $S$ in a real number larger than $1$. Note that $S$ is not restricted to
be an integer. We downsample the groundtruth images as the source images
$\mathbf{I}_{1}$ (i.e., the data) by bicubic interpolation, following the
traditional experimental setting in the literature of SISR, as the learned
weights are to be used for algorithm evaluation in comparison with state of
the art. We use the groundtruth images (taken from benchmark training dataset
BSD200 [21]) as target images $\mathbf{I}_{2}$ (i.e., the labels). For each
pixel $[i,j]$ in the an image $\mathbf{I}_{2}$, the RGB values of its
correspondence $(x,y)$ in $\mathbf{I}_{1}$ thus are known in the training
process. $\tilde{G}_{1}[x,y]$, $\tilde{R}_{1}[x,y]$ and $\tilde{B}_{1}[x,y]$,
as well as $\Delta G_{1}[x,y]$, $\Delta R_{1}[x,y]$ and $\Delta B_{1}[x,y]$
are computed on the pixels in each source image by interpolation function $f$.
Subsequently, Eq.9 can be rewritten as linear systems of the weights
$\omega_{(\cdot,\cdot)}$. For example,

$\displaystyle\begin{split}&\left[\begin{array}[]{cc}\Delta
R_{1}[x_{1},y_{1}]&\Delta B_{1}[x_{1},y_{1}]\\\ \vdots&\vdots\\\ \Delta
R_{1}[x_{K},y_{K}]&\Delta
B_{1}[x_{K},y_{K}]\end{array}\right]\left[\begin{array}[]{c}\omega_{g,r}\\\
\omega_{g,b}\end{array}\right]\\\
&=\left[\begin{array}[]{c}\tilde{G}_{1}[x_{1},y_{1}]-G_{2}[i_{1},j_{1}]\\\
\vdots\\\
\tilde{G}_{1}[x_{K},y_{K}]-G_{2}[i_{K},j_{K}]\end{array}\right]\end{split}$
(15)

for Eq.9a. Here integer $K=10000$ is the total number of randomly selected
training pixels, indexed by the subscripts of $x$, $y$, $i$ and $j$. Eq.9b and
Eq.9c can be re-organized similarly. The Mean Square Error (MSE) solutions or
Mean Absolute Error (MAE) solutions to the linear systems can be easily solved
in closed-form. This yields the weights $\omega_{g,r}$,
$\omega_{g,b}$,$\omega_{r,g}$,$\omega_{r,b}$,$\omega_{b,r}$ and
$\omega_{b,g}$.

Seemingly, the learned weight values should be specific to the scaling factor
$S$ and interpolation function $f$. However, we observe that weights learned
for large $S$ actually generalizes well to small $S$. Therefore, we suggest
conducting the training for a large scaling factor (e.g., $S=4$), and apply
the learned weights directly to applications whose scaling factor is smaller.

Table2 lists the weights we learned from Eq.15 by linear regression for $S=4$,
with the interpolation function $f$ being Bilinear, Bicubic and Lanczos (with
$5\times 5$ footprint).

Table 2: The integration weights learned by linear regression for $S=4$ and various interpolation functions $f$. | Bilinear | Bicubic | Lanczos
---|---|---|---
$\omega_{g,r}$ | 0.094 | 0.045 | 0.032
$\omega_{g,b}$ | 0.119 | 0.064 | 0.041
$\omega_{r,g}$ | 0.195 | 0.096 | 0.058
$\omega_{r,b}$ | 0.008 | 0.010 | 0.015
$\omega_{b,g}$ | 0.180 | 0.089 | 0.054
$\omega_{b,r}$ | -0.003 | 0.003 | 0.008

##  4 Experimental Results

We conduct two sets of experiments to evaluate the proposed color image
warping technique. The first experiment measures the SISR accuracy before and
after incorporating the cross-channel terms to intra-channel warping, at
various scaling factors and using different real time interpolation functions.
The second experiments compares our technique with the state of the art SISR
method Global Regression (GR) proposed by Timofte-De Smet-Van Gool in [23],
using the released Matlab source code. The GR method is the fastest SISR
method proposed in recent years to the best of our knowledge, and hence is
chosen for comparison in our experiment. It should be mentioned that, the
weight GR is significantly higher than our linear isotropic scheme used for
testing.

For comparison in equal conditions, the low resolution test images $I_{1}$ are
obtained by bicubicly downsampling the ground truth image $I_{2}$ at sampling
factor $S$. That is, the low resolution images have
$\left\lfloor\frac{M}{S}\right\rfloor$ rows and
$\left\lfloor\frac{N}{S}\right\rfloor$ columns. Laplacians are computed at
each pixel in each channel of the low resolution images. Both the Laplacian
maps and low resolution images are upsampled by the same interpolation
function $f$. In our experiments, $S$ is set to be 2, 3 or 4; and $f$ is set
to be Bilinear, Bicubic or Lanczos (of kernel size $5\times 5$).

Evaluation is carried out on benchmark datasets Set5, Set14, BSD100 and
Urban100. Beside the PSNR accuracy measure, we also report the Structure
Similarity (SSIM) measure [24] and running time (in seconds) of the SISR
algorithms. All experiments are conducted on an Intel Core i7-7700 3.60GHz CPU
with 8GB RAM, using Matlab.

###  4.1 Improvement on Popular Warping Methods

We compare the performance of popular real time warping methods and the
proposed method. Bilinear warping probably has the best trade-off between
performance and speed in the literature, as pointed out by Zitová-Flusser [8].
Hence it is recommended for video and animation applications, such as frame
registration for motion estimation in OpenCV and texture mapping in OpenGL.
Bicubic warping is another frequently adopted option, especially for still
image processing. For example, image editing softwares Photoshop and GIMP
employ bicubic interpolation for image resizing and perspective view
rectification.

Table3 presents the PSNR, SSIM and running time averaged over Set5 for various
$S$ and $f$. In this table, the plain intra-channel warping is referred to as
“Independent”, and the proposed inter-channel color warping is referred to as
“Correlated”. In all situations, the proposed color warping scheme achieves
evidential improvement, either in terms of PSNR or SSIM, in real time
computation. The most remarkable quantitative improvement is obtained for
bilinear interpolation, the PSNR values of which are increased by 1.58dB,
0.77dB and 0.98dB for 2, 3 and 4 times upscaling respectively. It can be
observed that the improvement decreases with the complexity of the
interpolation function. However, in the case of Lanczos interpolation (the
most complex function here), the accuracy of our color warping scheme for
$2\times$ upscaling is probably sufficient in many applications, and its
computation is merely in several milliseconds.

Table 3: PSNR, SSIM and computation time of plain intra-channel warping and the proposed color warping, averaged over Set5, for various scaling factor $S$ and interpolation function $f$. $f$ | $S$ | Independent | Correlated
---|---|---|---
PSNR | SSIM | Time | PSNR | SSIM | Time
Bilinear | $\times 2$ | 30.42 | 0.9566 | 0.0013 | 32.00 | 0.9673 | 0.0062
$\times 3$ | 27.78 | 0.9274 | 0.0012 | 28.55 | 0.9368 | 0.0052
$\times 4$ | 25.84 | 0.8966 | 0.0011 | 26.82 | 0.9096 | 0.0056
Bicubic | $\times 2$ | 31.81 | 0.9657 | 0.0020 | 32.37 | 0.9694 | 0.0075
$\times 3$ | 28.63 | 0.9365 | 0.0019 | 29.07 | 0.9413 | 0.0065
$\times 4$ | 26.70 | 0.9087 | 0.0018 | 27.13 | 0.9139 | 0.0068
Lanczos | $\times 2$ | 32.43 | 0.9687 | 0.0020 | 32.74 | 0.9706 | 0.0074
$\times 3$ | 29.05 | 0.9402 | 0.0018 | 29.29 | 0.9427 | 0.0063
$\times 4$ | 27.07 | 0.9129 | 0.0017 | 27.33 | 0.9157 | 0.0067

###  4.2 Comparison with State of the Art

To the best of our knowledge, among the state of the art SISR algorithms that
have Matlab source code released, GR [23] is the fastest. We compare the
accuracy and speed of our method with GR. In this experiment, we fix the
interpolation function $f$ to be Lanczos. Table 4 shows the comparison, and
Fig.2 shows a visual assessment. Relative to GR, the PSNR of our method is
generally around $0.2$dB lower for various scaling factors. However, on BSD100
the accuracy performance are very close. This might be attributed to the
integration weights learned from BSD200, and suggests the potential for higher
performance by using more sophisticated training strategy. It should be noted
that the proposed warping model has merely 6 parameters, and the computation
is tens of times faster.

Table 4: Comparison of PSNR and computation time of GR [23] the proposed color warping (with the interpolation function fixed to be Lanczos), averaged over four benchmark datasets, for various scaling factor $S$. Dataset | $S$ | GR | Ours
---|---|---|---
PSNR | Time | PSNR | Time
Set5 | $\times 2$ | 33.04 | 0.34 | 32.74 | 0.0074
$\times 3$ | 29.50 | 0.22 | 29.29 | 0.0063
$\times 4$ | 27.49 | 0.19 | 27.33 | 0.0067
Set14 | $\times 2$ | 29.85 | 0.67 | 29.62 | 0.0138
$\times 3$ | 26.92 | 0.46 | 26.75 | 0.0133
$\times 4$ | 25.21 | 0.38 | 25.08 | 0.0124
BSD100 | $\times 2$ | 28.90 | 0.45 | 28.97 | 0.0073
$\times 3$ | 26.38 | 0.32 | 26.30 | 0.0070
$\times 4$ | 25.03 | 0.26 | 24.97 | 0.0070
Urban100 | $\times 2$ | 26.30 | 2.31 | 26.20 | 0.0480
$\times 3$ | 23.69 | 1.46 | 23.50 | 0.0427
$\times 4$ | 22.28 | 1.24 | 22.08 | 0.0428

(a) Reference/PSNR

(b) GR/25.51 dB

(c) Our/25.37 dB

Fig. 2: ’Zebra’ image from Set14 with upscaling *3

##  5 Conclusion

In this paper, we have proposed a light weight color image warping technique,
which integrates intra-channel interpolation with cross-channel details.
Extensive experiments on four benchmark datasets have validated that, the
proposed technique substantially improves the most popular image warping
algorithms, in very simple form and at trivial computational cost. Our
algorithm has also been shown to push the performance of traditional Lanczos
upsampling scheme to be comparable with state of the art methods, while being
real time.

The proposed technique can be readily applied to interactive multimedia
applications that require both fast image warping and realistic visual
effects, such as Computer Gaming and Image Editing. As our integration of
cross-channel information is of light weight, it is suitable for multimedia
devices with limited computation power. It should be noted that our baseline
interpolation, integration formulation and weight training are all very basic.
Actually, each step can be extended for higher performance, if the computation
cost budget allows. Moreover, our method can also be used as an high quality
initialization to many CNN-based methods.

## References

  * [1] Christian Ledig, Lucas Theis, Ferenc Huszár, Jose Caballero, Andrew Cunningham, Alejandro Acosta, Andrew P Aitken, Alykhan Tejani, Johannes Totz, Zehan Wang, et al.,  “Photo-realistic single image super-resolution using a generative adversarial network.,”  in IEEE Conference on Computer Vision and Pattern Recognition, 2017\.
  * [2] Mason Woo, Jackie Neider, Tom Davis, and Dave Shreiner,  OpenGL programming guide: the official guide to learning OpenGL, version 1.2,  Addison-Wesley Longman Publishing Co., Inc., 1999.
  * [3] Leonard McMillan,  “Image-based rendering using image warping,”  in COM-PUTER GRAPHICS Proceedings, Annual Conference Sefies, ACM SIG-GRAP, New Orleans, Louisiana, 2009, pp. 123–125.
  * [4] Kaiming He, Huiwen Chang, and Jian Sun,  “Rectangling panoramic images via warping,”  ACM Transactions on Graphics (TOG), vol. 32, no. 4, pp. 79, 2013\.
  * [5] Michael Bleyer, Margrit Gelautz, Carsten Rother, and Christoph Rhemann,  “A stereo approach that handles the matting problem via image warping,”  in Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, 2009, pp. 501–508.
  * [6] Jean-Yves Bouguet,  “Pyramidal implementation of the affine lucas kanade feature tracker description of the algorithm,”  Intel Corporation, vol. 5, no. 1-10, pp. 4, 2001.
  * [7] Chris A Glasbey and Kantilal Vardichand Mardia,  “A review of image-warping methods,”  Journal of applied statistics, vol. 25, no. 2, pp. 155–171, 1998\.
  * [8] Barbara Zitova and Jan Flusser,  “Image registration methods: a survey,”  Image and vision computing, vol. 21, no. 11, pp. 977–1000, 2003\.
  * [9] Sebastiano Battiato, Giovanni Gallo, and Filippo Stanco,  “A locally adaptive zooming algorithm for digital images,”  Image and vision computing, vol. 20, no. 11, pp. 805–812, 2002\.
  * [10] Xin Li and Michael T Orchard,  “New edge-directed interpolation,”  IEEE transactions on image processing, vol. 10, no. 10, pp. 1521–1527, 2001.
  * [11] Hao Jiang and Cecilia Moloney,  “A new direction adaptive scheme for image interpolation,”  in International Conference on Image Processing, 2002, vol. 3.
  * [12] Jianchao Yang, John Wright, Thomas S Huang, and Yi Ma,  “Image super-resolution via sparse representation,”  IEEE transactions on image processing, vol. 19, no. 11, pp. 2861–2873, 2010.
  * [13] Jianchao Yang, Zhaowen Wang, Zhe Lin, Scott Cohen, and Thomas Huang,  “Coupled dictionary training for image super-resolution,”  IEEE transactions on image processing, vol. 21, no. 8, pp. 3467–3478, 2012.
  * [14] Dan Su and Philip Willis,  “Image interpolation by pixel-level data-dependent triangulation,”  Computer Graphics Forum, vol. 23, no. 2, pp. 189–201, 2004.
  * [15] Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang,  “Learning a deep convolutional network for image super-resolution,”  in European conference on computer vision, 2014, pp. 184–199.
  * [16] R. Szeliski,  Computer vision: algorithms and applications,  Springer Science & Business Media, 2010.
  * [17] Henrique S Malvar, Li-wei He, and Ross Cutler,  “High-quality linear interpolation for demosaicing of bayer-patterned color images,”  in Acoustics, Speech, and Signal Processing, 2004. Proceedings.(ICASSP’04). IEEE International Conference on. IEEE, 2004, vol. 3, pp. iii–485.
  * [18] Yan Niu, Jihong Ouyang, Wanli Zuo, and Fuxin Wang,  “Low cost edge sensing for high quality demosaicking,”  IEEE Transactions on Image Processing, 2018,  Early Access.
  * [19] Marco Bevilacqua, Aline Roumy, Christine Guillemot, and Marie-Line Alberi Morel,  “Low-complexity single-image super-resolution based on nonnegative neighbor embedding,”  in British Machine Vision Conference, 2012.
  * [20] Roman Zeyde, Michael Elad, and Matan Protter,  “On single image scale-up using sparse-representations,”  in International conference on curves and surfaces, 2010, pp. 711–730.
  * [21] David Martin, Charless Fowlkes, Doron Tal, and Jitendra Malik,  “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,”  in International Conference on Computer Vision, 2001, vol. 2, pp. 416–423.
  * [22] Jia-Bin Huang, Abhishek Singh, and Narendra Ahuja,  “Single image super-resolution from transformed self-exemplars,”  in IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 5197–5206.
  * [23] Radu Timofte, Vincent De Smet, and Luc Van Gool,  “Anchored neighborhood regression for fast example-based super-resolution,”  in International Conference on Computer Vision, 2013, pp. 1920–1927.
  * [24] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli,  “Image quality assessment: from error visibility to structural similarity,”  IEEE transactions on image processing, vol. 13, no. 4, pp. 600–612, 2004.


László Tóth1, Amin Honarmandi Shandiz1, Gábor Gosztolya1,2, Tamás Gábor Csapó3
---

# Adaptation of Tongue Ultrasound-Based Silent Speech Interfaces
Using Spatial Transformer Networks

###### Abstract

Thanks to the latest deep learning algorithms, silent speech interfaces (SSI)
are now able to synthesize intelligible speech from articulatory movement data
under certain conditions. However, the resulting models are rather speaker-
specific, making a quick switch between users troublesome. Even for the same
speaker, these models perform poorly cross-session, i.e. after dismounting and
re-mounting the recording equipment. To aid quick speaker and session
adaptation of ultrasound tongue imaging-based SSI models, we extend our deep
networks with a spatial transformer network (STN) module, capable of
performing an affine transformation on the input images. Although the STN part
takes up only about 10% of the network, our experiments show that adapting
just the STN module might allow to reduce MSE by 88% on the average, compared
to retraining the whole network. The improvement is even larger (around 92%)
when adapting the network to different recording sessions from the same
speaker.

Index Terms: silent speech interface, ultrasound tongue imaging, speaker
adaptation, spatial transformer network

##  1 Introduction

The goal of articulatory-to-acoustic mapping (AAM) is to synthesize speech
from data capturing the movement of the articulatory organs. This is the
purpose of the giving the background to the development of 'Silent Speech
Interface' systems (SSI [1, 2]). Ideally, these interfaces would record the
articulation and synthesize speech based on the movement of the organs –
without the user of the device actually producing any sound. The typical input
of AAM can be a video of the lip movements [3, 4, 5, 6, 7, 8], ultrasound
tongue imaging (UTI) [3, 9, 10, 11, 12, 13, 14, 15, 16, 17], or several other
modalities (e.g., MRI, EMA, PMA, EOS, radar, multimodal, etc.). All of the
articulatory tracking devices are highly sensitive to 1) the alignment of the
recording equipment across sessions, 2) the actual speaker's anatomy. For
example, in the case of ultrasound recordings, the probe fixing headset has to
be remounted onto the speaker for each recording session. This inevitably
causes the recorded ultrasound videos to become misaligned between each
recording session [18]. Moreover, there are large individual differences
across speakers, so even a system trained on the data of several speakers may
still perform poorly for a new speaker.

There have already been several cross-session and cross-speaker studies, of
which we mention only those related to imaging. To handle the session
dependency of UTI-based synthesis, Gosztolya et al. used data from different
sessions [19]. Ribeiro et al. reported that, for a speaker-independent system,
unsupervised model adaptation can improve the results for silent speech [20].
In a multi-speaker framework, Shandiz et al. experimented with x-vectors
features extracted from speakers, leading to a marginal improvement in the
spectral estimation step [17]. Zhang et al. evaluated UTI and lip video based
unconstrained multi-speaker voice recovery with a transfer learning strategy
and encoder-decoder architecture [21]. There have been more studies on multi-
speaker lip-to-speech synthesis [5, 6, 7, 8]. One of the first papers that
could produce intelligible speech for unseen speakers was based on WGAN with
new additional critics and losses [5]. Another study proposed speaker
disentanglement by inputting speaker identities or embeddings to the DNN [6].
An end-to-end ResNet-based model was claimed to outscore previous approaches
on unseen speakers [7], while Wang et al. applied cross-modal knowledge
transfer and voice conversion to generate good quality speech with high
naturalness [8].

Most of the above approaches hope to solve speaker sensitivity simply by
acquiring articulatory training data from a large quantity of speakers. In
this study, we experiment with a direct adaptation of an UTI-based SSI network
to the actual speaker or session. To avoid the need for a full retraining, we
extend our network with a spatial transformer network (STN) module and retrain
only this module during the adaptation step. The STN learns an affine
transformation on the input images, and our assumption is that this
transformation should mostly be able to compensate for the misalignment of the
recording device, and to a certain extent also for the inter-speaker
differences.

##  2 The UTI-to-Speech framework

Approaches for articulatory-to-acoustic mapping similar to ours were described
several times [15, 17, 22], so we just give a brief overview here (see also
Fig. 1). The input to our system is a sequence of ultrasound tongue images,
which record the movement of the articulators at a relatively high frame rate.
Our goal is to estimate the speech signal produced during articulation, but
instead of producing a speech signal, we estimate only a mel-spectrogram at
the network output, and we apply a neural vocoder to convert the spectrogram
to speech. The main advantage of this approach is that the mel-spectrum is a
very dense representation, which is easier to estimate than the speech signal
itself, and we can apply large pre-trained networks for the synthesis step,
such as WaveGlow [23]. Our training data consists of precisely time-aligned
pairs of ultrasound videos and speech signals, so the ultrasound-to-
spectrogram mapping can be performed on a simple frame-by-frame basis,
converting each ultrasound image to a spectral vector [11, 15]. While this
simple arrangement already performs reasonably well, significant improvement
can be achieved by involving the input context, that is, by using a block of
video frames as input instead of just one image. Several network architectures
have been proposed to process 3D blocks of input data, for video processing in
general [24, 25, 26], and for ultrasound input in particular [12, 14, 27, 16,
28]. In the experimental section we will experiment both with 2D and 3D
Convolutional Neural Networks (CNNs) for the mapping task. The problem could
also be addressed even in the lack of aligned training data using encoder-
decoder networks [20, 21] or video transformers [29, 30].

Figure 1: Illustration of the network architecture used, where FC denotes the
Fully-Connected layer.

##  3 The Spatial Transformer Network

The concept of the Spatial Transformer Network (STN) was motivated by the fact
that in image recognition the actual input is frequently shifted, rotated or
scaled compared to the training data, due to changes in camera angle,
distance, and other factors [31]. Although CNNs are somewhat invariant to
translations, their overall flexibility would greatly improve by introducing a
dynamic mechanism that can spatially transform the input image by an
appropriate transformation before classification. The STN is a network module
that can be inserted into a CNN architecture at any point, but typically it is
applied between the input and the first network layer. In this arrangement,
the STN performs an affine transformation on the input image, and returns a
transformed image of the same size. As it manipulates only the input, it can
be combined with practically any type of classification or regression network.
The affine transform defined by the STN is quite powerful and includes
translation, scaling, rotation, shearing and cropping as special cases.

The STN consists of three main parts, namely the localization network, the
grid generator and the sampler [31]. The grid generator and the sampler
together are responsible for executing the affine transformation defined by
the parameter $\theta$, which consists of 6 components in the case of a 2D
input (see Fig. 1).

The grid generator and the sampler are differentiable. This is vital for
propagating the error back to the third component, the localization network,
which learns to estimate the $\theta$ parameters, conditioned on the actual
input image. The localization network can take any form, but its uppermost
layer must be a regression layer to produce the $\theta$ values.

While we explained the STN concept assuming 2D images, the whole idea can be
naturally extended to 3D data blocks [31, 32]. The 3D variant has already been
used in visual speech recognition (e.g., lip reading) by Yu and Wang [33]. We
discuss the options for a 3D extension in Section 4.3.

Although the original concept uses the STN to decrease the variance of the
training data by transforming the input images to a canonical, expected pose,
it may also be useful in handling a domain mismatch between the training and
testing conditions. It was applied for domain adaptation in various image
processing situations [34, 35], and here we evaluate its efficiency for
speaker and session adaptation in UTI-to-speech conversion.

Figure 2: An example UTI image, before(left) and after(right) STN.

##  4 Experimental Set-Up

###  4.1 Data Acquisition and Preprocessing

Several Hungarian male and female subjects with normal speaking abilities were
recorded while reading sentences aloud (altogether 209 sentences each), of
which two female and two male speakers were selected for the speaker
adaptation experiments (with speaker IDs 048, 049, 102 and 103). The average
duration of the recordings was about 15 minutes per speaker. For the cross-
session experiments one female speaker from the 4 speakers (speaker 048) was
asked to record 4 additional sessions on different days (obviously, with
dismounting and remounting the UTI recording device between the sessions).
These 4 additional sessions are shorter than the main session of the given
speaker, with an average duration of 3.5 minutes. All sessions were divided
into training, validation and test partitions randomly, in a 85-10-5 ratio.
The speech signal and the ultrasound articulatory data were recorded in
parallel, and were synchronized using the software provided with the UTI
equipment, Figure 3. The tongue movement was recorded in a midsagittal
orientation using the “Micro” ultrasound system of Articulate Instruments Ltd.
at 81.67 fps. The speech signal was recorded with an omnidirectional condenser
microphone. For more details on the recording set-up see [10]. In the
experiments we used the raw scanline data of the ultrasound as the training
input, after being resized to 64×128 pixels. The intensity range of the images
was min-max normalized to the [-1, 1] interval. The speech signals were
converted to a mel-spectrogram with 80 spectral components, according to the
requirements of WaveGlow. These 80-component spectral vectors served as the
training targets for our neural networks, after standardization to zero mean
and unit variance.

###  4.2 Network configuration

In the first set of experiments we applied a simple 2D convolutional (2D-CNN)
network that transforms one ultrasound image to one spectral vector, and the
localization network of the STN also was a 2D-CNN. Although this approach is
suboptimal, it may already show whether the STN is suitable to perform domain
adaptation, while visualization and interpretation is easier in 2D. Fig. 1
illustrates this network arrangement, while Fig. 2 demonstrates the effect of
STN on an actual ultrasound image. The 2D-CNN network had a quite simple and
traditional structure. It consisted of 4 convolutional layers with
30-60-90-120 filters, with a MaxPooling layer after every second layer. The
convolutional processing was followed by a fully connected layer of 300
neurons, and the linear output layer of 80 neurons. All hidden layers applied
the Swish nonlinearity, and overfitting was minimized by placing dropout
layers (p=0.2) after each processing layer. We applied the mean squared error
(MSE) loss function, which was minimized using the Adam optimizer with a batch
size of 100 and a learning rate of $2\cdot 10^{-3}$. Training was halted using
early stopping on the validation set.

The input of the STN module is the same image as that of the spectral
regression module, and it also has to perform regression. Hence we used
exactly the same 2D-CNN architecture for the STN as for the main network.
However, as the STN has to estimate only 6 parameters (the $\theta$ vector)
instead of 80, we configured it to be much smaller than the regression
network. The number of free parameters in the STN module was only about 10% of
that of the whole network.

###  4.3 Extension to 3D

The simple framework presented above can be significantly improved by
extending the input from just one image to a sequence of video frames. Such 3D
blocks of input data can be processed by various network types, such as by
combining 2D-CNN and LSTM layers, by 3D-CNN models, or by using Convolutional
LSTM layers [12, 14, 28, 27, 36]. Here, we applied the 3D-CNN architecture
introduced by Tóth et al. [27, 17]. This network had the same global
architecture as the 2D-CNN presented above. However, the convolutions were
extended to the time axis as well, in a somewhat special arrangement. The
first Conv3D layer applied a large stride along the time axis, dividing the
sequence of 25 input images into 5 five blocks along time. The next two Conv3D
layers had a filter size of 1 along time, practically processing the 5 blocks
separately. Finally, the extracted information was fused along time by the
last Conv3D layer and by the topmost two fully connected layers.

As regards the STN, there are several options to extend it to 3D. First, we
could reformulate the affine transformation to operate on 3D data. However,
for our ultrasound data, translation, shearing and rotation along the time
axis seemed to be unnecessary. Hence, we chose to keep the transformation in
2D. The second question was whether the localization network should operate on
2D or 3D data. We chose the technically simpler solution, and used the same 2D
localization network as for the 2D-CNN. In this configuration, the STN applies
the same 2D transformation to all the 25 images of the actual input block.

##  5 Results and discussion

Table 1: MSE rates of the 2D-CNN on the dev set of the 4 speakers (top), and
for the 4 extra sessions of speaker 048 (down).

| s048 | s049 | s102 | s103
---|---|---|---|---
no STN | 0.361 | 0.387 | 0.390 | 0.343
with STN | 0.358 | 0.396 | 0.389 | 0.340

| s048-2 | s048-3 | s048-4 | s048-5
---|---|---|---|---
no STN | 0.505 | 0.462 | 0.464 | 0.504
with STN | 0.495 | 0.458 | 0.467 | 0.483

First of all, we examined whether the inclusion of the STN module has any
positive effect on the results when using only training data from one speaker
and one session. Table 1 shows the MSE scores of the 2D network, with and
without the STN module, for all the 4 speakers and the 4 extra sessions from
speaker 048. In this configuration, we expected no significant benefit from
the STN, and this is exactly what we see. The scores we obtained are quite
similar accross all speakers and the extra sessions of speaker 048. One
notable difference is that the results are consistently worse for the
additional sessions of speaker 048 than for her main session. The explanation
is that these extra sessions were much shorter than the multi-speaker
recordings (3.5 vs. 15 minutes).

In the adaptation experiments we considered the network trained on the data of
speaker 048 as our base model. First, we evaluated it on the samples from the
other speakers and the additional sessions of the same speaker without any
adaptation. These results will be considered as the base results (i.e. $0\%$
relative improvement due to adaptation). Then, we trained this baseline model
further on the data from the other speakers and sessions, using three
adaptation strategies. First, we allowed the full network to learn on the new
data, and the results obtained this way will be considered to have $100\%$
improvement (i.e. the best possible result achievable by adaptation). Next, we
froze the weights of the regression network, and allowed only the STN module
to adapt. In the final configuration, besides the STN we also allowed the re-
adjustment of the weights in the output layer. This layer simply calculates a
weighted linear combination of the hidden representations, so allowing it to
adapt might significantly decrease the regression error, while it contains
only about $1\%$ of all the weights in the network. In all the adaptation
experiments we used an initial learning rate that is an order of magnitude
smaller than that of the baseline, but otherwise there were no other changes
in the training process.

Figure 3: A person has been recorded his speech and tongue movements by the
Ultrasound Tongue Imaging device. Table 2: MSE rates of the 2D-CNN model of
spk048 for the other 3 speakers (top), and for the 4 extra sessions of the
same speaker (down), using various adaptation strategies.

| Adaptation method
---|---
Speaker | None | STN | STN + out | Full | Mean $\theta$
spk049 | 1.049 | 0.588 (-71%) | 0.517 (-82%) | 0.400 | 0.887 (-25%)
spk102 | 1.401 | 0.609 (-78%) | 0.449 (-94%) | 0.389 | 1.015 (-38%)
spk103 | 1.322 | 0.552 (-79%) | 0.469 (-88%) | 0.350 | 0.909 (-42%)
Average diff. | -0% | -76% | -88% | -100% | -35%

| Adaptation method
---|---
Session | None | STN | STN + out | Full | Mean $\theta$
ses-2 | 1.131 | 0.646 (-77%) | 0.547 (-93%) | 0.503 | 0.913 (-34%)
ses-3 | 0.998 | 0.619 (-69%) | 0.485 (-94%) | 0.451 | 0.934 (-11%)
ses-4 | 1.054 | 0.641 (-70%) | 0.522 (-90%) | 0.468 | 0.908 (-24%)
ses-5 | 1.174 | 0.604 (-85%) | 0.566 (-91%) | 0.506 | 0.955 (-32%)
Average diff. | -0% | -75% | -92% | -100% | -26%

Table 2 summarizes the results achieved with the different adaptation
strategies for the additional speaker (upper panel) and for the additional
sessions of the same speaker (lower panel). As the relative improvements are
more informative than the actual MSE values, in parentheses we display the
relative error reductions (considering the no-adaptation case as 0% and full
adaptation as 100%), and in the bottom line we summarize the average
improvements. The first column (`None', i.e. no adaptation) clearly shows that
the results are unacceptable without any adaptation – the error rates are
actually worse than those for a randomly initialization net without any
training. Moreover, the scores in the `Full' adaptation column are pretty
similar to those in the baseline table, so pre-training on speaker 048 seems
to be neither beneficial nor detrimental (interestingly even for the extra
sessions from the same speaker). By allowing only the STN module to adapt
(`STN' column) we can eliminate 75-76% of the performance gap between the non-
adapted and fully adapted models. This is pretty good, considering that the
regression network itself is not modified at all, we just simply allow the STN
to learn a more optimal affine transformation for the images of the given
speaker or session. The improvement is even larger if we also allow the linear
output layer of the regression network to adjust its weights to the new
speaker or session (`STN + out'). In this case, the average error reduction is
88% for cross-speaker and 92% for cross-session adaptation. The better cross-
session score is reasonable, as one would expect that the differences caused
by the misalignment of the device might be easier to compensate by an affine
transformation than the inherent anatomical differences between speakers.

Table 3: MSE rates for cross-speaker adaptation, using the 3D-CNN network. | Adaptation method
---|---
Speaker | None | STN | STN + out | Full
spk049 | 1.105 | 0.553 (-73%) | 0.497 (-80%) | 0.348
spk102 | 1.451 | 0.502 (-84%) | 0.416 (-91%) | 0.315
spk103 | 1.541 | 0.501 (-83%) | 0.418 (-90%) | 0.294
Avg.diff. | -0% | -80% | -87% | -100%

As in the baseline evaluation (Table 1) the STN had no effect on the results,
we speculated that maybe it is unnecessary to learn a separate transformation
for each image – perhaps learning one global $\theta$ for each speaker or
session would be enough. We examined the within-speaker variance of the 6
components of theta for the baseline models, and we indeed found that it was
very low, on the order of 0.01, while the between-speaker variance of the
$\theta$ vectors was 3-5 times larger. This fact seemed to underpin our
conjecture, so we preformed the following simple experiment. Instead of using
a unique theta for each image, we simply replaced the STN of the adapted
models with the mean $\theta$ vector over the samples of the given
speaker/session. The results are shown in the rightmost columns of Table 2
('mean $\theta$'), but the findings are negative: while there is a 25-35%
reduction in the error rates, it is very far from the best achievable. It
seems that even though the variance of $\theta$ is small, the minor nuances in
the learned transformations per image play an important role in the regression
accuracy.

Finally, we extended our experiments to 3D input blocks consisting of 25
subsequent images, using a 3D-CNN for the regression task and a 2D-CNN for the
STN module (cf. Section 4.3). In this case, the baseline model yielded an MSE
of 0.275 and 0.278 with and without the STN on the data set of speaker 048. As
speaker and session adaptation gave very similar results previously, we
repeated only the cross-speaker adaptation experiments with this network, and
the results are shown in Table 3. While all the error rates are typically
lower than they were for the 2D-CNN, the relative improvements with respect to
the various adaptation strategies are very similar.

##  6 Conclusions

Current tongue ultrasound-based SSI systems are sensitive to changing
speakers, or even to a slight displacement of the recording device. In this
study we examined whether an STN module is able to counterbalance these
factors by applying an affine transform on the input image. When applying a
simple 2D-CNN for spectral estimation, we found that allowing only the
adaptation of the STN module can reduce the error rate by about 75%, while
allowing also the linear output layer to adapt can compensate for 88-92% of
the error. We also extended the experiments to 3D input blocks, and we
observed similar tendencies, although the improvement was somewhat smaller.
Considering, however, that the STN module and the output layer altogether
contain only 10% of the weights of the full network, and that the improvements
were consistent for cross-session and cross-speaker set-ups, the proposed
method might allow for a quicker adaptation of the UTI-processing workflow. In
the future we plan to run experiments with a 3D localization network and with
smaller amounts of adaptation material.

##  7 Acknowledgements

The authors would like to express their gratitude to the following
institutions and funding sources for their valuable support in conducting this
research. This work was supported by the European Union project
RRF-2.3.1-21-2022-00004 within the framework of the Artificial Intelligence
National Laboratory and project TKP2021-NVA-09, implemented with the support
provided by the Ministry of Innovation and Technology of Hungary from the
National Research, Development, and Innovation Fund, financed under the
TKP2021-NVA funding scheme and FK 142163 grant. T.G.Cs. was supported by the
Bolyai János Research Fellowship of the Hungarian Academy of Sciences and by
the ÚNKP-22-5-BME-316 grant. The RTX A5000 GPU card used was donated by the
NVIDIA Corporation.

## References

  * [1] B. Denby, T. Schultz, K. Honda, T. Hueber, J. M. Gilbert, and J. S. Brumberg, ``Silent speech interfaces,'' _Speech Communication_ , vol. 52, no. 4, pp. 270–287, 2010.
  * [2] J. A. Gonzalez-Lopez, A. Gomez-Alanis, J. M. Martin Donas, J. L. Perez-Cordoba, and A. M. Gomez, ``Silent Speech Interfaces for Speech Restoration: A Review,'' _IEEE Access_ , vol. 8, pp. 177 995–178 021, 2020\.
  * [3] T. Hueber, E.-L. Benaroya, G. Chollet, G. Dreyfus, and M. Stone, ``Development of a silent speech interface driven by ultrasound and optical images of the tongue and lips,'' _Speech Communication_ , vol. 52, no. 4, pp. 288–300, 2010\.
  * [4] A. Ephrat and S. Peleg, ``Vid2speech: Speech reconstruction from silent video,'' in _Proc. ICASSP_ , 2017, pp. 5095–5099.
  * [5] R. Mira, K. Vougioukas, P. Ma, S. Petridis, B. W. Schuller, and M. Pantic, ``End-to-End Video-To-Speech Synthesis using Generative Adversarial Networks,'' _arXiv preprint arXiv:2104.13332_ , 2021.
  * [6] D. Oneaţă, A. Stan, and H. Cucu, ``Speaker disentanglement in video-to-speech conversion,'' in _Proc. EUSIPCO_ , 2021, pp. 46–50.
  * [7] N. Saleem, J. Gao, M. Irfan, E. Verdu, and J. P. Fuente, ``E2E-V2SResNet: Deep residual convolutional neural networks for end-to-end video driven speech synthesis,'' _Image and Vision Computing_ , vol. 119, p. 104389, 2022.
  * [8] D. Wang, S. Yang, D. Su, X. Liu, D. Yu, and H. Meng, ``VCVTS: Multi-speaker Video-to-Speech synthesis via cross-modal knowledge transfer from voice conversion,'' in _Proc. ICASSP_ , 2022.
  * [9] B. Denby and M. Stone, ``Speech synthesis from real time ultrasound images of the tongue,'' in _Proc. ICASSP_ , 2004, pp. 685–688.
  * [10] T. G. Csapó, T. Grósz, G. Gosztolya, L. Tóth, and A. Markó, ``DNN-based ultrasound-to-speech conversion for a silent speech interface,'' in _Proc. Interspeech_ , 2017, pp. 3672–3676.
  * [11] E. Tatulli and T. Hueber, ``Feature extraction using multimodal convolutional neural networks for visual speech recognition,'' in _Proc. ICASSP_ , 2017, pp. 2971–2975.
  * [12] E. Moliner and T. G. Csapó, ``Ultrasound-based silent speech interface using convolutional and recurrent neural networks,'' _Acta Acustica united with Acustica_ , vol. 105, no. 4, pp. 587–590, 2019.
  * [13] T. Csapó, M. Al-Radhi, G. Németh, G. Gosztolya, T. Grósz, L. Tóth, and A. Markó, ``Ultrasound-based silent speech interface built on a continuous vocoder,'' 2019.
  * [14] N. Kimura, M. Kono, and J. Rekimoto, ``Sottovoce: An ultrasound imaging-based silent speech interaction using deep neural networks,'' in _Proc. of CHI Conf. on Human Factors in Computing Systems_ , 2019.
  * [15] T. G. Csapó, C. Zainkó, L. Tóth, G. Gosztolya, and A. Markó, ``Ultrasound-based Articulatory-to-Acoustic Mapping with WaveGlow Speech Synthesis,'' in _Proc. Interspeech_ , 2020, pp. 2727–2731.
  * [16] C. Zainkó, L. Tóth, A. H. Shandiz, G. Gosztolya, A. Markó, G. Németh, and T. G. Csapó, ``Adaptation of Tacotron2-based Text-To-Speech for Articulatory-to-Acoustic Mapping using Ultrasound Tongue Imaging,'' in _Proc. ISCA SSW11_ , 2021, pp. 54–59.
  * [17] A. H. Shandiz, L. Tóth, G. Gosztolya, A. Markó, and T. G. Csapó, ``Neural Speaker Embeddings for Ultrasound-Based Silent Speech Interfaces,'' in _Proc. Interspeech 2021_ , 2021, pp. 1932–1936.
  * [18] T. G. Csapó and K. Xu, ``Quantification of Transducer Misalignment in Ultrasound Tongue Imaging,'' in _Proc. Interspeech_ , 2020, pp. 3735–3739.
  * [19] G. Gosztolya, T. Grósz, L. Tóth, A. Markó, and T. G. Csapó, ``Applying DNN Adaptation to Reduce the Session Dependency of Ultrasound Tongue Imaging-Based Silent Speech Interfaces,'' _Acta Polytechnica Hungarica_ , vol. 17, no. 7, pp. 109–124, 2020.
  * [20] M. S. Ribeiro, J. Sanger, J.-X. X. Zhang, A. Eshky, A. Wrench, K. Richmond, and S. Renals, ``TaL: a synchronised multi-speaker corpus of ultrasound tongue imaging, audio, and lip videos,'' in _Proc. SLT_ , 2021, pp. 1109–1116.
  * [21] J.-X. Zhang, K. Richmond, Zhen-Hua-Ling, and L.-R. Dai, ``TaLNet: Voice Reconstruction from Tongue and Lip Articulation with Transfer Learning from Text-to-Speech Synthesis,'' in _Proc. AAAI_ , 2021.
  * [22] Y. Yu, A. H. Shandiz, and L. Tóth, ``Reconstructing speech from real-time articulatory mri using neural vocoders,'' in _Proc. EUSIPCO_ , 2021, pp. 945–949.
  * [23] R. Prenger, R. Valle, and B. Catanzaro, ``Waveglow: A flow-based generative network for speech synthesis,'' in _Proc. ICASSP_ , 2019, pp. 3617–3621.
  * [24] S. Ji, W. Xu, M. Yang, and K. Yu, ``3d convolutional neural networks for human action recognition,'' _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 35, no. 1, pp. 221–231, 2013.
  * [25] D. Tran, H. Wang, L. Torresani, J. Ray, Y. LeCun, and M. Paluri, ``A closer look at spatiotemporal convolutions for action recognition,'' in _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2018, pp. 6450–6459.
  * [26] C. Feichtenhofer, H. Fan, J. Malik, and K. He, ``Slowfast networks for video recognition,'' in _Proceedings of the IEEE/CVF international conference on computer vision_ , 2019, pp. 6202–6211.
  * [27] L. Tóth and A. H. Shandiz, ``3D convolutional neural networks for ultrasound-based silent speech interfaces,'' in _Proc. ICAISC_. Springer, 2020, pp. 159–169.
  * [28] P. Saha, Y. Liu, B. Gick, and S. Fels, ``Ultra2speech – a deep learning framework for formant frequency estimation and tracking from ultrasound tongue images,'' in _Proc. MICCAI_ , 2020, pp. 473–482.
  * [29] G. Bertasius, H. Wang, and L. Torresani, ``Is space-time attention all you need for video understanding,'' _Proc. ICML_ , vol. 139, pp. 813–824, 2021.
  * [30] D. Serdyuk, O. Braga, and O. Siohan, ``Transformer-based video front-ends for audio-visual speech recognition,'' _arXiv preprint arXiv:2201.10439_ , 2022\.
  * [31] M. Jaderberg, K. Simonyan, A. Zisserman, and K. Kavukcuoglu, ``Spatial transformer networks,'' in _NIPS 28_ , 2015, pp. 2017–2025.
  * [32] A. Bas, P. Huber, W. A. P. Smith, M. Awais, and J. Kittler, ``3d morphable models as spatial transformer networks,'' in _Proc. Int. Conf. on Computer Vision (ICCV) Workshops_ , 2017, pp. 895–903.
  * [33] J. Yu and S. Wang, ``Visual speech recognition in natural scenes based on spatial transformer networks,'' in _Proc. Int. Conf. on Anti-counterfeiting, Security, and Identification (ASID)_ , 2020.
  * [34] W. Chu, W.-C. Hung, Y.-H. Tsai, D. Cai, and M.-H. Yang, ``Weakly-supervised caricature face parsing through domain adaptation,'' in _Proc. Int. Conf. on Image Processing (ICIP)_ , 2019, pp. 3282–3286.
  * [35] R. Robinson, Q. Dou, D. Coelho de Castro, K. Kamnitsas, M. de Groot, R. M. Summers, D. Rueckert, and B. Glocker, ``Image-level harmonization of multi-site data using image-and-spatial transformer networks,'' in _Proc. MICCAI_ , 2020, pp. 710–719.
  * [36] C. Zhao, J. Zhang, C. Wu, H. Wang, and K. Xu, ``Predicting tongue motion in unlabeled ultrasound video using convolutional LSTM neural networks,'' in _Proc. ICASSP_ , 2019, pp. 5926–5930.