# Compressible turbulent convection in highly stratified adiabatic background

John Panickacheril John  1 and Jörg Schumacher1,2  john.panickacheril-john@tu-
ilmenau.de 1 Institut für Thermo- und Fluiddynamik, Technische Universität
Ilmenau, Postfach 100565, D-98684 Ilmenau, Germany 2 Tandon School of
Engineering, New York University, New York City, New York 11201, USA

###### Abstract

Buoyancy-driven turbulent convection leads to a fully compressible flow with a
strong top-down asymmetry of first- and second-order statistics when the
adiabatic equilibrium profiles of temperature, density and pressure decay very
strongly across the convection layer. The growth of this asymmetry and the
formation of an increasingly thicker stabilized sublayer with a negative mean
convective heat flux at the top of the convection zone is reported here by a
series highly resolved three-dimensional direct numerical simulations beyond
the Oberbeck-Boussinesq and anelastic limits for dimensionless dissipation
numbers $0.1\leq D\leq 0.8$ at fixed Rayleigh number $Ra=10^{6}$ and
superadiabaticity. The highly stratified compressible convection regime
appears for $D>D_{\rm crit}\approx 0.65$, when density fluctuations collapse
to those of pressure; it is characterized by an up to nearly 50% reduced
global turbulent heat transfer and a sparse network of focused thin thermal
plumes falling through the top sublayer deep into the bulk.

###### keywords:

Compressible turbulent convection

##  1 Introduction

Solar convection close to the surface proceeds in the presence of an extremely
stratified adiabatic background equilibrium which is characterized by scale
heights, the distances across which temperature $T$, density $\rho$, or
pressure $p$ drop by an order of magnitude, as small as $\ell\sim 10^{3}$ km
amounting to only 5‰ of the total depth of the solar convection layer, $H\sim
2\times 10^{5}$ km (Christensen-Dalsgaard, 2002; Nordlund et al., 2009;
Schumacher & Sreenivasan, 2020). It is one prominent example for a fully
compressible convection flow which is highly asymmetric when up- and down-
welling thermal plumes are compared. Driven by the strong radiative cooling at
the top surface, cold plasma sinks into the interior with nearly the speed of
sound at the narrow edges of the granules, the convection cells in the solar
case. Plasma rises moderately to the surface across the granule interior with
a typical diameter of about $10^{3}$ km (Magic et al., 2013). It is yet open,
how deep the narrow plumes reach into the highly stratified interior (Cossette
& Rast, 2016; Brandenburg, 2016; Anders et al., 2019) and how much they get
focused by compressibility effects to overcome turbulent dispersion. Answers
to these questions would also provide insights on the vertical extent of
supergranules, the larger-scale convection cells in the background with an
extension of about 30 granule diameters (Rincon & Rieutord, 2018; Hanson et
al., 2020) and on the related solar conundrum of anomalously low velocity
fluctuation amplitudes from helioseismology in comparison to models (Hanasoge
et al., 2012). All this sets the physical motivation for the present study.

In this work, we want to study fully compressible convection in the presence
of highly stratified adiabatic equilibrium profiles for $T$, $\rho$, and $p$.
We thus analyse a series of direct numerical simulations (DNS) towards the
high stratification limit of $D\to 1$ for the dissipation number $D$ (defined
below). Our goal is to analyse the genuine compressibility effects isolated
from complex multi-physics, as in the solar example above, and from
temperature dependencies of dynamic viscosity and thermal conductivity. Our
three-dimensional DNS beyond the Oberbeck-Boussinesq (Chillà & Schumacher,
2012) and anelastic regimes (Ogura & Phillips, 1962; Verhoeven et al., 2015;
Jones et al., 2022) lead to strongly differing boundary layer dynamics at the
top and bottom, such as a sublayer formation at the top with mean negative
convective heat flux which causes a reduction of the global turbulent heat
transfer by almost 50% for $D=0.8$. In this layer, narrow focused thermal
plumes form at the top and fall deeply into the highly stratified bulk of the
layer. Our study thus demonstrates that some typical phenomena in natural
convection flows can be understood in significantly simpler compressible flow
configurations even though being far away from the natural flow case in terms
of Rayleigh and Prandtl numbers.

##  2 Simulation model and parameters of compressible convection

The equations of motion for compressible convection are given by

$\displaystyle\partial_{t}\rho+\partial_{i}(\rho u_{i})$ $\displaystyle=0\,,$
(1) $\displaystyle\partial_{t}(\rho u_{i})+\partial_{j}(\rho u_{i}u_{j})$
$\displaystyle=-\partial_{i}p+\partial_{j}\sigma_{ij}+\rho g\delta_{i,3}\,,$
(2) $\displaystyle\partial_{t}(\rho e)+\partial_{j}(\rho eu_{j})$
$\displaystyle=-p\partial_{i}u_{i}+\partial_{i}(k\partial_{i}T)+\sigma_{ij}S_{ij}\,,$
(3) $\displaystyle p$ $\displaystyle=\rho RT\quad\textrm{ where }\quad
R=C_{p}-C_{v}.$ (4)

These equations correspond to mass, momentum and energy conservation laws
along with the ideal gas equation of state. Here, $\rho$, $\rho u_{i}$, $p$,
$\rho e$, $T$ are the mass density, momentum density components, pressure,
internal energy density, and temperature, respectively. The viscous stress
tensor is
${\bm{\sigma}}=2\mu{\bm{S}}+2\mu{\bm{I}}({\bm{\nabla}}\cdot{\bm{u}})/3$ with
the Kronecker tensor ${\bm{I}}$ and the rate of strain tensor
${\bm{S}}=({\bm{\nabla}}{\bm{u}}+{\bm{\nabla}}{\bm{u}}^{T})/2$. The dynamic
viscosity $\mu$ is assumed to be constant in these simulations. The thermal
conductivity $k$ is related to the viscosity through the Prandtl number,
$k=\mu C_{p}/Pr$. In our DNS, $Pr=0.7$. $C_{p}$ and $C_{v}$ correspond to
specific heat at constant pressure and volume, respectively. Their ratio,
$\gamma=C_{p}/C_{v}=1.4$ for a diatomic gas.

A uniform grid is used in $x$– and $y$–directions along with periodic boundary
conditions. In wall-normal $z$–direction, a non-uniform grid with a point
clustering near the walls which follows a hyperbolic tangent stretching
function. Spatial derivatives are calculated by a 6th-order compact scheme
(Lele, 1992) for all points except near the walls (Baranwal et al., 2022);
there 4th- and 3rd–order compact schemes are used at the last two grid points
near the wall. No-slip, isothermal boundary conditions are applied at the top
and bottom. The boundary condition for $p$ is evaluated using the
$z$-component of the momentum equation at $z=0,H$, $\partial p/\partial
z=\partial\sigma_{iz}/\partial x_{i}+\rho g$. The fields are advanced in time
by a low storage 3rd-order Runge-Kutta with a Courant number of ${\rm
CFL}=0.5$.

Case | $D$ | $T^{B}_{a}/T^{T}_{a}$ | $\rho^{B}_{a}/\rho^{T}_{a}$ | $M_{t}^{\rm max}$ | $Nu$ | $Re$ | $\lambda_{\rho_{\rm min}}$ | $\delta_{u_{z}T}$ | $\left(\lambda^{\rm max}_{\rho^{\prime}T^{\prime}},\delta_{\rho^{\prime}T^{\prime}}\right)$ | Style
---|---|---|---|---|---|---|---|---|---|---
1 | $0.10$ | $1.10$ | $1.3$ | $0.14$ | $7.94\pm 0.05$ | $424\pm 2$ | 0.002 | 0.001 | $\left(0.001,0\right)$ |
2 | $0.34$ | $1.52$ | $2.8$ | $0.32$ | $7.64\pm 0.01$ | $414\pm 1$ | 0.012 | 0.005 | $\left(0.006,0.012\right)$ |
3 | $0.50$ | $2.00$ | $5.7$ | $0.41$ | $6.93\pm 0.02$ | $407\pm 2$ | 0.024 | 0.010 | $\left(0.012,0.024\right)$ |
4 | $0.60$ | $2.50$ | $9.9$ | $0.49$ | $6.24\pm 0.01$ | $370\pm 8$ | 0.037 | 0.016 | $\left(0.018,0.037\right)$ |
5 | $0.65$ | $2.86$ | $13.8$ | $0.55$ | $5.90\pm 0.01$ | $364\pm 2$ | 0.042 | 0.019 | $\left(0.021,0.042\right)$ |
6 | $0.70$ | $3.33$ | $20.3$ | $0.54$ | $5.80\pm 0.04$ | $425\pm 2$ | 0.051 | 0.023 | $\left(0.023,0.051\right)$ |
7 | $0.75$ | $4.00$ | $32.0$ | $0.60$ | $5.14\pm 0.01$ | $399\pm 0$ | 0.068 | 0.033 | $\left(0.028,0.068\right)$ |
8 | $0.80$ | $5.00$ | $55.9$ | $0.61$ | $4.46\pm 0.03$ | $356\pm 1$ | 0.098 | 0.051 | $\left(0.035,0.098\right)$ |
Table 1:  List of the direct numerical simulations. All cases have $Pr=0.7$,
$Ra\approx 10^{6}$, $\epsilon=0.1$, $\gamma=C_{p}/C_{v}=1.4$, and an aspect
ratio, $\Gamma=L/H=4$ resolved by $N_{x}:N_{y}:N_{z}=512:512:256$ grid points.
We list dissipation number $D$, ratios of adiabatic temperatures and densities
from bottom to top, the global maximum turbulent Mach number $M_{t}^{\rm
max}$, the Nusselt and Reynolds numbers, $Nu$ and $Re$. The statistical
error/difference for the Nusselt number and Reynolds number are less than
$2\%$ for all cases. Furthermore, we display boundary layer scales for
density, convective heat flux, and density-temperature correlations. In all
figures that follow, coloured lines are given as in this table.

