# EOG Artifact Removal from Single and Multi-channel EEG Recordings through
the combination of Long Short-Term Memory Networks and Independent Component
Analysis

Behrad TaghiBeyglou1111 B. TaghiBeyglou was with School of Electrical
Engineering, Sharif University of Technology, Tehran, Iran.
<EMAIL_ADDRESS>Fatemeh Bagheri2222 F. Bagheri was with
Biomedical Engineering Department, Amirkabir University of Technology, Tehran,
Iran. Institute of Biomedical Engineering, University of Toronto, Toronto,
Canada  KITE- Toronto Rehabilitation Institute, University Health Network,
Toronto, Canada  Deparment of Medical Biophysics, University of Toronto,
Toronto, Canada

###### Abstract

Introduction: Electroencephalogram (EEG) signals have gained significant
popularity in various applications due to their rich information content.
However, these signals are prone to contamination from various sources of
artifacts, notably the electrooculogram (EOG) artifacts caused by eye
movements. The most effective approach to mitigate EOG artifacts involves
recording EOG signals simultaneously with EEG and employing blind source
separation techniques, such as independent component analysis (ICA).
Nevertheless, the availability of EOG recordings is not always feasible,
particularly in pre-recorded datasets. Objective: In this paper, we present a
novel methodology that combines a long short-term memory (LSTM)-based neural
network with ICA to address the challenge of EOG artifact removal from
contaminated EEG signals. Approach: Our approach aims to accomplish two
primary objectives: 1) estimate the horizontal and vertical EOG signals from
the contaminated EEG data, and 2) employ ICA to eliminate the estimated EOG
signals from the EEG, thereby producing an artifact-free EEG signal. Main
results: To evaluate the performance of our proposed method, we conducted
experiments on a publicly available dataset comprising recordings from 27
participants. We employed well-established metrics such as mean squared error,
mean absolute error, and mean error to assess the quality of our artifact
removal technique. Significance: Furthermore, we compared the performance of
our approach with two state-of-the-art deep learning-based methods reported in
the literature, demonstrating the superior performance of our proposed
methodology.

###### keywords:

EEG , EOG removal , Long short-term memory , Independent component analysis ,
artifact removal, deep neural networks, signal processing.

††journal: arXiV

##  1 Introduction

Brain-computer interfaces (BCIs) are systems that transform raw brain
activities of humans into interpretable information that can control other
devices or reflect the state of rehabilitation or disease status [1, 2, 3].
One of the most prominent modalities to record brain activation is the
Electroencephalogram (EEG), which primarily records scalp voltages non-
invasively [4, 5, 6]. These voltages represent the superposition of many
neural activities and contribute to different tasks [7]. However, as EEG
recording is performed on the scalp, physiological activities other than brain
signals can interfere and introduce artifacts [8]. These unwanted
interferences, known as artifacts, need to be removed to ensure accurate EEG
analysis [9, 10].

One of the most significant physiological artifacts is the Electrooculogram
(EOG), which directly relates to eye movements and blinks. Due to the close
proximity between the eyes and the scalp, the effect of the EOG artifact is
undeniable [11, 12]. Eye blinking and movements generate spike-like signal
waveforms, with peak amplitudes reaching up to $800\mu$V and occurring within
a short period of $200$–$400$ ms [13, 14]. Additionally, EOG artifacts often
overlap with EEG signals in the frequency domain, particularly at low
frequencies, as well as in the time domain [1, 15]. Since the extraction of
features and information from EEG signals for classification and analysis
requires a clean signal, the removal of EOG artifacts is essential for
subsequent steps [9, 10].

Since the introduction of the EEG signal by Berger [16], researchers have
faced the challenge of EOG artifacts [17]. With the development of
mathematical tools, significant attention has been devoted to removing EOG
artifacts from EEG signals [18, 19, 20]. Most of these methods rely on blind
source separation (BSS), particularly independent component analysis (ICA).
However, these techniques are most effective when EOG signals are available
alongside EEG recordings. In cases where EOG signals are not available, which
is common in pre-recorded datasets, researchers have sought methods to
directly estimate EOG signals from contaminated or artifact-included EEG
signals. To provide a comprehensive review of the literature, we categorize
different EOG artifact removal methods into three categories: Single Channel,
Multi-channel EEG recording, and methods effective for both scenarios, as some
approaches claim to remove EOG artifacts using multi-channel EEG recordings,
while others rely on a single channel.

##  2 Literature review

###  2.1 Single Channel

In 1973, Girton et al. introduced a simple hardware-based method for online
removal of EOG artifacts from EEG recordings [20]. Their approach involved
using a circuit to record EOG signals from two electrodes (related to
Horizontal and Vertical movements of the eyes) and another circuit to record
EEG signals. By utilizing a scaling and subtraction circuit, they successfully
removed Horizontal EOG (HEOG) and Vertical EOG (VEOG) from the raw EEG.
However, this method relied on the availability of hardware and separate EOG
electrodes, making it limited in its applicability. Nonetheless, it served as
the foundation for subsequent advancements in the field of EOG artifact
removal.

Nearly 20 years later, in 2004, He et al. proposed an adaptive filtering-based
method for EOG removal [21]. Their approach required an EOG electrode, which,
when combined with another EEG electrode, facilitated the construction of an
adaptive filter with a filter order of M=3. While effective, the main
limitation of this method was its dependence on the presence of a separate EOG
recording.

In 2014, Hu et al. introduced a novel approach combining an adaptive neural
fuzzy inference system (ANFIS) and a functional link neural network (FLNN) to
remove EOG and electromyogram (EMG) artifacts from EEG signals [22]. The
method involved an adaptive filtering algorithm that adjusted the parameters
of the fuzzy inference and neural network. Although successful in artifact
removal, this technique still relied on the availability of raw artifact data
in the early stages of the removal process.

Another approach proposed in 2014 by Maddirala et al. utilized singular
spectrum analysis (SSA) and adaptive noise canceler (ANC) to remove EOG
artifacts from contaminated EEG signals [23]. The technique involved grouping
the SSA components to construct an EOG reference signal for ANC. Using this
estimated reference signal, an adaptive filter was employed to remove the EOG
artifact. Performance evaluation using RRMSE (relative root mean square error)
and MAE (Mean Absolute Error) metrics demonstrated the effectiveness of the
proposed algorithms. Notably, this method only required EEG electrodes and did
not rely on the presence of dedicated EOG channels, enhancing its reliability
in situations where EOG electrodes were not available.

In 2020, Noorbasha et al. presented a research study introducing the overlap
segmented adaptive singular spectrum analysis (Ov-ASSA) combined with adaptive
noise canceler (ANC) technique for EOG artifact removal [24]. In this
approach, the first one or two reconstructed components of the Ov-SSA
technique were adaptively grouped and used as a reference EOG signal for ANC
in a single-channel EEG recording system. The algorithm’s performance was
evaluated using simulated and real data, with RRMSE and MAE metrics employed
as performance measures.

###  2.2 Multi Channel

In 2000, Jung et al. presented a method for removing various artifacts from
EEG records using BSS through ICA [25]. Their study demonstrated the
effectiveness of ICA in detecting, separating, and removing contamination from
different artifactual sources in EEG records. However, their method relied on
the availability of EOG electrodes during EEG recording, which is a
significant limitation. Nonetheless, their contribution in the field of
artifact removal has been influential and serves as inspiration for our
proposed method.

Shahabi et al. introduced a method for removing eye blink artifacts from EEG
signals using EOG reference electrodes [26]. They employed an autoregressive
(AR) process to model the EEG activity and an output-error model for eye
blinks. By utilizing a Kalman filter, they estimated the actual EEG by merging
the two models. The performance of their method was evaluated using the BCI
competition 2008, dataset II-a, which is also one of the datasets used for
performance evaluation in our study.

Nguyen et al. published a framework in 2012 that combined artificial neural
networks and wavelet transformation for EEG artifact removal [1]. Their
algorithm utilized the universal approximation characteristics of neural
networks and the time-frequency properties of the wavelet transform. The
neural network was trained on a simulated dataset with known ground truths.
Notably, their framework did not require EOG electrodes during the artifact
removal procedure, distinguishing it from many other EEG artifact removal
algorithms. They also compared their results with the ICA method using a time-
consuming procedure for artifact rejection on simulated and real datasets.

Zeng et al. proposed a mixed-method approach in which they used singular
spectrum analysis (SSA) as a blind source separation technique to extract
different components, including artifactual and neural components, from EEG
signals [18]. They then employed empirical mode decomposition (EMD) to denoise
the components affected by EOG artifacts. The artifactual components were
projected back to the electrode space and subtracted from the EEG signals to
obtain clean EEG. The experimental results on artificially contaminated EEG
data and publicly available real EEG data were reported. However, their method
still required EOG recording electrodes during EEG recording.

In 2015, Yang et al. proposed the ICA-based multivariate empirical mode
decomposition (IMEMD) method for removing EOG artifacts from multichannel EEG
signals [27]. They decomposed the EEG signals using multivariate empirical
mode decomposition (MEMD) into multiple multivariate intrinsic mode functions
(MIMFs). The artifactual components were extracted by reconstructing the MIMFs
corresponding to EOG artifacts. By performing ICA on the contaminated signals,
the EOG-related independent components (ICs) were identified and removed.
Finally, the denoised EEG signals were reconstructed by applying the inverse
transform of ICA and MEMD. They evaluated the performance of their method
using signal-to-noise ratio (SNR) and mean squared error (MSE). However, this
technique still required EOG electrodes during EEG recording.

###  2.3 Both single channel and multi channel

In 2016, Yang et al. proposed a method that combines a sparse autoencoder
(SAE) with recursive least square adaptive (RLS) filtering to remove EOG
artifacts without the need for reference electrodes [14]. During the offline
step, the SAE model learns information from EOG signals, and in the online
stage, the trained SAE model is used to extract preliminary EOG artifacts from
the raw EEG signal. The recursive least square adaptive filter then uses the
identified EOG artifacts as a reference signal to remove interference without
parallel EOG recordings. Their method was evaluated using a classification
accuracy metric.

In 2018, Yang et al. published another research in which they utilized a deep
learning network (DLN) to remove EOG artifacts from EEG signals [4]. The DLN
consisted of a multi-layer perceptron (MLP) with an SAE-based approach to
estimate non-artifactual EEG from the raw EEG signal. The proposed method
involved two stages. In the offline stage, training samples without artifacts
were used to train the DLN, enabling it to reconstruct denoised EEG signals
and learn the high-order statistical moments of the EEG. In the online stage,
the learned DLN was applied as a filter to remove EOG artifacts from
contaminated EEG signals.

Considering the advantages and limitations of previous methods, our study aims
to present a method for EOG artifact removal without requiring an EOG
reference electrode during the procedure, applicable to both single-channel
and multi-channel EEG signals. Additionally, we aim to estimate both
horizontal (HEOG) and vertical (VEOG) EOG signals and remove these artifactual
components from the EEG signal. Our approach is based on long short-term
memory (LSTM) neural networks and ICA. The study involves offline and online
stages. In the offline stage, the LSTM network learns features of EOG signals
from EEG recordings after signal normalization. In the online stage, the EEG
signals are processed by the LSTM network to estimate the EOG signal, which is
then used as the EOG reference electrode. By performing ICA, all components,
including artifactual and non-artifactual ones, are extracted. Components
similar to the estimated EOG signal are identified as artifactual and removed
from the sources. Finally, the denoised EEG signal is estimated by back-
projecting to the electrode subspace.

