# Joint magnitude estimation and phase recovery using Cycle-in-Cycle GAN for
non-parallel speech enhancement

###### Abstract

For the lack of adequate paired noisy-clean speech corpus in many real
scenarios, non-parallel training is a promising task for DNN-based speech
enhancement methods. However, because of the severe mismatch between input and
target speeches, many previous studies only focus on the magnitude spectrum
estimation and remain the phase unaltered, resulting in the degraded speech
quality under low signal-to-noise ratio conditions. To tackle this problem, we
decouple the difficult target w.r.t. original spectrum optimization into
spectral magnitude and phase, and a novel Cycle-in-Cycle generative
adversarial network (dubbed CinCGAN) is proposed to jointly estimate the
spectral magnitude and phase information stage by stage under unpaired data.
In the first stage, we pretrain a magnitude CycleGAN to coarsely estimate the
spectral magnitude of clean speech. In the second stage, we incorporate the
pretrained CycleGAN with a complex-valued CycleGAN as a cycle-in-cycle
structure to simultaneously recover phase information and refine the overall
spectrum. Experimental results demonstrate that the proposed approach
significantly outperforms previous baselines under non-parallel training. The
evaluation on training the models with standard paired data also shows that
CinCGAN achieves remarkable performance especially in reducing background
noise and speech distortion.

Index Terms—  speech enhancement, non-parallel, Cycle-in-Cycle GAN, magnitude
spectrum estimation, phase recovery

##  1 Introduction

Speech is often corrupted by various types of background noise in real-world
environments. These interferences severely degrade the performance of many
speech-related tasks such as automatic speech recognition (ASR) systems or
hearing aids [1]. To cope with these degradations, speech enhancement (SE) can
be applied to remove or attenuate the background noise for better speech
quality and intelligibility. Recent SE approaches based on deep neural
networks (DNNs) have demonstrated their superior capability in dealing with
non-stationary noise under low signal-to-noise ratio (SNR) conditions [2].
These DNN-based methods can be roughly categorized into masking-based
approaches [3] and mapping-based approaches [4]. More recently, generative
adversarial networks (GANs) [5] begin to gain development in the single-
channel SE area [6, 7, 8, 9, 10, 11, 12], in which a generator ($G$) aims to
conduct the enhancement process and a discriminator ($D$) aims to distinguish
between the generated speech features and real clean ones. For these standard
supervised methods, the networks are trained to minimize the discrepancy
between the enhanced output and the clean target.

Despite the impressive performance of the above DNN-based SE approaches, they
always need numerous paired clean-noisy samples to conduct supervised training
and improve the generalization. However, in some real scenarios, it is
troublesome to record parallel clean-noisy pairs, and we can only obtain clean
speech that mismatches the source noisy speech. In this case, standard
supervised SE methods usually cannot produce satisfactory performance. To
tackle this problem, cycle-consistent GAN (CycleGAN) has been developed to
conduct unsupervised SE, which was originally used for unpaired image-to-image
translation [13]. In the SE area, CycleGAN-based approaches demonstrate the
remarkable ability in preserving the speech structure and improving speech
quality in both paired and unpaired cases [14, 15, 16, 17, 18]. Nevertheless,
these conventional CycleGAN-based approaches have several limitations for non-
parallel SE. First, to ensure cycle consistency of the original noisy speech
and target clean speech for unpaired data, the enhanced signal always contains
the original noise information. In other words, cycle-consistency-based
methods remain audible residual noise, which is challenging to be eliminated
completely. Second, previous CycleGAN-based SE methods only focus on magnitude
spectrum estimation and remain the noisy phase unaltered. This is because
simultaneously restoring the clean spectral magnitude and phase information
still remains to be a tough task, let alone when noisy-clean pairs are
mismatched. However, recent studies revealed the importance of the accurate
phase for improving perceptual speech quality, especially in low SNR cases
[19].

Fig. 1:  System flowchart of CinCGAN. The magnitude and complex cycle are
shown in the yellow and blue dotted boxes, respectively. (a) The diagram of
the forward cycle in CinCGAN. (b) The diagram of the backward cycle in
CinCGAN.

Multi-stage algorithms based on curriculum learning begin to thrive in the SE
area [20, 21], where the original difficult task is decomposed into multiple
simpler sub-tasks and a better result can be induced progressively. Motivated
by these studies, we couple a magnitude spectrum estimation CycleGAN (MCGAN)
and a complex spectrum refined CycleGAN (CCGAN) as a Cycle-in-Cycle GAN
(CinCGAN) paradigm to estimate the clean spectral magnitude and phase
information step-by-step under non-parallel training. In the first step, the
complex original spectrum is decoupled into the spectral magnitude and phase,
and only the magnitude of the clean spectrum is estimated by MCGAN.
Subsequently, a complex-valued CCGAN is integrated with the pretrained MCGAN
to jointly recover the clean spectral magnitude and phase information. To be
specific, we couple the estimated magnitude spectrum by MCGAN and original
phase as a coarse-estimated complex spectrum, and feed it into CCGAN to
estimate both real and imaginary (RI) components of the clean spectrum.

The remainder of the paper is organized as follows. In Section 2, the overall
diagram is described in detail. The experimental setup is presented in Section
3, while the results and analysis are given in Section 4. Some conclusions are
drawn in Section 5.

##  2 Algorithm description

###  2.1  Cycle-in-Cycle Scheme

To simultaneously estimate magnitude and phase information of clean spectrum
under non-parallel training conditions, we integrate a magnitude spectrum
estimation CycleGAN (MCGAN) with a complex spectrum restoration CycleGAN
(CCGAN) as a Cycle-in-Cycle GAN (CinCGAN). The proposed CinCGAN consists of a
forward noisy-clean-noisy cycle and a backward clean-noisy-clean cycle, as
shown in Fig. 1. In the forward cycle, the enhancement procedure can be
divided into two steps. First, we decouple the complex spectrum into spectral
magnitude and phase, and only the amplitude is processed. A magnitude
generator $G^{mc}_{X\rightarrow Y}$ is utilized to convert the compressed
noisy spectral magnitude $\left|{X}_{t,f}\right|$ to the ”clean” one
$|\widetilde{X}^{mc}|$, while a magnitude discriminator $D^{mc}_{Y}$ is used
to distinguish the generated spectral magnitude (_i.e._ ,
$|\widetilde{X}_{t,f}|$) from the clean one (_i.e._ ,
$\left|{Y}_{t,f}\right|$), where $t$ and $f$ denote the time and frequency
indices. To constrain the cycle consistency under unpaired data, an inverse
generator $F^{mc}_{Y\rightarrow X}$ is employed to reconstruct the noisy
spectral magnitude from the generated one. Subsequently, we couple the
estimated spectral magnitude and the original phase (_i.e._ ,
$e^{j\theta_{X_{t,f}}}$) into a coarsely estimated complex spectrum and feed
it into CCGAN. By estimating both RI parts of the clean spectrum, CCGAN can
refine the spectral magnitude and recover the clean phase. Specifically, a
complex-valued generator $G^{cc}_{X\rightarrow Y}$ attempts to generate the
enhanced complex spectrum $\widetilde{X}^{cc}$ as clean one, and a complex
discriminator $D^{cc}_{Y}$ attempts to discriminate between the generated and
clean complex spectrum (_i.e._ , ${Y}_{t,f}$). Similarly, an inverse complex
generator $F^{cc}_{Y\rightarrow X}$ brings the enhanced complex spectrum
$\widetilde{X}^{cc}$ back to the original noisy domain. To better constrain
the mapping optimization space, a backward clean-noisy-clean cycle is also
designed and jointly trained with the forward cycle. In the testing stage, the
whole enhancement procedure can be formulated as:

$\displaystyle|\widetilde{X}^{mc}|=G^{mc}_{X\rightarrow
Y}\left(|X_{t,f}|;\Phi_{1}\right),$ (1)
$\displaystyle\widetilde{X}^{mc}=|\widetilde{X}^{mc}|\exp(j\theta_{{X}_{t,f}}),$
(2) $\displaystyle\widetilde{X}^{cc}=G^{cc}_{X\rightarrow
Y}\left(\widetilde{X}^{mc};\Phi_{2}\right),\vspace{-3mm}$ (3)

where $|X_{t,f}|$ and $\theta_{{X}_{t,f}}$ denote the noisy spectral magnitude
and phase, respectively. $|\widetilde{X}^{mc}|$ and $\widetilde{X}^{cc}$
denote the outputs of the noisy-to-clean magnitude generator
$G^{mc}_{X\rightarrow Y}$ and complex generator $G^{cc}_{X\rightarrow Y}$ with
parameter set $\Phi_{\left(\cdot\right)}$.

###  2.2 Network Architecture

Following the same configuration as [17], the magnitude generator in MCGAN
consists of three major parts, including three downsampling layers as the
encoder, an adaptive attention-in-attention (AIA) module and three upsampling
layers as the decoder, as shown in Fig. 2 (a). Each downsampling/upsampling
layer block is composed of a 2D (de)convolution layer, followed by instance
normalization (IN), Parametric Relu activation function (PReLU) and gated
liner units (GLUs). The proposed AIA module consists of six adaptive time-
frequency attention (ATFA) modules and an adaptive hierarchical attention
(AHA) module, as proposed in [17]. More specifically, ATFA consists of two
branches: an adaptive time attention branch (ATAB) and an adaptive frequency
attention branch (AFAB). The two branches cooperate to capture the global
dependencies along temporal and spectral axes in parallel. The AHA module
integrates different hierarchical feature maps given all outputs of ATFA
modules and fuses the global context of all intermediate outputs with
different weights, thus helping to guide the feature learning progressively.
The output channels of each layer in the encoder are (16, 32, 64), while the
kernel size and strides are set to (3, 5) and (1, 2) along the time and
frequency axes, respectively.

Fig. 2:  (a) The diagram of the magnitude generators in MCGAN. (b) The diagram
of the complex generators in CCGAN.

We also design the pure complex-valued generators of CCGAN similar as [18],
which consists of eight complex-valued encoder/decoder blocks, and an complex-
valued AIA module, as shown in Fig. 2(b). The complex encoder/decoder blocks
include complex-valued 2D (de)convolutions, complex IN and real-valued PReLU.
The complex-valued convolution simulates the correlation between real and
imaginary parts of feature maps by complex multiplication, which is composed
of four traditional convolution operations [22]. To further capture long-term
dependency with complex-valued features, we utilize the above attention-in-
attention mechanism in the complex-valued manner. The number of the channels
for each complex encoder layer is (32, 32, 64, 64, 128, 128, 256, 256), with
kernel size and strides set to (3, 3) and (1, 2), respectively.

We utilize multi-scale discriminators that can force the generator to produce
feature maps with global consistency and finer details, as pointed out in
[17]. Both the magnitude and complex discriminators are composed of six real-
valued 2D convolutions, spectral normalization (SN) and PReLU. The kernel size
for each 2D convolutions layer is (3, 5) except (1, 1) for the last layer in
the temporal and frequency axis, respectively. The number of channels
throughout the 2D convolutions is (32, 32, 64, 64, 128, 1), with stride set to
(1, 2). Note that we concatenate RI components along the channel axis and feed
them into the complex discriminator.

###  2.3 Loss Function

In both MCGAN and CCGAN, we use the following three losses to jointly optimize
the forward and backward cycle, namely relativistic adversarial loss, cycle-
consistency loss, and an identity mapping loss. Taking the magnitude CycleGAN
(MCGAN) as an example, for the forward noisy-to-clean mapping, the
relativistic average least-square (RaLS) adversarial loss [23] is used to
enforce the enhanced magnitude spectra indistinguishable from the clean ones,
which is calculated as follows.

$\displaystyle\mathcal{L}_{Radv}(D^{mc}_{Y})=\mathbb{E}_{y}\left[(D^{mc}_{Y}(y)-\mathbb{E}_{x}D^{mc}_{Y}(G_{X\rightarrow
Y}(x))-1)^{2}\right]+$
$\displaystyle\mathbb{E}_{x}\left[(D^{mc}_{Y}(G_{X\rightarrow
Y}(x))-\mathbb{E}_{y}(D^{mc}_{Y}(y))+1)^{2}\right],$
$\displaystyle\mathcal{L}_{Radv}(G^{mc}_{X\rightarrow
Y})=\mathbb{E}_{x}\left[(D^{mc}_{Y}(G_{X\rightarrow
Y}(x)))-\mathbb{E}_{y}D^{mc}_{Y}(y)-1)^{2}\right]$
$\displaystyle+\mathbb{E}_{y}\left[(D^{mc}_{Y}(y)-\mathbb{E}_{x}(D^{mc}_{Y}(G_{X\rightarrow
Y}(x)))+1)^{2}\right],\vspace{-2mm}$ (4)