Incompressible Rayleigh-Bénard convection is determined by the Rayleigh and
Prandtl numbers, $Ra$ and $Pr$. In compressible convection, additional
dimensionless parameters need to be introduced. The first one is the
dissipation number $D$ given by

$D=\frac{gH}{C_{p}T^{B}}\,.$ (5)

The parameter $D$ can be considered as a rescaled dry adiabatic lapse rate as
it determines differences in the thermodynamic properties across the layer
depth in a purely isentropic process. The adiabatic temperature ratio between
top (T) and bottom (B) is $T^{T}_{a}/T^{B}_{a}=(1-D)$, the corresponding
density ratio follows to
$\rho^{T}_{a}/\rho^{B}_{a}=\left[T^{T}_{a}/T^{B}_{a}\right]^{1/\left(\gamma-1\right)}$.
The subscript ”a” corresponds to the adiabatic equilibrium without convection.
In order to initiate convective motion, the actual temperature gradient across
the layer must be greater than due to a purely adiabatic process (Jeffreys,
1930), i.e., $T^{B}-T^{T}>T_{a}^{B}-T_{a}^{T}$. In our study, the bottom plate
is taken as the reference, thus $T^{B}=T_{a}^{B}$ is a constant for all cases;
the instability criterion reduces to $T^{T}<T_{a}^{T}$.

The second compressible parameter is the superadiabaticity $\epsilon$, which
is given by

$\epsilon=\frac{T_{a}^{T}-T^{T}}{T_{B}}\,.$ (6)

Superadiabaticity is the excess relative temperature gradient with respect to
the adiabatic equilibrium state, or in other words, a measure of the driving
of convection. If one takes the characteristic free-fall velocity,
$U_{f}=\sqrt{\epsilon gH}$, then a free-fall Mach number,

$M_{f}(z)=\frac{U_{f}}{\gamma R\langle T(z)\rangle_{A,t}}=\sqrt{\frac{\epsilon
D}{\gamma-1}}\sqrt{\frac{T_{B}}{\langle T(z)\rangle_{A,t}}}\,,$ (7)

follows. Thus strength of compressibility is thus a function of both,
$\epsilon$ and $D$. From (7), it is expected that compressibility at the top
is higher than at the bottom. For $\epsilon\rightarrow 0$, we get the
anelastic approximation where the acoustic waves are filtered out from the
equations. Rayleigh Bérnard convection in the Oberbeck-Boussinesq (OB) regime
follows for $\epsilon,D\rightarrow 0$. On the basis of the new parameters, we
will use $Ra=\epsilon gH^{3}/\left(\nu\kappa\right)$ and Prandtl number,
$Pr=\nu/\kappa$, where $\kappa$ is the temperature diffusivity and
$\nu=\mu/\rho$ the kinematic viscosity. In this work, we focus on a constant
and low but finite superadiabaticity, $\epsilon=0.1$ for a range of $D$ from
$0.1$ to $0.8$. All other governing parameters such as the $Ra$, $Pr$ and
aspect ratio, $\Gamma$ are kept constant, see Table. 1. In all figures that
follow, unless specified otherwise, we use the plotting style given in the
last column.

Figure 1:  Snapshots of the normalized superadiabatic temperature fields.
Isosurfaces of normalized $T_{sa}$ for case 1 with $D=0.1$ (left) and case 8
with $D=0.8$ (right) are shown. The right panel shows a ridge network with
thin thermal plumes emanating from the top deeply into the convection layer.
The superadiabatic temperature is normalized in both panels to $[-1,0]$.
(Left) Contour levels at the bottom and top correspond to -0.27 and -0.73,
respectively. The bottom and top contours to the right stand for -0.09 and
-0.61, respectively, see the color bars.

\begin{picture}(0.0,0.0)\end{picture}

##  3 Growing asymmetry with increased stratification

The adiabatic temperature profile, $T_{a}(z)=T^{B}_{a}(1-Dz)$, is the
equilibrium profile. Thus the part of temperature corresponding to convection
is the superadiabatic temperature defined as

$T_{sa}({\bm{x}},t)=T({\bm{x}},t)-T_{a}(z)\,.$ (8)

In figure 1, we compare the contour visualization of $T_{sa}$ for the two
extreme cases in our series with $D=0.1$ (low) and $D=0.8$ (high). The left
panel corresponds to case 1, closest to the OB limit. Not surprisingly, the
plumes from the top and bottom are more or less symmetric in shape and
frequency. However, this top-bottom-symmetry disappears completely for case 8
in the right panel. One observes broader thick plumes from the bottom compared
to a network of slender thinner plumes originating from the top. Recently,
John & Schumacher (2023) showed that thinner plumes from the top boundary are
a common characteristic of compressible convection regardless of the relative
boundary layer thickness, but here the differences are very pronounced.

Similar asymmetry is observed for the mean profile of the superadiabatic
temperature, $\langle T_{sa}(z)\rangle_{A,t}$, where the notation
$\langle\cdot\rangle_{A,t}$ corresponds to an average over the horizontal
directions and time. In figure 2(a), we plot the normalized superadiabatic
temperature, $\langle T_{sa}(z)\rangle_{A,t}/\Delta\langle
T_{sa}\rangle_{A,t}$ as a function of depth $z$ where $\Delta\langle
T_{sa}\rangle_{A,t}=\langle T_{sa}(z=0)\rangle_{A,t}-\langle
T_{sa}(z=H)\rangle_{A,t}$. Clearly one observes that the asymmetry increases
considerably with growing dissipation number. As $D$ increases, the bulk
temperature is increasingly closer to the bottom temperature. It is also
evident from the figure that the thermal boundary layer thickness at the top
is large in comparison to that at the bottom as $D$ increases.

Figure 2:  Asymmetry of the (a) normalized superadiabatic temperature and (b)
the turbulent Mach number $M_{t}$. The inset in (a) corresponds to mean
thermal boundary layer thickness $\lambda^{B,T}(D)$. Symbols $+$ and $\rm o$
are for the top and bottom boundary layer thickness, respectively. The inset
in (b) shows $M_{f}(D)$ from (7). The red, blue and green lines are at
$z/H=0,1$, and 0.5. The black squares are the maximum Mach numbers as given
Table 1, which are observed in the simulations over the whole space and time.
The legend in (b) corresponds to the inset figure.

\begin{picture}(0.0,0.0)\end{picture}

We can quantify this thickness as the mean of a local thermal boundary layer
thickness (Scheel & Schumacher, 2014) which is given by

$\tilde{\lambda}^{B,T}(x,y)=\frac{1}{H}\left|\langle
T_{sa}\left(0,H\right)-T_{sa}\left(H/2\right)\rangle_{A,t}\right|\left|\frac{\partial\langle
T_{sa}\rangle_{A,t}}{\partial z}\right|^{-1}_{z=0,H}.$ (9)

We plot $\lambda^{B,T}=\langle\tilde{\lambda}^{B,T}\rangle_{A,t}$ as a
function of $D$ in the inset of figure 2(a). The difference between the
boundary layer thickness at the bottom and top plate increases with $D$ in
agreement with the differences observed in the distributions of the
instantaneous, local thermal boundary layer thicknesses (John & Schumacher,
2023). The top boundary layer thickness increases drastically with $D$ whereas
bottom boundary layer thickness is almost invariant. This along with increased
asymmetry in the bulk towards the bottom plate, as seen in the normalized
$T_{sa}$-profiles, suggests that at high $D$, the mixing is not as efficient
as in OB and compressibility effects are mainly localized near the top
boundary.

An asymmetry similar to that of $\langle T_{sa}(z)\rangle_{A,t}$ is observed
for the turbulent Mach number too with significant localized compressibility
seen near the top boundary. The turbulent Mach number is defined as

$M_{t}(z)=\frac{u_{\rm rms}(z)}{\sqrt{\gamma R\langle
T(z)\rangle_{A,t}}}\quad\mbox{with}\quad u_{\rm
rms}(z)=\sqrt{\sum_{i=1}^{3}\left[\langle u_{i}^{2}(z)\rangle_{A,t}-\langle
u_{i}(z)\rangle^{2}_{A,t}\right]}\,.$ (10)

The turbulent Mach number at the top increases monotonically with $D$ whereas
the behaviour at the bottom boundary is non-monotonic with $D$. At the bottom,
$M_{t}$ increases up to $D=0.5$ and then decreases with $D$ as a result
increasing the asymmetry between the boundaries further as seen in figure
2(b).

In the inset of figure 2(b), we compare the resulting maximum turbulent Mach
number in the flow with the estimate from (7). The black squares are the
maximum local turbulent Mach number for each $D$ encompassing all time steps
and spatial locations. The red and blue lines correspond to $M_{f}$-estimates
using the bottom and top temperatures, $T^{B}$ and $T^{T}$, respectively. One
finds that the resulting maxima of $M_{t}$ are in between these estimates.
This is due to no-slip boundary conditions at the wall; thus the velocity
peaks only at the outer layer of the boundary layer and in the bulk. Since the
(mean) bulk temperature is closer to that at the bottom plate, the values of
$M_{t}$ are closer to $M_{f}\left(T^{B}\right)$. A really close agreement for
the simulations is seen, when the Mach number is estimated using the bulk
temperature, green line in the inset figure. This suggests that a simple
estimate such as in (7) can give us a good idea about the level of
compressibility in the flow. We would like to point out that in compressible
turbulence $M_{t}$ is not sufficient to characterize the system; an additional
parameter $\delta$, the ratio of dilatational to solenoidal root mean square
velocities, needs to be included as discussed in Donzis & John (2020); John et
al. (2019, 2020, 2021).

##  4 Fluctuations of state variables and regime transition