Throughout this paper, vectors are represented by bold small letters, e.g.,
$\mathbf{x}$, matrices by capital bolded letters, e.g., $\mathbf{X}$, and
scalars by normal italic fonts, e.g., $x$. Element-wise multiplication is
denoted by ”$\circ$”, matrix multiplication by ”$\times$”, and concatenation
by ”$[,]$”. The $l^{2}$-norm of a vector is denoted as ”$|\lvert.|\lvert$”,
and the transpose and inverse of a matrix are denoted by superscript $\top$
and $-1$, respectively.

The remainder of the paper is organized as follows: Section 3 introduces the
proposed method and the employed dataset and provides a detailed description
of the methodology. Section 4 presents the experimental setup and results.
Finally, Sections 5 and 6 conclude the paper with a discussion of the
findings.

##  3 Method

###  3.1 LSTM

LSTM networks were first introduced by Hochreiter et al. in 1997 as a specific
type of recurrent neural network capable of learning long-term dependencies
[28]. These networks were designed to address the challenge of capturing and
retaining information over extended periods. Notably, LSTM networks excel at
memorizing information for prolonged durations by leveraging their unique
architecture [29, 30]. The structure of LSTM networks enables them to learn
and retain information effectively, which distinguishes them from other
recurrent neural networks. In general, recurrent neural networks consist of a
repetitive sequence of simple neural network modules or units.

The core component of LSTM networks is the LSTM cell, illustrated in Figure 1.
The annotations used in the figure and in this section are inspired by [31].

Figure 1: Diagram of an LSTM cell.

In this context, $\mathbf{x}_{t}$ represents the input vector at the current
time step, while $\mathbf{s}_{t}$ and $\mathbf{h}_{t}$ are the cell state and
the output of the hidden layer at time $t$, respectively. Both
$\mathbf{s}_{t}$ and $\mathbf{h}_{t}$ are initialized as zeros (i.e.,
$\mathbf{s}_{0}=\mathbf{h}_{0}=\mathbf{0}$). The LSTM cell computes
$\mathbf{s}_{t}$ and $\mathbf{h}_{t}$ as follows:

$\displaystyle\mathbf{s}_{t}$
$\displaystyle=\mathbf{s}_{t-1}\circ\mathbf{f}_{t}+\mathbf{\tilde{i}}_{t},$
(1) $\displaystyle\mathbf{h}_{t}$
$\displaystyle=\mathbf{o}_{t}\circ\textrm{tanh}(\mathbf{s}_{t}),$ (2)

where $\mathbf{f}_{t}$, $\mathbf{i}_{t}$, and $\mathbf{o}_{t}$ denote the
outputs of the forget, input, and output gates, respectively.

The forget gate determines the amount of information from $\mathbf{s}_{t-1}$
that should be discarded. Its output is computed as:

$\mathbf{f}_{t}=\sigma(\mathbf{W}_{f}\times[\mathbf{h}_{t-1},\mathbf{x}_{t}]+\mathbf{b}_{f}),$
(3)

where $\mathbf{W}_{f}$ and $\mathbf{b}_{f}$ are trainable coefficients, and
$\sigma(\cdot)$ represents the sigmoid function defined as:

$\sigma(z)=\frac{1}{1+e^{-z}}.$ (4)

The input gate determines the amount of information from $\mathbf{x}_{t}$ and
$\mathbf{h}_{t-1}$ that should be retained in the current cell state
$\mathbf{s}_{t}$. Its output is computed as:

$\mathbf{i}_{t}=\mathbf{\tilde{i}}_{t}\circ\mathbf{\tilde{s}}_{t},$ (5)

where

$\displaystyle\mathbf{\tilde{i}}_{t}$
$\displaystyle=\sigma(\mathbf{W}_{i}\times[\mathbf{h}_{t-1},\mathbf{x}_{t}]+\mathbf{b}_{i}),$
(6) $\displaystyle\mathbf{\tilde{s}}_{t}$
$\displaystyle=\textrm{tanh}(\mathbf{W}_{s}\times[\mathbf{h}_{t-1},\mathbf{x}_{t}]+\mathbf{b}_{s}),$
(7)

with $\mathbf{W}_{i}$, $\mathbf{W}_{s}$, $\mathbf{b}_{i}$, and
$\mathbf{b}_{s}$ representing weighting matrices and biases, respectively. The
function $\textrm{tanh}(\cdot)$ denotes the hyperbolic tangent function
defined as:

$\textrm{tanh}(z)=\frac{e^{z}-e^{-z}}{e^{z}+e^{-z}}.$ (8)

The output gate determines the amount of information from $\mathbf{s}_{t}$
that should be included in the output of the hidden layer $\mathbf{h}_{t}$.
Its output is computed as:

$\mathbf{o}_{t}=\sigma(\mathbf{W}_{o}\times[\mathbf{h}_{t-1},\mathbf{x}_{t}]+\mathbf{b}_{o}),$
(9)

where $\mathbf{W}_{o}$ and $\mathbf{b}_{o}$ represent the corresponding
weighting matrix and bias vector.

In this study, we employ a deep LSTM network consisting of four layers, each
with 64 hidden units and accompanied by dropout rates of 0.1, 0.3, 0.3, and
0.1, respectively. The last layer is connected to a fully connected network
with an output size of 2, representing the VEOG and HEOG signals.

###  3.2 ICA

EEG signals recorded from scalp electrodes represent a combination of neural
activations occurring in different areas of the brain [32, 33]. The scalp
electrode’s location determines the weights assigned to each neural activity.
The primary objective is to extract the dominant activity from these combined
activations to gain insights into the underlying brain sources. While invasive
recording allows for direct measurement of specific neural activity, it
requires surgery and is generally not feasible or desirable. Instead, scalp
EEG recordings provide a convenient means to capture neural activity, albeit
with a mixture of different activations weighted differently across electrode
channels.

If we had knowledge of these weights, we could compute the potentials in the
brain sources using a sufficient number of electrodes [32]. However, the exact
weights and the contribution of each neural activity are typically unknown.
Consequently, we rely on BSS methods to overcome this challenge.

ICA is a widely recognized BSS algorithm [34] and aims to extract sources by
maximizing non-gaussianity and minimizing mutual information [35]. It assumes
that the observed signals are linear mixtures of these independent sources.
Various implementations of ICA, such as Infomax, SOBI, JADE, FastICA, and
kernel-independent component analysis, have been proposed [36]. ICA is well-
suited for EEG signals due to two key assumptions about the recordings:

  1. 1.

The observations from different electrodes on the scalp correspond to neural
activations from distinct brain areas, assuming independence between areas.

  2. 2.

The generation and recording speed of EEG signals are much slower than the
speed of electromagnetic wave propagation, resulting in negligible delays
between the origin of neural activity and the recording electrodes.

These assumptions allow ICA to effectively separate the contributions of
different sources from scalp EEG recordings, providing valuable insights into
the underlying neural activity.

Consider $\mathbf{x}j$ as the vector signal of the j-th channel, and
$\mathbf{a}_{1}$ to $\mathbf{a}_{N_{s}}$ as the projection vectors.
Additionally, let $\mathbf{s}_{1}$ to $\mathbf{s}_{N_{s}}$ represent the
vector sources to be estimated. The number of channels, sources, and time
samples are denoted as $N_{c}$, $N_{s}$, and $T$, respectively. The overall
problem can be expressed as follows (Equation (10)):

$\mathbf{x}_{j}=a_{j,1}\mathbf{s}_{1}+\cdots+a_{j,N_{s}}\mathbf{s}_{N_{s}}$
(10)

To generalize Equation (10) to all channels and sources using the projection
vectors, we can write Equation (11). Each vector is defined as follows:
$\mathbf{x}_{j}=[x_{j,1},\cdots,x_{j,T}],j=1,...,N_{c}$,
$\mathbf{a}_{i}=[a_{i,1},\cdots,a_{i,N_{c}}]^{T},i=1,...,N_{s}$, and
$\mathbf{s}_{i}=[s_{i,1},\cdots,s_{i,T}],i=1,...,N_{s}$.

$\begin{bmatrix}\mathbf{x}_{1}\\\ \vdots\\\ \mathbf{x}_{N_{c}}\\\
\end{bmatrix}=\begin{bmatrix}\mathbf{a}_{1}&\cdots&\mathbf{a}_{N_{s}}\end{bmatrix}\times\begin{bmatrix}\mathbf{s}_{1}\\\
\vdots\\\ \mathbf{s}_{N_{s}}\\\
\end{bmatrix}\qquad\Rightarrow\qquad\mathbf{X}=\mathbf{A}\times\mathbf{S}=\sum_{i=1}^{N_{s}}\mathbf{a}_{i}^{T}\mathbf{s}_{i}$
(11)

The above equation represents the problem of finding the sources matrix
($\mathbf{S}$) and the mixing matrix ($\mathbf{A}$). To reduce the dependency
on the center and variance of the observations, the observed data
($\mathbf{X}$) should be whitened. The centering and whitening procedure can
be performed using Algorithm 1.

Algorithm 1 Centering and Whitening algorithm

0: $\mathbf{X}$ as the observation

0: $\mathbf{\tilde{X}}$ as the whitened matrix

1: $\mathbf{\hat{X}}\leftarrow\mathbf{X}-\mathbb{E}\\{\mathbf{X}\\}$

2: Eigenvalue decomposition of the covariance matrix of $\mathbf{\hat{X}}$
$\rightarrow$
$\mathbf{V}\mathbf{D}\mathbf{V}^{T}=\mathbb{E}\\{\mathbf{\hat{X}}\mathbf{\hat{X}}^{T}\\}$

3: $\mathbf{P}\leftarrow\mathbf{D}^{-\frac{1}{2}}\mathbf{V}^{T}$

4: $\mathbf{\tilde{X}}\leftarrow\mathbf{P}\mathbf{\hat{X}}$

5: return $\mathbf{\tilde{X}}$

After performing the whitening procedure, the next step is to solve the
problem of finding the sources and mixing matrix arrays. In this study, we
utilized FastICA, which is based on a fixed-point iteration scheme that aims
to maximize the non-Gaussianity of $\mathbf{W}^{T}\mathbf{X}$ [34]. The
iteration process is achieved using the Newton iterative method [37]. To
measure non-Gaussianity, FastICA employs a non-quadratic nonlinear function
$f(u)$, its first derivative $g(u)$, and its second derivative
$g^{\prime}(u)$. These functions are defined in Equation (12). The algorithm
for extracting multiple components using the FastICA method is outlined in
Algorithm 2, where $\mathbf{1}$ represents the column vector of size $T\times
1$ filled with the value 1.

$f(u)=-e^{-\frac{u^{2}}{2}}\xrightarrow{\frac{\partial}{\partial
u}}g(u)=ue^{-\frac{u^{2}}{2}}\xrightarrow{\frac{\partial}{\partial
u}}g^{\prime}(u)=(1-u^{2})e^{-\frac{u^{2}}{2}}$ (12)

Algorithm 2 FastICA algorithm

0: $\mathbf{\tilde{X}}\in{\rm I\\!R}^{N_{c}\times T}$ as the whitened
observation, and $N_{s}$ which is the number of proposed components ($\leq
N_{s}$)

0: $\mathbf{W}\in{\rm I\\!R}^{N_{c}\times N_{s}}$ as the unmixing matrix,
$\mathbf{S}\in{\rm I\\!R}^{N_{s}\times T}$ independent component matrix

1: Initialize $\mathbf{W}$ with random values
($\sim\mathcal{N}(\mu=0,\,\sigma^{2}=1)$)

2: for $i$ in $1$ to $N_{s}$ do

3: while $\mathbf{w}_{i}$ changes do

4:
$\mathbf{w}_{i}\leftarrow\frac{1}{T}\mathbf{\tilde{X}}g(\mathbf{w}_{i}^{T}\mathbf{\tilde{X}})^{T}-\frac{1}{T}g^{\prime}(\mathbf{w}_{i}^{T}\mathbf{\tilde{X}})\mathbf{1}\mathbf{w}_{i}$