Analogously, $\mathcal{L}_{Radv}(D^{mc}_{X})$ and
$\mathcal{L}_{Radv}(F^{mc}_{Y\rightarrow X})$ are employed in the inverse
magnitude clean-to-noisy mapping. The cycle-consistency loss and identity
mapping loss are utilized to constrain the mapping functions under non-
parallel training conditions and to further preserve the linguistic
composition of the original speech, which can be defined as:

$\displaystyle\mathcal{L}_{cycle}(G^{mc}_{X\rightarrow Y},F^{mc}_{Y\rightarrow
X})=\mathbb{E}_{x}\left[\left\|F_{Y\rightarrow X}(G^{mc}_{X\rightarrow
Y}(x))-x\right\|_{1}\right]$
$\displaystyle+\mathbb{E}_{y}\left[\left\|G^{mc}_{X\rightarrow
Y}(F^{mc}_{Y\rightarrow X}(y))-y\right\|_{1}\right],$ (5)
$\displaystyle\mathcal{L}_{id}(G^{mc}_{X\rightarrow Y},F^{mc}_{Y\rightarrow
X})=\mathbb{E}_{x}\left[\left\|F^{mc}_{Y\rightarrow
X}(x)-x\right\|_{1}\right]+$
$\displaystyle\mathbb{E}_{y}\left[\left\|G^{mc}_{X\rightarrow
Y}(y)-y\right\|_{1}\right],\vspace{-2mm}$ (6)

In the first step, we pretrain MCGAN alone with $\mathcal{L}_{MCGAN}$ until
convergence, which is calculated as:

$\begin{gathered}\mathcal{L}_{MCGAN}=\mathcal{L}_{Radv}(G^{mc}_{X\rightarrow
Y},D^{mc}_{Y})+\mathcal{L}_{Radv}(F^{mc}_{Y\rightarrow X},D^{mc}_{X})\\\
+\lambda_{cycle}\mathcal{L}_{cycle}(G^{mc}_{X\rightarrow
Y},F^{mc}_{Y\rightarrow X})+\lambda_{id}\mathcal{L}_{id}(G^{mc}_{X\rightarrow
Y},F^{mc}_{Y\rightarrow X}),\\\ \end{gathered}$ (7)

where $\lambda_{cycle}$ and $\lambda_{id}$ are empirically set to be 5 and 10,
respectively. Then, we jointly fine-tune MCGAN and CCGAN with the same
relativistic adversarial loss, cycle-consistency loss, and identity mapping
loss, where the parameters of the first sub-network are initialized with the
pretrained MCGAN. The overall loss functions of CinCGAN are expressed as:

$\begin{gathered}\mathcal{L}_{CinCGAN}=\gamma\mathcal{L}_{MCGAN}+\mathcal{L}_{CCGAN}\\\
\end{gathered}$ (8)

where $\mathcal{L}_{CCGAN}$ denotes the loss function for the complex cycle.
With some trials, we empirically set $\gamma=0.1$ in all the following
experiments .

##  3 Experiments

###  3.1 Datasets

The dataset used in this work is publicly available as proposed in [24], which
is a selection of the Voice Bank corpus [25] with 28 speakers for training and
another 2 unseen speakers for testing. The training set consists of 11,572
mono audio samples, while the test set contains 824 utterances by 2 speakers
(one male and one female). For the training set, audio samples are mixed
together with one of the 10 noise types (_i.e._ , two artificial (babble and
speech shaped) and eight real noise from the DEMAND database [26]) at four
SNRs of $\left\\{0\rm{dB},5\rm{dB},10\rm{dB},15\rm{dB}\right\\}$. The testing
utterances are created with 5 unseen test-noise types (all from the DEMAND
database) at SNRs of
$\left\\{2.5\rm{dB},7.5\rm{dB},12.5\rm{dB},17.5\rm{dB}\right\\}$. The original
raw waveforms are downsampled from 48kHz to 16kHz for our task.

###  3.2 Implementation Setup

The Hanning window of length 32ms is applied, with 75% overlap between
adjacent frames. The 512-point STFT is utilized, leading to a 257-D spectral
feature. Due to the efficacy of power compression in dereverberation and
denoising tasks [21], we conduct the power compression toward the spectral
magnitude while leaving the phase unaltered, and the optimal compression
coefficient is set to 0.5. For the non-parallel training strategy, we randomly
crop a fixed-length segment (_i.e._ , 108 frames) from a randomly selected
noisy audio file as the input, while the target is a randomly selected clean
audio file whose content is totally different from the input audio. We adopt
the Adam optimizer [27] with the momentum term $\beta_{1}=0.9$,
$\beta_{2}=0.999$. We divide our training process into two steps. Firstly, we
pretrain the MCGAN with an initial learning rate of 0.0002 for discriminators
and 0.0005 for generators, respectively. Subsequently, we jointly fine-tune
the pretrained MCGAN and CCGAN as the proposed CinCGAN with a learning rate of
0.0002 for all generators and discriminators, while the learning rate is
linearly decayed after 40 epochs. We use $\mathcal{L}_{id}$ in each cycle only
for the first 20 epochs. The sound demo samples are available online
111https://github.com/yuguochencuc/CinCGAN-SE.

##  4 Results and Analysis

Table 1: Experimental results among different models including GAN-based
methods and the proposed CycleGAN-based methods under unpaired data. Best
scores are highlighted in bold.

Methods | Feature type | Magnitude | Complex | PESQ | STOI(%) | CSIG | CBAK | COVL | SegSNR | DNSMOS
---|---|---|---|---|---|---|---|---|---|---
fc | bc | fc | bc
Unprocessed |  – |  – |  – |  – |  – | 1.97 | 92.1 | 3.35 | 2.44 | 2.63 | 1.68 | 3.02
GAN-based methods |
MGAN | Magnitude | $\times$ | $\times$ | $\times$ | $\times$ | 2.03 | 91.6 | 3.54 | 2.78 | 2.72 | 5.28 | 2.72
MGAN+fc | Magnitude | $\checkmark$ | $\times$ | $\times$ | $\times$ | 2.58 | 92.8 | 3.81 | 3.03 | 3.19 | 5.28 | 3.26
CGAN | RI components | $\times$ | $\times$ | $\times$ | $\times$ | 1.86 | 88.9 | 3.17 | 2.62 | 2.64 | 2.98 | 2.63
CGAN+fc | RI components | $\times$ | $\times$ | $\checkmark$ | $\times$ | 2.32 | 91.2 | 3.48 | 2.74 | 3.18 | 4.67 | 3.04
Proposed CycleGAN-based Systems |
MCGAN | Magnitude | $\checkmark$ | $\checkmark$ | $\times$ | $\times$ | 2.67 | 93.2 | 3.86 | 3.20 | 3.21 | 7.23 | 3.47
CCGAN | RI components | $\times$ | $\times$ | $\checkmark$ | $\checkmark$ | 2.56 | 92.1 | 3.67 | 3.10 | 3.16 | 5.38 | 3.42
CinCGAN (I) | Magnitude + RI | $\checkmark$ | $\times$ | $\checkmark$ | $\times$ | 2.70 | 93.4 | 3.93 | 3.24 | 3.25 | 7.34 | 3.44
CinCGAN (II) | Magnitude + RI | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$ | 2.77 | 93.6 | 3.96 | 3.02 | 3.30 | 4.49 | 3.49
CinCGAN (III) | Magnitude + RI | $\checkmark$ | $\times$ | $\checkmark$ | $\checkmark$ | 2.73 | 93.5 | 3.94 | 3.27 | 3.29 | 7.98 | 3.51
CinCGAN (IV) | Magnitude + RI | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | 2.84 | 94.1 | 4.10 | 3.36 | 3.37 | 8.91 | 3.53

Table 2: Comparison with other GAN and Non-GAN based systems under standard
paired data. ”$-$” denotes that the result is not provided in the original
paper.

Methods | Feature type | PESQ | STOI(%) | CSIG | CBAK | COVL
---|---|---|---|---|---|---
Unprocessed |  – | 1.97 | 92.1 | 3.35 | 2.44 | 2.63
GAN-based Systems |
SEGAN [6] | Waveform | 2.16 | 92.5 | 3.48 | 2.94 | 2.80
MMSEGAN [9] | Gammatone | 2.53 | 93.0 | 3.80 | 3.12 | 3.14
SERGAN [7] | Waveform | 2.51 | 93.7 | 3.78 | 3.23 | 3.16
CP-GAN [10] | Waveform | 2.64 | 94.0 | 3.93 | 3.29 | 3.28
MetricGAN [8] | Magnitude | 2.86 |  – | 3.99 | 3.18 | 3.42
CRGAN [11] | Magnitude | 2.92 | 94.0 | 4.16 | 3.24 | 3.54
SASEGAN [12] | Waveform | 2.36 | 93.5 | 3.54 | 3.08 | 2.93
Non-GAN based Systems |
Wave-U-net [28] | Waveform | 2.64 |  – | 3.56 | 3.08 | 3.09
DFL-SE [29] | Waveform |  – |  – | 3.86 | 3.33 | 3.22
CRN-MSE [30] | Magnitude | 2.61 | 93.8 | 3.78 | 3.11 | 3.24
GCRN [31] | RI components | 2.51 | 94.0 | 3.71 | 3.24 | 3.09
DCCRN [32] | RI components | 2.68 | 93.9 | 3.88 | 3.18 | 3.27
TFSNN [33] | Waveform | 2.79 |  – | 4.17 | 3.27 | 3.49
Proposed CycleGAN-based approaches |
MCGAN | Magnitude | 2.74 | 93.6 | 3.96 | 3.25 | 3.29
CCGAN | RI components | 2.60 | 92.8 | 3.82 | 3.12 | 3.20
CinCGAN | Magnitude+ RI | 2.86 | 94.4 | 4.18 | 3.38 | 3.42

We use the following objective metrics to evaluate speech enhancement
performance: the perceptual evaluation of speech quality (PESQ) [34], short-
time objective intelligibility (STOI) [35], segmental signal-to-noise ratio
(SegSNR), and three mean opinion score (MOS) prediction (_i.e._ , signal
distortion evaluation (CISG), the intrusiveness of background noise (CBAK) and
overall effect (COVL)) [36]. We also evaluate the subjective quality by DNSMOS
[37], which is a robust non-intrusive perceptual speech quality metric. Higher
values of all metrics indicate better performance.

###  4.1 Comparison with GAN-based and CycleGAN-based methods under non-
parallel Training

To investigate the effectiveness of our proposed method under a non-parallel
training strategy, we re-implement two GAN based SE approaches with spectral
magnitude and complex spectrum as input (_i.e._ , MGAN and CGAN), and their
corresponding models with a forward cycle (fc) (_i.e._ , MGAN+fc and CGAN+fc).
We then investigate the performance of different cycle structures, including
MCGAN, CCGAN fed with the original noisy complex spectrum, and CinCGANs with
different cycles (_i.e._ , CinCGAN (I)-(IV)). To be specific, CinCGAN (I)
employs the magnitude and complex fc, while CinCGAN (II) and CinCGAN (III) add
the backward cycle in MCGAN and CCGAN, respectively. Furthermore, CinCGAN (IV)
employs all the forward and backward cycle in the magnitude and complex
domain. From Table 1, one can have the following observations. First, MGAN+fc
outperforms other GAN-based methods in terms of all metrics. This fully
validates the efficacy of fc under a non-parallel training strategy. Second,
when adding a backward cycle (bc) as standard CycleGANs, MCGAN and CCGAN
consistently surpass MGAN+fc and CGAN+fc, respectively. Even only estimating
the spectral magnitude and remaining the noisy phase unaltered, MCGAN achieves
better performance than CCGAN, which aims to estimate both RI components of
the clean spectrum. The possible reason is that it seems quite difficult to
simultaneously restore the magnitude spectrum and phase information under
unpaired data by a one-stage strategy, due to the severe mismatch between
input and target. Third, by integrating CCGAN with pretrained MCGAN, CinCGANs
with different cycles improve speech performance at different degrees. When
only using spectral magnitude fc and complex fc as shown in Figure.1 (a),
CinCGAN (I) outperforms MCGAN and CCGAN in terms of most metrics. This
phenomenon reveals the significance of simultaneously recovering spectral
magnitude and phase information under non-parallel training conditions. That
is to say, when the phase is incorporated in the complex cycle, better speech
quality can be derived. When adding the spectral magnitude bc, CinCGAN (II)
further achieves average 0.07 PESQ, 0.2% STOI, and 0.05 DNSMOS improvements
than CinCGAN (I), which demonstrates the merit of bc in improving perceptual
speech quality and intelligibility. However, the spectral magnitude bc forces
the noisy-to-clean generators in MCGAN to preserve original speech
information, thus leading to residual noise. We also notice that using both
spectral magnitude bc and complex bc can facilitate all generators, achieving
notable improvements. For example, going from CinCGAN (I) to CinCGAN (IV),
average 0.14, 0.7%, 1.57dB and 0.11 improvements are observed in terms of
PESQ, STOI, SegSNR and DNSMOS, respectively. This indicates that CinCGAN with
both cycles in MCGAN and CCGAN consistently improves speech quality, while
effectively reducing residual noise and speech distortion.