A major consequence of compressibility are the fluctuations of $p$, $\rho$,
and $T$. In figure 3, we compare these fluctuations as a function of depth for
the two extreme cases of our DNS series, $D=0.1$ and $0.8$. John & Schumacher
(2023) showed that $D\leq 1-\epsilon$, which would result in a maximally
possible theoretical value of $D=0.9$ in our DNS series. The relative
fluctuations are defined as $X_{\rm rms}/\langle X\rangle_{A,t}$ with $X_{\rm
rms}=(\langle X^{2}\rangle_{A,t}-\langle X\rangle_{A,t}^{2})^{1/2}$. Panel (a)
at $D=0.1$ shows that pressure fluctuations are negligible. The magnitude of
the density and temperature fluctuations are of the same order of magnitude
and the profiles of both quantities are practically identical. These findings
are consistent with the OB limit at $\epsilon,D\rightarrow 0$.

In figure 3(b), the relative fluctuations of $p$, $\rho$, and $T$ at $D=0.8$
are shown. As result of high compressibility at $D=0.8$ at the top, see figure
2(b) for $M_{t}(z)$, the thermodynamic fluctuations are considerably higher
compared to the bottom one. Close to the top boundary, the pressure and
density fluctuations dominate the temperature fluctuations. This is due to the
constant temperature imposed at the boundaries which results in a isothermal
condition. For an assumed ideal gas, $p\propto\rho$. Indeed this is evident
for $D=0.8$ at the top boundary as the relative fluctuations of density and
pressure perfectly match at the top boundary. This would be also true at
bottom boundary and even for $D=0.1$, but it is not evident due to low
$M_{t}$. Compressibility effects imply in the following discussion that we
face a large magnitude of pressure fluctuations.

As we move away from the wall, fluctuations of both, $p$ and $\rho$, decrease,
but are not perfectly correlated anymore due to the growing temperature
fluctuations. After a critical distance from the wall, pressure fluctuations
continue to decrease; the density fluctuations reach a minimum at
$z=z_{\ast}$. This local minimum of $\rho_{\rm rms}/\langle\rho\rangle_{A,t}$
exists when pressure and temperature fluctuations become comparable in
magnitude. For the highest dissipation number $D=0.8$, density fluctuations at
$z<z_{\ast}$ are driven by temperature fluctuations in contrast to
$z>z_{\ast}$. One also notes that the location of the secondary maximum of
$\rho_{\rm rms}/\langle\rho\rangle_{A,t}$ does not match the local maxima of
$T_{\rm rms}/\langle T\rangle_{A,t}$. This implies that pressure fluctuations,
although not dominant, remain significant and thus compressibility effects too
(Jones et al., 2022).

Figure 3:  Relative fluctuations of thermodynamic quantities for (a) $D=0.1$
and (b) $D=0.8$. (c) Maximum of the relative thermodynamic fluctuations versus
dissipation number, $D$. For all figures, the red, blue and green correspond
to temperature, density and pressure respectively. Inset of (c) shows four
different distances $\lambda_{X}$ as a function of $D$. For the definition of
the specific $\lambda_{X}$ we refer to the text in section 4 and the legend of
the inset.

\begin{picture}(0.0,0.0)\end{picture}

These features are absent for $D=0.1$ in panel (a). Thus there seems to exist
a critical $D$ where a transition to the compressibility-dominated convection
regime exists that is qualitatively different to the OB or anelastic limits.
In order to identify the critical $D$, we plot the maximum relative
thermodynamic fluctuations with respect to $D$ in figure 3(c). We find a
$D_{\rm crit}\approx 0.65$ for the present $(Ra,\epsilon)$. For $D<D_{\rm
crit}$, maxima of density and temperature fluctuations dominate over those of
the pressure fluctuations. As $D$ increases, the maximum pressure fluctuations
start to increase as well which is connected with increasing differences in
the maxima of density and temperature fluctuations. For $D\geq D_{\rm crit}$,
the maxima of pressure and density fluctuations collapse and are larger than
the one for the temperature fluctuations.

Clearly, there are multiple physical processes and thus scales at play near
the top boundary at high $D$. This seems to be similar to the turbulent
boundary layer with its inner and outer scales (Smits et al., 2011). The
behaviour of density can be analyzed to demarcate various regimes inside the
boundary layers in compressible convection. We thus introduce the following
boundary layer scales: (1) $\lambda_{T_{\rm max}}$ as the distance between the
top wall and the location of the local maximum of the relative temperature
fluctuations. (2) $\lambda_{\rho_{\rm max}}$ as the distance between the top
wall and the maximum of the relative density fluctuations. Note that in OB and
even anelastic convection, $\lambda_{T_{\rm max}}=\lambda_{\rho_{\rm max}}$
(Jones et al., 2022). Furthermore, (3) $\lambda_{\rho_{\rm min}}$ as the
distance between the top wall and the location of the local minimum of the
relative density fluctuations. This distance corresponds to the location where
we start to observe significant compressibility effects. Finally, (4)
$\lambda_{\rho_{{\rm max},T}}$ as the distance between the top wall and the
local maximum of the relative density fluctuations corresponding to the
maximum of local temperature fluctuations.

Next, we display all 4 scales $\lambda_{X}$ versus $D$ in the inset of figure
3(c). $\lambda_{T_{\rm max}}$ increases with $D$ which is consistent with the
growth of the thermal boundary layer thickness at the top, $\lambda^{T}$, see
figure 2(a). We see a sudden transition of $\lambda_{\rho_{\rm max}}$ at
$D_{crit}=0.65$ indicating that the maximum for $D>D_{\rm crit}$ is at $z=H$.
Note also that for $D<0.65$, $\lambda_{\rho_{\rm max}}=\lambda_{\rho_{{\rm
max},T}}$. Both scales increase with $D$, but start to deviate from
$\lambda^{T}$ at $D<D_{\rm crit}=0.65$. This deviation starts for $D\geq 0.5$;
it suggests that pressure fluctuations are important already, but not as
dominant as discussed before. It might thus indicate that anelastic
approximation is strictly speaking not applicable even for $D\geq 0.5$.
Finally, $\lambda_{\rho_{\rm min}}$ starts from nearly zero at $D=0.1$ with
negligible compressibility. With increasing dissipation number
$\lambda_{\rho_{\rm min}}$ increases. This is a new top boundary layer regime
induced by the increasing stratification. It is in line with an increase of
$\lambda_{T_{\rm max}}$ and $\lambda^{T}$, see again the inset of figure 2(a).
We discuss the implications of this sublayer for the turbulent heat transfer
in the next section.

##  5 Stratification impact on turbulent heat transfer

A systematic increase of the level of stratification, controlled by $D$, at
relatively small superadibaticity of $\epsilon=0.1$ has a significant impact
on the normalized stress profiles as seen in figure 4. We show the variation
of the normalized convective heat flux profiles $J_{c}(z)$ in (a) and the
normalized density-temperature correlation profiles $J_{\rho}(z)$ in (b) for
all 8 runs. They are qualitatively similar, but sign-flipped behaviour is
observed for both second-order correlations in (a,b). An increasing sublayer
with an average “anti-convection”-behaviour, i.e., negative convective heat
flux and positive density-temperature correlation can be detected. For
$z\lesssim 0.9$ both mean profiles correspond to unstable convection.

Next, we relate the depth of these stabilized regions to $\lambda_{\rho_{\rm
min}}$ and thus to the regime transition at $D=D_{\rm crit}$. We therefore
define three further distances or scales: (1) $\delta_{u_{z}T}$ as the width
of negative normalized heat flux with respect to the top boundary. (2)
$\lambda^{\rm max}_{\rho^{\prime}T^{\prime}}$ as the distance from the top
boundary to the location of the positive maximum of the normalized density-
temperature correlation. (3) $\delta_{\rho^{\prime}T^{\prime}}$ as the width
of the region of positive density-temperature correlations.

Figure 4:  Profiles of plane-time averaged stresses versus $D$. (a) Normalized
convective heat flux $J_{c}(z)=\langle
u_{z}^{\prime}T^{\prime}(z)\rangle_{A,t}/[u_{z,{\rm rms}}(z)\langle
T(z)\rangle_{A,t}]$. (b) Normalized density-temperature correlation
$J_{\rho}(z)=\langle\rho^{\prime}T^{\prime}(z)\rangle_{A,t}/[\langle\rho(z)\rangle_{A,t}\langle
T(z)\rangle_{A,t}]$. (c) Nusselt number $Nu$ versus dissipation number $D$,
see eq. (11). Inset of (c) displays $1/Nu$ versus $\lambda_{\rho_{min}}$. The
dashed line corresponds to a fit $A+B\lambda_{\rho_{\rm min}}$ with
coefficients $A=0.12$ and $B=1.05$. Line colouring corresponds to table 1.

\begin{picture}(0.0,0.0)\end{picture}

These length scales along with $\lambda_{\rho_{\rm min}}$ are listed in table
1 for all DNS. We observe that for all runs, except for $D=0.1$, the region of
positive correlation exactly matches with the compressibility-dominated
region, i.e., $\delta_{\rho^{\prime}T^{\prime}}=\lambda_{\rho_{\rm min}}$.
Thus, there is also a direct proportionality between regions dominated by
compressibility and anti-convection with $\delta_{{u_{z}}T}<\lambda_{\rho_{\rm
min}}$. For $D\leq 0.70$, we also find that
$\delta_{{u_{z}}T}\approx\lambda^{\rm max}_{\rho^{\prime}T^{\prime}}$, the
location of the maximum positive correlation. For $D>0.7$, we get
$\delta_{{u_{z}}T}>\lambda^{\rm max}_{\rho^{\prime}T^{\prime}}$, suggesting an
increasing influence of compressibility. Furthermore,
$\delta_{{u_{z}}T}<\lambda_{\rho_{\rm min}}$ implies that on average weak
convective motions can occur at heights, $\lambda_{\rho_{\rm min}}\gtrsim
z\gtrsim\delta_{{u_{z}}T}$, even though there are already significant
compressibility effects. We recall for this discussion that even for the
highest $D$, we still observe focused cold thermal plumes which detach from
the top thus indicating locally positive convective flux. All these processes
take places in the region where density fluctuations are tightly correlated to
pressure fluctuations as we discussed in figure 3(c). We also report that that
there is no trend of the location of maximum turbulent Mach number,
$M^{\ast}_{t}=\max_{z}M_{t}(z)$, see figure 2(b), with $D$ even though
$M_{t}^{\ast}$ increases with $D$. For all the cases considered, the distance
of this maximum to the top wall is $H-z_{\ast}\in[0.16,0.174]$, significantly
larger than $\lambda_{\rho_{\rm min}}$. This again suggests the inadequacy of
$M_{t}$ alone to characterize general compressibility conditions (Donzis &
John, 2020). Particularly for the strongly asymmetric cases at highest $D$,
the behaviour of compressible convection is completely different for same
$M_{t}$-magnitude at either sides of $z_{\ast}$.