5:
$\mathbf{w}_{i}\leftarrow\mathbf{w}_{i}-{\sum}_{j=1}^{i-1}(\mathbf{w}_{i}^{T}\mathbf{w}_{j})\mathbf{w}_{j}$

6:
$\mathbf{w}_{i}\leftarrow\frac{\mathbf{w}_{i}}{\left\lVert\mathbf{w}_{i}\right\rVert}$

7: end while

8: end for

9:
$\mathbf{W}\leftarrow[\mathbf{w}_{1},\mathbf{w}_{2},\cdots,\mathbf{w}_{N_{s}}]$

10: $\mathbf{S}\leftarrow\mathbf{W}^{T}\mathbf{\tilde{X}}$

11: return $\mathbf{W},\mathbf{S}$

###  3.3 Combining LSTM and ICA

This paper introduces two different approaches: a single-channel approach and
a multi-channel approach. Each approach consists of two stages: an offline
stage and an online stage. In the offline stage, a Deep LSTM network is
trained to learn the underlying time series, sequences, and features in both
contaminated EEG signals and EOG recordings. In the online stage, the
estimated EOGs are treated as external EOG channels and are combined with
other EEG recording electrodes. Subsequently, ICA is applied to extract both
clean and artifactual sources. Finally, the artifactual sources are removed
during back projection, resulting in a cleaned recording. The estimated EOG
recordings are also provided as part of the outputs.

####  3.3.1 Training LSTM network in offline stage

To train the proposed neural network, the number of input and output channels
needs to be specified. Since our network is designed for both single-channel
and multi-channel EOG artifact removal, we consider both single and multiple
inputs to account for the removal of EOG artifacts in each dataset.

Our framework is also flexible in terms of the number of input samples. This
means that when we feed the network with varying sample sizes, the output can
be generated without any limitations on data segmentation or other sampling
techniques. During the training scheme, we used a batch size of 250 samples
per iteration at each epoch. However, the segment size is not fixed, and
during testing, we can feed the network with different input sample sizes.

For the loss function during the training of the Deep LSTM network, we
utilized Mean Squared Error (MSE). This loss function, as formulated in
Equation (13), is chosen due to its differentiability and sensitivity to
outliers.

$MSE=\frac{1}{N}\sum_{i=1}^{N}(\hat{y}_{i}-y_{i})^{2}$ (13)

To optimize the Deep LSTM network weights and parameters, we employed the Adam
algorithm [38]. This algorithm utilizes first-order gradient optimization
based on adaptive momentum estimation. It is easy to implement,
computationally efficient, requires minimal memory, is invariant to diagonal
rescaling of gradients, and is well-suited for problems involving large
datasets and/or parameters [38].

We trained the network for 50 epochs and implemented an early stopping scheme
with a patience of 2 epochs based on validation loss to avoid over-fitting on
training data [39].

Once the training is completed, the network is ready to estimate EOG
recordings (with 2 channels for VEOG and HEOG) and the next step involves
using the ICA algorithm to extract non-artifactual sources and reconstruct the
clean EEG recordings.

###  3.4 Investigating ICA for EOG artifact removal in online stage

The estimated EOG recordings are combined with the normalized EEG recordings
(without EOG references) to form the EOG electrodes. The resulting observation
is then subjected to whitening, followed by the application of the FastICA
algorithm to extract sources ($\mathbf{S}$) and the unmixing matrix
($\mathbf{W}$). At this stage, it is necessary to identify and remove EOG
artifactual sources in subsequent steps. To accomplish this, we compute the
absolute correlation index (due to the scale uncertainty in the ICA method)
between each source and the estimated EOG recordings from the previous step.
Sources with a high correlation (above 0.8) are considered potential
artifactual sources, and their corresponding columns in the mixing matrix
($\mathbf{A}$) are set to zero to eliminate their effect during the
reconstruction stage. The complete procedure is outlined in Algorithm 3. It is
important to note that while our formulation assumes two EOG estimations
(generated by the Deep LSTM network, mainly related to VEOG and HEOG), our
framework can accommodate different number of EOG channels. The ”;” command
represents vertical concatenation, and corrcoeff refers to the temporal
correlation measurement function.

Algorithm 3 Application of ICA in EOG artifact removal

0: $\mathbf{X}_{N}$ as the normalized EEG recording, $\mathbf{\hat{Y}}$ as the
estimated EOG from Deep LSTM network, and $N_{c}$, $N_{E}$ as the number of
EEG and EOG channels, respectively

0: $\mathbf{X}_{r}$ as the reconstructed EEG

1: $\mathbf{Z}\leftarrow[\mathbf{X_{N}};\mathbf{\hat{Y}}]$

2: $\mathbf{Z}\leftarrow\texttt{Whitening}(\mathbf{Z})$ using Algorithm 1

3: $N\leftarrow N_{c}+N_{E}$

4: $[\mathbf{S},\mathbf{W}]\leftarrow\texttt{FastICA}(\mathbf{Z},N)$ using
Algorithm 2

5: for $i$ in $1$ to $N$ do

6:
$corr_{1}=|\texttt{corrcoeff}\left(\mathbf{\hat{Y}}(1,:),\mathbf{S}(i,:)\right)|$

7:
$corr_{2}=|\texttt{corrcoeff}\left(\mathbf{\hat{Y}}(2,:),\mathbf{S}(i,:)\right)|$

8: end for

9: $id_{1}=\texttt{find}\left(corr_{1}\geq 0.8\right)$

10: $id_{2}=\texttt{find}\left(corr_{2}\geq 0.8\right)$

11: $\mathbf{A}\leftarrow\mathbf{W}^{-1}$

12: $\mathbf{A}(:,id_{1})\leftarrow\mathbf{0}$

13: $\mathbf{A}(:,id_{2})\leftarrow\mathbf{0}$

14: $\mathbf{X_{r}}\leftarrow\mathbf{A}\mathbf{S}$

15: $\mathbf{X}_{r}\leftarrow\mathbf{X}_{r}(1:N_{c},:)$

16: return $\mathbf{X_{r}}$

To remove EOG artifacts from single-channel EEG recordings, Algorithm 3 can
still be applied with a slight modification in the notation. Instead of
representing the observation as a matrix $\mathbf{X}$, it is now represented
as a vector $\mathbf{x}$. The overall process of model development is
illustrated in Figure 2, showcasing the different steps involved in EOG
artifact removal.

Figure 2: The flowchart of the proposed model.

###  3.5 Evaluation metrics

The performance of the EOG estimation and the effectiveness of the algorithm
in reconstructing clean EEG signal are evaluated using three metrics: Mean
Squared Error (MSE), Mean Error (ME), and Mean Absolute Error (MAE). Assuming
the non-contaminated original EEG signal is denoted as
$\mathbf{y}\in\mathbb{R}^{N}$ and the estimated artifact-free signal is
$\hat{\mathbf{y}}$, the metrics are defined as follows:

$\displaystyle MSE$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}i)^{2}$ (14)
$\displaystyle MAE$ $\displaystyle=\frac{1}{N}\sum{i=1}^{N}\lvert
y_{i}-\hat{y}i\rvert$ (15) $\displaystyle ME$
$\displaystyle=\frac{1}{N}\sum{i=1}^{N}(y_{i}-\hat{y}_{i}).$ (16)

###  3.6 Dataset

One of the crucial aspects in developing an artifact removal algorithm is the
availability of suitable data for evaluating its effectiveness. Since real
datasets often lack a direct benchmark for artifact removal and only provide
measurements of accuracy or performance in specific tasks, it becomes
necessary to generate synthetic data that includes both EOG recordings and
pure EEG recordings. In this project, we utilized a semi-simulated EEG/EOG
dataset consisting of fifty-four datasets recorded from twenty-seven healthy
subjects (age: $27.17\pm 5.2$).

Each dataset includes nineteen EEG electrodes (FP1, FP2, F3, F4, C3, C4, P3,
P4, O1, O2, F7, F8, T3, T4, T5, T6, Fz, Cz, Pz) placed according to the 10-20
International System (as shown in Fig.3), along with four EOG electrodes used
to record VEOG and HEOG eye movements in a bipolar scheme [40].

The recorded datasets consist of 30 seconds of closed eyes without any eye
movements to capture clean and pure EEG signals. The signals were sampled at a
frequency of 200 Hz and filtered using a bandpass filter (0.5 to 40 Hz). In
parallel, EOG recordings were acquired from the same subjects and during the
same time duration. These EOG signals were then filtered using a 0.5 to 5 Hz
bandpass filter and added to the pure EEG signals using Equation (17) to
contaminate them.

The parameters $a$ and $b$ in Equation (17) represent subject- and dataset-
specific coefficients that were computed using linear regression to introduce
blink and eye movement artifacts into the pure EOG signals.

$\mathbf{X}_{\textrm{Contaminated}}=\mathbf{X}_{\textrm{Pure}}+a\mathbf{X}_{\textrm{VEOG}}+b\mathbf{X}_{\textrm{HEOG}}$
(17)

Figure 3: Channel location for the dataset used in this study.

###  3.7 Pre-processing

To mitigate the impact of varying signal values and potential outliers in EEG
and EOG signals during the optimization algorithm’s training phase, it is
crucial to normalize these signals. This normalization involves applying a
zero mean unit variance normalization to both the EEG and EOG channels.
Additionally, the means and standard deviations of the signals are stored in
vectors for reconstruction purposes after network training. The normalization
procedure is outlined in Algorithm 4.

Algorithm 4 Normalization Algorithm

0: $\mathbf{X}\in{\rm I\\!R}^{N_{c}\times T}$ as input EEG observation,
$\mathbf{Y}\in{\rm I\\!R}^{N_{E}\times T}$ as the output EOG recording

0: $\mathbf{X}_{N},\mathbf{Y}_{N}$ as the normalized signals,
$\bm{\mu},\bm{\sigma}$ as the mean, and standard deviation parameters

1: $\bm{\mu}\leftarrow[\hskip 4.26773pt]$

2: $\bm{\sigma}\leftarrow[\hskip 4.26773pt]$

3: for $i$ in $1$ to $N_{c}$ do

4: $\mu_{x}\leftarrow\frac{1}{T}\sum_{j=1}^{T}\mathbf{X}(i,j)$

5:
$\sigma_{x}\leftarrow\sqrt{\frac{\sum_{j=1}^{T}(\mathbf{X}(i,j)-\mu_{x})^{2}}{T}}$

6: $\mathbf{X}_{N}(i,:)\leftarrow\frac{\mathbf{X}(i,:)-\mu_{x}}{\sigma_{x}}$

7: $\bm{\mu}\leftarrow[\bm{\mu},\mu_{x}]$

8: $\bm{\sigma}\leftarrow[\bm{\sigma},\sigma_{x}]$

9: end for

10: for $i$ in $1$ to $N_{E}$ do

11: $\mu_{y}\leftarrow\frac{1}{T}\sum_{j=1}^{T}\mathbf{Y}(i,j)$

12:
$\sigma_{y}\leftarrow\sqrt{\frac{\sum_{j=1}^{T}(\mathbf{Y}(i,j)-\mu_{y})^{2}}{T}}$

13: $\mathbf{Y}_{N}(i,:)\leftarrow\frac{\mathbf{Y}(i,:)-\mu_{y}}{\sigma_{y}}$

14: $\bm{\mu}\leftarrow[\bm{\mu},\mu_{y}]$

15: $\bm{\sigma}\leftarrow[\bm{\sigma},\sigma_{y}]$

16: end for