###  4.2 Comparison with other competitive methods under standard parallel
training

As presented in Table 2, our proposed methods are compared with other state-
of-the-art (SOTA) SE approaches under standard paired data, including GAN-
based methods (_e.g._ , SEGAN [6], MMSEGAN [9], SERGAN [7], MetricGAN [8],
SASEGAN [12] and CRGAN [11]) and Non-GAN based methods (_i.e._ , Wave-U-net
[28], DFL-SE [29], CRN-MSE [30], GCRN [31], DCCRN [32], and TFSNN [33]). Note
that we re-implement CRN-MSE, GCRN and DCCRN on Voice Bank + DEMAND dataset,
and directly use the reported scores of other methods in their original
papers. One can have the following observations. Compared with the existing
competitive GAN-based methods, CinCGAN provides better performance on STOI,
CSIG and CBAK scores, which verifies the merit of CinCGAN in reducing speech
distortion and suppressing residual noise components. Although CRGAN achieves
a better PESQ score than CinCGAN, which is because that PESQ is directly used
as the metric loss for CRGAN, CinCGAN provides average 0.4%, 0.02 and 0.14
improvements than CRGAN with respect to STOI, CSIG and CBAK, respectively.
When compared with Non-GAN based methods, CinCGAN also provides competitive
performance in terms of all metrics. For example, CinCGAN achieves average
0.18, 0.7% ,0.30 0.20 and 0.15 improvements than DCCRN in PESQ, STOI, CSIG,
CBAK and COVL, respectively.

##  5 Conclusions

In this paper, we propose a novel Cycle-in-Cycle GAN dubbed CinCGAN to jointly
recover the spectral magnitude and phase information of clean speech for non-
parallel speech enhancement. In the first stage, a standard magnitude spectrum
estimation CycleGAN is pretained to only estimate the spectral magnitude under
unpaired data. In the second stage, we couple the pretrained CycleGAN with a
complex-valued CycleGAN as a cycle-in-cycle structure to further refine the
magnitude spectrum and recover phase information. More specifically, the
CinCGAN aims at estimating both real and imaginary components of the clean
speech complex spectrum. Comprehensive experiments on Voice Bank + DEMAND
dataset demonstrate that the proposed approach consistently outperforms
original GAN-based and CycleGAN-based baselines under mismatched noisy-clean
pairs. When experiments are conducted on standard parallel data, the proposed
approach also shows its effectiveness in improving speech quality and reducing
speech distortion.

## References

  * [1] P. C. Loizou,  Speech enhancement: theory and practice,  CRC press, 2013.
  * [2] D. L. Wang and J. Chen,  “Supervised speech separation based on deep learning: An overview,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 26, no. 10, pp. 1702–1726, 2018.
  * [3] Y. Wang, A. Narayanan, and D. L. Wang,  “On training targets for supervised speech separation,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 22, no. 12, pp. 1849–1858, 2014.
  * [4] Y. Xu, J. Du, L-R. Dai, and C-H. Lee,  “A regression approach to speech enhancement based on deep neural networks,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 23, no. 1, pp. 7–19, 2014.
  * [5] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio,  “Generative adversarial nets,”  in Proc. NIPS, 2014, pp. 2672–2680.
  * [6] S. Pascual, A. Bonafonte, and J. Serra,  “Segan: Speech enhancement generative adversarial network,”  in Proc. Interspeech, 2017, pp. 3642–3646.
  * [7] D. Baby and S. Verhulst,  “Sergan: Speech enhancement using relativistic generative adversarial networks with gradient penalty,”  in Proc. ICASSP. IEEE, 2019, pp. 106–110.
  * [8] S. W. Fu, C. F. Liao, Y. Tsao, and S. D. Lin,  “Metricgan: Generative adversarial networks based black-box metric scores optimization for speech enhancement,”  in Proc. ICML. PMLR, 2019, pp. 2031–2041.
  * [9] M. H. Soni, N. Shah, and H. A. Patil,  “Time-frequency masking-based speech enhancement using generative adversarial network,”  in Proc. ICASSP. IEEE, 2018, pp. 5039–5043.
  * [10] G. Liu, K. Gong, X. Liang, and Z. Chen,  “Cp-gan: Context pyramid generative adversarial network for speech enhancement,”  in Proc. ICASSP. IEEE, 2020, pp. 6624–6628.
  * [11] Z. Zhang, C. Deng, Y. Shen, D. S. Williamson, Y. Sha, Y. Zhang, H. Song, and X. Li,  “On loss functions and recurrency training for gan-based speech enhancement systems,”  in Proc. Interspeech, 2020, pp. 3266–3270.
  * [12] H. Phan, H. Le Nguyen, O. Y. Chén, P. Koch, N. Q. Duong, I. McLoughlin, and A. Mertins,  “Self-attention generative adversarial network for speech enhancement,”  in Proc. ICASSP. IEEE, 2021, pp. 7103–7107.
  * [13] J. Y. Zhu, T. Park, P. Isola, and A. A Efros,  “Unpaired image-to-image translation using cycle-consistent adversarial networks,”  in Proc. ICCV, 2017, pp. 2223–2232.
  * [14] Z. Meng, J. Li, and Y. Gong,  “Cycle-consistent speech enhancement,”  in Proc. Interspeech 2018, 2018, pp. 1165–1169.
  * [15] Y. Xiang and C. Bao,  “A parallel-data-free speech enhancement method using multi-objective learning cycle-consistent generative adversarial network,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 28, pp. 1826–1838, 2020.
  * [16] Y. Wang, G. Yu, J. Wang, H. Wang, and Q. Zhang,  “Improved relativistic cycle-consistent gan with dilated residual network and multi-attention for speech enhancement,”  IEEE Access, vol. 8, pp. 183272–183285, 2020.
  * [17] G. Yu, Y. Wang, C. Zheng, H. Wang, and Q. Zhang,  “Cyclegan-based non-parallel speech enhancement with an adaptive attention-in-attention mechanism,”  in Proc. Asia-Pacific Signal Inf. Process. Assoc. (APSIPA), 2021, pp. 523–529.
  * [18] G. Yu, Y. Wang, H. Wang, Q. Zhang, and C. Zheng,  “A two-stage complex network using cycle-consistent generative adversarial networks for speech enhancement,”  Speech Commun., vol. 134, pp. 42–54, 2021.
  * [19] K. Paliwal, K. Wójcicki, and B. Shannon,  “The importance of phase in speech enhancement,”  Speech Commun., vol. 53, no. 4, pp. 465–494, 2011.
  * [20] A. Li, W. Liu, C. Zheng, C. Fan, and X. Li,  “Two heads are better than one: A two-stage complex spectral mapping approach for monaural speech enhancement,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 29, pp. 1829–1843, 2021.
  * [21] A. Li, W. Liu, X. Luo, G. Yu, C. Zheng, and X. Li,  “A simultaneous denoising and dereverberation framework with target decoupling,”  in Proc. Interspeech, 2021, pp. 2801–2805.
  * [22] C. Trabelsi, O. Bilaniuk, Y. Zhang, D. Serdyuk, S. Subramanian, J. F. Santos, S. Mehri, N. Rostamzadeh, Y. Bengio, and C. J. Pal,  “Deep complex networks,”  arXiv preprint arXiv:1705.09792, 2017.
  * [23] A. Jolicoeur-Martineau,  “The relativistic discriminator: a key element missing from standard gan,”  arXiv preprint arXiv:1807.00734, 2018.
  * [24] C. Valentini-Botinhao, X. Wang, S. Takaki, and J. Yamagishi,  “Investigating rnn-based speech enhancement methods for noise-robust text-to-speech,”  in Proc. SSW, 2016, pp. 146–152.
  * [25] C. Veaux, J. Yamagishi, and S. King,  “The voice bank corpus: Design, collection and data analysis of a large regional accent speech database,”  in Proc. O-COCOSDA/CASLRE. IEEE, 2013, pp. 1–4.
  * [26] J. Thiemann, N. Ito, and E. Vincent,  “The diverse environments multi-channel acoustic noise database: A database of multichannel environmental noise recordings,”  JASA, vol. 133, no. 5, pp. 3591–3591, 2013.
  * [27] D. Kingma and J. Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint arXiv:1412.6980, 2014.
  * [28] D. Stoller, S. Ewert, and S. Dixon,  “Wave-u-net: A multi-scale neural network for end-to-end audio source separation,”  arXiv preprint arXiv:1806.03185, 2018.
  * [29] F. Germain, Q. Chen, and V. Koltun,  “Speech denoising with deep feature losses,”  in Proc. Interspeech, 2019, pp. 2723–2727.
  * [30] K. Tan and D. L. Wang,  “A convolutional recurrent neural network for real-time speech enhancement.,”  in Proc. Interspeech, 2018, pp. 3229–3233.
  * [31] K. Tan and D. L. Wang,  “Learning complex spectral mapping with gated convolutional recurrent networks for monaural speech enhancement,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 28, pp. 380–390, 2019.
  * [32] Y. Hu, Y. Liu, S. Lv, M. Xing, and L. Xie,  “Dccrn: Deep complex convolution recurrent network for phase-aware speech enhancement,”  in Proc. Interspeech, 2020, pp. 2472–2476.
  * [33] W. Yuan,  “A time–frequency smoothing neural network for speech enhancement,”  Speech Commun, vol. 124, pp. 75–84, 2020.
  * [34] A. Rix, J. Beerends, M. Hollier, and A. Hekstra,  “Perceptual evaluation of speech quality (PESQ)-a new method for speech quality assessment of telephone networks and codecs,”  in Proc. ICASSP. IEEE, 2001, vol. 2, pp. 749–752.
  * [35] C. H. Taal, R. C. Hendriks, R. Heusdens, and J. Jensen,  “A short-time objective intelligibility measure for time-frequency weighted noisy speech,”  in Proc. ICASSP. IEEE, 2010, pp. 4214–4217.
  * [36] Y. Hu and P. C. Loizou,  “Evaluation of objective quality measures for speech enhancement,”  IEEE/ACM Trans. Audio. Speech, Lang. Process., vol. 16, no. 1, pp. 229–238, 2007.
  * [37] C. Reddy, V. Gopal, and R. Cutler,  “Dnsmos: A non-intrusive perceptual objective speech quality metric to evaluate noise suppressors,”  arXiv preprint arXiv:2010.15258, 2020.


# Poiseuille flow for a simplified pseduoplastic rheology

Chris Reese
Department of Physics
Lewis and Clark Community College
Godfrey, IL 62035
<EMAIL_ADDRESS>
###### Abstract

Poiseuille flow in cylindrical and planar geometries with a simplified,
pseudoplastic (shear thinning) rheology characterized by constant viscosity
plateaus above and below a transition strain rate is considered. Analytical,
steady state solutions for velocity profile and volume flux are formulated.
Transient flow development is addressed numerically and compared to the theory
in the steady state limit. Stationary flow is approached after the momentum
diffusion timescale based on the spatially dominant kinematic viscosity. For
large viscosity ratio and shear thinning region confined near the domain
boundary, velocity distributions are quasi-plug like with large boundary to
interior strain rate ratio.

##  1 Introduction

Due to analytical complexity, instances of closed form solutions for
Poiseuille flow with non-Newtonian rheologies are limited. While the
Ostwald–de Waele power-law fluid admits solutions in planar and cylindrical
geometries, the rheological model suffers from lack of effective viscosity
limits at extreme strain rates. A truncated power law model alleviates this
drawback [Lavrov, 2015] and can be adopted to fit Carreau fluid
characteristics [Wrobel, 2020]. Nouar and Frigaard [2009] demonstrate the
existence of stable solutions for combined plane Couette-Poiseuille flow of a
shear thinning Carreau fluid for certain rheological parameters. Analytical
Poiseuille flow solutions have also been obtained for a number of rheological
models in the context of the effects of boundary slip condition [Ferras,
2012]. Volume flow rate in layer and pipe geometries can be calculated for
several generalized Newtonian rheologies utilizing both a variational
principle and the so-called Weissenberg, Rabinowitsch, Mooney and Schofield
method [Sochi, 2014, 2015a, 2015b, 2016]. Griffiths [2020] addresses channel
flow for both shear-thinning and shear-thickening Carreau fluids finding exact
solutions for particular rheological parameters.