Finally, we plot the Nusselt number defined as

$Nu=\frac{1}{2}\left[\left.-\frac{H}{\epsilon T_{B}}\frac{d\langle
T_{sa}\rangle_{A,t}}{dz}\right|_{z=0}+\left.-\frac{H}{\epsilon
T_{B}}\frac{d\langle T_{sa}\rangle_{A,t}}{dz}\right|_{z=H}\right]$ (11)

with dissipation as function of $D$ in figure 4(c). The Nusselt number is
calculated as the average of the heat flux at the top and bottom. Since the
flows are statistically steady, the Nusselt number at the top and bottom are
the same and the maximum mean Nusselt number difference between the plate is
found to be less than $3\%$ for all cases. Clearly visible is that the Nusselt
number $Nu$ monotonically decreases with $D$. At critical $D=0.65$, there
seems to be slight change in behaviour but is not significant. To demonstrate
conclusively that inefficient heat transfer is due to compressibility, we plot
the inverse of the Nusselt number versus $\lambda_{\rho_{\rm min}}$ in the
inset of (c), we observe that the trend can be reasonably approximated by a
linear fit.

##  6 Final discussion

The objective of the present study was to take the least complex, but fully
compressible turbulent convection system to study the genuine effects of
compressibility isolated from other major contributions to a non-Boussinesq
behavior, such as temperature-dependent material properties, beyond the
Oberbeck-Boussinesq and anelastic limits. We aimed at a configuration that is
operated in the limit of strongly decreasing mean profiles of the
thermodynamic state variables, the high stratification regime. It corresponds
to the dissipation number of $D\to 1$ for small superadiabaticity $\epsilon$.
One prominent example for such a configuration might be the convection
dynamics close to the solar surface. Even though the aim of our work is not to
provide a realistic modeling of particularly this complex multi-physics flow
and transport process, that involves the coupling to magnetic fields and
radiation transfer, the present simulations provide an insight into generic
physical mechanisms which lead to highly asymmetric thermal plume motions, the
fundamental local process in any convection flow that triggers fluid
turbulence and transports heat. The limit cases of our simplified
configuration leads already to a sparse network of thin plume structures which
resembles similarities to the granule network at the Sun, see again figure
1(b). These coherent structures “leak“ through the stratified top region and
thus maintain a significantly reduced global turbulent heat transport. Future
studies could aim at an increase of the aspect ratio and Rayleigh number
together with a decrease of the Prandtl number and $T$-dependent material
parameters.

Acknowledgements