17: return $\mathbf{X}_{N},\mathbf{Y}_{N},\bm{\mu},\bm{\sigma}$

##  4 Results

In this section, various metrics are evaluated for single and multi channel
EOG artifact removal from EEG recordings.

The entire procedure is divided into three stages: 1) Training the Deep LSTM
network using contaminated EEG recordings and EOG reference channels with
training segment samples. 2) Estimating EOG recordings for the testing data.
3) Combining EEG recordings with the estimated EOG activity, performing ICA,
and reconstructing the artifact-free recording. Since both clean and
contaminated recordings are provided in the dataset, we can measure the
performance of both EOG estimation and EEG denoising using the proposed
algorithm and compare it with previous literature. The overall network
architecture for this dataset consists of nineteen input channels and two
output channels for EOG estimation. The MSE loss values during the 50 epochs
of training are depicted in Figure 4. It is important to note that training
was stopped early due to early stopping criteria with a patience of two
epochs.

Figure 4: Changes in MSE Loss During Training Scheme. The red plot represents
the validation set, while the blue plot corresponds to the training data.
Epoch 0 denotes the first epoch.

At this stage, we can obtain the EOG estimations from the trained Deep LSTM
network. To evaluate the performance, we randomly selected 30 % of the semi-
simulated datasets using a cross-subject paradigm. This ensures a fair
evaluation, as training on segments from one subject and testing on segments
from the same subject could introduce bias. The remaining test data was
divided into 250 different segments with variable lengths and overlapping
samples. The results are presented in Figure 5 and Table 1. Additionally, for
a visual comparison, a sample output of the Deep LSTM network for both HEOG
and VEOG estimation is shown in Figure 6, alongside the original EOG
activation in the dataset.

Figure 5: The performance of VEOG and HEOG estimation over 250 iterations. Table 1: Quantified EOG estimation performance on test set over 250 iterations. EOG channel | Error (Mean $\pm$ STD)
---|---
VEOG | $0.051\pm 0.011$
HEOG | $0.040\pm 0.007$
Average | $\mathbf{0.046\pm 0.005}$

(a) HEOG

(b) VEOG

Figure 6: Comparison between estimated EOG and its original waveform in a
sample segment.

To perform the final analysis, we evaluate the cleaning of the artifactual EEG
recordings on the test data, which was not seen during training. The estimated
EOG signals are concatenated with the observed contaminated data, and then the
FastICA algorithm is applied to extract the artifactual sources and clean the
contaminated EEG signals. For each of the nineteen channels, we compare the
estimated clean EEG with the pure signals provided in the semi-simulated
dataset. The results for each channel are shown in Figure 7.

In order to visualize the time series data, Figure 8 presents the
reconstructed EEG signals for all nineteen channels alongside the
corresponding pure EEG signals. To facilitate comparison and eliminate the
effect of scale uncertainty, the signals are min-max normalized between 0 and
1 in each channel before plotting.

Figure 7: Multi-channel EEG de-artifactualization performance with proposed
method. The arrows represent the range of standard deviation. Figure 8: Sample
of Cleaned EEG Using the Proposed Method: Zoom-In for Details. The red curves
represent the original pure EEG, while the blue curves depict the output of
the proposed algorithm.

The results clearly indicate that the overall error between the reconstructed
clean EEG and the pure EEG is relatively low (approximately 0.05 MSE, 0.16
MAE, and 0.02 ME) across all channels in the main dataset. To further evaluate
the effectiveness of our proposed method, we compared it with two recent deep
learning-based approaches: SAE-RLS by [14], and DLN-SAE by [4]. The comparison
results are presented in Table 2, highlighting the superior performance of our
suggested method over the other deep learning-based approaches.

Table 2: Comparison of the proposed method with the existing similar literature. The reported values are average errors between estimated artifact-free EEG and the original EEG recordings.  Method | MSE | MAE | ME
---|---|---|---
SAE-RLS [14] | $0.15$ | $0.39$ | $\mathbf{0}$
DLN-SAE [4] | $0.08$ | $0.29$ | $0.01$
Proposed method | $\mathbf{0.05}$ | $\mathbf{0.16}$ | $0.02$
Figure 9: Detailed visualization of the performance of the proposed algorithm
in Fp1 channel.

##  5 Discussion

This paper presents a novel method for estimating EOG signals using a deep
LSTM network and subsequently removing EOG artifacts from contaminated EEG
recordings through ICA. The results in Figures 8 and 9 demonstrate that the
proposed algorithm performs effectively even in frontal channels such as Fp1,
successfully detecting EOG peaks in EEG recordings and eliminating their
effects using ICA. This is attributed to the strong performance of ICA in
removing EOG artifacts and the potential of LSTM networks in various EEG-
related tasks [31, 41, 42, 43], indicating their ability to learn underlying
relationships in EEG data.

The results in Table 2 indicate that the proposed method achieves superior
performance in terms of MSE and MAE compared to the DLN-SAE method [4], both
of which are metrics that are unaffected by the sign of the EOG amplitude.
However, for the ME metric, which is highly sensitive to the polarity of the
EOG amplitude, the DLN-SAE method demonstrates better performance. The
enhanced performance in terms of MSE and MAE can be attributed to the
utilization of a more intricate network architecture, specifically LSTM, in
contrast to the combination of MLP and SAE, as previously shown by [44].

While the proposed pipeline primarily focuses on EOG artifact removal from EEG
signals, it can also be employed for EOG estimation in applications such as
gaze detection using EEG signals when EOG recordings are not available.

However, it is important to acknowledge certain limitations of this approach.
Firstly, the fine-tuning and adjustment of the LSTM network were carried out
through trial and error, suggesting the need for comprehensive approaches to
optimize hyperparameters, including the number of hidden layers and LSTM
units, in different applications and datasets. Secondly, the availability of a
gold standard non-contaminated dataset, as demonstrated in [40], was crucial
for this study. Therefore, further investigations are required to ensure the
applicability of this approach in real-world scenarios. For instance, transfer
learning could be explored to train models on one dataset and evaluate their
performance on other datasets. Given the high subject-to-subject variability,
the generalizability of the approach remains a key question to be addressed in
future research. Lastly, while FastICA was utilized as the core method in the
ICA stage of this paper, other BSS techniques may outperform FastICA in
different applications. Exploring alternative BSS methods could be an avenue
for future directions in this study.

These limitations highlight important aspects for future research in order to
address them and further enhance the efficacy and applicability of the
proposed method.

##  6 Conclusion

In this paper, we introduced an LSTM-ICA methodology for effectively removing
EOG artifacts from EEG recordings in the absence of external EOG recordings.
The proposed approach was evaluated using a dataset that included both
contaminated and non-contaminated EEG recordings. The performance of the model
was assessed in two scenarios: single-channel and multi-channel EEG data. Our
methodology demonstrated superior performance compared to existing deep
learning-based approaches, emphasizing the potential of combining ICA and LSTM
in future EEG studies.

## Conflict of Interest

The authors have no conflicts of interest to disclose.

## Ethics statement

This research did not involve the collection of any data from human subjects.

## Funding statement

This research was not supported by any funding resources.