In this work, a simplified pseudoplastic rheology [James et al., 2021] is
adopted to investigate planar and cylindrical Poiseuille flow. The generalized
Newtonian viscosity (§2) is a shear thinning, three parameter function
characterized by small and large strain rate viscosity plateaus and a
rheological transition stress. Analytical, steady state solutions for velocity
profile and volume rate are developed in §3. §4 presents numerical results for
flow development from a quiescent initial condition and approach to steady
state. A discussion of results is presented in §5. The appendix (§6)
benchmarks the numerical method utilized in the study against analytical
solutions for isoviscous, transient flow.

##  2 Rheological model

For a generalized Newtonian rheology, effective viscosity is a function of
strain rate. The deviatoric stress tensor

$\boldsymbol{\tau}=2\,\eta_{\rm
eff}(\dot{\varepsilon})\,\boldsymbol{\dot{\varepsilon}}\;,$ (2.1)

where

$\boldsymbol{\dot{\varepsilon}}=\frac{1}{2}\left[\gradient\boldsymbol{u}+\gradient\boldsymbol{u}^{\rm
T}\right]\;,$ (2.2)

is the strain rate tensor and $\boldsymbol{u}$ is fluid velocity. The strain
rate and stress invariants

$\dot{\varepsilon}=\left[\frac{1}{2}\;\Tr\left(\boldsymbol{\dot{\varepsilon}}\boldsymbol{\dot{\varepsilon}}^{\rm
T}\right)\right]^{1/2}\,,\qquad\tau=\left[\frac{1}{2}\;\Tr\left(\boldsymbol{\tau}\boldsymbol{\tau}^{\rm
T}\right)\right]^{1/2}\;.$ (2.3)

The piecewise linear approximation to pseudoplastic rheology [James et al.,
2021] adopted in this study is

$\eta_{\rm
eff}(\dot{\varepsilon})=\begin{cases}\eta&\dot{\varepsilon}\leq\dot{\varepsilon}_{\rm
c}\\\
\dfrac{\tau_{c}}{2\dot{\varepsilon}}\left(1-\dfrac{\mu}{\eta}\right)+\mu&\dot{\varepsilon}\geq\dot{\varepsilon}_{\rm
c}\\\ \end{cases}$ (2.4)

where $\dot{\varepsilon}_{\rm c}$ is the transition strain rate for onset of
shear thinning, $\tau_{c}=2\,\eta\,\dot{\varepsilon}_{\rm c}$ is the
transition stress invariant, $\eta$ is the low strain rate viscosity, and
$\mu$ is the high strain rate viscosity. This rheological model is compared to
the isoviscous case and a Cross fluid in Fig. 1.

Figure 1: A qualitative comparison of various rheologies. (left) The stress
invariant versus strain invariant for isoviscous (blue), simplified
pseudoplastic (red), and Cross fluids (green). All fluids have the same small
strain rate viscosity. The simplified pseduoplastic and Cross fluids approach
the same asymptotic large strain rate viscosity. The parameter $\chi=\eta/\mu$
is the ratio of small strain rate to large strain rate viscosities. Solid and
dashed curves correspond to $\chi$ = 3 and 10, respectively. (right) Effective
viscosities of the rheological models.

##  3 Steady state

For incompressible fluid flow driven solely by a pressure gradient, the Cauchy
momentum equation is

$-\gradient
p+\gradient\cdot\boldsymbol{\tau}=\rho\;\left[\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{u}\cdot\gradient\boldsymbol{u}\right]$ (3.1)

where $p$ is pressure and $\rho$ is fluid density. In this section, planar and
cylindrical steady state solutions are developed for the generalized Newtonian
viscosity described in §2.

###  3.1 Planar geometry

Consider a plane layer perpendicular to the $y$-axis and bounded by no slip
surfaces at $y=0$ and $y=a$. Flow is driven by a uniform pressure gradient,
$p(x)=-\nabla p\,x$ with $\nabla p>0$. In this case, the steady state velocity
$\boldsymbol{u}=u(y)\,\hat{\boldsymbol{i}}$ and is symmetric about the layer
mid-plane. That is, $u(y)$ is maximum at the stationary point $y=a/2$ and
decreases to zero at the boundaries.

The non-zero components of the stress tensor

$\tau_{xy}=2\,\eta_{\rm
eff}(\dot{\varepsilon})\,\dot{\varepsilon}_{xy}=2\,\eta_{\rm
eff}(\dot{\varepsilon})\,\left(\frac{1}{2}\frac{du}{dy}\right)=\eta_{\rm
eff}(\dot{\varepsilon})\,\frac{du}{dy}$ (3.2)

The momentum equation reduces to

$\frac{dp}{dx}=-\nabla p=\frac{d\tau_{xy}}{dy\hfill}$ (3.3)

Integrating and applying the symmetry condition $\tau_{xy}=0$ at $y=a/2$
yields the shear stress

$\tau_{xy}=-\nabla p\left(y-\frac{a}{2}\right)$ (3.4)

The stress invariant $\tau(y)=\nabla p\left|y-a/2\right|$ increases linearly
with distance from the mid-plane with maximum stress invariant at the no-slip
layer boundaries, $\tau(0)=\tau(a)=\nabla p\;a/2$ (Fig. 2). For transition
stress $\tau_{c}\geq\tau(a)$, the entire layer is in the high viscosity
rheological regime. For $\tau_{c}<\tau(a)$, there are two $y$ positions where
$\tau=\tau_{c}$. Denoting the upper half-plane position where this occurs as
$\xi$ gives

$\xi=\frac{a}{2}+\frac{\tau_{c}}{\nabla p}\;.$ (3.5)

In the lower half-plane, the position below which the stress invariant exceeds
$\tau_{c}$, is $y=a-\xi$ (Fig. 2).

Figure 2: An example planar solution with $\chi=3$ and $\xi=0.75$. (left)
Stress invariant (blue) is symmetric about the mid-layer (black dashed)
increasing linearly from zero to maxima $\nabla p\,a/2$ at the boundaries. For
$\tau_{c}<\tau(a)$, shear thinning regions exist for $y>\xi$ and $y<a-\xi$
(green dashed). (middle) Effective viscosity (blue). (right) Velocity profile
(blue). In the upper half layer the velocities above and below the shear
thinning transition are denoted by $u_{+}$ and $u_{-}$, respectively.

The velocity profile $u(y)$ is formulated for the upper half-plane, $a/2\leq
y\leq a$. Reflection symmetry about the mid-plane gives the velocity in the
lower half-plane. Let $u_{-}$ be the velocity for $a/2\leq y<\xi$ where
$\tau<\tau_{c}$ and $\eta_{\rm eff}=\eta$. Then, Eq. (3.4) becomes

$\eta\frac{du_{-}}{dy\hfill}=-\nabla p\left(y-\frac{a}{2}\right)$ (3.6)

Integrating

$u_{-}=\frac{\nabla
p\,a^{2}}{2\eta}\,\frac{y}{a}\left(1-\frac{y}{a}\right)+k_{1}$ (3.7)

where $k_{1}$ is a constant. Let $u_{+}$ be the velocity above the viscosity
transition level, i.e., for $\xi\leq y\leq a$. Upon substituting for the
strain rate invariant, $\dot{\varepsilon}=1/2\,|du_{+}/dy|$ in the effective
viscosity, Eq. (3.4) yields,

$\frac{\eta}{\chi}\left[\frac{\tau_{c}}{\eta\left|\dfrac{du_{+}}{dy\hfill}\right|}\left(\chi-1\right)+1\right]\dfrac{du_{+}}{dy\hfill}=-\nabla
p\left(y-\frac{a}{2}\right)$ (3.8)

where $\chi=\eta/\mu$ is the high to low viscosity ratio. In the upper half
plane, the velocity gradient is monotonic and less than zero, i.e.,
$|du_{+}/dy|=-du_{+}/dy$. Substituting, eliminating $\tau_{c}$ using Eq.
(3.5), integrating, and applying the no slip boundary condition $u_{+}(a)=0$
yields

$u_{+}=\frac{\nabla
p\,a^{2}}{2\eta}\left[\chi\frac{y}{a}\left(1-\frac{y}{a}\right)-2(\chi-1)\left(\frac{\xi}{a}-\frac{1}{2}\right)\left(1-\frac{y}{a}\right)\right]$
(3.9)

Matching solutions, $u_{-}(\xi)=u_{+}(\xi)$ yields the constant $k_{1}$ and

$u_{-}=\frac{\nabla
p\,a^{2}}{2\eta}\left[\frac{y}{a}\left(1-\frac{y}{a}\right)+(\chi-1)\left(1-\frac{\xi}{a}\right)^{2}\right]$
(3.10)

Velocity profiles for some values of $\chi$ and $\xi$ are shown in Fig. (3).

Various limits can be checked against isoviscous cases. For $\eta=\mu$,
$\chi=1$ and $u_{+}=u_{-}=u_{\rm iso}$ where the isoviscous solution

$u_{\rm iso}=\frac{\nabla
p\,a^{2}}{2\eta}\,\frac{y}{a}\left(1-\frac{y}{a}\right)$ (3.11)

When $\tau_{c}=\tau(a)$, $\xi=a$, the entire layer is in the high viscosity
regime, and $u_{-}=u_{\rm iso}$. In the large strain rate limit,
$\tau_{c}\ll\tau(a)$ and $\xi\rightarrow a/2^{+}$. Apart from a narrow, small
strain rate sublayer around $y=a/2$, the rest of the layer is in the low
viscosity regime. In this case, $u_{+}\approx\chi\,u_{\rm iso}$ which
corresponds to isoviscous flow with viscosity $\mu$.

Figure 3: Velocity profile dependence on viscosity ratio $\chi$ and
rheological transition level $\xi$. In this graph, the coordinate $y$ is non-
dimensionalized by layer thickness $a$ and the velocity has been scaled by the
maximum isoviscous velocity $u_{\rm iso,max}=(\nabla p\,a^{2})/(8\eta)$.
(left) The isoviscous case (solid black line) is compared to pseudoplastic
solutions (blue lines) for $\chi=3$ and $\xi$ = (0.825, 0.750, 0.625, 0.500).
Layer regions where $\tau>\tau_{c}$ are denoted by dashed blue lines. (center)
Velocity profiles for $\chi$ = 10 and the same values of $\xi$. (right) Volume
rate (per unit length in the $z$ direction) as a function of rheological
transition position $\xi$. The isoviscous case (dashed black line) is compared
to pseudoplastic solutions for $\chi$ = (3,5,10) shown in (blue, green, red)
respectively.

The volume rate (per unit length in the $z$ direction)

$Q=2\left[\int_{a/2}^{\xi}u_{-}(y)\,dy+\int_{\xi}^{a}u_{+}(y)\,dy\right]$
(3.12)

Substituting for the velocities and evaluating yields

$\frac{Q\hfill}{Q_{\rm
iso}}=\left(2\,\chi-1\right)-2\,(\chi-1)\,\frac{\xi^{2}}{a^{2}}\left(3-2\,\frac{\xi}{a}\right)$
(3.13)

where

$Q_{\rm iso}=\frac{\nabla p\,a^{3}}{12\eta}$ (3.14)

Checking the isoviscous limit $\eta=\mu$, $\chi=1$, yields $Q=Q_{\rm iso}$.
Likewise, for $\tau_{c}=\tau(a)$, $\xi=a$, and $Q=Q_{\rm iso}$. Finally, in
the large strain rate limit $\tau_{c}\ll\tau(a)$, $\xi\rightarrow a/2^{+}$,
and $Q\approx\chi\,Q_{\rm iso}$.

###  3.2 Cylindrical geometry

Consider cylindrical conduit geometry with axial coordinate $z$ and no slip
surface at radius $r=R$. Flow is driven by a uniform pressure gradient,
$p(z)=-\nabla p\,z$ with $\nabla p>0$. In this case, the steady state velocity
$\boldsymbol{u}=u(r)\,\hat{\boldsymbol{z}}$.

The non-zero components of the stress tensor

$\tau_{rz}=2\,\eta_{\rm eff}\,\dot{\varepsilon}_{rz}=2\,\eta_{\rm
eff}\,\left(\frac{1}{2}\frac{du}{dr}\right)=\eta_{\rm eff}\,\frac{du}{dr}$
(3.15)

The momentum equation reduces to

$\frac{dp}{dz}=-\nabla p=\frac{1}{r}\frac{d\hfill}{dr}\left(r\tau_{rz}\right)$
(3.16)

Integrating and requiring that $\tau_{rz}$ be finite at $r=0$ yields