J.P.J. is supported by grant no. SCHU 1410/31-1 of the Deutsche
Forschungsgemeinschaft and a Postdoctoral Fellowship of the Alexander von
Humboldt Foundation. The authors gratefully acknowledge the Gauss Centre for
Supercomputing e.V. (https:// www.gauss-centre.eu) for funding this project by
providing computing time through the John von Neumann Institute for Computing
(NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).

The authors report no conflict of interest.

## References

  * Anders et al. (2019) Anders, E. H., Lecoanet, D. & Brown, B. P. 2019 Entropy rain: Dilution and compression of thermals in stratified domains. Astrophys. J. 884 (1), 65.
  * Baranwal et al. (2022) Baranwal, A., Donzis, D. A. & Bowersox, R. D. W. 2022 Asymptotic behaviour at the wall in compressible turbulent channels. J. Fluid. Mech. 933.
  * Brandenburg (2016) Brandenburg, A. 2016 Stellar mixing length theory with entropy rain. Astrophys. J. 832 (1), 6.
  * Chillà & Schumacher (2012) Chillà, F. & Schumacher, J. 2012 New perspectives in turbulent Rayleigh-Bénard convection. Euro. Phys. J. E 35 (7), 1–25.
  * Christensen-Dalsgaard (2002) Christensen-Dalsgaard, J. 2002 Helioseismology. Rev. Mod. Phys. 74 (4), 1073.
  * Cossette & Rast (2016) Cossette, J.-F. & Rast, M. P. 2016 Supergranulation as the largest buoyantly driven convective scale of the Sun. Astrophys. J. Lett. 829 (1), L17.
  * Donzis & John (2020) Donzis, D. A. & John, J. P. 2020 Universality and scaling in homogeneous compressible turbulence. Phys. Rev. Fluids 5, 084609.
  * Hanasoge et al. (2012) Hanasoge, S. M., Duvall Jr, T. L. & Sreenivasan, K. R 2012 Anomalously weak solar convection. Proc. Nat. Acad. Sci. 109 (30), 11928–11932.
  * Hanson et al. (2020) Hanson, C. S., Duvall, T. L., Birch, A. C., Gizon, L. & Sreenivasan, K. R 2020 Solar east-west flow correlations that persist for months at low latitudes are dominated by active region inflows. Astro. & Astrophys. 644, A103.
  * Jeffreys (1930) Jeffreys, H. 1930 The instability of a compressible fluid heated below. In Math. Proc. Camb. Phil. Soc., , vol. 26, pp. 170–172. Cambridge University Press.
  * John et al. (2019) John, J. P., Donzis, D. A & Sreenivasan, K. R 2019 Solenoidal scaling laws for compressible mixing. Phys. Rev. Lett. 123 (22), 224501.
  * John et al. (2020) John, J. P., Donzis, D. A. & Sreenivasan, K. R 2020 Compressibility effects on the scalar dissipation rate. Combust. Sci. Technol. 192 (7), 1320–1333.
  * John et al. (2021) John, J. P., Donzis, D. A. & Sreenivasan, K. R. 2021 Does dissipative anomaly hold for compressible turbulence? J. Fluid. Mech. 920, A20.
  * John & Schumacher (2023) John, J. P. & Schumacher, J. 2023 Strongly superadiabatic and stratified limits of compressible convection. arXiv preprint arXiv:2302.03621 .
  * Jones et al. (2022) Jones, C. A., Mizerski, K. A & Kessar, M. 2022 Fully developed anelastic convection with no-slip boundaries. J. Fluid Mech. 930, A13.
  * Lele (1992) Lele, S. K. 1992 Compact finite-difference schemes with spectral-like resolution. J. Comp. Phys. 103, 16–42.
  * Magic et al. (2013) Magic, Z., Collet, R., Asplund, M., Trampedach, R., Hayek, W., Chiavassa, A., Stein, R. F. & Nordlund, Å. 2013 The Stagger-grid: A grid of 3D stellar atmosphere models–I. Methods and general properties. Astro. & Astrophys. 557, A26.
  * Nordlund et al. (2009) Nordlund, Å., Stein, R. F. & Asplund, M. 2009 Solar surface convection. Living Rev. Solar Phys. 6 (1), 1–117.
  * Ogura & Phillips (1962) Ogura, Y. & Phillips, N. A. 1962 Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci 19 (2), 173–179.
  * Rincon & Rieutord (2018) Rincon, F. & Rieutord, M. 2018 The Sun’s supergranulation. Living Rev. Solar Phys. 15, 1–74.
  * Scheel & Schumacher (2014) Scheel, J. D. & Schumacher, J. 2014 Local boundary layer scales in turbulent Rayleigh–Bénard convection. J. Fluid. Mech. 758, 344–373.
  * Schumacher & Sreenivasan (2020) Schumacher, J. & Sreenivasan, K. R. 2020 Colloquium: Unusual dynamics of convection in the Sun. Rev. Mod. Phys. 92 (4), 041001.
  * Smits et al. (2011) Smits, A. J., McKeon, B. J. & Marusic, I. 2011 High–Reynolds number wall turbulence. Annu. Rev. Fluid Mech. 43, 353–375.
  * Verhoeven et al. (2015) Verhoeven, J., Wiesehöfer, T. & Stellmach, S. 2015 Anelastic versus fully compressible turbulent Rayleigh–Bénard convection. Astrophys.J. 805 (1), 62\.


# A Robust Lane Detection Associated with Quaternion Hardy Filter

Wenshan Bi, Dong Cheng and Kit Ian Kou This work was supported in part by the
Science and Technology Development Fund, Macau SAR FDCT/085/2018/A2 and the
Guangdong Basic and Applied Basic Research Foundation (No. 2019A1515111185).
W. Bi and K.I. Kou are with the Department of Mathematics, Faculty of Science
and Technology, University of Macau, Macau, China (e-mail:
<EMAIL_ADDRESS>kikou@um.edu.mo).D. Cheng is with the Research Center for
Mathematics and Mathematics Education, Beijing Normal University at Zhuhai,
519087, China (e-mail: chengdong720@163.com).

###### Abstract

In this article, a robust color-edge feature extraction method based on the
Quaternion Hardy filter is proposed. The Quaternion Hardy filter is an
emerging edge detection theory. It is along with the Poisson and conjugate
Poisson smoothing kernels to handle various types of noise. Combining with the
Quaternion Hardy filter, Jin’s color gradient operator and Hough transform,
the color-edge feature detection algorithm is proposed and applied to the lane
marking detection. Experiments are presented to demonstrate the validity of
the proposed algorithm. The results are accurate and robust with respect to
the complex environment lane markings.

Keywords:Lane line detection, Hilbert transform, Quaternion Hardy filter,
Color image.  AMS Class Classification: 30E25, 11E88, 42B10

A great deal of research has recently been conducted on intelligent vehicle
systems in color image processing [1, 2, 3]. It consists of road lane
detection, pedestrian avoidance, and obstacle detection. Among them, road lane
detection is regarded as the most indispensable. The edge detector is
considered to be the basic module for road lane detection. Sangwine et al [4,
5] introduced a well established color image edge detection algorithm.

In order to derive a robust and discriminative performance on lane marking,
challenging scenarios such as shadows caused by illumination or weather
changes are necessary. In 2014, Wu et al [6] introduced the lane-mark detector
which is combining with Kalman filter and multi-adaptive thresholds method.
Later on, Li et al [7] in 2016 proposed the lane-mark detector based on Canny
filter and Hough transform to extract lane boundaries or markings in complex
scenarios. In order to provide robust and discriminative performance on lane
marking, edge detectors which can handle various type of noises are necessary.
Quaternion Hardy function [8] and Quaternion analytic signal [8] consist of
Poisson and conjugate Poisson kernels. In 2016, Hu et al analyzed the
Quaternion Hardy function [8] and Quaternion analytic signal [8] for edge
detection problem. It was proved that these edge detectors has significant
properties associated with Quaternion Fourier transform. They can handle
various types of noises and may be favorable for lane markings in complex
scenarios.

Quaternion analysis recently draws wide attention in color image processing
[8, 9, 10, 11, 12]. The recently developed concepts of quaternion Fourier
transform [13, 14, 15, 16, 17], quaternion representation [10] and quaternion
analytic signal [8, 9, 18], which are based on quaternion algebra, have been
found to be an indispensable tool in the representation of the
multidimensional signal. In earlier 1992, Ell [19] studied the two-sided
quaternion Fourier transform and apply it to the color image analysis. In
[20], Pei et al. gave a detailed interpretation of the quaternion Fourier
transform. In [13, 14, 15], Hitzer proposed the generalized quaternion Fourier
transform by follow-up the works of Ell, Sanguine, and Ernst. As the
generalization of the quaternion Fourier transform, quaternion linear
canonical transform was firstly studied in Kou er al. [11]. In 2017, Kou er
al. [9] studied the properties of quaternion linear canonical transform and
applied it to envelop detection. In [16], Dong et al. give the Plancherel and
inversion theorems of quaternion Fourier transform in the square-integrable
signals space. In 2016, Zou er al. [10] proposed two quaternion
representation-based classification algorithms and applied them to color face
recognition. Quaternion representation processes the color image in a holistic
way and the structural information between color channels is preserved.
Obviously, the experiments show the superiority of the proposed algorithm for
color reconstruction and face recognition. In 2017, the quaternion analytic
signal are studied and applied to edge detection [8].

The goal of this study is introducing a novel lane-mark extraction algorithm
which can handle complex scenarios. Firstly, edge enhancement based on
Quaternion Hardy filter [18] was used to enhance the pre-processing lane
image’s edges. Then, the edges were detected by Jin’s color detector [21]
which is based on Di Zenzo’s multichannel image gradient operator and combined
with the straight lines which were detected by Hough transform [22]. Moreover,
an inside lane line extraction algorithm was proposed on the basis of slope
constraint. Finally, the aim of marking the inside lane was realized. An
overview flowchart of the proposed lane-mark extraction algorithm is shown in
Fig. 1. By experimenting with Canny lane markers, the proposed algorithm can
realize the detection. Moreover, the approach can overcome the influence of
uneven light and be able to eliminate the interference from the side lane line
and the falling leaf, etc. Lane line detection is conducive to the vehicle
running on its road. The contributions of this paper are summarized as
follows.

  * 1.

A novel lane detector associated with Quaternion Hardy filter for color-edge
road lane extraction is proposed and analyzed. Firstly, input the color road
image, then edge enhancement based on Quaternion Hardy filter was applied to
obtain the ”smooth” lane images. Second, the edges were recorded by Jin’s
color gradient detector. Finally, combine with the straight lines obtained by
Hough transform.

  * 2.

We have carried out experimental research on the identification of road
markings in various complex environments. Experimental results show the
effectiveness and flexibility of the proposed method.

The key notations and acronyms used in this article are summarized in Table 1.
The rest of this paper is organized as follows. Section 1 introduces some
preliminaries of the quaternions, quaternion Fourier transform, and quaternion
Hardy space. The definitions of quaternion analytic signal and quaternion
Hardy filter are recalled in Section 2. Section 3 introduces the proposed lane
line detector, and the experimental results are presented in Section 4.
Finally, concluding remarks are drawn in Section 5.

Table 1: Key notations and acronyms used in this paper. Notation | Description
---|---
$\mathbb{R}^{2}$ | 2-dimensional Euclidean space
$\mathbb{H}$ | quaternion space
$\mathbf{Q}^{2}(\mathbb{C}_{\mathbf{i}\mathbf{j}}^{+})$ | quaternion Hardy space
$\mathcal{F}_{D}[\cdot]$ | discrete quaternion Fourier transform
${\mathcal{F}}_{D}^{-1}[\cdot]$ | discrete inverse quaternion Fourier transform
$\mathcal{F}[\cdot]$ | quaternion Fourier transform
${\mathcal{F}}^{-1}[\cdot]$ | inverse quaternion Fourier transform
$f_{q}$ | quaternion analytic signal
QHF | Quaternion Hardy filter

##  1 Preliminaries

In order to mention the proposed method, some basic knowledge about quaternion
algebra [11], quaternion Fourier transform [16, 19], quaternion Hardy space
[8], Jin’s gradient operator [21] and Hough transform [23] are reviewed as
follows.

###  1.1 Quaternion algebra and quaternion Fourier transform

Let $\mathbb{H}$ denote the Hamiltonian skew field of quaternion

$\displaystyle\mathbb{H}:=\\{q=q_{0}+q_{1}\mathbf{i}+q_{2}\mathbf{j}+q_{3}\mathbf{k}|\
q_{0},q_{1},q_{2},q_{3}\in\mathbb{R}\\},$ (1.1)

which is an associative anti-commutative four-dimensional algebra. A
quaternion-valued function $f:\mathbb{R}^{2}\rightarrow\mathbb{H}$ can be
represented by

$\displaystyle
f(x_{1},x_{2})=f_{0}(x_{1},x_{2})+f_{1}(x_{1},x_{2})\mathbf{i}+f_{2}(x_{1},x_{2})\mathbf{j}+f_{3}(x_{1},x_{2})\mathbf{k},$
(1.2)

where $f_{n}:\mathbb{R}^{2}\rightarrow\mathbb{R}(n=0,1,2,3)$.

###  1.2 Quaternion Hardy space

Let $f\in L^{1}(\mathbb{R}^{2},\mathbb{H})$ . Then the quaternion Fourier
transform (QFT) [16, 24, 25] of $f$ is defined by

$\mathcal{F}[f](w_{1},w_{2}):=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}e^{-\mathbf{i}w_{1}x_{1}}f(x_{1},x_{2})e^{-\mathbf{j}w_{2}x_{2}}d{x_{1}}d{x_{2}},$
(1.3)

where $w_{l}(l=1,2)$ denote the 2D angular frequencies. Furthermore, if $f$ is
an integrable $\mathbb{H}$-valued function defined on $\mathbb{R}^{2}$, the
inverse quaternion Fourier transform (IQFT) [16, 24, 25] of $f$ is defined by

${\mathcal{F}}^{-1}[f](x_{1},x_{2}):=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}e^{\mathbf{i}w_{1}x_{1}}f(w_{1},w_{2})e^{\mathbf{j}w_{2}x_{2}}d{w_{1}}d{w_{2}}.$
(1.4)

Let ${\mathbb{C}}^{+}:=\\{z|z=x+si,x,x\in\mathbb{R},s>0\\}$, namely upper half
complex plane and ${\mathbf{H}}^{2}(\mathbb{C}^{+})$ be the Hardy space on the
upper half complex plane. Let
$\mathbb{C}_{\mathbf{i}\mathbf{j}}:=\\{(z_{1},z_{2})|z_{1}=x_{1}+s_{1}\mathbf{i},z_{2}=x_{2}+s_{2}\mathbf{j},x_{l},s_{l}\in\mathbb{R},l=1,2\\}$.
And a subset of $\mathbb{C}_{\mathbf{i}\mathbf{j}}$ is defined by
$\mathbb{C}_{\mathbf{i}\mathbf{j}}^{+}:=\\{(z_{1},z_{2})|z_{1}=x_{1}+s_{1}\mathbf{i},z_{2}=x_{2}+s_{2}\mathbf{j},x_{l},s_{l}\in\mathbb{R},s_{l}>0,l=1,2\\}.$
The quaternion Hardy space
$\mathbf{Q}^{2}(\mathbb{C}_{\mathbf{i}\mathbf{j}}^{+})$ consists of all
functions $h$ satisfying [Hu1],