## References

  * [1] H.-A. T. Nguyen, J. Musson, F. Li, W. Wang, G. Zhang, R. Xu, C. Richey, T. Schnell, F. D. McKenzie, J. Li, Eog artifact removal using a wavelet neural network, Neurocomputing 97 (2012) 374–389.
  * [2] R. A. Ramadan, A. V. Vasilakos, Brain computer interface: control signals review, Neurocomputing 223 (2017) 26–44.
  * [3] Z. Qiu, J. Jin, H.-K. Lam, Y. Zhang, X. Wang, A. Cichocki, Improved sffs method for channel selection in motor imagery based bci, Neurocomputing 207 (2016) 519–527.
  * [4] B. Yang, K. Duan, C. Fan, C. Hu, J. Wang, Automatic ocular artifacts removal in eeg using deep learning, Biomedical Signal Processing and Control 43 (2018) 148–158.
  * [5] S. Ferdowsi, V. Abolghasemi, Multi layer spectral decomposition technique for erd estimation in eeg $\mu$ rhythms: An eeg–fmri study, Neurocomputing 275 (2018) 1836–1845.
  * [6] N. Cudlenco, N. Popescu, M. Leordeanu, Reading into the mind’s eye: Boosting automatic visual recognition with eeg signals, Neurocomputing 386 (2020) 281–292.
  * [7] R. Grech, T. Cassar, J. Muscat, K. P. Camilleri, S. G. Fabri, M. Zervakis, P. Xanthopoulos, V. Sakkalis, B. Vanrumste, Review on solving the inverse problem in eeg source analysis, Journal of neuroengineering and rehabilitation 5 (1) (2008) 25.
  * [8] R. J. Croft, R. J. Barry, Removal of ocular artifact from the eeg: a review, Neurophysiologie Clinique/Clinical Neurophysiology 30 (1) (2000) 5–19.
  * [9] V. Mehdizadehfar, F. Ghassemi, A. Fallah, H. Pouretemad, Eeg study of facial emotion recognition in the fathers of autistic children, Biomedical Signal Processing and Control 56 (2020) 101721.
  * [10] P. Shooshtari, G. Mohammadi, B. Molaee Ardekani, M. B. Shamsollahi, Removing ocular artifacts from eeg signals using adaptive filtering and armax modeling, in: Proceeding of World Academy of Science, Engineering and Technology, Vol. 11, 2006, pp. 277–280.
  * [11] J. W. Scott, P. E. Scott-Johnson, The electroolfactogram: a review of its history and uses, Microscopy research and technique 58 (3) (2002) 152–160.
  * [12] S. Karimi, B. Molaee-Ardekani, M. B. Shamsollahi, C. Leroy, P. Derambure, Automatic ocular correction in eeg recordings using maximum likelihood estimation, in: IEEE International Symposium on Signal Processing and Information Technology, IEEE, 2013, pp. 000164–000169.
  * [13] D. Hagemann, E. Naumann, The effects of ocular artifacts on (lateralized) broadband power in the eeg, Clinical Neurophysiology 112 (2) (2001) 215–231.
  * [14] B. Yang, K. Duan, T. Zhang, Removal of eog artifacts from eeg using a cascade of sparse autoencoder and recursive least squares adaptive filter, Neurocomputing 214 (2016) 1053–1060.
  * [15] M. Fatourechi, A. Bashashati, R. K. Ward, G. E. Birch, Emg and eog artifacts in brain computer interface systems: A survey, Clinical neurophysiology 118 (3) (2007) 480–494.
  * [16] H. Berger, Über das elektrenkephalogramm des menschen, Archiv für Psychiatrie und Nervenkrankheiten 94 (1) (1931) 16–60.
  * [17] K. Imaizumi, F. Takahashi, G. Yoshida, K. Ogawa, The origin of the electrooculogram (eog), in: VI th ISCERG Symposium, Erfurt, 1967, pp. 29–46.
  * [18] T. Elbert, W. Lutzenberger, B. Rockstroh, N. Birbaumer, Removal of ocular artifacts from the eeg—a biophysical approach to the eog, Electroencephalography and clinical neurophysiology 60 (5) (1985) 455–463.
  * [19] A. Gevins, C. L. Yeager, S. Diamond, J. Spire, G. Zeitlin, A. Gevins, Automated analysis of the electrical activity of the human brain (eeg): A progress report, Proceedings of the IEEE 63 (10) (1975) 1382–1399.
  * [20] D. Girton, J. Kamiya, A simple on-line technique for removing eye movement artifacts from the eeg, Electroencephalography and clinical Neurophysiology 34 (2) (1973) 212–216.
  * [21] P. He, G. Wilson, C. Russell, Removal of ocular artifacts from electro-encephalogram by adaptive filtering, Medical and biological engineering and computing 42 (3) (2004) 407–412.
  * [22] J. Hu, C.-s. Wang, M. Wu, Y.-x. Du, Y. He, J. She, Removal of eog and emg artifacts from eeg using combination of functional link neural network and adaptive neural fuzzy inference system, Neurocomputing 151 (2015) 278–287.
  * [23] A. K. Maddirala, R. A. Shaik, Removal of eog artifacts from single channel eeg signals using combined singular spectrum analysis and adaptive noise canceler, IEEE Sensors Journal 16 (23) (2016) 8279–8287.
  * [24] S. K. Noorbasha, G. F. Sudha, Removal of eog artifacts from single channel eeg–an efficient model combining overlap segmented assa and anc, Biomedical Signal Processing and Control 60 (2020) 101987.
  * [25] T.-P. Jung, S. Makeig, C. Humphries, T.-W. Lee, M. J. Mckeown, V. Iragui, T. J. Sejnowski, Removing electroencephalographic artifacts by blind source separation, Psychophysiology 37 (2) (2000) 163–178.
  * [26] H. Shahabi, S. Moghimi, H. Zamiri-Jafarian, Eeg eye blink artifact removal by eog modeling and kalman filter, in: 2012 5th International Conference on BioMedical Engineering and Informatics, IEEE, 2012, pp. 496–500.
  * [27] G. Wang, C. Teng, K. Li, Z. Zhang, X. Yan, The removal of eog artifacts from eeg signals using independent component analysis and multivariate empirical mode decomposition, IEEE journal of biomedical and health informatics 20 (5) (2015) 1301–1308.
  * [28] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural computation 9 (8) (1997) 1735–1780.
  * [29] A. Gensler, J. Henze, B. Sick, N. Raabe, Deep learning for solar power forecasting—an approach using autoencoder and lstm neural networks, in: 2016 IEEE international conference on systems, man, and cybernetics (SMC), IEEE, 2016, pp. 002858–002865.
  * [30] M. Hafezi, N. Montazeri, S. Saha, K. Zhu, B. Gavrilovic, A. Yadollahi, B. Taati, Sleep apnea severity estimation from tracheal movements using a deep learning model, IEEE Access 8 (2020) 22641–22649.
  * [31] P. Wang, A. Jiang, X. Liu, J. Shang, L. Zhang, Lstm-based eeg classification in motor imagery tasks, IEEE transactions on neural systems and rehabilitation engineering 26 (11) (2018) 2086–2095.
  * [32] M. Ungureanu, C. Bigan, R. Strungaru, V. Lazarescu, Independent component analysis applied in biomedical signal processing, Measurement Science Review 4 (2) (2004) 18.
  * [33] T.-P. Jung, C. Humphries, T.-W. Lee, S. Makeig, M. J. McKeown, V. Iragui, T. J. Sejnowski, Extended ica removes artifacts from electroencephalographic recordings, in: Advances in neural information processing systems, 1998, pp. 894–900.
  * [34] A. Hyvärinen, E. Oja, Independent component analysis: algorithms and applications, Neural networks 13 (4-5) (2000) 411–430.
  * [35] A. Taymourtash, F. Ghassemi, Independent component analysis of sparse-transformed eeg signals for adhd/normal adults’ classification, in: 2015 23rd Iranian Conference on Electrical Engineering, IEEE, 2015, pp. 151–155.
  * [36] V. Matic, W. Deburchgraeve, S. Van Huffel, Comparison of ica algorithms for ecg artifact removal from eeg signals, in: IEEE-EMBS Benelux Chapter Symposium, 2009\.
  * [37] A. Hyvarinen, Fast ica for noisy data using gaussian moments, in: 1999 IEEE International Symposium on Circuits and Systems (ISCAS), Vol. 5, IEEE, 1999, pp. 57–61.
  * [38] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [39] R. Caruana, S. Lawrence, C. L. Giles, Overfitting in neural nets: Backpropagation, conjugate gradient, and early stopping, in: Advances in neural information processing systems, 2001, pp. 402–408.
  * [40] M. A. Klados, P. D. Bamidis, A semi-simulated eeg/eog dataset for the comparison of eog artifact rejection techniques, Data in brief 8 (2016) 1004–1006.
  * [41] P. Nagabushanam, S. Thomas George, S. Radha, Eeg signal classification using lstm and improved neural network algorithms, Soft Computing 24 (2020) 9981–10003.
  * [42] G. Zhang, V. Davoodnia, A. Sepas-Moghaddam, Y. Zhang, A. Etemad, Classification of hand movements from eeg using a deep attention-based lstm network, IEEE Sensors Journal 20 (6) (2019) 3113–3122.
  * [43] C. Lee, J. An, Lstm-cnn model of drowsiness detection from multiple consciousness states acquired by eeg, Expert Systems with Applications 213 (2023) 119032.
  * [44] Q. Li, Y. Liu, Y. Shang, Q. Zhang, F. Yan, Deep sparse autoencoder and recursive neural network for eeg emotion recognition, Entropy 24 (9) (2022) 1187\.


Our method is capable of taking a 2D scan of a real wood specimen (left,
photograph) and then automatically matching corresponding parameters for a
procedural solid wood appearance model based on [LDHM16]. Such a model is
fully three-dimensional and allows realistic rendering of carved wooden solids
(right, Monte Carlo simulated render) based on the input photographs.

# Automatic inference of a anatomically meaningful solid wood texture from a
single photograph

Thomas K. Nindel1,20000-0002-8003-2672, Mohcen Hafidi10000-0002-1931-0493,
Tomáš Iser10000-0003-3013-8994, and Alexander Wilkie10000-0002-3536-6577
1Charles University, Faculty of Mathematics and Physics, Czech Republic
2Berufsakademie Sachsen, Germany Corresponding author

###### Abstract

Wood is a volumetric material with a very large appearance gamut that is
further enlarged by numerous finishing techniques. Computer graphics has made
considerable progress in creating sophisticated and flexible appearance models
that allow convincing renderings of wooden materials. However, these do not
yet allow fully automatic appearance matching to a concrete exemplar piece of
wood, and have to be fine-tuned by hand. More general appearance matching
strategies are incapable of reconstructing anatomically meaningful volumetric
information. This is essential for applications where the internal structure
of wood is significant, such as non-planar furniture parts machined from a
solid block of wood, translucent appearance of thin wooden layers, or in the
field of dendrochronology.

In this paper, we provide the two key ingredients for automatic matching of a
procedural wood appearance model to exemplar photographs: a good
initialization, built on detecting and modelling the ring structure, and a
phase-based loss function that allows to accurately recover growth ring
deformations and gives anatomically meaningful results. Our ring-detection
technique is based on curved Gabor filters, and robustly works for a
considerable range of wood types. <ccs2012> <concept>
<concept_id>10010147.10010371.10010372.10010376</concept_id>
<concept_desc>Computing methodologies Reflectance modeling</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10010147.10010371.10010382.10010383</concept_id>
<concept_desc>Computing methodologies Image processing</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10010405.10010432.10010437.10010438</concept_id>
<concept_desc>Applied computing Environmental sciences</concept_desc>
<concept_significance>300</concept_significance> </concept> </ccs2012>

[500]Computing methodologies Reflectance modeling [500]Computing methodologies
Image processing [300]Applied computing Environmental sciences

††volume: 1††issue: 1

##  1 Introduction

In _wood rendering_ and _wood appearance modeling_ , the current state of the
art is based either on highly precise BSSRDF measurements (BTFs, appearance
scans), or procedural models. Measurement-based approaches yield pixel-perfect
matches at the cost of high storage and memory-bandwidth requirements,
acquisition complexity, and inability to edit the data after acquisition. On
the other hand, procedural wood models allow artistic control and editing, but
they are difficult to match to given wood samples, both if an exact match to a
given piece of wood is needed, or only with regard to a general wood type.
Recent work on using optimisation to match procedural material models to
observations [SLH∗20] has been quite successful in a broad range of settings,
including wood. However, as figure 1 shows, current approaches still fail to
properly match a given wood sample _down to its internal 3D structure_ – a
feature that is needed to make wood grain wrap correctly around a solid 3D
object. Specifically, as figure 1 shows, extant techniques are capable of
generating an internal ring structure - but not one that really convincingly
matches the wood grain pattern seen on the top surface, like in the results we
show for our technique in figure 7.

Our main contribution, described in section 3, is that we propose a robust,
deterministic method for the automatic inference of a locally fit, procedural
3D material model for wood. Our approach allows one to obtain realistic solid
wood textures that can be carved to any geometry (such as the elephants in
Fig. Automatic inference of a anatomically meaningful solid wood texture from
a single photograph). Our approach is not based on machine learning and does
not rely on any training datasets.

Our work builds on the assumption that by accurately identifying the basic
ring structure and its deformations, we can satisfy not only the requirements
from wood rendering, that is to attain closely matching structural appearance
(Fig. 2, bottom left), but also the requirements of _dendrochronology_ for a
precisely aligned identification of ring boundaries (Fig. 2, bottom right). In
the dendrochronology field (section 2.3), core samples of living trees are
inspected for their growth rings variations, and the resulting data is
interesting for climate research, archaeology, and art history. Automatic
approaches have to be robust to the wide range of wood anatomical features, as
certain wood species (e.g., diffuse-porous hardwoods) are very difficult to
work with. Our approach performs well even in these settings, and in fact, we
used it in dendrochronological contexts as a means of verifying its
performance.

Figure 1: False-color volume rendering of a solid texture that was fit using a
differentiable material graph, (left and middle) compared to our approach. The
density of the volume is set to $\sigma_{t}=15$ for growth-lateness of greater
$0.9$, $\sigma_{t}=0.5$ otherwise, both with a homogeneous single-scattering
albedo $\alpha=(0.805,0.578,0.359)^{T}$ (RGB). The model in the first image
was optimized with our board orientation estimate, but without an intial guess
for the deformation field as a prior. The second image additionally uses our
initial deformation field. The third image shows a visualization of our
result; insets are the respective frontal, non-volume renderings. Notice how
in spite of the very good match in terms of frontal appearance, the left and
middle volumes show blurry ring boundaries. Figure 2: Method pipeline
overview. Our method starts with a real wood specimen (top left) that is
scanned or photographed from one side. The 2D image is then processed to
detect the wood growth rings (see also Fig. 3). This process is already useful
for automatic ring labeling in dendrochronology (bottom right), where the
input image is a 2D scan of a flattened drilled cylinder of a tree. For use in
solid wood rendering (bottom left), the 2D ring pattern has to be further
processed and fitted in a 3D space, since we are interested in parameters of a
solid 3D appearance model based on [LDHM16]. As only a 2D scan is available as
input, the 3D fitting can only be expected to approximate the actual specimen
in its local neighborhood. The final path-traced render is using the fitted
parameters and model on an arbitrary user-specified 3D mesh.

##  2 Related work

###  2.1 Modelling wood appearance

Wood is a complex material with an intricate volumetric structure. Hartmann et
al. [HKRFM17] dive into the complex process that controls the formation of
wood cells, whose variance is ultimately responsible for the emergence of the
characteristic growth patterns. From a technical standpoint, wood is a bio-
composite formed from three different kinds of molecules: cellulose,
hemicellulose, and lignin. Lignin is the main chemical responsible for the
wood color. It is preferentially found in areas of thick cell walls [LVS∗18],
which explains why late growth shows a higher color saturation than early
growth in coniferous species. Chemically replacing lignin with a transparent
polymer results in transparent wood, a material with interesting optical
properties [VCL∗18].