$\tau_{rz}=-\frac{\nabla p}{2}\,r$ (3.17)

The stress invariant $\tau(r)=\nabla p/2\,r$ is an increasing linear function
of radial position with maximum at the conduit boundary $\tau(R)=\nabla
p/2\,R$. For $\tau(R)<\tau_{c}$, there is a radial position $\xi<R$ where the
stress invariant equals the rheological transition stress,

$\xi=\frac{2\,\tau_{c}}{\nabla p}$ (3.18)

For $r<\xi$, the effective viscosity $\eta_{\rm eff}=\eta$. As in the planar
geometry case, denote the velocity in this region as $u_{-}$. Then,
substituting for $\tau_{rz}$ in Eq. (3.17) gives

$\eta\,\frac{du_{-}}{dr\hfill}=-\frac{\nabla p}{2}\,r$ (3.19)

Integrating,

$u_{-}=-\frac{\nabla p\,R^{2}}{4\,\eta}\,\left(\frac{r}{R}\right)^{2}+c_{1}$
(3.20)

where $c_{1}$ is a constant. For $r>\xi$, the effective viscosity is given by
Eq. (2.4) with strain rate invariant $\frac{1}{2}\,|du_{+}/dr|$. Letting
$u_{+}$ be the velocity in the region $\xi\leq r\leq R$, Eq. (3.17) becomes

$\frac{\eta}{\chi}\left[\frac{\tau_{c}}{\eta\left|\dfrac{du_{+}}{dr\hfill}\right|}\left(\chi-1\right)+1\right]\dfrac{du_{+}}{dr\hfill}=-\frac{\nabla
p}{2}\,r$ (3.21)

The velocity gradient is a monotonically decreasing function of $r$, i.e.,
$|du_{+}/dr|=-du_{+}/dr$. Making this substitution, eliminating $\tau_{c}$
using Eq. (3.18), integrating, and applying the no slip boundary condition
$u_{+}(R)=0$ yields

$u_{+}=\frac{\nabla
p\,R^{2}}{4\,\eta}\,\left[\chi\left(1-\frac{r^{2}}{R^{2}}\right)-2\,(\chi-1)\,\frac{\xi}{R}\,\left(1-\frac{r}{R}\right)\right]$
(3.22)

Matching solutions $u_{-}(\xi)=u_{+}(\xi)$ gives the constant $c_{1}$,

$u_{-}=\frac{\nabla
p\,R^{2}}{4\,\eta}\,\left[\chi+(\chi-1)\left(\frac{\xi^{2}}{R^{2}}-2\frac{\xi}{R}\right)-\frac{r^{2}}{R^{2}}\right]$
(3.23)

In the isoviscous limit $\eta=\mu$ and $\chi=1$, $u_{+}=u_{-}=u_{\rm iso}$
where the isoviscous velocity

$u_{\rm iso}=\frac{\nabla
p\,R^{2}}{4\,\eta}\left(1-\frac{r^{2}}{R^{2}}\right)$ (3.24)

When $\tau_{c}=\tau(R)$, $\xi=R$, the entire conduit is in the high viscosity
regime, and $u_{-}=u_{\rm iso}$. In the large strain rate limit,
$\tau_{c}\ll\tau(R)$ and $\xi\ll R$. Apart from a small strain rate region
around $r=0$, the rest of the fluid is in the low viscosity regime. In this
case, $u_{+}\approx\chi\,u_{\rm iso}$ which corresponds to isoviscous flow
with viscosity $\mu$.

Figure 4: Velocity profile dependence on viscosity ratio $\chi$ and
rheological transition level $\xi$. The velocity is scaled by the maximum
isoviscous velocity $u_{\rm iso,max}=(\nabla p\,R^{2})/(4\eta)$. (left) The
isoviscous case (solid black line) is compared to pseudoplastic solutions
(blue lines) for $\chi=3$ and $\xi$ = (0.75, 0.50, 0.25, 0.00). Layer regions
where $\tau>\tau_{c}$ are denoted by dashed blue lines. (center) Velocity
profiles for $\chi$ = 10 and the same values of $\xi$. (right) Volume rate as
a function of rheological transition radius $\xi$. The isoviscous case (dashed
black line) is compared to pseudoplastic solutions for $\chi$ = (3,5,10) shown
in (blue, green, red) respectively.

The volume rate

$Q=2\pi\left[\int_{0}^{\xi}r\,u_{-}(r)\,dr+\int_{\xi}^{R}r\,u_{+}(r)\,dr\right]$
(3.25)

Substituting for the velocities and integrating yields

$\frac{Q\hfill}{Q_{\rm
iso}}=\chi-\frac{4}{3}\,(\chi-1)\,\frac{\xi}{R}\left(1-\frac{1}{4}\,\frac{\xi^{3}}{R^{3}}\right)$
(3.26)

where

$Q_{\rm iso}=\frac{\pi}{8}\,\frac{\nabla p\,R^{4}}{\eta}$ (3.27)

Checking the isoviscous limit $\eta=\mu$, $\chi=1$, yields $Q=Q_{\rm iso}$.
Likewise, for $\tau_{c}=\tau(R)$, $\xi=R$, and $Q=Q_{\rm iso}$. Finally, in
the large strain rate limit $\tau_{c}\ll\tau(R)$, $\xi\rightarrow 0$, and
$Q\approx\chi\,Q_{\rm iso}$.

Figure 5: Velocity evolution (blue) and steady state velocities (black) in
planar geometry. Velocities are scaled by the maximum steady state value
$u_{\rm max}=1/8\,[1+4(\chi-1)(1-\xi)^{2}]$. The location of the rheological
transition $\dot{\varepsilon}=\dot{\varepsilon}_{c}$ is shown in green (bold
line/dots). (top) Velocity profiles for $\chi=3$, steady state values $\xi$ =
0.825, 0.750, 0.625 (left, middle, right), and times $t$ = 0.1, 0.2, 0.3, 0.5.
(bottom) As for the top with $\chi=10$ and times $t$ = 0.3, 0,5, 1.0, 1.5.

3

##  4 Transient flow

In this section, numerical simulations of flow development from a static
initial state are considered. The numerical method is benchmarked against
isoviscous theory for planar and cylindrical geometries (see §6). The approach
to steady state is compared to theoretical results developed in §3. In
particular, velocity distribution, extent of shear thinning, and volume rate
are considered.

Figure 6: (top) Evolution of the rheological transition location (left) and
volume rate (right) for $\chi$ = 3. The three cases correspond to steady state
values $\xi$ = 0.825, 0.750, 0.625 (blue, green, red). The volume rate is
normalized by the isoviscous value based on $\eta$ (Eq. 3.14). Steady values
are shown with dashed lines. (bottom) As for the top with $\chi$ = 10.

###  4.1 Planar geometry

For the unidirectional flow considered here, the velocity gradient is
perpendicular to the velocity resulting in a vanishing advective derivative
$\boldsymbol{u}\cdot\gradient\boldsymbol{u}$. Substituing for $\tau_{xy}$ and
pressure $p(x)$ in the momentum equation Eq. (3.1),

$\frac{\partial u}{\partial t}=\frac{\partial}{\partial y}\left(\nu_{\rm
eff}(\dot{\varepsilon})\,\frac{\partial u}{\partial
y}\right)+\frac{1}{\rho}\,\nabla p$ (4.1)

where $\nu_{\rm eff}(\dot{\varepsilon})=\eta_{\rm eff}/\rho$ is the effective
kinematic viscosity. Rewriting $\eta_{\rm eff}$ in terms of the viscosity
ratio $\chi$,

$\nu_{\rm
eff}(\dot{\varepsilon})=\begin{cases}\nu&\dot{\varepsilon}\leq\dot{\varepsilon}_{\rm
c}\\\
\dfrac{\nu}{\chi}\left[\dfrac{\tau_{c}}{2\,\eta\,\dot{\varepsilon}}\left(\chi-1\right)+1\right]&\dot{\varepsilon}\geq\dot{\varepsilon}_{\rm
c}\\\ \end{cases}$ (4.2)

where $\nu=\eta/\rho$. Non-dimensionalization with the following scales

$[y]=a\qquad[t]=\frac{a^{2}}{\nu}\qquad[u]=\frac{\nabla
p\,a^{2}}{\rho\,\nu}\qquad[\dot{\varepsilon}]=\frac{\nabla
p\,a}{\rho\,\nu}\qquad[\tau]=\nabla p\,a$ (4.3)

yields the momentum equation

$\frac{\partial u}{\partial t}=\frac{\partial}{\partial y}\left(\nu_{\rm
eff}^{\prime}(\dot{\varepsilon})\,\frac{\partial u}{\partial y}\right)+1$
(4.4)

where

$\nu_{\rm
eff}^{\prime}(\dot{\varepsilon})=\begin{cases}1&\dot{\varepsilon}\leq\dot{\varepsilon}_{\rm
c}\\\
\dfrac{1}{\chi}\left[\dfrac{\tau_{c}}{2\,\dot{\varepsilon}}\left(\chi-1\right)+1\right]&\dot{\varepsilon}\geq\dot{\varepsilon}_{\rm
c}\\\ \end{cases}$ (4.5)

The non-dimensional, transition stress invariant $\tau_{c}=\xi-1/2$. The
numerical computation was initialized with a small parabolic velocity profile
$u(y,0)=\epsilon(y-y^{2})$ with $\epsilon$ = 10-4. The spatial discretization
interval $\delta y$ = 0.05.

Channel velocity evolution for $\chi$ = 3,10 and $\xi$ = 0.625, 0.750, 0.825
is shown in Fig. (5). Depending on $\chi$ and $\xi$, steady state is
approached after non-dimensional times $\sim$ 0.5-2. For large rheological
transition stress, the fluid is essentially isoviscous with $\eta_{\rm
eff}\sim\eta$ and steady state is reached after non-dimensional time $\sim$
0.5. When the bulk of the fluid is in the low viscosity regime (small $\xi$)
and $\chi$ is large, the flow adjustment time is longer. In the limit of fully
developed shear thinning ($\xi\rightarrow$ 0), the momentum diffusion time is
governed by $\eta_{\rm eff}\sim\mu$ and time to reach steady state scales
linearly with $\chi$. Pseudoplastic transition location and volume rate as a
function of time are shown in Fig. (6).

Figure 7: Velocity evolution (blue) and steady state velocities (black).
Velocities are scaled by the maximum steady state value $u_{\rm
max}=1/4\,[\chi-(\chi-1)\,\xi\,(2-\xi)]$. The location of the rheological
transition $\dot{\varepsilon}=\dot{\varepsilon}_{c}$ is shown in green (bold
line/dots). (top) Velocity profiles for $\chi=3$, steady state values $\xi$ =
0.75, 0.50, 0.25 (left, middle, right), and times $t$ = 0.1, 0.3, 0.5, 1.0.
(bottom) As for the top with $\chi=10$ and times $t$ = 0.3, 1.0, 2.0, 3.0.

###  4.2 Cylindrical geometry

In cylindrical geometry, the transient momentum equation is

$\frac{\partial u}{\partial t}=\frac{1}{r}\,\frac{\partial}{\partial
r}\left(r\,\nu_{\rm eff}(\dot{\varepsilon})\,\frac{\partial u}{\partial
r}\right)+\frac{1}{\rho}\,\nabla p$ (4.6)

with $\nu_{\rm eff}$ given by Eq. (4.2). Non-dimensionalizing with the
following scales

$[r]=R\qquad[t]=\frac{R^{2}}{\nu}\qquad[u]=\frac{\nabla
p\,R^{2}}{\rho\,\nu}\qquad[\dot{\varepsilon}]=\frac{\nabla
p\,R}{\rho\,\nu}\qquad[\tau]=\nabla p\,R$ (4.7)

yields

$\frac{\partial u}{\partial t}=\frac{1}{r}\,\frac{\partial}{\partial
r}\left(r\,\nu_{\rm eff}^{\prime}(\dot{\varepsilon})\,\frac{\partial
u}{\partial r}\right)+1$ (4.8)

where $\nu_{\rm eff}^{\prime}$ is given by Eq. (4.5). The non-dimensional,
rheology transition stress invariant $\tau_{c}=\xi/2$. The computation was
initialized with a small parabolic velocity profile satisfying the no slip
boundary condition, $u(r,0)=\epsilon(1-r^{2})$ with $\epsilon$ = 10-4. The
spatial discretization interval $\delta r$ = 0.05.