$\displaystyle\left\\{\begin{aligned}
&\frac{\partial}{\partial\overline{z_{1}}}h(z_{1},z_{2})=0;\\\
&h(z_{1},z_{2})\frac{\partial}{\partial\overline{z_{2}}}=0;\\\
&\left(\sup\limits_{\begin{subarray}{c}s_{1}>0\\\
s_{2}>0\end{subarray}}\int_{\mathbb{R}^{2}}|h(x_{1}+s_{1}\mathbf{i},x_{2}+s_{2}\mathbf{j})|^{2}dx_{1}dx_{2}\right)^{\frac{1}{2}}<\infty,\end{aligned}\right.$
(1.5)

where
$\frac{\partial}{\partial\overline{z_{1}}}:=\frac{\partial}{\partial{x_{1}}}+\mathbf{i}\frac{\partial}{\partial{s_{1}}}$,
$\frac{\partial}{\partial\overline{z_{2}}}:=\frac{\partial}{\partial{x_{2}}}+\mathbf{j}\frac{\partial}{\partial{s_{2}}}$.

###  1.3 Jin’s gradient operator

Let $g$ be an $M\times{N}$ color image that maps a point $(x,y)$ to a vector
$(g_{1}(x,y)$, $g_{2}(x,y)$, $g_{3}(x,y))$. Specifically, for color image $g$,
let $(g_{1}(x,y)$, $g_{2}(x,y)$, $g_{3}(x,y))$ denotes the red, the green, and
the blue channel, respectively. Then the square of the variation of $f$ at the
position $(x,y)$ with the distance $b$ in the direction $\theta$ is given by

$\displaystyle\begin{aligned}
dg^{2}&:=\|g(x+b{\cos\theta,y+b{\sin\theta})-g(x,y)}\|_{2}^{2}\\\
&\approx\sum\limits_{i=1}^{3}\left(\frac{\partial{g_{i}}}{\partial{x}}b\cos\theta+\frac{\partial{g_{i}}}{\partial{y}}b\sin\theta\right)^{2}\\\
&=b^{2}f(\theta),\end{aligned}$ (1.6)

where

$\displaystyle\begin{aligned}
g(\theta):=&2\sum\limits_{i=1}^{3}\frac{\partial{g_{i}}}{\partial{x}}\frac{\partial{g_{i}}}{\partial{y}}\cos\theta\sin\theta\\\
&+\sum\limits_{i=1}^{3}\left(\frac{\partial{g_{i}}}{\partial{x}}\right)^{2}\cos^{2}\theta+\sum\limits_{i=1}^{3}\left(\frac{\partial{g_{i}}}{\partial{y}}\right)^{2}\sin^{2}\theta.\end{aligned}$
(1.7)

Define

$\displaystyle\left\\{\begin{aligned}
&H:=\sum\limits_{i=1}^{3}\left(\frac{\partial{g_{i}}}{\partial{x}}\right)^{2};\\\
&J:=\sum\limits_{i=1}^{3}\left(\frac{\partial{g_{i}}}{\partial{y}}\right)^{2};\\\
&K:=\sum\limits_{i=1}^{3}\frac{\partial{g_{i}}}{\partial{x}}\frac{\partial{g_{i}}}{\partial{y}}.\end{aligned}\right.$
(1.8)

Then the gradient magnitude $f_{\max}$ of the Jin’s gradient operator is given
by

$\displaystyle\begin{aligned} g_{\max}(\theta_{\max}):&=\max_{0\leq\theta\leq
2\pi}{g(\theta)}\\\
=&\frac{1}{2}\bigg{(}H+K+\sqrt{(H-K)^{2}+(2J)^{2}}\bigg{)}.\end{aligned}$
(1.9)

The gradient direction is defined as the value $\theta_{\max}$ that maximizes
$g(\theta)$ over $0\leq\theta\leq 2\pi$

$\displaystyle\theta_{\max}:=$
$\displaystyle\mbox{sgn}(J)\arcsin\bigg{(}\frac{g_{\max}-H}{2g_{\max}-H-K}\bigg{)},$
(1.10)

where $(H-K)^{2}+J^{2}\neq 0$,
$\mbox{sgn}(J)=\left\\{\begin{array}[]{ll}1,&{J\geq 0;}\\\
-1,&{J<0.}\end{array}\right.$ While $(H-K)^{2}+J^{2}=0$, $\theta_{\max}$ is
undefined.

###  1.4 Hough transform

The Hough transform is done on a binary image, obtained after processing the
original image by an edge detector [26, 23, 27]. This is an ingenious method
that converts such global curve detection problem into an efficient peak
detection problem in parameter space. In addition, Dorst et al [28, 29, 30,
31] studied the application of plane-based geometric algebra in the
representation of straight lines. The basic idea of Hough transform is to use
the duality relationship between the points and lines in the image domain and
the parameter domain. In particular, the slope-intercept form is a classical
parametric form for straight line detection, as is shown in Eq. (1.11):

$\displaystyle f(x,y)=y-kx-b=0,$ (1.11)

where $(x,y)$ is the coordinate of the pixel being mapped, and the parameters
$k$ and $b$ are the slop and y-intercept of the line respectively. It’s worth
noting that this method is sensitive to the choice of coordinate axes on the
image plane because both the parameters become unbounded when $k$ is infinite.

Duda and Hart [32] lucidly resolved the issue of unboundedness by mapping a
point in image space to a sinusoidal curve in $\rho-\theta$ parameter space.
Here we use HoughlinesP method to detect lines, which is based on statistics
not only could detect the two endpoints of the straight lines but also have
the advantage of high efficiency. $\rho$ and $\theta$ satisfy the following
relationship:

$\displaystyle\rho(\theta)=x\cos{\theta}+y\sin{\theta},$ (1.12)

where $\theta$ is the angle between the line and the $Y$-axis, $\rho$ is the
distance from the origin to the line in the image coordinate, $\rho>0$ and
$\rho(\theta+\pi)=-\rho(\theta)$.

The first step of Hough transform is initialize accumulator to all zeros,
followed by the parameter space is quantized in intervals of $\Delta\rho$ and
$\Delta\theta$ and corresponding accumulators are created to collect the
evidence (or vote) of object pixels satisfying Eq.(1.12). Consider the ($m$,
$n$)th accumulator corresponding to intervals [($m-1)\Delta\rho$,
$m\Delta\rho$) and [($n-1)\Delta\theta$, $n\Delta\theta$ ) where
$m,n\in\mathcal{Z}^{+}$ . For each object pixel satisfying Eq.(1.12) with
parameters in the above range, the vote in the ($m$, $n$)th accumulator is
incremented by 1. In the peak detection phase, the accumulators having number
of votes above a critical threshold correspond to straight lines in the image.
Therefore, the line detection problem is converted to peak statistics
problems.

##  2 Quaternion Hardy filter

The analytic signal [33] of a given real signal $f$ is defined by

$f_{a}(x):=f(x)+i\mathcal{H}[f](x),\;\;x\in\mathbb{R},$ (2.1)

where $\mathcal{H}[f]$ denotes the Hilbert transform [34] of $f$. By a direct
computation, the Fourier transform of $f_{a}$ is given by

$\widehat{f_{a}}(w)=\left(1+\mbox{sgn}(w)\right)\widehat{f}(w),\;\;w\in\mathbb{R}.$
(2.2)

Figure 1: The flowchart of the proposed lane detection method.

The analytic signal of $f$ can be regarded as the output signal of a filter
with input $f$. The system function of this filter is given by

$H_{1}(w):=1+\mbox{sgn}(w).$ (2.3)

Motivated by the success of analytic signal in the real signal processing, we
try to use the alternative tool in the 2D signal processing.

Give a 2D quaternion-valued signal $f$. The quaternion analytic signal [35]
$f_{q}$ is defined by

$f_{q}(x_{1},x_{2}):=f(x_{1},x_{2})+\mathbf{i}\mathcal{H}_{x_{1}}[f](x_{1})+\mathcal{H}_{x_{2}}[f](x_{2})\mathbf{j}+\mathbf{i}\mathcal{H}_{{x_{1}}{x_{2}}}[f](x_{1},x_{2})\mathbf{j},$
(2.4)

where $\mathcal{H}_{x_{1}}$, $\mathcal{H}_{x_{2}}$ are the partial Hilbert
transforms with respect to $x_{1},x_{2}$ respectively, while
$\mathcal{H}_{{x_{1}}{x_{2}}}$ the 2D Hilbert transform. By a direct
computation [Kou2], the quaternion Fourier transform of $f_{q}$ is written as

$\mathcal{F}[{f_{q}}](w_{1},w_{2})=[1+\mbox{sgn}(w_{1})][1+\mbox{sgn}(w_{2})]\mathcal{F}[{f}](w_{1},w_{2}).$
(2.5)

The quaternion analytic signal of $f$ can be regarded as the output signal of
a filter with input $f$. The system function of this filter is given by

$H_{2}(w_{1},w_{2}):=[1+\mbox{sgn}(w_{1})][1+\mbox{sgn}(w_{2})].$ (2.6)

By the inverse quaternion Fourier transform on the both side of (2.5), we
obtain

$f_{q}(x_{1},x_{2})=\mathcal{F}^{-1}[H_{2}\mathcal{F}[{f}]](x_{1},x_{2}).$
(2.7)

In this paper, we apply quaternion Hardy filter as a lane detector of color
image. The system function of quaternion Hardy filter is defined by

$\displaystyle H_{3}(w_{1},w_{2},s_{1},s_{2}):=$ $\displaystyle
e^{-\mid{w_{1}}\mid{s_{1}}}e^{-\mid{w_{2}}\mid{s_{2}}}H_{2}(w_{1},w_{2})$
(2.8) $\displaystyle=$ $\displaystyle
e^{-\mid{w_{1}}\mid{s_{1}}}e^{-\mid{w_{2}}\mid{s_{2}}}[1+\mbox{sgn}(w_{1})][1+\mbox{sgn}(w_{2})].$

Based on the above analysis, we have the following theorem.

Figure 2: The first row is the original image in rainy night driving
condition. The second row and the third row are the lane detection results
which are captured by Jin’s algorithm and the the proposed algorithm,
respectively.

###### Theorem 2.1

Let $f\in L^{2}(\mathbb{R}^{2},\mathbb{H})$, and

$\displaystyle
f_{Q}(x_{1},x_{2},s_{1},s_{2})={\mathcal{F}}^{-1}[H_{3}\mathcal{F}[f]](x_{1},x_{2},s_{1},s_{2}).$
(2.9)

Where $H_{3}$ is given by (2.8). Then
$f_{Q}\in\mathbf{Q}^{2}(\mathbb{C}_{\mathbf{i}\mathbf{j}}^{+})$.

The proof of Theorem 2.1 can be found in [QHF]. It states that the output
signal of the quaternion Hardy filter belongs to the quaternion Hardy space
$\mathbf{Q}^{2}(\mathbb{C}_{\mathbf{i}\mathbf{j}}^{+})$.

###### Remark 2.2

The quaternion analytic signal is a special case of quaternion Hardy filter.
That is, when $s_{1}=s_{2}=0$ in Theorem 2.1, the system function
$H_{3}(x_{1},x_{2},0,0)$ reduces to $H_{2}(w_{1},w_{2})$, and the signal
$f_{Q}(x_{1},x_{2},0,0)$ reduces to $f_{q}(x_{1},x_{2})$ which is defined by
equation (2.7).

###### Remark 2.3

In image processing, quaternion Hardy filter can not only suppress the high
frequency of image, but also enhance the edge feature. It’s worth noting that,
the system function $H_{3}(x_{1},x_{2},s_{1},s_{2})$ plays a key role in
quaternion Hardy filter. It has the advantage of both Hilbert transform and
low-pass filter. On the one hand, the factor
$[1+\mbox{sgn}(w_{1})][1+\mbox{sgn}(w_{2})]$ in
$H_{3}(x_{1},x_{2},s_{1},s_{2})$ performs Hilbert transform on the input
signal. It can enhance the edge feature. On the other hand,
$e^{-\mid{w_{1}}\mid{s_{1}}}e^{-\mid{w_{2}}\mid{s_{2}}}$ could improve the
ability of noise immunity for the quaternion Hardy filter.

##  3 The proposed lane detection algorithm

The flowchart is shown in Fig. 1. In the following, let’s consider the details
of the proposed algorithm. They are divided by the following steps.

  * Step 1. Input the road image $f(x_{1},x_{2})$. By quaternion representation, it can be written as

$\displaystyle
f(x_{1},x_{2})=f_{R}(x_{1},x_{2})\mathbf{i}+f_{G}(x_{1},x_{2})\mathbf{j}+f_{B}(x_{1},x_{2})\mathbf{k},$
(3.1)

where $f_{R},f_{G}$ and $f_{B}$ represent the red, green and blue components
of color image $f$, respectively.

  * Step 2. Compute the discrete quaternion Fourier transform (DQFT) [36], $\mathcal{F}_{D}[f](w_{1},w_{2})$, of the image $f$ from Step 1, we have

$\displaystyle\mathcal{F}_{D}[f](w_{1},w_{2})=\sum_{x_{1}=0}^{M-1}\sum_{x_{2}=0}^{M-1}e^{-\mathbf{i}2{\pi}(w_{1}x_{1}/M)}f(x_{1},x_{2})e^{-\mathbf{j}2{\pi}(w_{2}x_{2}/M)}.$
(3.2)

Figure 3: The first row is the original image in snowy night driving
condition. The second row and the third row are the lane detection results
which are captured by Jin’s algorithm and the the proposed algorithm,
respectively.

  * Step 3. For the fixed positive $s_{1}$ and $s_{2}$, calculate $H_{3}$ by Eq. (2.8), we have

$H_{3}(w_{1},w_{2},s_{1},s_{2}):=e^{-\mid{w_{1}}\mid{s_{1}}}e^{-\mid{w_{2}}\mid{s_{2}}}[1+\mbox{sgn}(w_{1})][1+\mbox{sgn}(w_{2})].$

Multiplying $\mathcal{F}_{D}[f]$ with $H_{3}$, from Step 2 and 3,
respectively. By Theorem 2.1, we have

$\displaystyle\mathcal{F}_{D}[f_{Q}](w_{1},w_{2})=H_{3}(w_{1},w_{2},s_{1},s_{2})\mathcal{F}_{D}[f](w_{1},w_{2}).$
(3.3)

  * Step 4. Compute the inverse DQFT of $\mathcal{F}_{D}[f_{Q}]$, from Step 3,

$\displaystyle\begin{aligned}
f_{Q}(x_{1},x_{2})=&\sum_{w_{1}=0}^{M-1}\sum_{w_{2}=0}^{M-1}e^{\mathbf{i}2{\pi}(w_{1}x_{1}/M)}\mathcal{F}_{D}[f_{Q}](w_{1},w_{2})e^{\mathbf{j}2{\pi}(w_{2}x_{2}/M)}\\\
=&\sum_{w_{1}=0}^{M-1}\sum_{w_{2}=0}^{M-1}e^{\mathbf{i}2{\pi}(w_{1}x_{1}/M)}[1+\mbox{sgn}(w_{1})][1+\mbox{sgn}(w_{2})]e^{-{3}\mid{w_{1}}\mid}e^{-{20}\mid{w_{2}}\mid}\\\
&\mathcal{F}_{D}[f](w_{1},w_{2})e^{\mathbf{j}2{\pi}(w_{2}x_{2}/M)}.\end{aligned}$
(3.4)

We obtain $f_{Q}$.

  * Step 5. Extract the vector part of $f_{Q}$, we obtain

$\displaystyle\mathbf{Vec}(f_{Q}):=h_{1}\mathbf{i}+h_{2}\mathbf{j}+h_{3}\mathbf{k},$
(3.5)

where $h_{k}$, $k=1,2,3$ are real-valued functions.

  * Step 6. Applying the Jin’s gradient operator [21] to $\mathbf{Vec}(f_{Q})$.

  * Step 7. Furthermore, perform the Hough operator [22].

  * Step 8. Finally, we show the final result by placing the resulting lane lines on the original image.

##  4 Experiment

To verify the effectiveness of quaternion hardy filter-based lane detection
algorithm, we conduct experiments on 27 video clips. The experiments are
programmed in Matlab R2016b. The proposed algorithm is compared with previous
work [1], [37], [38], [39], [40], [41], and [21]. Here, Jin [21] means that
the proposed method without steps 2-4 and provides a baseline reference
results. The previous work involved six different approaches that are typical
of the field of lane line detection.
The video clips used in the experiment were from public dataset [42], which
were shot using dashboard cameras in the US and Korea. This data set contains
a total of 49 video clips, and 27 are randomly selected for this experiment.
The 27 video clips consist of 44159 frame images which divided into 6
different weather conditions: day-clear (5), day-rainy (5), day-snowy (5),
night-clear (5), night-rainy (5) and night-snowy (2). The number in
parentheses after the weather label represent the number of video clips used
in the experiment.

###  4.1 Visual results

Visual results including rain night and snow night are illustrate in Fig. 2
and Fig. 3, respectively. When car lights form bright lines on the road on a
rainy night, the comparison method may fail. It affects the detection of lane
lines. For example, in the second row and fourth column of Fig. 2, Jin’s
method is obviously affected by the light, thus making a misjudgment of the
lane line. In addition, Jin’s method (e.g., the second row and first column in
Fig. 3) detects the boundary line of the navigation screen and ignores the
actual lane line on the right. Gradient-based approaches [39] may preserve the
boundary portion of the navigation screen while ignoring the actual lane lines
when they are reflective navigation screens are present. In general, the
proposed method can achieve satisfactory results under complex weather
conditions and illumination conditions. The proposed method cannot, however,
detect lane markings when affected by surrounding buildings or windshield
wipers, as shown in Fig. 4.

###  4.2 Quantitative comparisons

For each image, two boundary lines are drawn on the long edge of the lane mark
in the presence of the lane mark, and the real value of the ground is obtained
manually. The detection rate is obtained by the following

$\displaystyle DR=\frac{correct}{total}\times 100\%$ (4.1)

where $correct$ represents the number of images correctly detected and $total$
represents the total number of images participated in the detection.

Figure 4: The first row is the original image in complex driving condition.
The second row and the third row are the false lane detection results which
are captured by Jin’s algorithm and the the proposed algorithm, respectively.

We present a quantitative evaluation for seven methods in Table 2. It can be
seen from Table 2 that all the algorithms can generate good results when lane
images are captured with good conditions, as shown in clear days and nights.
For clear day sets, QHF-based method achieves a 99.4% accuracy among the
compared methods. For rain day and snow day sets, the proposed method keeps
steady in the accuracy. On clear night sets, the accuracy of the proposed
method is still in the Top three. In the rainy night data set, the accuracy of
the proposed method was still 88.5%, despite the fact that raindrops falling
on the windshield caused considerable interference with the lane lines in the
image. For the snowy night data set, the other six methods were not
considered, and the proposed method was accurate to 95.4%. We have to say that
there are some situations that are hard to dealt with. Fig. 4 gives several
failure samples. For those images affected by wipers, illuminations and
weather conditions, our method faces challenges. A potential solution to solve
this problem is to use $\Lambda$ ROI (region of interest) [41]. However, ROI-
based lane detection method require much computation cost. Therefore, more
optimization works needs to be done in the future.

Table 2: Comparison of detection rate of different algorithms for color image
lane detection on 44159 images .

DR (%) | Son[1] | Aly[37] | Jung[38] | Yoo[39] | Liu[40] | Lee[41] | Jin[21] | Ours
---|---|---|---|---|---|---|---|---
Clear-Day | 93.1 | 96.3 | 88$\sim$98 | 96.36 | - | 99 | 99.3 | 99.4
Rain-Day | 89 | - | - | 96.19 | 98.67 | 96.9 | 96.8 | 97
Snow-Day | - | - | - | - | - | 93 | 94.2 | 95.5
Clear-Night | 94.3 | - | 87.7 | 96.5 | 94.06 | 98 | 95.7 | 96.1
Rain-Night | 93 | - | - | - | - | 93.6 | 79.4 | 88.5
Snow-Night | - | - | - | - | - | - | 89.9 | 95.4
Processing |  |  |  |  |  |  |  |
time/ | 33 | 20 | 3.9 | 50 | 200 | 35.4 | 300 | 800
frame(ms) |  |  |  |  |  |  |  |
Platform | PC | PC | PC | PC | PC | ARM | PC | PC
CPU | Unknow | Core 2 | 2.8GHz | Unknow | 2.3GHz | A9 GHz | 2.3GHz | 2.3GHz
|  | @2.4GHz |  |  |  | @800MHz |  |

##  5 Conclusion and future work

In this paper, we proposed a QHF-based lane detection method for color image.
It exploits quaternion Hardy filter to generate clear feature image, and it
use Jin’s gradient operator and Hough transform to provided lane lines, which
ensures accuracy lane lines. Our method could handle with multiple situations
such as various weather conditions, illuminations and lane markings.
Experiments conducted on the dataset indicate that the proposed method
achieves satisfactory results. The results were accurate and robust with
respect to the complex environment lane markings. Finally, we are interested
in extending our proposed quaternion lane line detection idea to color video
for lane line detection in the future. In order to optimize the algorithm, we
will also consider fast algorithms [43] and algorithm optimization software
[44].

## Acknowledgment

This research was supported by Macao Science and Technology Development Fund
(No. FDCT/085/2018/A2).

## References

  * [1] Jongin Son, Hunjae Yoo, Sanghoon Kim, and Kwanghoon Sohn.  Real-time illumination invariant lane detection for lane departure warning system.  Expert Systems with Applications, 42(4):1816–1824, 2015.
  * [2] Robert E Ralston and Justin D Bloom.  Dynamic routing intelligent vehicle enhancement system, September 6 2016\.  US Patent 9,435,652.
  * [3] Yun Luo and Dieter Hoetzer.  Video based intelligent vehicle control system, October 20 2015.  US Patent 9,165,468.
  * [4] Todd A Ell, Nicolas Le Bihan, and Stephen J Sangwine.  Quaternion Fourier transforms for signal and image processing.  John Wiley & Sons, 2014.
  * [5] Stephen John Sangwine.  Colour image edge detector based on quaternion convolution.  Electronics Letters, 34(10):969–971, 1998.
  * [6] Pei-Chen Wu, Chin-Yu Chang, and Chang Hong Lin.  Lane-mark extraction for automobiles under complex conditions.  Pattern Recognition, 47(8):2756–2767, 2014.
  * [7] Yadi Li, Liguo Chen, Haibo Huang, Xiangpeng Li, Wenkui Xu, Liang Zheng, and Jiaqi Huang.  Nighttime lane markings recognition based on canny detection and hough transform.  In 2016 IEEE International Conference on Real-time Computing and Robotics (RCAR), pages 411–415. IEEE, 2016.
  * [8] Xiao-Xiao Hu and Kit Ian Kou.  Phase-based edge detection algorithms.  Mathematical Methods in the Applied Sciences, 41(11):4148–4169, 2018.
  * [9] Kit Ian Kou, Ming-Sheng Liu, João Pedro Morais, and Cuiming Zou.  Envelope detection using generalized analytic signal in 2d qlct domains.  Multidimensional Systems and Signal Processing, 28(4):1343–1366, 2017.
  * [10] Cuiming Zou, Kit Ian Kou, and Yulong Wang.  Quaternion collaborative and sparse representation with application to color face recognition.  IEEE Transactions on image processing, 25(7):3287–3302, 2016.
  * [11] Kit Ian Kou, Jianyu Ou, and Joao Morais.  Uncertainty principles associated with quaternionic linear canonical transforms.  Mathematical Methods in the Applied Sciences, 39(10):2722–2736, 2016.
  * [12] Yan Yang, Kit Ian Kou, and Cuiming Zou.  Edge detection methods based on modified differential phase congruency of monogenic signal.  Multidimensional Systems and Signal Processing, 29(1):339–359, 2018\.
  * [13] Eckhard Hitzer.  The quaternion domain fourier transform and its properties.  Advances in Applied Clifford Algebras, 26(3):969–984, 2016.
  * [14] Eckhard Hitzer and Stephen J Sangwine.  The orthogonal 2d planes split of quaternions and steerable quaternion fourier transformations.  In Quaternion and Clifford Fourier transforms and wavelets, pages 15–39. Springer, 2013.
  * [15] Eckhard Hitzer.  General two-sided quaternion fourier transform, convolution and mustard convolution.  Advances in Applied Clifford Algebras, 27(1):381–395, 2017.
  * [16] Cheng Dong and Kou Kit Ian.  Plancherel theorem and quaternion fourier transform for square integrable functions.  complex variables and elliptic equations, 64(2):223–242, 2019.
  * [17] Xiao-Xiao Hu and Kit Ian Kou.  Quaternion fourier and linear canonical inversion theorems.  Mathematical Methods in the Applied Sciences, 40(7):2421–2440, 2017\.
  * [18] Wenshan Bi, Dong Cheng, and Kit Ian Kou.  A robust color edge detection algorithm based on quaternion hardy filter.  arXiv preprint arXiv:2001.01800, 2020.
  * [19] Todd Anthony Ell.  Hypercomplex spectral transformations.  PhD thesis, University of Minnesota, 1992.
  * [20] Soo-Chang Pei, Jian-Jiun Ding, and Ja-Han Chang.  Efficient implementation of quaternion fourier transform, convolution, and correlation by 2-d complex fft.  IEEE Transactions on Signal Processing, 49(11):2783–2797, 2001\.
  * [21] Lianghai Jin, Hong Liu, Xiangyang Xu, and Enmin Song.  Improved direction estimation for di zenzo’s multichannel image gradient operator.  Pattern recognition, 45(12):4300–4311, 2012.
  * [22] Dana H Ballard.  Generalizing the hough transform to detect arbitrary shapes.  Pattern recognition, 13(2):111–122, 1981.
  * [23] Zezhong Xu, Bok-Suk Shin, and Reinhard Klette.  Accurate and robust line segment extraction using minimum entropy with hough transform.  IEEE Transactions on Image Processing, 24(3):813–822, 2014.
  * [24] Eckhard MS Hitzer.  Quaternion fourier transform on quaternion fields and generalizations.  Advances in Applied Clifford Algebras, 17(3):497–517, 2007.
  * [25] Eckhard MS Hitzer.  Directional uncertainty principle for quaternion fourier transform.  Advances in Applied Clifford Algebras, 20(2):271–284, 2010.
  * [26] John Illingworth and Josef Kittler.  A survey of the hough transform.  Computer vision, graphics, and image processing, 44(1):87–116, 1988\.
  * [27] Khary Popplewell, Kaushik Roy, Foysal Ahmad, and Joseph Shelton.  Multispectral iris recognition utilizing hough transform and modified lbp.  In 2014 IEEE International Conference on Systems, Man, and Cybernetics (SMC), pages 1396–1399. IEEE, 2014.
  * [28] Leo Dorst.  A guided tour to the plane-based geometric algebra pga.
  * [29] Leo Dorst and Arnold WM Smeulders.  Best linear unbiased estimators for properties of digitized straight lines.  IEEE transactions on pattern analysis and machine intelligence, (2):276–282, 1986.
  * [30] Leo Dorst and Arnold WM Smeulders.  Discrete representation of straight lines.  IEEE Transactions on Pattern Analysis and Machine Intelligence, (4):450–463, 1984.
  * [31] Leo Dorst and Robert PW Duin.  Spirograph theory: A framework for calculations on digitized straight lines.  IEEE transactions on pattern analysis and machine intelligence, (5):632–639, 1984.
  * [32] Richard O Duda and Peter E Hart.  Use of the hough transformation to detect lines and curves in pictures.  Communications of the ACM, 15(1):11–15, 1972.
  * [33] John Garnett.  Bounded analytic functions, volume 236.  Springer Science & Business Media, 2007.
  * [34] Frederick W King.  Hilbert transforms, volume 1.  Cambridge University Press Cambridge, 2009.
  * [35] Thomas Bulow and Gerald Sommer.  Hypercomplex signals-a novel extension of the analytic signal to the multidimensional case.  IEEE Transactions on signal processing, 49(11):2844–2852, 2001\.
  * [36] Stephen John Sangwine.  The discrete quaternion fourier transform.  1997\.
  * [37] Mohamed Aly.  Real time detection of lane markers in urban streets.  In 2008 IEEE Intelligent Vehicles Symposium, pages 7–12. IEEE, 2008\.
  * [38] Soonhong Jung, Junsic Youn, and Sanghoon Sull.  Efficient lane detection based on spatiotemporal images.  IEEE Transactions on Intelligent Transportation Systems, 17(1):289–295, 2015.
  * [39] Hunjae Yoo, Ukil Yang, and Kwanghoon Sohn.  Gradient-enhancing conversion for illumination-robust lane detection.  IEEE Transactions on Intelligent Transportation Systems, 14(3):1083–1094, 2013.
  * [40] Guorong Liu, Shutao Li, and Weirong Liu.  Lane detection algorithm based on local feature extraction.  In 2013 Chinese Automation Congress, pages 59–64. IEEE, 2013.
  * [41] Lee Chanho and Ji-Hyun Moon.  Robust lane detection and tracking for real-time applications.  IEEE Transactions on Intelligent Transportation Systems, 19(12):4043–4048, 2018.
  * [42] https://drive.google.com/file/d/1315Ry7isciL3nRvU5SCXM_-4meR2MyI/view?usp=sharing.
  * [43] Eckhard Hitzer and J. Sangwine Stephen.  The orthogonal 2d planes split of quaternions and steerable quaternion fourier transformations.  Quaternion and Clifford Fourier transforms and wavelets, pages 15–39, 2013.
  * [44] Gaalop.  http://www.gaalop.de/dhilden/.