The basic approaches for modelling wood in photorealistic rendering can be
divided into _material scanning_ ones, which essentially treat the material as
a surface of certain optical properties, and _procedural_ , which typically
model the whole three-dimensional wood interior. Our method is based on
finding the corresponding parameters for a procedural model, such that a
scanned surface of a wood can be transformed into a fully procedural and
editable 3D model.

#### Material scanning and BTFs

The appearance of a specific specimen of wood can be measured and encapsulated
in various ways. Bidirectional Texture Functions (BTFs) [DVGNK99] store the
full 6-dimensional surface reflectance using, essentially, lookup tables.
Material scanners are available that capture diffuse reflectance, spatially
varying roughness and a normal map, which can be fed into an appropriate BRDF
implementation, such as the Disney Principled BRDF [BS12]. Both approaches can
yield impressive realism for the specimen they were applied to, but suffer
from immense storage requirements. Henzler et al. [HMR20] propose an
interesting data-driven approach to generate solid textures from 2D specimen,
wood being one of their use-cases. Their approach has the advantage of
providing an infinite solid texture, but represents the input texture only
qualitatively.

#### General solid noise functions

Procedural noise can be a versatile tool to capture the texture of a wide
range of materials, Perlin Noise [Per85] and Gabor Noise (e.g. [LLDD09]) being
prominent examples. Matching the noise function parameters to obtain a desired
input texture is also possible: For example, [GLLD12] show a method to match
the power spectrum of Gabor Noise to exemplars, while [BMM∗21] show an example
of fitting the parameters of a Perlin Noise even in the presence of
discontinuities. Since our goal is to closely match the ring structure, using
power-spectrum based noise would require additional and hard to control steps,
since matching the rings requires a precise match of both the noise signals’
amplitude and _phase_.

#### Domain-specific procedural models

A state of the art appearance model of wood is the work of Liu et al.
[LDHM16], on whose work we based our results (Fig. 2). Their model uses
various distortion functions to modulate a cylindrical coordinate system.
Together with a sawtooth-like function that describes the color saturation
changes in the earlywood / latewood transitions, they also provide expressions
for diffuse reflectance and the fibre direction. These can be used as base
functions to define inputs to a BRDF shader. By the explicit way they model
wood “age” on a sub-ring scale, their model can be seen as anatomically
informed. Once the model parameters are tuned, a process that requires manual
work by an expert, it is very expressive and can be used to model many species
of wood. We are directly building on this model, and our contribution is that
we are able to automatically recover its main parameters. We also recover a
color map that describes coarse-level variations as a function of ring growth
time.

Larsson et al. [LIY∗22] extend the expressive power of wood-specific
procedural solid textures by the integration of wood knots, based on a
skeletal description of the branches and a distance-field based formulation of
their respective area of influence to the growth rings formation.

###  2.2 Procedural appearance matching

Differentiable material graphs allow for the use of gradient-based
optimization to polish material parameters, after a good initial guess has
been found [GHYZ20, SLH∗20]. Since the success of appearance matching depends
on (and can be limited by) the expressive power of the underlying procedural
model, neural loss functions [AAL16] seem to be best suited in finding
plausible parameters for a wide range of material graphs, because they
evaluate appearance differences in qualitative, perceptual terms. This allows
these algorithms to converge to satisfactory results even if the
expressiveness of the underlying graph is limited. However, none of these
works provide a 3D solid texture that can be used to "carve" objects from or
accurately replicate the curvature of the growth rings present in the target.

[LP00] show a half-automatic system that can extract parameters for a
procedural 3D solid wood texture. From a thresholded input image, they are
able to estimate ring frequency and spatial orientation of a given wood plank,
based on image statistics. Ring variations are modelled statistically in terms
of turbulence intensity and frequency. Their orientation estimates are based
on statistical properties of the thresholded image, while ours is based on
accurately identified rings. Our method also accurately recovers the
deformation field.

###  2.3 Dendrochronology

Automated dendrochronological measurements rely on image processing techniques
that use either photos [WQZ∗10, FDBP16, FDBP17, Len20], or computed-tomography
(CT) data [MGSS∗21] as their input. Salient points of these methods is the
application of 2D image processing techniques to the problem of growth ring
identification. They achieve reasonable accuracy for coniferous species, but
their accuracy greatly degrades in the presence of pores.

Notable data driven approaches include [FD18], who use a U-Net that was
trained on manually labeled data, while [PAH∗22] use a Mask R-CNN
architecture. Neural networks seem to improve detection rates compared to the
approaches based on traditional 2D image processing especially for ring porous
wood, a category that our ring detection method shows on-par performance. Data
driven approaches create likelihood maps for ring boundaries, that are
subsequently thresholded to get exact locations. In contrast, our approach
directly pinpoints the ring boundaries without the ambiguity of a likelihood.

Martinez et al. [MGSS∗21] extend the idea of tree ring extraction to 3D data
and reconstruct tree ring isosurfaces from X-ray computed tomography data
based on edge detection. Basing the method on CT data makes the approach very
robust, because of strong correlation between X-ray intensities and material
density. The limited availability of CT scanners hinders the practical
applicability of their approach.

In a more broader context, systems for the identification of tree species from
either microscopic or macroscopic images build on segmenting anatomical
features of wood, including its rings. Martins [Mar18] contains a survey of
the proposed methods. Datasets of dendrochronological measurements are
available through the work of [FDBP17], the WIAD database [RSB∗21], both
containing labels. The DendroElevator [oM22] database accumulates data from
several sources, containing some very detailed tree ring images (up to
microscopic scale) obtained from wood core samples. Some of their data are
labelled.

###  2.4 Finding rings in images

Another field that considers the enhancement, identification, and segmentation
of ring-like structures is _fingerprint identification_. The key differences
to our domain are the bandwidth of ring frequencies (fingerprint ridges are
quite evenly spaced), and the presence of branching (fingerprint ridges can
split and merge, which usually is not the case for tree rings). This bandwidth
limitation makes transferring techniques that are informed by or based on
frequency-domain information (e.g., [THG16, LWL∗20]) somewhat more
problematic, but not impossible. To our knowledge, the connection between the
fingerprint enhancement domain and growth ring identification was first
discovered in [Jon08]. We base our approach on _curved Gabor filters_ [Got11],
which perform very well on fingerprint images. Gabor filters [Gab46] have been
used across many disciplines for signal analysis and processing, and are often
amongst the top performers from the suitable candidates. A very interesting
property is their capability for retrieving phase information [Kon09], further
[GR19, AG21].

Figure 3: Ring detection process. We start with an input 2D texture (left) and
compute its local orientation. We then compute a curved region around each
point in the image and apply a Gabor filter on it. Accumulating over all
pixels gives the Gabor filter phase response. We then trace the rings in the
phase image to detect their positions (right).

##  3 Method

###  3.1 Overview

Our overall goal (Figs. Automatic inference of a anatomically meaningful solid
wood texture from a single photograph, 2) is to reconstruct a subset of
parameters that we can plug into the procedural wood appearance model proposed
by Liu et al. [LDHM16]. Their starting point is a cylindrical coordinate
system that has the z-axis aligned with the stem of the tree, and they apply
distortion functions that work on both the radial distance and the azimuth.
They further propose several versions of the spatially varying distortion
field, like procedural noise, or using a lookup table (“distortion texture”).
This can be either applied in a radially symmetric way, where the texture-
space y-coordinates correspond to radial distances and the x-coordinates to
values of z in tree space, or using a more elaborate helical wrapping to allow
distortions that vary with the azimuth.

Our approach aims at recovering this distortion texture, and a basic color
map, from a single photograph of the tangential plane of a wood cut. Please
also refer to Figure 3 for a graphical overview.

We recover this data using a three-stage process. First, we calculate both an
estimate of the board orientation, and an initial estimate of the distortion
field from the locations of growth rings (Sec. 3.3, 3.4). The growth rings are
found using curved Gabor filters (Sec. 3.2). Second, we polish the distortion
field using a signal-phase informed optimization scheme (Sec. 3.6). Finally, a
color map is calculated.

###  3.2 Curved Gabor filters

The key idea to curved Gabor filters is to calculate them in an image over
curved regions that follow local orientation (Fig. 3, left to middle). The
initial local orientation is obtained from the 2D image gradient: we use the
Scharr operator [Sch00] after smoothing the image with a small Gaussian
kernel. The gradient is averaged over a rectangular neighborhood, rotated
$90$\mathrm{\SIUnitSymbolDegree}$$ and normalized to obtain the initial local
orientation field.

Once the local orientation is known, we compute a curved region around each
point in the image. Each curved region consists of $2p+1$ parallel contours
with $2q+1$ points each, in our implementation $p=q=80$. They are constructed
by first sampling _perpendicular_ to the local orientation for $p$ steps, into
both positive and negative directions. From each of the $2p+1$ points, a
sampling walk _along_ the local orientation in positive and negative
directions for $q$ steps, gives the pixel values of the curved region. As the
step size is 1 pixel for both, we can look at the result as a pixelated patch
of $2p\times 2q$ pixels that follows the local curvature.

Because the patch is rotated so that it follows the local orientation, growth
rings will always be aligned horizontally when sampled isometrically. This
allows estimating the local ring-frequency of a patch from the number of peaks
is contains.

Using the (fixed) ring orientation inside each patch $P$ and the detected
frequency, we then calculate the convolution of the patch with a complex Gabor
kernel $K=\\{K_{\mathrm{re}},K_{\mathrm{im}}\\}$ of an appropriate size to get
the filtered patch $P\textquoteright$:

$P^{\prime}=P\ast K,$ (1) $\displaystyle
K_{\mathrm{re}}(x,y,\theta,f,\sigma_{x},\sigma_{y})$
$\displaystyle=\exp\bigg{(}-\frac{1}{2}\bigg{[}\frac{x_{\theta}^{2}}{\sigma_{x}^{2}}+\frac{y_{\theta}^{2}}{\sigma_{y}^{2}}\bigg{]}\bigg{)}\cdot\cos(2\pi
fx_{\theta}),$ (2) $\displaystyle
K_{\mathrm{im}}(x,y,\theta,f,\sigma_{x},\sigma_{y})$
$\displaystyle=\exp\bigg{(}-\frac{1}{2}\bigg{[}\frac{x_{\theta}^{2}}{\sigma_{x}^{2}}+\frac{y_{\theta}^{2}}{\sigma_{y}^{2}}\bigg{]}\bigg{)}\cdot\sin(2\pi
fx_{\theta}),$ (3) $\displaystyle x_{\theta}$
$\displaystyle=x\cdot\cos\theta+y\cdot\sin\theta,$ (4) $\displaystyle
y_{\theta}$ $\displaystyle=-x\cdot\sin\theta+y\cdot\cos\theta.$

We accumulate the resulting signals $P^{\prime}$ over all patches, giving a
complex valued filtered image $I_{f}=(I_{\mathrm{re}},I_{\mathrm{im}})$.
Transforming to polar coordinates gives a magnitude-image
$I_{\mathrm{mag}}=\sqrt{I_{\mathrm{re}}^{2}+I_{\mathrm{im}}^{2}}$, and a
phase-image $I_{\phi}=\arctan_{2}({I_{\mathrm{im}},I_{\mathrm{re}}})$. From
the phase information, we can estimate the location of tree ring boundaries at
where $I_{\phi}(x,y)=\frac{\pi}{2}$, the location of the signals negative zero
crossing. The phase image contains a sign ambiguity that we resolve later
(Sections 3.6, 4.3).

As as side-effect to calculating the Gabor filter on curved regions, it is
also possible to estimate local curvature and local ring frequency. The local
curvature can be computed from the absolute value of the differences of the
local orientation of the contours’ middle points and their end points. This
can be used to identify anomalies such as knots.

The final local orientation field is obtained from the phase image using the
same procedure that calculated the initial local orientation estimate, which
gives a smooth field that follows along the growth ring’s trajectory.