Velocity evolution for $\chi$ = 3,10 and $\xi$ = 0.75, 0.50, 0.25 is shown in
Fig. (7) for pipe flow. Velocities approach steady state after non-dimensional
times $\sim$ 1-3 depending on spatial extent of shear thinning and viscosity
contrast. For sufficiently large $\tau_{c}$, the fluid is essentially
isoviscous with $\eta_{\rm eff}\sim\eta$ and the momentum diffusion time is
$\sim$ 1\. When more of the fluid is in the low viscosity regime (small $\xi$)
and the viscosity contrast $\chi$ is large, the flow adjustment time is
longer. In the limit of fully developed shear thinning ($\xi\rightarrow$ 0),
momentum diffusion is governed by $\eta_{\rm eff}\sim\mu$ and the flow
adjustment timescale $\sim\chi$. Evolution of pseudoplastic transition
location and volume rate are shown in Fig. (8).

Figure 8: (top) Evolution of the rheological transition location (left) and
volume rate (right) for $\chi$ = 3. The three cases correspond to steady state
values $\xi$ = 0.75, 0.50, 0.25 (blue, green, red). The volume rate is
normalized by the isoviscous value based on $\eta$ (Eq. 3.26). Steady values
are shown with dashed lines. (bottom) As for the top with $\chi$ = 10.

##  5 Discussion

In this paper, pseudoplastic fluid flow in flat, channel and pipe geometries
is considered. A simplified model which captures fundamental characteristics
of more complicated shear thinning fluid rheologies is adopted. The model
admits closed form solutions for steady state velocity distributions expressed
as piecewise polynomial functions in spatial coordinates. Likewise, volume
rates are simple analytical functions of rheological parameters and driving
pressure gradient. Analytical simplicity of the results may prove useful for
parameter space exploration prior to adoption of more sophisticated rheologies
and/or appeal to finite element simulations.

Transient flow evolution from a static state is investigated via numerical
computation. It is found that model fluid flow approaches theoretical steady
state on timescales governed by the spatially dominant effective kinematic
viscosity. That is to say, the reduced effective viscosity governs stationary
flow development timescale when shear thinning occurs throughout the majority
of the fluid. In the steady state limit, numerical results are in good
agreement with theory.

Figure 9: Planar, non-dimensional velocity profiles for isoviscous (black) and
pseduoplastic (blue) cases with equal average (plug) velocity $\bar{u}$. The
pseudoplastic case has $\chi$ = 100 and $\xi$ = 0.9 and shear thinning regions
are denoted by dashed lines.

For large viscosity ratio $\chi$ and shear thinning regions confined near no-
slip boundaries the flow can exhibit plug-like characteristics. This is
illustrated in Fig. 9 for a plane layer with shear thinning case $\chi$ = 100
and $\xi$ = 0.9 such that narrow, large strain-rate layers exist adjacent to
domain boundaries. The velocity profile is compared to the isoviscous case for
equal average velocities.

The simplified rheology and results considered in this work may find
application in models of dike emplacement and/or magma flow in volcanic
conduits [Gonnermann and Manga, 2007] as silicic magmas exhibit shear thinning
behavior [Lavallée et al., 2007; Jones et al.,, 2020; Vasseur et al., 2023].
The model may also prove useful in industrial applications when volume rate
control of generalized Newtonian fluids is of interest (e.g., in the case of
polymeric liquids [Bird et al., 1987]).

##  6 Appendix: Numerical methods

In this appendix, the numerical method is benchmarked against isoviscous
analytical solutions for transient development of Poiseuille flow in planar
and cylindrical geometries from zero velocity to steady state. The transient
momentum equations (Eqs. 4.4, 4.8) are of the form

$\frac{\partial u(r,t)}{\partial t}=\mathcal{D}\left[u(r,t)\right]$ (6.1)

where $\mathcal{D}$ is a non-linear differential operator. The Python package
py-pde [Zwicker, 2020] provides methods for solving partial differential
equations of this form. The py-pde, method-of-lines scheme utilizes implicit,
Adams backward differentiation [Hindmarsh, 1983; Petzold, 1983].

###  6.1 Planar geometry

Consider flow start-up from an initial zero velocity state $u(y,0)=0$ subject
to no slip conditions at the channel boundaries, $u(0,t)=u(1,t)=0$. The non-
dimensional, analytical solution of the inhomogeneous diffusion equation (Eq.
4.4) with $\nu^{\prime}_{\rm eff}=1$ subject to these initial and boundary
conditions is

$u(y,t)=4\,\sum_{n\;{\rm odd}}\frac{\sin n\pi
y}{(n\pi)^{3}}\left(1-e^{-n^{2}\pi^{2}t}\right)\;.$ (6.2)

Integrating the velocity profile over the layer thickness gives the volume
rate

$Q(t)=8\,\sum_{n\;{\rm odd}}\frac{1-e^{-n^{2}\pi^{2}t}}{(n\pi)^{4}}\;.$ (6.3)

For comparison to computational results, sums were truncated at 11 non-zero
terms. For non-dimensional time $t$ = 1, the next series term for the volume
rate is $\sim$ 10-5 times the partial sum.

Discretization and initial condition are the same as for the pseudoplastic
computations (§4.1). Numerical and analytical results are in good agreement
during development of steady state flow (Fig. 10). Velocity distribution and
volume flux solutions approach steady state solutions at time $t\sim$ 0.5. For
time $t$ = 1.0 the absolute error between analytical and numerical velocity
profiles integrated over the layer is $\sim$ 10-5.

Figure 10: Benchmark for startup of isoviscous planar Poiseuille flow. (left)
Numerical (green solid line) and analytical (blue dashed line) velocity
profiles at non-dimensional times 0.05, 0.15, 0.25, and 0.50. Solutions
approach the steady state analytical profile (red dashed line) at time $t\sim$
0.5. (right) Numerical (green squares/line) and analytical (blue circles)
transient volume rate. The non-dimensional steady state value $Q=$ 1/12 is
shown with a dashed black horizontal line.

###  6.2 Cylindrical geometry

For initial condition $u(r,0)=0$ and no slip boundary condition $u(1,t)=0$,
the non-dimensional, analytical velocity profile

$u(r,t)=\frac{1}{4}\left(1-r^{2}\right)-2\,\sum_{n}\frac{J_{0}(\lambda_{n}r)}{\lambda_{n}^{3}\;J_{1}(\lambda_{n})}e^{-\lambda_{n}^{2}t}\;,$
(6.4)

where $J_{0}$ and $J_{1}$ are Bessel functions of the first kind of order zero
and one respectively and $\lambda_{n}$ are the positive roots of $J_{0}$. The
associated volume rate

$Q(t)=\frac{\pi}{8}-4\pi\sum_{n}\frac{e^{-\lambda_{n}^{2}t}}{\lambda_{n}^{2}}$
(6.5)

For comparison to computational results, sums were truncated at 10 terms. For
non-dimensional time $t$ = 1, the next series term for the volume rate is
vanishingly small.

Discretization and initial condition are the same as for the pseudoplastic
computations (§4.2). Numerical and analytical results are in good agreement
during development of steady state flow (Fig. 11). Velocity distribution and
volume flux solutions approach steady state solutions at time $t\sim$ 1.0. For
time $t$ = 1.0 the absolute error between analytical and numerical velocity
profiles integrated over the layer is $\sim$ 3 $\times$ 10-4.

Figure 11: Benchmark for startup of isoviscous cylindrical Poiseuille flow.
(left) Numerical (green solid line) and analytical (blue dashed line) velocity
profiles at non-dimensional times 0.05, 0.15, 0.25, 0.50, and 1.0. Solutions
approach the steady state analytical profile (red dashed line) at time $t\sim$
1.0. (right) Numerical (green squares/line) and analytical (blue circles)
transient volume rate. The non-dimensional steady state value $Q=\pi/8$ is
shown with a dashed black horizontal line.

## References

  * Bird et al. [1987] Bird, R., Armstrong, R., and Hassager, O. (1987), Dynamics of Polymeric Liquids. Wiley.
  * Ferras [2012] Ferrás, L.L., Nobréga, J.M., and Pinho, F.T., Analytical solutions for Newtonian and inelastic non-Newtonian flows with wall slip, J. Non-Newt. Fluid Mech., 175-176, 76-88.
  * Gonnermann and Manga [2007] Gonnermann, H.M. and Manga, M. (2007), The fluid mechanics inside a volcano, Ann. Rev. Fluid Mech., 39, 321-356.
  * Griffiths [2020] Griffiths, P.T. (2020), Non-Newtonian channel flow—exact solutions, IMA J. Appl. Math., 85, 263-279.
  * Hindmarsh [1983] Hindmarsh, A.C. (1983), ODEPACK, A systemized collection of ODE solvers, Scientific Computing (ed. R.S. Stepleman et al.), IMACS Trans. Sci. Comp. 1, 55-64.
  * James et al. [2021] James, F., M’ Baye, M., Msheik, K. and Nguyen, D. (2021), A lubrication equation for a simplified model of shear thinning fluid, ESAIM: Proceedings and Surveys, 70, pp.147-165.
  * Jones et al., [2020] Jones, T.J., Llewellini, E.W., and Mader, H.M., (2020), The use of a shear-thinning polymer as a bubbly magma analogue for scaled laboratory experiments, J. Volcanol. Geotherm. Res., 392, 106768.
  * Lavallée et al. [2007] Lavallée, Y., Hess, K.U., Cordonnier, B., and Dingwill, D.B. (2007), Non-Newtonian rheological law for highly crystalline dome lavas, Geology, 35, 843-846.
  * Lavrov [2015] Lavrov A., Flow of truncated power-law fluid between parallel walls for hydraulic fracturing applications, J. Non-Newt. Fluid Mech., 223, 141-146.
  * Nouar and Frigaard [2009] Nouar, C. and Frigaard, I. (2009), Stability of plane Couette–Poiseuille flow of shear-thinning fluid, Phys. Fluids, 21, 064104.
  * Petzold [1983] Petzold, L. (1983), Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations, SIAM J. Sci. Stat. Comput., 4, 136-148.
  * Sochi [2014] Sochi, T. (2014), Using the Euler–Lagrange variational principle to obtain flow relations for generalized Newtonian fluids, Rheol. Acta, 53, 15-22.
  * Sochi [2015a] Sochi, T. (2015a), Further validation to the variational method to obtain flow relations for generalized Newtonian fluids, Korea-Aust. Rheol. J., 27, 113–124.
  * Sochi [2015b] Sochi, T. (2015b), Analytical solutions for the flow of Carreau and Cross fluids in circular pipes and thin slits, Rheol. Acta, 54, 745-756.
  * Sochi [2016] Sochi, T. (2016), Variational approach for the flow of Ree–Eyring and Casson fluids in pipes, Int. J. Model. Simul. Sci. Comput., 7, 1650007.
  * Vasseur et al. [2023] Vasseur, J., Wadsworth, F.B., and Dingwell, D.B. (2023), Shear thinning and brittle failure in crystal-bearing magmas arise from local non-Newtonian effects in the melt, Earth Planet. Sci., 603, 117988.
  * Wrobel [2020] Wrobel, M. (2020), An efficient algorithm of solution for the flow of generalized Newtonian fluid in channels of simple geometries, Rheol. Acta, 59, 651-663.
  * Zwicker [2020] Zwicker, D. (2020), py-pde: A Python package for solving partial differential equations, J. Open Source Softw., 5(48), 2158.


# VP-SLAM: A Monocular Real-time Visual SLAM with Points, Lines and Vanishing
Points

Andreas Georgis  School of ECE, NTUA
Athens, Greece
<EMAIL_ADDRESS>Panagiotis Mermigkas  School of ECE, NTUA
Athens, Greece
<EMAIL_ADDRESS>Prof. Petros Maragos  School of ECE, NTUA
Athens, Greece
<EMAIL_ADDRESS>
###### Abstract

Traditional monocular Visual Simultaneous Localization and Mapping (vSLAM)
systems can be divided into three categories: those that use features, those
that rely on the image itself, and hybrid models. In the case of feature-based
methods, new research has evolved to incorporate more information from their
environment using geometric primitives beyond points, such as lines and
planes. This is because in many environments, which are man-made environments,
characterized as ”Manhattan world”, geometric primitives such as lines and
planes occupy most of the space in the environment. The exploitation of these
schemes can lead to the introduction of algorithms capable of optimizing the
trajectory of a Visual SLAM system and also helping to construct an exuberant
map. Thus, we present a real-time monocular Visual SLAM system that
incorporates real-time methods for line and VP extraction, as well as two
strategies that exploit vanishing points to estimate the robot’s translation
and improve its rotation.Particularly, we build on ORB-SLAM2, which is
considered the current state-of-the-art solution in terms of both accuracy and
efficiency, and extend its formulation to handle lines and VPs to create two
strategies the first optimize the rotation and the second refine the
translation part from the known rotation. First, we extract VPs using a real-
time method and use them for a global rotation optimization strategy. Second,
we present a translation estimation method that takes advantage of last-stage
rotation optimization to model a linear system. Finally, we evaluate our
system on the TUM RGB-D benchmark and demonstrate that the proposed system
achieves state-of-the-art results and runs in real time, and its performance
remains close to the original ORB-SLAM2 system.