###  3.3 Ring detection

To detect a growth ring, we seed the algorithm at any point $p=(x,y)$ where
$I_{\phi}(x,y)\approx{}\frac{\pi}{2}$, and walk along the local orientation
field into both positive and negative directions. By following the orientation
field instead of just thresholding the phase image, our algorithm can trace
rings even across small anomalies.

###  3.4 Estimating board location

With the detected visible growth rings it becomes possible to reason about the
boards orientation within the 3D tree coordinate system. Using this positional
information then allows us to calculate an estimate of the distortion field
that would lead to the observed pattern. We show this for tangential cut
planes, but the approach can easily be extended to other cuts as well.

Our position estimate starts with automatically finding the projection of the
tree center (the $z$ axis) onto the plane. It is usually located between the
two rings with the largest mean distance from one another, with the remaining
inter-ring distances being monotonically decreasing to both sides. This
projection of the tree-center also gives the translation of the tangential
board along its y-axis. From the median average inter-ring distance, we can
induce the average scale of the rings in image-space coordinates. To find the
x-translation of the observed plane, we fit the radial scaling factor and
translation on the XY-plane using a brute force search in the neighbourhood of
the initial estimate, which only takes seconds to compute. The resulting sign
of the x-translation is ambiguous since we are only observing a plane.

###  3.5 Modeling growth distortions and initial guess

In their model, Liu et. al. [LDHM16] model distortions of points in tree space
$q$ in the radial and/or tangential direction, magnitudes $m_{r}$ and $m_{t}$.
The magnitudes are spatially varying and can be any function $f\colon R^{3}\to
R$, such as a Perlin noise, or a texture. The location of the distorted point
$q^{\prime}$ is given by

$q^{\prime}=d(p)=q+m_{r}(q)\hat{r}+m_{t}(q)\hat{t}.$ (5)

$\hat{r}$ and $\hat{t}$ are unit vectors into the radial and tangential
direction, respectively. As the initial guess for the distortion field, we use
the difference between the observed and the ideal, undistorted rings in image
space and map this displacement to radial distortions $m_{r}(q)\hat{(}r)$
using the board location estimate. Tangential distortions $m_{t}(q)$ are not
used here due to the ambiguity of radial and tangential distortions when only
observing a planar projection of the deformed tree rings.

###  3.6 Polishing the initial deformation field

A sufficiently accurate initial guess of the deformation field is important
for accurate convergence of the final polishing step. The periodic nature of
the texture function constrains gradient-based optimization to the interval of
the period the starting point is contained in. More precisely, the initial
deformation field needs to have a local phase error of less than
$\frac{\pi}{2}$ for best results.

Formally, we want to minimize an energy $E$ that is a function of the
reference image $J$ and a rendered image $I$.

$\operatorname*{arg\,min}_{m_{r},c,s_{r}}E(J,I)$ (6)

$I$ is obtained from a rendering operator, $\operatorname*{\mathfrak{R}}$,
that takes the boards orientation, given by the linear transformation $T$, a
ring scaling factor $s_{r}$, the distortion field $d(q)$ and a color map
$c(q_{r})$ to an rendered image $I$:

$I=\operatorname*{\mathfrak{R}}(T,s_{r},d,c)$ (7)

$s_{r}$ encodes the mean width of growth rings in radial units, the color map
$c$ describes, for each ring individually, the color-changes observable as the
tree growth forms early and late growth. The distortion field encodes the ring
growth variations. In its simplest form, the rendering function is just a
projection of the solid texture to image space by a cutting plane. Note that
all functions need to be differentiable in order for gradient descent
optimization to be applicable.

#### Loss function

Optimizing using a per-pixel loss, or a feature loss, will result in a good
appearance matched _image_ , even good estimates for BRDF parameters other
than color, through the use of a differentiable renderer as
$\operatorname*{\mathfrak{R}}$. Problems arise due to the co-optimization of
both distortion field and color map and will lead to a fit that lacks
anatomical meaning (Fig. 1). This is due to the non-orthogonality of the
parameters: From a reasonable gamut of colors in the color map, individual
surface pixels can be modulated by using the distortion field, making it a
proxy to achieve a certain surface color. This does not move the rings as a
whole, but rather individual points. The result is a non-smooth, "fuzzy"
deformation field that can even lead to ring fold-over.

To overcome this, the optimization needs to be constrained, or the metric
improved. We have experimented implementing a monotonic constraint on the
distorted radii to discourage foldover, and we also tried enforcing smoothness
by convolving the distortion field with a Gaussian kernel in exponentially
increasing intervals during optimization. However, we were not able to attain
significant improvements. Another potential issue arises from the periodic
discontinuity the procedural model contains, which is a direct result of the
underlying sawtooth-like function that describes the growth periods of non-
tropical wood grain (cf. figure 8). It has been shown that not accounting for
discontinuities can lead to convergence issues (shown for example by Loubet
et. al. in [LHJ19]) . We implemented an experimental reparameterizing renderer
on top of the Mitsuba 3 renderer [JSR∗22] to correctly handle texture
discontinuities, and also tried smoothing the falling edge of the yearly
growth boundary by convolving with a small Gaussian kernel. Neither lead to
satisfactory convergence behaviour in our experiments.

We finally opted for improving the optimization metric by introducing using a
loss based on signal phase. We also separate the optimization of the
deformation field from the color map into a two stage process. Our loss
encodes the difference in signal phase between points on the reference image,
and points on the rendered image of the current iteration.

$E=\big{(}J_{\phi}-I_{\phi}\big{)}^{2}$ (8)

The phase-based loss function lets the optimizer fit the locations of entire
rings holistically. One could implement a phase-based loss on top of a bank of
regular Gabor filters that recover the phase information for both the
reference image $I$ and the rendering of the current iteration. Propagating
loss gradients back through this filter bank is very costly, however. Due to
the constrained nature of the underlying solid texture (tree rings cannot be
arbitrarily placed, but follow certain rules), the phase information can be
directly obtained: rendering the fractional part of the radial coordinate in
the distorted tree’s cylindrical coordinate system, multiplied by $2\pi$
circumvents this bottleneck. The phase response of the curved Gabor filter
calculated from the input image earlier then serves as the optimization
target.

The phase ambiguity on the reference image can be resolved by assuming a fixed
orientation, and flipping the phase values by 180° beyond the projected tree
center, given it is visible in the image. We also note that the loss function
of the optimization needs to respect the cyclic nature of phase information
when taking differences.

#### Optimization algorithm and hyperparameters

We implemented the fitting procedure using a custom differentiable rasterizer.
The rasterizer is built on top of TensorFlow [AAB∗15], which is used as a
framework for reverse-mode automatic differentiation. We use their
implementation of the Adam optimizer [KB14] to drive the optimization loop.
The learning rate was set to $0.03$, and we run the optimization for $500$
epochs. A single run takes less than a minute on a NVIDIA RTX 2080 GPU for a
target texture with a resolution of $750x735$.

Figure 4: The procedural wood model is an infinite 3D texture centered around
a virtual tree (left). An arbitrary 3D mesh, such as a cube (middle), can be
positioned within the 3D texture coordinate system accordingly to where the
object was carved from the tree. This gives a textured solid wood object
(right) that can be rendered. Figure 5: Volumetric appearance is highly
dependent on grain orientation - shown here are renderings of 1.8mm thick
veneer sheets backlit by a spherical light source. The volumetric properties
are derived from our fit (right), and extruded along the surface normal
(left). Insets show side-views to illustrate grain Figure 6: Volumetric
rendering of a lampshade that was cut from a fit model using our method (Fig.
7, 2nd column). Both single scattering albedo and volumetric density are
heterogeneous. Figure 7: Results obtained with our approach. Rows from top to
bottom: Input texture, detected rings, frontal view of fit model (diffuse
reflectance only), volumetric visualization (false color), fit color map.

##  4 Results

###  4.1 Structurally matching wood appearance

We apply our method to several pieces of coniferous wood, where most rings are
identified correctly. Difficulties arise where the plane grazes a growth ring
tangentially, where tiny variations of the estimated board location and the
ring distortion can create huge changes in appearance. This can be seen in
Figure 7, 1st and 3rd column, where mistakenly identified rings add to the
problem, leading to a noisy, implausible distortion estimate in that area. In
contrast, the specimen in the third column does not suffer from this problem,
since it does not contain the projection of the tree center.

Through the fit, we obtain a solid texture (see also Figure 4) that also
supplies _semantic information_ that we use to derive other BRDF parameters.
In Figure Automatic inference of a anatomically meaningful solid wood texture
from a single photograph, we show a rendering that modulates the roughness and
displacement as a function to the lateness of growth of each ring, which
results in a very realistic appearance with specular highlights that are in
alignment with the late growth.

###  4.2 Volumetric appearance

For small material thicknesses, volumetric light transport becomes significant
for the appearance of wood. Examples are decorative use-cases such as veneer
lamp shades, or functional usages such as translucent wood touch panels, or
architectural use of delignified wood. To illustrate the importance of
subsurface reconstruction of the wood grain, we simplify the cellular anatomy
into an approximate participating medium. Using our fit, we set the spatially
varying density proportionally to the lateness of growth, latewood being more
dense than earlywood. The volumetric properties were estimated to a scattering
coefficient of $\mu_{t}\approx 16.0mm^{-1}$ for latewood, $\mu_{t}\approx
6.5mm^{-1}$ for earlywood, and a (wavelength dependent) single scattering
albedo derived from the surface reflectance.

To come up with the density parameters, we applied the Inverse Adding-Doubling
Method [Pra11] on data provided in [SKSG18], who performed measurements on
Sugi wood, a coniferous subspecies, using micro spectrometer hardware. The
heterogeneous single scattering albedo is calculated using the surface albedo
mapping function proposed by Elek et. al [ESZ∗17], their Eqn. 4. Fig. 5
compares the translucent appearance of both a model that was fit using our
pipeline, and a simple extrusion of the surface ring pattern into the depth of
the volume. The appearance difference is explained by the orientation of the
latewood shells, which is perpendicular in the extruded case, but follows wood
grain in the fit case. In Fig. 6 we show a volumetric rendering of a
lampshade. The characteristic warm-colored translucent appearance wood shows
when it is backlit is captured very well.

###  4.3 Dendrochonology

We confirm the robustness of our ring detection method by running it on a
dendrochonological dataset. There, accuracy of detection is very important,
since it is used as primary data for various other fields, including climate
research.

To detect rings in dendrochonolgy images, the images are first being scaled to
uniform width, and then converted to $HSV$ colorspace, from which only the
V-channel is used. The detection procedure is executed as detailed above. The
phase sign ambiguity is resolved by assuming an image orientation with the
youngest rings at the bottom of the image. The recovered rings are then
compared to ground-truth by measuring the closest distance between the ground-
truth label (one x/y coordinate pair for each ring). If the distance is
smaller than $3$ pixels, conforming with the evaluation shown in [FDBP17], the
ring is counted as a match. Table 1 summarizes the performance of our
approach. Performance for coniferous woods is good. For ring-porous woods, an
anatomical variety that has shown to be problematic in previous work, our
results are exellent. The high accuracy of our fit means that it can also
supply additional data for dendrochronological evaluations. Fig. 8 shows a
partial color map, superimposed with the inverse lightness, which correlates
with the time the latewood transition occurred in the respective year.