##  I Introduction

Visual SLAM (vSLAM) systems try to estimate a robot’s location based on the
multi-view geometry of the scene, combined with computer vision algorithms,
while generating a 3D map of the environment. It is a critical tool for 3D
reconstruction, image refinement, 3D holographic application, visual place
recognition, AR/VR reality, and autonomous vehicles like micro air vehicles
(MAVs). Various vSLAM approaches have been created based on various sensors
such as single camera, stereo camera, RGB-D camera and event-camera. Moreover,
feature-based approaches have traditionally attracted the most attention,
since, they rely on well-suited computer vision algorithms to extract features
and are more resistant to changes in light than direct methods. However, in
low-textured or man-made environments where the extracted point features are
not well-distributed or sufficient, incorporating other geometric primitives
from multi-view geometry into the system, such as lines planes, or VPs
(Vanishing Points) can boost the robustness of these systems. [1], [2], [3],
[4].

On the other hand, most applications in practice work on certain scenarios,
such as man-made environments.In these environments, to boost the performance
of the vSLAM system, the hypothesis of the Manhattan World (_MW_) is used [5].
The _MW_ , is a man-made environment with significant structural regularity,
and with most areas of surrounding environment being described as a box world
with three mutual orthogonal dominant directions. As a result, each _MW_ is
associated with a frame, which is denoted as Manhattan Frame (_MF_) and can be
inferred from the VPs, which is the intersection of the image projections of
the 3D parallel lines in _MW_. As a result, employing the VPs can lead to a
reduction in pose drift [6], [7].

Motivated by the insights made above, we present a vSLAM system that
integrates simple computer vision algorithms for extracting lines, and VPs to
reduce the drift in pose and optimize it. More specifically, the main
contributions of this paper are summarized below::

  * •

Employ real-time computer vision algorithms for extract- ing lines and VPs..

  * •

An optimization strategy of absolute rotation based on VPs.

  * •

A simple linear system for estimate translation.

##  II Related Work

Next we briefly review the related works on vSLAM system with the focus on
features and on leveraging the structural regularity to improve system
performance. First, we have ORB-SLAM2 [8] which is a popular feature-based
monocular vSLAM system, that extends the multi-threaded and keyframe-based
architecture of PTAM [9]. It uses ORB features, builds a co-visibility graph
and performs loop closing and localization tasks.To improve the robustness of
point-based methods, the authors in [9] extracted lines from the environment
and propose an algorithm to integrate them into a monocular Extended Kalman
Filter SLAM system (EKF-SLAM). Finally, in PL-SLAM [1] points and lines
extracted concurrently into a point-based system.

On the other hand, we have vSLAM systems based on the _MW_ assumption, like
[6] which proposes a 2 stage _MF_ tracking to estimate the rotation of the
pose based on the VPs and parallel lines and a refinement strategy to pose
optimization. [10] proposes a mean-shift algorithm to track the rotation of
_MF_ across scenes, while using 1-D density alignments for translation
estimation. OPVO [11] improves the translation estimation by using the KLT
tracker. Both methods require two planes to be visible in the frame at all
times. LPVO [12] eases this requirement by incorporating lines into the
system. Structural lines are aligned with the axes of the dominant _MF_ and
can be integrated into the mean shift algorithm, improving robustness. Hence,
for LPVO, only a single plane is required to be visible in the scene, given
the presence of lines. Drift can still occur in translation estimation, as it
relies on frame-to-frame tracking.

##  III System Overview

The proposed vSLAM system, addresses three major challenges: real-time
complexity, optimize global rotation utilizing orthogonality and parallelism
constraints, and a simple translation refinement using the information of
global rotation. So, in the following paragraphs, each module of the proposed
VP-SLAM system is presented briefly.

The proposed VP-SLAM system is illustrated in “Fig. 1”, which is based on the
ORB-SLAM2 [8] architecture. The system is composed of two modules: the front-
end and the back-end. The front-end is responsible for real-time VO, providing
ego-motion estimation, with local optimization for the current frame and
keyframe decision, while the back-end is responsible for the map
representation, local map optimization, and the strategy for inserting and
culling keyframes.

Figure 1: VP-SLAM pipeline.

Similar to [1], in the front-end section, points and lines are output in
parallel in each RGB frame. Then, similar to ORB-SLAM2 [8], we use the
constant velocity motion model to obtain an initial pose estimation, and then
use the points and lines to optimize it. Then, to further optimize the
rotation, we present an optimization approach that incorporates the extracted
VP and the information about the parallelism of the lines. Then, knowing the
optimized rotation, we can refine the translation by solving a linear system,
as shown in [6]. Finally, note that throughout the text, tables and vectors
are in bold capital letters, while scalar gradients are represented as simple
letters. Also, to represent the camera pose of the $i$-frame relative to the
global coordinate system, we use the notation (
$\mathbf{R_{iw}},\mathbf{t_{iw}}$ ) where $\mathbf{R_{iw}}\in\emph{SO(3)}$ is
the rotation matrix and $\mathbf{t_{iw}}\in\mathbb{R}^{3xn}$ the translation
vector.

We also use $\mathbf{\delta_{k}^{i}}$ to represent the set of directions of
VPs in the $C_{i}$ frame, and with $\mathbf{d_{k}}$ as the set of three
orthogonal dominant directions in _MW_. Additionally, we define
“$[\cdot]_{x}$“ as a 3 x 3 skew symmetric matrix, and we use ”$\sim$” to
represent an equality up to a scale.

###  III-A Real time Complexity

The computational complexity is one of the most important concerns to address
in vSLAM algorithms. To preserve the real-time properties of ORB-SLAM2 [8], we
carefully selected, and utilized, fast methods for operating line extraction
and VPs estimation.

####  III-A1 Lines

To detect line segments in the image, we use the robust LSD (Line Segment
Detector) detector [13] which is an _O(n)_ line segment detector where $n$ is
the number of pixels in the image. To describe, the line segments we use the
LBD (Line Band Descriptor) descriptor [14] which relies on local appearance of
lines and geometric constraints while preserving the real time complexity and
robustness against image artifacts.

####  III-A2 VPs Estimation

To extract the VPs we adopt the [15] method, which is based on 2 lines and
guarantees that the execution time is _O(n)_ and the extracted vanishing
points are optimal and orthogonal VPs in _MW_ . The main idea of the method is
to exploit the Gaussian sphere as the parameter space of the rotation, with
the principal point $(x_{0},y_{0})$T being the center of the sphere. Thus, two
parallel lines in 3D are projected onto the Gaussian sphere as two great
circles that intersect at a point. The direction of this point from the origin
of the sphere is considered as a candidate vanishing point direction ( $v_{1}$
). To achieve real-time performance, a polar grid is created that intersects
the image plane and spans the latitude and longitude of half of the Gaussian
sphere with a size of 90x360 and a precision of 1o. Thus, a pair of lines
intersecting at a point in the image plane contributes to the corresponding
cell of the polar grid by weight:

$\emph{score}=\emph{score}+||l_{0}||\times||l_{1}||\times\sin(2\theta)$

where $||l_{0}||$ and $||l_{1}||$ represent the length of two line segments in
pixel and $\theta$ is the angle between them and _score_ is the accumulated
score on each polar grid cell. After creating the first candidates VPs
($v_{1}$) directions and considering the orthogonal constraint, the second VP
($v_{2}$) must lie on the great circle that is perpendicular to the great
circle of ($v_{1}$) generating 360 evenly distributed candidates VPs
$({v_{i}\textsuperscript{2}}_{i=1}^{360})$ directions. Given the first
($v_{1}$) and the second ($v_{2}$) VPs directions, the third candidate VPs
directions can be obtained by the cross product of each pair of $v_{1}$ and
$({v_{i}\textsuperscript{2}}_{i=1}^{360})$.Thus, given all candidate
directions of VPs, the best estimate of VPs are calculated from the set (VPs)
with the highest score, where the score of a (VPs) hypothesis set is the sum
of the scores of three polar grid cells that belong to the three associated VP
directions.

###  III-B Absolute Rotation Optimization

Inspired by [6], [7] we propose a strategy that further optimizes the camera
rotation $C_{i}$ by exploiting the extracted VPs and the set of three mutually
orthogonal dominant directions $\mathbf{d_{k}}$ = [$\mathbf{d_{i}}^{1}$ $|$
$\mathbf{d_{i}}^{2}$ $|$ $\mathbf{d_{i}}^{3}$] of the Manhattan world (
emphMW) where $\mathbf{d_{i}}^{j}$ is a column vector. Specifically, given an
image with a cluster of 3D parallel lines in the scene, these lines must be
aligned with a dominant direction $\mathbf{d_{k}^{i}}$ in _MW_. Thus, given at
least two clusters of lines in the image, the normal vector of the great
circle on the Gaussian sphere of the corresponding line in the associated
cluster $\mathbf{s_{i}}$ must be perpendicular to the dominant direction of
the cluster. Therefore, for the $i$-th camera frame $C_{i}$ , we can formulate
a least-square problem to calculate the set of dominant directions
$\mathbf{d_{k}^{i}}$ in the current frame as below:

$\mathbf{S}^{T}\mathbf{d_{k}^{i}}=0$ (1)

where $\mathbb{S}\in\mathbb{R}^{3xn}$ and n is the number of lines in the
cluster. The normals $\mathbf{s_{i}}$ represent the columns of $\mathbf{S}$
can be obtained from $\mathbf{s_{i}}$ = $\mathbf{K^{T}}$ $\mathbf{l_{i}}$
where $\mathbf{K}$ is the camera matrix and $\mathbf{l_{i}}$ =
$(\mathbf{sp_{i}}\times\mathbf{ep_{i}}/||\mathbf{sp_{i}}\times\mathbf{ep_{i}}||)$
where $\mathbf{sp_{i}},\mathbf{ep_{i}}\in\mathbb{R}^{2}$ are the endpoints of
the line on the image plane. Finally, solving this problem with SVD we obtain
the set of mutual orthogonal dominant directions $\mathbf{d_{k}^{i}}$ in the
current $C_{i}$ frame. Also, since, the direction of the VPs reflect the
orientation of the current $C_{i}$ frame w.r.t _MW_ , and the initial set
$\mathbf{d_{k}}$ computed from the Eq. (1) (instead of computed from VPs on
the inital frame as in [7]) present the orientation of the initial frame
$C_{0}$ w.r.t _MW_ we can construct the following equation that connects the
estimated VP with the set of initial mutual orthogonal directions
$\mathbf{d_{k}}$ as follow:

$\mathbf{vp_{k}^{i}}\sim\mathbf{K}\mathbf{R_{iw}}\mathbf{d_{k}}$

Thus, as a result, we conclude that:

$\mathbf{\delta_{k}^{i}}\sim\mathbf{d_{k}^{i}}\sim\mathbf{R_{iw}}\mathbf{d_{k}}$
(2)

where, $\mathbf{\delta_{k}^{i}}$ is the Vanishing Direction (VD) of the VP
$\mathbf{vp_{k}^{i}}$.

$\mathbf{\delta_{k}^{i}}=\mathbf{K^{-1}}\mathbf{vp_{k}^{i}}$

Therefore, to optimize further the absolute rotation $\mathbf{R_{iw}}$ of the
current $C_{i}$ frame with respect to condition Eq. (2) we define the
following cost function to minimize:

$\mathbf{E(\omega_{i})}=\sum_{k=1}^{3}\mathbf{E_{k}(\omega_{i})}=\sum_{k=1}^{3}\arccos{(\mathbf{\delta_{k}^{i}}\cdot\mathbf{R_{iw}}\mathbf{d_{k}})}$
(3)

Note that if the initial frame $C_{0}$ has not at least two clusters with
enough lines, we continue with the next frame until we find a frame $C_{i}$
that meets the condition. Additionally, the $\omega_{i}$ in Lie algebra is
mapping from $\mathbf{R_{iw}}$ in Lie group. We employ the Levenberg–Marquardt
(LM) optimization to minimize the cost function $\mathbf{E(\omega_{i})}$. The
Jacobians of the cost function $\mathbf{E(\omega_{i})}$ are:

$\begin{split}\mathbf{J_{k}}&=-\frac{1}{\sqrt{1-(\delta_{k}^{i}\cdot\mathbf{d_{k}})^{2}}}\cdot\delta_{k}^{i}\frac{\partial\mathbf{d_{k}}}{\partial
w_{i}}\\\
&=\frac{1}{\sqrt{1-(\delta_{k}^{i}\cdot\mathbf{d_{k}})^{2}}}\cdot\delta_{k}^{i}\cdot([\mathbf{d_{k}}]_{x}),\end{split}$

The initial value of $\mathbf{R_{iw}}$ is obtained by optimizing both the
reprojection error of lines and reprojection error of points.

###  III-C Translation Refinement

Having optimized the global rotation $\mathbf{R_{iw}}$ from the previous step
we can use this information to construct a linear system to solve for refine
the translation. Moreover, the availability of global rotation
$\mathbf{R_{iw}}$ can lead to simplify the initial non-linear reprojection
error to a linear one. Specifically, from the:

$\mathbf{E_{p}}=\sum_{i=1}^{N}\rho(||\mathbf{x_{i}}-\pi(\mathbf{X_{i}},\mathbf{R_{iw}},\mathbf{t_{iw}},\mathbf{K})||^{2})$

with known the absolute rotation, we lead to:

$\mathbf{t_{iw}}=((\mathbf{R_{iw}}\cdot\mathbf{X_{i}})^{(3)}+t^{(3)})\cdot\begin{bmatrix}\frac{u_{i}-c_{x}}{f_{x}}\\\
\frac{v_{i}-c_{y}}{f_{y}}\end{bmatrix}-\begin{bmatrix}\mathbf{(R_{iw}}\cdot\mathbf{X_{i}})^{(1)}+t^{(1)}\\\
\mathbf{(R_{iw}}\cdot\mathbf{X_{i}})^{(2)}+t^{(2)}\end{bmatrix},$ (4)

where $[]^{(j)}$ is $j$-th row of a vector. Re-arranging (4) we lead to
minimizing the following system:

$\mathbf{t_{iw}}=\min(||\mathbf{A}\cdot\mathbf{t}-\mathbf{b}||^{2})$ (5)

where:

$\mathbf{A}=\begin{bmatrix}-1&&0&&\frac{u_{i}-c_{x}}{f_{x}}\\\
0&&-1&&\frac{v_{i}-c_{y}}{f_{y}}\end{bmatrix}$

and

$\mathbf{b}=\begin{bmatrix}(\mathbf{R_{iw}}\cdot\mathbf{X_{i}})^{(1)}-(\mathbf{R_{iw}}\cdot\mathbf{X_{i}})^{(3)}\cdot\frac{u_{i}-c_{x}}{f_{x}}\\\
(\mathbf{R_{iw}}\cdot\mathbf{X_{i}})^{(2)}-(\mathbf{R_{iw}}\cdot\mathbf{X_{i}})^{(3)}\cdot\frac{v_{i}-c_{y}}{f_{y}}\end{bmatrix}$

The system in Eq. (5) is a least square problem that can be solved under a
RANSAC framework in order to find the maximum inlier set from all pairs of
observations of the current ${C_{i}}$ frame to solve the following normal
equation:

$\mathbf{A^{T}}\cdot\mathbf{A}\cdot{t}=\mathbf{A}\cdot\mathbf{b}$

##  IV Experiment Evaluation

To validate the performance of the proposed SLAM system in terms of
computation time and accuracy, we compare our SLAM system with the state-of-
the art ORB-SLAM2 [8] in real-world scenarios using the TUM RGB-Benchmark
[10]. All experiments were carried out with an AMD Ryzen 5 2600 6-core CPU. In
this paper, we aim at demonstrating the effectiveness of exploiting VPs and
dominant directions in a _MW_ , for pose optimization. For, a fair comparison,
we disable the loop closing function from both systems. This is due to the
fact that when the loop closing module is enabled, the two systems will
converge to the same trajectory and have the same absolute pose error, so we
will not see the results of our method.

###  IV-A Localization Accuracy in the TUM RGB-D Benchmark

We test our method on TUM-RGB-D dataset [16] which consists of several real-
world camera sequences which contains a variety of scenes, like cluttered
areas and scenes containing varying degrees of structure and texture, with
ground truth trajectories and images with (640x480) resolution, recorded at
full frame rate (30 Hz). Table.I shows the absolute translation RMSE (Root
Mean Square Error) of the compared VO methods and SLAM methods, and the
estimated trajectories of different sequences in TUM-RGBD benchmark [16].

TABLE I: TUM-RGB-D results RMSE ATE (m) TUM-RGB-D Sequence | Ours | ORB-SLAM2 | PL-SLAM
---|---|---|---
fr1/desk | 0.032 | $\mathbf{0.020}$ | 0.032
fr1/plant | $\mathbf{0.020}$ | $\mathbf{0.020}$ | 0.097
fr2/desk | 0.012 | $\mathbf{0.009}$ | 0.012
fr3/large-cabinet | $\mathbf{0.100}$ | 1.420 | $\mathbf{0.100}$
fr3/long-office-house-validation | 0.019 | $\mathbf{0.010}$ | 0.019
fr3/structure-texture-far | 0.010 | 0.012 | $\mathbf{0.008}$
fr3/structure-texture-far-near | $\mathbf{0.014}$ | 0.015 | 0.014
fr3/nostr-text-near-withLoop | 0.177 | $\mathbf{0.015}$ | 0.177
fr3/nostr-text-far | $\mathbf{0.029}$ | 0.048 | 0.029

All trajectories were aligned with 7DoF with the ground truth before computing
the ATE (Absolute Translation Error) error with the script provided by the
benchmark [16]. As shown in Table.I even with disable the Loop Closure module,
the accuracy of the proposed VP-SLAM, is near the state of the art. This is
due to our concern of keeping the performance of VP-SLAM as real time as
possible and concurrently utilizing additional information from the geometry
of the environment. Moreover, we propose two simple programmable modules for
optimization and refinement, compared to ORB-SLAM2 that uses non-linear pose
graph optimizations at the back end. Additional, disabling the Loop Closure
module the accumulating error of rotation can be alleviated to some extent,
but as shown the trajectory optimized by our method is close to the ORB-SLAM2
[8]. Additionally, in Fig. 2, (left) we see the trajectory that produces our
VP-SLAM, in contrast with the ORB-SLAM2 Fig. 2 (right). Finally, in Fig. 3 we
present the line clusters that computed after extracting the VPs.

Figure 2: Estimated absolute trajectories and GT on TUM fr2/desk dataset for
VP-SLAM (left) and ORB-SLAM2 (right) respectively. TABLE II: Time Execution of
Tracking and Mapping Mean execution time in (ms)
---
Thread |  Operation |  Ours |  ORB-SLAM2 |  PL-SLAM
|  Feature Extraction |  96,94 |  15,79 |  96,64
|  VPs Extraction |  46,56 |  —- |  —-
Tracking |  Initial Pose Estimation |  71,1 |  49,06 |  71,1
|  Track Local Map |  29,93 |  26,61 |  29,93
|  Total |  244,53 |  91,46 |  197,97
|  Keyframe Insertion |  47,90 |  45,82 |  47,90
|  Map Feature Culling |  2,48 |  1,89 |  2,48
Local Mapping |  Map feature Creation |  135,93 |  21,93 |  135,93
|  Keyframe Culling |  72,52 |  70,07 |  72,52
|  Local BA |  1499,37 |  1434,57 |  1499,37
|  Total |  1758,2 |  $\mathbf{1574,28}$ |  1711,64

Figure 3: The three line clusters are shown in red, green and blue after VPs
extraction with ORB-features in green rectangles.

###  IV-B Time Complexity

We also evaluate the runtime complexity of our VP-SLAM against the state of
the art ORB-SLAM2 [8]. We summarize In Table.II the time required for each
subtask in “Tracking“ and “Local Mapping“ in our VP-SLAM and ORB-SLAM2 [8].
Note that the most time-consuming task is the local BA and the map feature
creation in both systems. In conclusion, despite the line feature detection
and VPs extraction combined with the optimization strategy increase the cost
time, the whole system can still perform in real time.

##  V Conclusion

In this paper, we propose a real time monocular visual SLAM system that
exploits the structure of man-made environments to further optimize the pose.
More specifically, it is especially suitable for environments with more
geometric structures, because it can detect VPs and lines from a single image.
We propose two methods, one of which is leverage the information of VPs and
the parallelism constraint of lines to construct an optimization strategy for
rotation, while the second one uses the information from the former
optimization to build a distinct refinement strategy for translation. Finally,
the experiments on benchmark datasets with real world scenes show that the
accuracy of the proposed system is close to the state of the art ORB-SLAM2
[8]. Additionally, the performance maintains real-time and indicate that drift
could be further reduced.

In the future, we intend to strengthen our system with a stronger optimization
technique and test it in more indoor environments.

## References

  * [1] A. Pumarola, A. Vakhitov, A. Agudo, A. Sanfeliu and F. Moreno-Noguer, ”PL-SLAM: Real-time monocular visual SLAM with points and lines,” 2017 IEEE International Conference on Robotics and Automation (ICRA), 2017, pp. 4503-4508, doi: 10.1109/ICRA.2017.7989522.
  * [2] Gomez-Ojeda, Ruben et al. “PL-SLAM: A Stereo SLAM System Through the Combination of Points and Line Segments.” IEEE Transactions on Robotics 35 (2019): 734-746.
  * [3] Hosseinzadeh, Mehdi et al. “Structure Aware SLAM Using Quadrics and Planes.” ACCV (2018).
  * [4] H. Li, J. Yao, J. -C. Bazin, X. Lu, Y. Xing and K. Liu, ”A Monocular SLAM System Leveraging Structural Regularity in Manhattan World,” 2018 IEEE International Conference on Robotics and Automation (ICRA), 2018, pp. 2518-2525, doi: 10.1109/ICRA.2018.8463165.
  * [5] J. Straub, N. Bhandari, J. J. Leonard and J. W. Fisher, ”Real-time manhattan world rotation estimation in 3D,” 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2015, pp. 1913-1920, doi: 10.1109/IROS.2015.7353628.
  * [6] J. Liu and Z. Meng, ”Visual SLAM With Drift-Free Rotation Estimation in Manhattan World,” in IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 6512-6519, Oct. 2020, doi: 10.1109/LRA.2020.3014648.
  * [7] H. Li, J. Yao, J. -C. Bazin, X. Lu, Y. Xing and K. Liu, ”A Monocular SLAM System Leveraging Structural Regularity in Manhattan World,” 2018 IEEE International Conference on Robotics and Automation (ICRA), 2018, pp. 2518-2525, doi: 10.1109/ICRA.2018.8463165.
  * [8] R. Mur-Artal, J. M. M. Montiel and J. D. Tardós, ”ORB-SLAM: A Versatile and Accurate Monocular SLAM System,” in IEEE Transactions on Robotics, vol. 31, no. 5, pp. 1147-1163, Oct. 2015, doi: 10.1109/TRO.2015.2463671.
  * [9] G. Klein and D. Murray, “Parallel tracking and mapping for small AR workspaces,” in 2007 6th IEEE and ACM international symposium on mixed and augmented reality. IEEE, 2007, pp. 225–234.
  * [10] Zhou, Yi & Kneip, Laurent & Rodríguez, Cristian & li, Hongdong. (2017). Divide and Conquer: Efficient Density-Based Tracking of 3D Sensors in Manhattan Worlds. 3-19. 10.1007/978-3-319-54193-81.
  * [11] P. Kim, B. Coltin, and H. J. Kim, “Visual Odometry with Drift-free Rotation Estimation Using Indoor Scene Regularities,” in British Machine Vision Conference (BMVC), 2017
  * [12] P. Kim, B. Coltin, and H. J. Kim, “Low-drift visual odometry in structured environments by decoupling rotational and translational motion,” in 2018 IEEE international conference on Robotics and automation (ICRA). IEEE, 2018, pp. 7247–7253.
  * [13] R. Grompone von Gioi, J. Jakubowicz, J. -M. Morel and G. Randall, ”LSD: A Fast Line Segment Detector with a False Detection Control,” in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 32, no. 4, pp. 722-732, April 2010, doi: 10.1109/TPAMI.2008.300.
  * [14] L. Zhang and R. Koch. An efficient and robust line segment matching approach based on LBD descriptor and pairwise geometric consistency. JVCIR, 24(7):794–805, 2013.
  * [15] X. Lu, J. Yaoy, H. Li, Y. Liu and X. Zhang, ”2-Line Exhaustive Searching for Real-Time Vanishing Point Estimation in Manhattan World,” 2017 IEEE Winter Conference on Applications of Computer Vision (WACV), 2017, pp. 345-353, doi: 10.1109/WACV.2017.45.
  * [16] J. Sturm, N. Engelhard, F. Endres, W. Burgard and D. Cremers, ”A benchmark for the evaluation of RGB-D SLAM systems,” 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2012, pp. 573-580, doi: 10.1109/IROS.2012.6385773.