Figure 8: From the extracted color map, it becomes possible to gain insight into the specifics of each growth cycle, here the latewood-transition. The plot shows the normalized and inverted lightness value of the HSL-transformed colorbar, smoothed by a gaussian kernel. The lightness is proportional to the cell-wall thickness and lignin content. |  | Ours |  | [FD18] | [FDBP17]
---|---|---|---|---|---
species | rings | sen | pre | sen | pre | sen | pre
larch | 99 | 0.77 | 0.83 |  |  | 0.90 | 1.00
spruce | 115 | 0.93 | 0.96 |  |  | 0.99 | 1.00
pine | 108 | 0.79 | 0.91 |  |  | 0.93 | 0.97
ash | 114 | 0.98 | 0.97 | 0.97 | 0.98 | 0.45 | 0.85
birch | 30 | 0.17 | 0.06 |  |  | 0.86 | 0.74
Table 1: Comparison of sensitivity (sen) and precision (pre) of our ring-
detection method and two state of the art methods on a dendrochronological
dataset. [FD18] was developed especially with ringporous wood in mind. Our
method performs especially well on ringporous wood (ash) and gives decent
results on confierous species (first three rows). Diffuseporous woods such as
birch are problematic, where our approach produces many false positives.

##  5 Conclusion

In this paper, we demonstrated curved Gabor filters to be an efficient
technique for the extraction of tree ring information from single images of
planar wooden boards. This information can then be used to build a realistic,
anatomically meaningful 3D procedural texture that closely matches the
structure of the log the original board was cut from.

The main focus of our work was to establish an automatic pipeline that is fast
and accurate. The technology we ported and adapted from fingerprint detection
not only serves as a tree ring detector with great performance, but also
enables the efficient optimization of the ring deformations based on signal
phase. The reconstruction of additional BSDF parameters by using
differentiable rendering could further improve the model: but as the results
of this paper show, the approach we currently use is already suitable for
production work. Our pipeline further enables volumetric rendering of wood’s
subsurface light transport, and can pave a way to finding a faster,
approximate BSSRDF, which we see as a interesting direction for further work.

We also see a potential of our robust ring extraction technique to be adopted
by the dendrochronology community, where ring labelling is an important task
that our method solves as a side effect.

## References

  * [AAB∗15] Abadi M., Agarwal A., Barham P., Brevdo E., Chen Z., Citro C., Corrado G. S., Davis A., Dean J., Devin M., Ghemawat S., Goodfellow I., Harp A., Irving G., Isard M., Jia Y., Jozefowicz R., Kaiser L., Kudlur M., Levenberg J., Mané D., Monga R., Moore S., Murray D., Olah C., Schuster M., Shlens J., Steiner B., Sutskever I., Talwar K., Tucker P., Vanhoucke V., Vasudevan V., Viégas F., Vinyals O., Warden P., Wattenberg M., Wicke M., Yu Y., Zheng X.:  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.  URL: http://tensorflow.org/.
  * [AAL16] Aittala M., Aila T., Lehtinen J.:  Reflectance modeling by neural texture synthesis.  _ACM Transactions on Graphics (ToG) 35_ , 4 (2016), 1–13.
  * [AG21] Alaifari R., Grohs P.:  Gabor phase retrieval is severely ill-posed.  _Applied and Computational Harmonic Analysis 50_ (2021), 401–419.
  * [BMM∗21] Bangaru S., Michel J., Mu K., Bernstein G., Li T.-M., Ragan-Kelley J.:  Systematically differentiating parametric discontinuities.  _ACM Trans. Graph. 40_ , 107 (2021), 107:1–107:17.
  * [BS12] Burley B., Studios W. D. A.:  Physically-based shading at disney.  In _ACM SIGGRAPH_ (2012), vol. 2012, vol. 2012, pp. 1–7.
  * [DVGNK99] Dana K. J., Van Ginneken B., Nayar S. K., Koenderink J. J.:  Reflectance and texture of real-world surfaces.  _ACM Transactions On Graphics (TOG) 18_ , 1 (1999), 1–34.
  * [ESZ∗17] Elek O., Sumin D., Zhang R., Weyrich T., Myszkowski K., Bickel B., Wilkie A., Křivánek J.:  Scattering-aware texture reproduction for 3D printing.  _ACM Transactions on Graphics (Proceedings of SIGGRAPH Asia) 36_ , 6 (2017), 241:1–241:15.
  * [FD18] Fabijańska A., Danek M.:  Deepdendro – a tree rings detector based on a deep convolutional neural network.  _Computers and Electronics in Agriculture 150_ (2018), 353–363.  doi:https://doi.org/10.1016/j.compag.2018.05.005.
  * [FDBP16] Fabijańska A., Danek M., Barniak J., Piórkowski A.:  A comparative study of image enhancement methods in tree-ring analysis.  In _International Conference on Image Processing and Communications_ (2016), Springer, pp. 69–78.
  * [FDBP17] Fabijańska A., Danek M., Barniak J., Piórkowski A.:  Towards automatic tree rings detection in images of scanned wood samples.  _Computers and Electronics in Agriculture 140_ (2017), 279–289.
  * [Gab46] Gabor D.:  Theory of communication. part 1: The analysis of information.  _Journal of the Institution of Electrical Engineers-Part III: Radio and Communication Engineering 93_ , 26 (1946), 429–441.
  * [GHYZ20] Guo Y., Hašan M., Yan L., Zhao S.:  A bayesian inference framework for procedural material parameter estimation.  In _Computer Graphics Forum_ (2020), vol. 39, Wiley Online Library, pp. 255–266.
  * [GLLD12] Galerne B., Lagae A., Lefebvre S., Drettakis G.:  Gabor noise by example.  _ACM Transactions on Graphics (ToG) 31_ , 4 (2012), 1–9.
  * [Got11] Gottschlich C.:  Curved-region-based ridge frequency estimation and curved gabor filters for fingerprint image enhancement.  _IEEE Transactions on Image Processing 21_ , 4 (2011), 2220–2227.
  * [GR19] Grohs P., Rathmair M.:  Stable gabor phase retrieval and spectral clustering.  _Communications on Pure and Applied Mathematics 72_ , 5 (2019), 981–1043.
  * [HKRFM17] Hartmann F. P., K Rathgeber C. B., Fournier M., Moulia B.:  Modelling wood formation and structure: power and limits of a morphogenetic gradient in controlling xylem cell proliferation and growth.  _Annals of forest science 74_ , 1 (2017), 1–15.
  * [HMR20] Henzler P., Mitra N. J., Ritschel T.:  Learning a neural 3d texture space from 2d exemplars.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2020), pp. 8356–8364.
  * [Jon08] Jonsson C.:  _Detection of annual rings in wood_.  Master’s thesis, Linköping University, SE-601 74 Norrköping, Sweden, 11 2008.
  * [JSR∗22] Jakob W., Speierer S., Roussel N., Nimier-David M., Vicini D., Zeltner T., Nicolet B., Crespo M., Leroy V., Zhang Z.:  Mitsuba 3 renderer, 2022.  https://mitsuba-renderer.org.
  * [KB14] Kingma D. P., Ba J.:  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [Kon09] Kong A. W.-K.:  An analysis of gabor detection.  In _International Conference Image Analysis and Recognition_ (2009), Springer, pp. 64–72.
  * [LDHM16] Liu A. J., Dong Z., Hašan M., Marschner S.:  Simulating the structure and texture of solid wood.  _ACM Transactions on Graphics (TOG) 35_ , 6 (2016), 1–11.
  * [Len20] Lenty B.:  Tree-ring growth measurements automation using machine vision.  In _Photonics Applications in Astronomy, Communications, Industry, and High Energy Physics Experiments 2020_ (2020), vol. 11581, International Society for Optics and Photonics, p. 115810V.
  * [LHJ19] Loubet G., Holzschuch N., Jakob W.:  Reparameterizing discontinuous integrands for differentiable rendering.  _Transactions on Graphics (Proceedings of SIGGRAPH Asia) 38_ , 6 (Dec. 2019).  doi:10.1145/3355089.3356510.
  * [LIY∗22] Larsson M., Ijiri T., YoshidaA H., Huber J. A., Fredriksson M., Broman O., Igarashi T.:  Procedural texturing of solid wood with knots.  _ACM Transactions on Graphics (TOG)_ (2022).
  * [LLDD09] Lagae A., Lefebvre S., Drettakis G., Dutré P.:  Procedural noise using sparse gabor convolution.  _ACM Transactions on Graphics (TOG) 28_ , 3 (2009), 1–10.
  * [LP00] Lefebvre L., Poulin P.:  Analysis and synthesis of structural textures.  In _Graphics Interface_ (2000), vol. 2000, pp. 77–86.
  * [LVS∗18] Li Y., Vasileva E., Sychugov I., Popov S., Berglund L.:  Optically transparent wood: Recent progress, opportunities, and challenges.  _Advanced Optical Materials 6_ , 14 (2018), 1800059.
  * [LWL∗20] Le N. T., Wang J.-W., Le D. H., Wang C.-C., Nguyen T. N.:  Fingerprint enhancement based on tensor of wavelet subbands for classification.  _IEEE Access 8_ (2020), 6602–6615.
  * [Mar18] Martins A. L. R.:  Towards automatic identification of woods from microscopic images.
  * [MGSS∗21] Martinez-Garcia J., Stelzner I., Stelzner J., Gwerder D., Schuetz P.:  Automated 3d tree-ring detection and measurement from x-ray computed tomography.  _Dendrochronologia 69_ (2021), 125877.
  * [oM22] of Minnesota U.:  DendroElevator website, 2022.  URL: https://dendro.elevator.umn.edu.
  * [PAH∗22] Poláček M., Arizpe A., Hüther P., Weidlich L., Steindl S., Swarts K.:  Automation of tree-ring detection and measurements using deep learning.  _bioRxiv_ (2022).  doi:10.1101/2022.01.10.475709.
  * [Per85] Perlin K.:  An image synthesizer.  _ACM Siggraph Computer Graphics 19_ , 3 (1985), 287–296.
  * [Pra11] Prahl S. A.:  Everything i think you should know about inverse adding-doubling.  _Oregon Medical Laser Center, St. Vincent Hospital 1344_ (2011), 1–74.
  * [RSB∗21] Rademacher T., Seyednasrollah B., Basler D., Cheng J., Mandra T., Miller E., Lin Z., Orwig D. A., Pederson N., Pfister H., et al.:  The wood image analysis and dataset (wiad): Open-access visual analysis tools to advance the ecological data revolution.  _Methods in Ecology and Evolution 12_ , 12 (2021), 2379–2387.
  * [Sch00] Scharr H.:  _Optimal operators in digital image processing_.  PhD thesis, 2000.
  * [SKSG18] Sugimoto H., Kawabuchi S., Sugimori M., Gril J.:  Reflection and transmission of visible light by sugi wood: effects of cellular structure and densification.  _Journal of Wood Science 64_ , 6 (2018), 738–744.
  * [SLH∗20] Shi L., Li B., Hašan M., Sunkavalli K., Boubekeur T., Mech R., Matusik W.:  Match: Differentiable material graphs for procedural material capture.  _ACM Trans. Graph. 39_ , 6 (Dec. 2020), 1–15.
  * [THG16] Thai D. H., Huckemann S., Gottschlich C.:  Filter design and performance evaluation for fingerprint image segmentation.  _PloS one 11_ , 5 (2016), e0154160.
  * [VCL∗18] Vasileva E., Chen H., Li Y., Sychugov I., Yan M., Berglund L., Popov S.:  Light scattering by structurally anisotropic media: A benchmark with transparent wood.  _Advanced Optical Materials 6_ , 23 (2018), 1800999.
  * [WQZ∗10] Wang H.-j., Qi H.-n., Zhang G.-q., Li W.-z., Wang B.-h.:  An automatic method of tree-rings boundary detection on wood micro-images.  In _2010 International Conference on Computer Application and System Modeling (ICCASM 2010)_ (2010), vol. 2, IEEE, pp. V2–477.