# Identification and classification of exfoliated graphene flakes from
microscopy images using a hierarchical deep convolutional neural network

Soroush Mahjoubi
Department of Civil, Environmental and Ocean Engineering
Stevens Institute of Technology
Hoboken, New Jersey 07030, USA
Fan Ye
Department of Electrical and Computer Engineering
University of Massachusetts Amherst
Amherst, Massachusetts 01003, USA
Yi Bao
Department of Civil, Environmental and Ocean Engineering
Stevens Institute of Technology
Hoboken, New Jersey 07030, USA
<EMAIL_ADDRESS>
Weina Meng
Department of Civil, Environmental and Ocean Engineering
Stevens Institute of Technology
Hoboken, New Jersey 07030, USA
<EMAIL_ADDRESS>
Xian Zhang
Department of Mechanical Engineering
Stevens Institute of Technology
Hoboken, New Jersey 07030, USA
<EMAIL_ADDRESS>

###### Abstract

Identification of the mechanically exfoliated graphene flakes and
classification of the thickness is important in the nanomanufacturing of next-
generation materials and devices that overcome the bottleneck of Moore’s Law.
Currently, identification and classification of exfoliated graphene flakes are
conducted by human via inspecting the optical microscope images. The existing
state-of-the-art automatic identification by machine learning is not able to
accommodate images with different backgrounds while different backgrounds are
unavoidable in experiments. This paper presents a deep learning method to
automatically identify and classify the thickness of exfoliated graphene
flakes on Si/SiO2 substrates from optical microscope images with various
settings and background colors. The presented method uses a hierarchical deep
convolutional neural network that is capable of learning new images while
preserving the knowledge from previous images. The deep learning model was
trained and used to classify exfoliated graphene flakes into monolayer (1L),
bi-layer (2L), tri-layer (3L), four-to-six-layer (4-6L), seven-to-ten-layer
(7-10L), and bulk categories. Compared with existing machine learning methods,
the presented method possesses high accuracy and efficiency as well as
robustness to the backgrounds and resolutions of images. The results indicated
that our deep learning model has accuracy as high as 99% in identifying and
classifying exfoliated graphene flakes. This research will shed light on
scaled-up manufacturing and characterization of graphene for advanced
materials and devices.

##  1 Introduction

Two-dimensional (2D) materials, endowed with superior electrical [1],
ferromagnetic [2, 3], semiconducting [4], superconducting [5], thermal [6-10],
and optical properties [11, 12], are the next-generation materials in the
areas of mechanical, photonic, and electronic. 2D materials are widely
obtained by mechanical exfoliation using scotch tapes [13, 14], which ensures
the most pristine properties and conditions of 2D materials. Mechanical
exfoliation methods yield flakes of 2D materials with random thickness,
shapes, and morphology on the substrate, followed by a manual search under an
optical microscope by human intuitions. For scientific and device development
purposes, the ideal flakes have a uniform thickness. Many interesting
phenomena only occur at samples of one to a few atomic layers. In current
practices, flake thickness is mainly evaluated by human operators through
examination of the color and contrast of flakes from the microscope images.
The flakes are then further evaluated by material characterization equipment,
such as Raman microscope and atomic force microscope [15-20]. There are
several disadvantages in this process: manual examination inevitably yields
errors due to different colors and background conditions regarding color and
brightness; and characterization of 2D flakes is costly, inefficient, and
time-consuming. Therefore, an efficient and accurate method for identifying
the thickness is highly desired [21-24].

A few automatic methods to identify 2D flakes from optical microscope images
were reported in references [17, 25-29]. Five major limitations were
identified from the existing automatic methods: (1) The existing methods were
designed for specific microscope conditions such as background color and light
brightness, compromising the applicability for different conditions; (2) the
existing machine learning models are not capable of learning new micrographs,
so re-training the existing models using new data leads to catastrophic loss
of the knowledge learned from old data [30]; (3) optical microscope images
often suffer from artifacts such as poor contrast, vignetting, and
overexposure. Multiple methods were developed to enhance microscope images for
human vision [31, 32], but not applicable for machine vision; (4) imbalanced
class distribution in training images compromises the evaluation of under-
represented classes [33]. The areas of the thin flakes are relatively small
compared with the areas of the substrate (background) and bulk flakes. (5) The
classification precision of graphene thickness was rough. The flakes were
categorized into three classes [25], which are insufficient in many
applications.

In this work, we present a novel hierarchical deep learning method based on an
unsupervised classification model and multiple semantic segmentation models
and evaluate the performance of the method in automatically identifying
exfoliated graphene flakes and determining thicknesses from optical microscope
images. The flakes are categorized into six classes: monolayer (1L), bi-layer
(2L), tri-layer (3L), four-to-six-layer (4-6L), seven-to-ten-layer (7-10L),
and bulk. Novel computer vision techniques are presented to improve microscope
images using an image quality metric and improved adaptive gamma correction
[34]. The proposed method leverages weak learning, data augmentation,
iterative stratification [35], and weighted cross-entropy loss to improve the
accuracy and the generalization performance. The effects of the resolution and
the background condition of images are considered to improve the robustness of
the deep learning model. This method promotes the capability of processing and
manufacturing 2D materials and devices by improving efficiency and accuracy
while minimizing human intervention.

##  2 Methodology

We produced 2D flakes using the scotch tape mechanical exfoliation method [13,
14] from bulk crystals (Graphene Supermarket). After bulk crystal was folded
on scotch tape for multiple times, 2D flakes were directly transferred onto
the SiO2/Si substrate by a 2-minute scratch. The SiO2 thickness was 285 nm for
the best optical contrast to identify the thickness of 2D flakes. An optical
microscope (Nikon Eclipse LV150N) was used to capture RGB images of graphene
flakes from 1L up to 10 L on the SiO2/Si substrate using different microscope
settings. A total of 273 high-resolution optical microscope images were
captured. An atomic force microscope (Bruker BioScope Resolve) and Raman
spectroscope (Renishaw inVia Confocal) were used to evaluate the thickness of
each flake. Pixel-wise ground truth labels of MATLAB were annotated with the
thickness evaluation results. The pixels were categorized into seven classes:
background, 1L, 2L, 3L, 4-6L, 7-10L, and bulk.

Figure 1 shows the flowchart with five main steps: (1) Dataset development:
Optical microscopy images are captured to develop a dataset for training a
machine learning model (Section 2.1). (2) Dataset enhancement: The quality of
images is improved using multiple methods developed in this study (Section
2.2). (3) Dataset standardization and categorization: The images are
standardized and categorized using an unsupervised machine learning model
(Section 2.3). (4) Semantic segmentation: A deep convolutional neural network
is trained to identify exfoliated graphene flakes and quantify layer numbers
(Section 2.4). (5) Performance evaluation: The performance of the framework is
evaluated using performance metrics (Section 2.5).

Figure 1: Flowchart of the robust machine intelligence framework for
classification of 2D flakes.

###  2.1 Statistical Analysis

Statistical analysis was conducted using the ground truth labels to evaluate
the distributions of the seven classes. The weight of the i-th class was
defined as the number of pixels of the i-th class divided by the total number
of pixels. Table 1 lists the results of the mean, median, and maximum class
weights of the seven classes, and the percentage of images with a class zero-
class weight. The mean class weights of classes 1, 3, 4, and 5 were less than
1%. The median class weights of classes 1 to 6 were zero, while the median
class weight of background was higher than 95%. The percentages of images with
zero class weight for classes 1 (1L) and 5 (7-10L) were higher than 75%. The
statistical results indicated that the different classes in the dataset were
highly imbalanced. The interested classes 1 to 6 were under-represented. The
background class was over-represented. The statistical analysis was extended
to the colors of images. It was found that different images had different
background colors and brightness (see Figure 2). Due to the variation of pixel
color and intensity, it is unsuitable to use color-based segmentation methods
[17]. In summary, the dataset of images introduces challenges due to the
imbalanced classes and scattering of background colors and brightness.

In addition, microscopy images introduce challenges due to overexposure,
specular reflection, vignetting, dirt, and out-of-focus errors: (1)
Overexposure: When excessive light is received by the microscope, overexposure
may occur and cause loss of detail in images. (2) Specular reflection:
Overhead lighting may cause specular reflection of graphene flakes and lead to
over-saturation and loss of detail. (3) Vignetting: Vignetting appears as a
radial darkening towards the corners of an image. (4) Dirt: Dirt and debris
cause blur images. (5) Out-of-focus: Improper focus of the microscope causes
blurred images and loss of details.

(a)

(b)

Figure 2: Scatter plots of pixel intensity levels: (a) the pixel intensity of
the background, and (b) the average pixel intensity of images. Each dot
represents a microscope image. The color of dot represents the calculated
color value. Table 1: Statistics of class weights

Designation | Description | Mean class | Median class | Maximum class | Images with zero
---|---|---|---|---|---
of class | of class | weight | weight | weight | class weight
0 | Background | 91.99% | 95.87% | 99.83% | 0.00%
1 | $1L^{1}$ | 0.59% | 0.00% | 19.72% | 75.63%
2 | 2L | 1.77% | 0.00% | 62.62% | 55.51%
3 | 3L | 0.70% | 0.00% | 15.46% | 55.15%
4 | 4-6L | 0.68% | 0.00% | 12.17% | 61.76%
5 | 7-10L | 0.47% | 0.00% | 27.51% | 79.04%
6 | Bulk | 3.81% | 0.00% | 40.61% | 51.84%

  * 1

“L” stands for layer(s).

###  2.2 Optimized adaptive gamma correction

To improve the quality of images, we present an optimized adaptive gamma
correction method (refer Figure 1), through integrating the improved adaptive
gamma correction method [34], noise-aware image quality metrics [38], and
particle swarm optimization [37]. The development and integration of the
methods are elaborated in the following subsections.

####  2.2.1 Improved adaptive gamma correction

The improved adaptive gamma correction method modifies light intensity while
preserving the color properties of images with overexposure problems in the
following eight steps [34]:

  1. 1.

RGB images are converted to Y $C_{b}$ $C_{r}$ color space, where Y is the luma
component for light intensity; $C_{b}$ and $C_{r}$ are the blue-difference and
red-difference chroma component, respectively.

  2. 2.

The luma component is used to convert an image into a negative image:

${Y\left(x,y\right)}_{N}=255-Y\left(x,y\right)$ (1)

where $Y\left(x,y\right)_{N}$ is the negative image; $Y\left(x,y\right)$ is
the luma component of the input image; $x$ and $y$ are the coordinates of
images: $x=1,2,\ldots,M$, and $y=1,2,\ldots,N$, where $M$ and $N$ are the
width and height of an image in pixel.

  3. 3.

The probability density function of pixels with an intensity level $l$ is
expressed as:

$PDF\left(l\right)=\frac{n_{l}}{MN}$ (2)

where $n_{l}$ is the number of pixels with an intensity level $l$.

  4. 4.

The weighting distribution function is expressed as:

$W\left(l\right)={PDF}_{max}\left[\frac{PDF\left(l\right)-{PDF}_{min}}{{PDF}_{max}-{PDF}_{min}}\right]^{\alpha}$
(3)

where ${PDF}_{max}$ and ${PDF}_{min}$ are the maximum and minimum probability
density functions of pixels with an intensity level l, respectively; $\alpha$
is the adjusted parameter; and $l_{max}$ is the maximum intensity level of the
luma component.

  5. 5.

The cumulative distribution function is expressed as:

$CDF\left(l\right)=\frac{\sum_{i=1}^{l}W\left(i\right)}{\sum_{i=1}^{l_{max}}W\left(i\right)}$
(4)

  6. 6.

The modified luma component of the negative image is obtained from the
following transformation function:

${I^{\prime}\left(x,y\right)}_{N}=l_{max}\left(\frac{l}{l_{max}}\right)^{1-CDF\left(l\right)}$
(5)

  7. 7.

The modified luma component of the input image is determined as:

$Y^{\prime}\left(x,y\right)=255-{Y^{\prime}\left(x,y\right)}_{N}$ (6)

  8. 8.

The modified luma component and chroma components of the image are converted
to a RGB image, as the improved image.

####  2.2.2 Noise-aware image quality metric

A noise-aware image quality metric was used to assess image quality based on
the gradient, entropy, and noise of an image [38]. We introduced three metrics
based on gradient, entropy, and noise. The gradient-based metric
($M_{gradient}$) was used to evaluate the edge information.

${\widetilde{g}}_{i}=\frac{1}{N_{g}}[\log{\lambda(}g_{i}-\gamma)+1]$ (7)

$N_{g}=\log{\left[\lambda\left(1-\gamma\right)+1\right]}$ (8)

$G_{j}=\sum_{i\in c_{j}}{\widetilde{g}}_{i}\ ,\ j=1,2,\ldots,N_{c}$ (9)

$M_{gradient}=\frac{K_{G}E\left(G\right)}{S\left(G\right)}$ (10)

where ${\widetilde{g}}_{i}$ is the amount of gradient information at pixel
$i$; $g_{i}$ is the gradient magnitude at pixel $i$, which is estimated using
Sobel operator; $\lambda$ is the control parameter; $\gamma$ is the activation
threshold; $N_{g}$ and $K_{G}$ are normalization factors; $G_{j}$ is the grid
cell; $N_{c}$ is the number of grid cells; $E$ and $S$ are the mean and
standard deviation operators, respectively.

The entropy-based metric was used to estimate the information of an image:

$M_{entropy}=-k_{e}\sum_{i=0}^{255}{P(i)\log_{2}{P(i)}}$ (11)

where $P(i)$ is the probability of intensity level $i$ in a grayscale image;
$k_{e}$ is a normalization factor.

The noise-based metric is defined as follows:

$M_{noise}=\sum_{j=1}^{3}{\sqrt{\frac{\pi}{2}}\frac{1}{N_{p}}\sum_{i}H\left(i\right)U\left(i\right)\left|I_{j}\ast
K\right|(i)}$ (12)

where $N_{p}$ is the number of pixels whose pixel value is 1 in $H.U$, where
“.” is the elementwise multiplication operator; $H(i)$ is the homogeneous
region mask for pixel $i$; $U(i)$ is the overexposed-underexposed region mask
for pixel $i$; $I_{j}$ is the pixel value of the input image in channel $j$;
and “$\ast$” is the convolution operator; $K$ is a noise estimation kernel
that was used to evaluate the noise level of an image with Gaussian noise, as
defined in Equation (13); the inhomogeneous, overexposed, and underexposed
regions of images were removed by using a homogeneous mask and overexposed-
underexposed mask defined in Equation (14) and Equation (15), since the noise
estimation kernel overestimates noises in those regions.

$K=\left[\begin{array}[]{ccc}1&-2&1\\\ -2&4&-2\\\ 1&-2&1\end{array}\right]$
(13) $H(i)=\left\\{\begin{array}[]{c}1\ \ \ g_{i}\leq\delta\\\ 0\ \ \
g_{i}>\delta\end{array}\right.$ (14) $U(i)=\left\\{\begin{array}[]{c}1\ \ \
{\tau}_{l}\leq I_{i}\leq{\tau}_{u}\\\ 0\ \ \ \ \ otherwise\end{array}\right.$
(15)

where $\delta$ is the adaptive threshold; $\tau_{l}$ and $\tau_{u}$ are the
lower and upper bounds, respectively; and $I_{i}$ is the pixel value of pixel
$i$. Finally, the image quality metric is derived by combining gradient,
entropy, and noise-based metrics:

$f\left(I\right)=AM_{gradient}+BM_{entropy}-CM_{noise}$ (16)

where $A$, $B$, and $C$ are user parameters, which are set to 0.4, 0.6, and
0.6 according to [34]. A high $f\left(I\right)$ implies higher edge and
texture details and lower noise. Therefore, it is promising to maximize the
quality metric to obtain a well-exposed microscopy image.

####  2.2.3 Particle swarm optimization

With the noise-aware image quality metric, the particle swarm optimization
algorithm was used to optimize the adjusted parameter, which is the only
control parameter of the improved adaptive gamma correction method. The image
quality metrics were used to define the objective function for optimization.
We used 20 search agents that moved in the search space following Equation
(17) to seek the optimal solutions in 30 iterations for five independent runs:

$\phi_{i}^{\prime}=\phi_{i}+V_{i}^{\prime}$ (17) $V_{i}^{\prime}=\omega
V_{i}+c_{1}r_{1}\left({pbest}_{i}-\phi_{i}\right)+c_{2}r_{2}\left(gbest-\phi_{i}\right)$
(18) $\omega=\left(\frac{itermax-
iter}{itermax-1}\right)\left(\omega_{max}-\omega_{min}\right)+\omega_{min}$
(19)

where $\phi_{i}^{\prime}$ and $\phi_{i}$ are the new and current positions of
the $i$-th search agent; $V_{i}^{\prime}$ is the current velocity of the
$i$-th agent; $V_{i}$ is the previous velocity of the $i$-th agent;
${pbest}_{i}$ is the best solution obtained by the $i$-th agent; $gbest$ is
the global best solution; $\omega$ is the inertia weight; $\omega_{max}$ and
$\omega_{min}$ are the maximum and minimum inertia weights, respectively;
$c_{1}$ and $c_{2}$ are the acceleration coefficients; $r_{1}$ and $r_{2}$ are
random numbers in the range of 0 to 1.

###  2.3 Image standardization and categorization

The images with different resolutions were resized to a standard size (256 ×
256 pixels). Next, the resized images were normalized by dividing the pixel
intensities in each channel by 255, so that the pixel intensity values were in
the range of 0 to 1. With the different distributions of colors in the images,
an unsupervised classification model was developed based k-means++ clustering
with squared Euclidean distance to categorize the images into different groups
according to the chroma channels ($C_{b}$ and $C_{r}$, see Figure 3) [39].

Figure 3: The geography of groups in the $C_{b}$ and $C_{r}$ space. Each dot
represents a microscope image. The color of each dot represents the mean pixel
intensity. The microscope images are grouped by k-means clustering based on
the blue-difference and red-difference chroma components.

###  2.4 Semantic segmentation models

A semantic segmentation model was developed based on an object-contextual
representation network with high-resolution network [40, 41]. The model
consisted of multiple convolutional neural networks for pixel-wise
classification, and a convolutional neural network was trained for each class.
A weak learning method was presented to enable the neural network to learn
from new images while preserving the knowledge learned from previous images.
Three methods were used for the imbalanced dataset: (1) A multi-class cross-
entropy loss function was presented to mitigate the extent of imbalance. (2)
An iterative stratification method was proposed to split the highly imbalanced
dataset into train and test sets. (3) Image augmentation was performed to
artificially enlarge the dataset. These methods are elaborated in the
following subsections.

####  2.4.1 Object-contextual representation

The object-contextual representation network has six main steps (Figure 1):

  1. 1.

Feature extraction: HRNetV2-W18 is used as the feature extractor by the neural
network with 20 million parameters and pretrained using ImageNet for faster
convergence [40].

  2. 2.

Coarse segmentation: With the extracted features, a coarse segmentation map is
generated for each class. The pixel value indicates the class of a pixel. The
map is obtained according to the output of HRNetV2-W18 and a 1×1 convolution
unit [40].

  3. 3.

Object region representation: The object region representation for class $k$
is defined as:

$f_{k}=\sum_{i\in I}{{\overline{m}}_{k,i}x_{i}}$ (20)
${\overline{m}}_{k,i}=softmax\left(m_{k,i}\right)=\frac{e^{m_{k,i}}}{\sum_{j\in
I}e^{m_{k,j}}}$ (21)

where $f_{k}$ is the object region representation for class $k$; $x_{i}$ is
the output of the backbone network for pixel $i$; ${\overline{m}}_{k,i}$ is
the normalized value of $m_{k,i}$; $m_{k,i}$ is the value of pixel $i$ in
coarse segmentation map for class $k$; the $softmax$ function is used to
normalize the coarse segmentation maps.

  4. 4.

Pixel-region relation: The relation between pixels with each class is defined
as:

$\omega_{k,i}=\frac{e^{\psi\left(x_{i}\right).\psi(f_{k})}}{\sum_{j=1}^{K}{\psi\left(x_{i}\right).\psi(f_{j})}}$
(22)

where $K$ is the number of classes; $\psi$ is the transformation function
consisting of 1×1 convolution, batch normalization, and $ReLU$ activation
function.

  5. 5.

Object contextual representation: The object contextual representation is
computed by aggregating the pixel-region relation with all coarse segmentation
maps:

$y_{i}=\psi\left[\sum_{k=1}^{K}\omega_{k,i}\psi\left(f_{k}\right)\right]$ (23)

where $y_{i}$ is the object contextual representation for pixel $i$; $f_{j}$
is the object region representation for class $k$; $K$ is the number of
classes; and $\omega_{k,i}$ is the pixel-region relation for class $k$ and
pixel $i$.

  6. 6.

Augmented representation: The final semantic segmentation map is obtained
using the transformation function, the output of the backbone network, and
object contextual representation:

$z_{i}=\psi([x_{i}\ y_{i}]T)$ (24)

where $z_{i}$ is the final segmentation solution for pixel $i$; $x_{i}$ is the
output of backbone network for pixel $i$; and $y_{i}$ is the object contextual
representation for pixel $i$.

####  2.4.2 Weak learning

Weak learning was proposed to address two challenges in training of the
convolutional neural networks: (1) The dataset is limited and becomes smaller
when it is divided into different groups. (2) It is time-consuming to
repeatedly train the convolutional neural network when new images are added to
the dataset.

Weak learning was performed in two steps: (1) A convolutional neural network
was trained using the whole dataset. (2) Transfer learning was performed. The
neural network was retained for each group of images while the learning rate
was near zero. The learning rate is a hyperparameter of convolutional neural
network and controls the changes of the weights in training. Since the
learning rate was low, the change of weights was small, so the knowledge
learned from previous images was reserved. Meanwhile, the convolutional neural
network learned new images and slightly modified the weights to achieve high
performance for the new images.

####  2.4.3 Multi-class cross-entropy loss

Cross-entropy was applied to obtain the discrepancy between the predicted
distribution mask $P\left(x\right)$ and the true distribution mask
$D\left(x\right)$:

$CE\left(P,D\right)=\sum_{x\epsilon\mathrm{\Omega}}{D\left(x\right)\mathrm{log}\mathrm{}(P\left(x\right))}$
(25)

where $\Omega$ is the mask region. The multi-class cross-entropy loss function
is expressed as:

${Loss}_{C}=\sum_{i=1}^{K}\sum_{x\epsilon\Omega}{d_{x,i}\mathrm{log}\mathrm{}(p_{x,i})}$
(26)

where $K$ is the number of classes; $d_{x,i}$ is the true one-hot distribution
probability of pixel $x$ on class $i$; and $p_{x,i}$ is the predicted
probability for pixel $x$ and class $i$.

Since the annotated pixels of under-represented classes were scarce, the
effect of the under-represented classes on the loss function was limited.
Thus, the neural network tends to maximize the prediction accuracy for over-
represented classes while ignoring the under-presented classes. This research
assigned sample weights to mitigate the representation problem:

${Loss}_{W}=\sum_{i=1}^{K}\sum_{x\epsilon\Omega}{{w_{i}d}_{x,i}\mathrm{log}\mathrm{}(p_{x,i})}$
(27)

where $w_{i}$ is the sample weight for class $i$.

There was no definitive advice on how to set the sample weights. We proposed
to calculate the sample weights using the following equation:

$w_{i}={(\frac{1}{\mu_{i}})}^{\beta}$ (28)

where $\mu_{i}$ is the class weight for class $i$; and $\beta$ is an
adjustment factor. The adjustment factor was optimized by the particle swarm
optimization to minimize the image quality metric.

####  2.4.4 Iterative stratification

When the dataset is small and highly imbalanced, it is inappropriate to
randomly split the dataset into train and test sets, because it is possible
that under-represented classes are missing in one of the sets and thus the
trained model fails to represent the whole dataset. We proposed to employ an
iterative stratification method for multi-label data [35]. The whole dataset
was divided into more than two subsets while the class weights of the sets
were the same.

####  2.4.5 Image augmentation

Augmentation was performed to increase the dataset size and improve the
generalizability of the semantic segmentation model [44]. The training dataset
was enlarged using four strategies: (1) Resizing: The images were resized from
256×256 pixels to 320×256 pixels. (2) Random cropping: The images were
randomly cropped to 256×256 pixels. (3) Random flipping: The images were
randomly flipped horizontally and vertically with a probability of 50%. (4)
Photometric distortion: The brightness, contrast, saturation, and hue of
images were randomly modified with a probability of 50%.

###  2.5 Performance metrics

Six performance metrics were used to evaluate the semantic segmentation model:

$Pixel\
accuracy=\frac{\sum_{k=1}^{K}{(TP_{k}+TN_{k})}}{\sum_{k=1}^{K}{(TP_{k}+FP_{k}+TN_{k}+FN_{k})}}$
(29) $Mean\
accuracy=\frac{1}{K}\sum_{k=1}^{K}\left(\frac{TP_{k}+TN_{k}}{TP_{k}+FP_{k}+TN_{k}+FN_{k}}\right)$
(30)
$mIoU=\frac{1}{K}\sum_{k=1}^{K}\left(\frac{TP_{k}}{TP_{k}+FP_{k}+FN_{k}}\right)$
(31)
$Precision=\frac{1}{K}\sum_{k=1}^{K}\left(\frac{TP_{k}}{TP_{k}+FP_{k}}\right)$
(32)
$Recall=\frac{1}{K}\sum_{k=1}^{K}\left(\frac{TP_{k}}{TP_{k}+FN_{k}}\right)$
(33) $F1=2\times\frac{Precision}{Precision+Recall}$ (34)

where $TP_{k}$ and $TN_{k}$ are respectively the true positives and true
negatives corresponding to class $k$; $FP_{k}$ and $FN_{k}$ are respectively
the false positives and false negatives corresponding to class $k$; $K$ is the
number of classes; and $mIoU$ is mean intersection-over-union.

##  3 Results

Four pairs of representative images before and after applying the optimized
adaptive gamma correction are compared in Figure 4. The comparison shows that
the optimized adaptive gamma correction method is capable of improving the
visibility of the microscope images of graphene regardless of the background
color and brightness. The texture and edges of the improved images become
easier to distinguish from the background. The optimal adjusted parameter
obtained by the particle swarm optimization performed to minimize the image
quality metric was 0.561. The lowest image quality metric value was 1.48. The
pixel intensity levels of the fourth pair of images are compared in Figure 5.
The comparison shows that the oversaturation and oversaturated regions are
greatly reduced by impeding the optimized adaptive gamma correction. An index
defined by Equation (35) is proposed to indicate the effect of the optimized
adaptive gamma correction for improvement of the images:

$\delta=\frac{A_{o}}{A_{T}}$ (35)

where $\delta$ is the oversaturated index; $A_{o}$ is the number of pixels
with an intensity level over 253; and $A_{T}$ is the total number of pixels.
The oversaturated indices of original and modified images for the red, green,
and blue channels are highly reduced as shown in Figure 5 (e) to Figure 5 (g).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4: The effect of optimized adaptive gamma correction on the microscope
images. The upper row are the original images, and the lower row of images are
the enhanced images by the optimized adaptive gamma correction.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 5: An investigation on an overexposed image: (a) the original image,
(b) the improved image, (c) the pixel intensity histogram of the red channel
of the original image, (d) the pixel intensity histogram of the red channel of
the improved image, (e) oversaturated index of the red channel, (f)
oversaturated index of the green channel, and (g) oversaturated index of the
blue channel. In (a) and (b), pixels with intensity levels over 253 represent
oversaturation and are shown in red color.

The proposed method (HRNet+OCR) was compared with three state-of-the-art
convolutional neural networks, which are UNet++ [45], PSPNet [46], and
DeepLabv3+ [47]. The proposed models were trained by stochastic gradient
descent. The learning rate was 0.1. The momentum was 0.9. The weight decay was
0.0005. The batch size was eight. The maximum number of iterations was 10,000.
The cross-entropy loss was used as the loss function. The performance metrics
of the different methods are compared in Table 2. The proposed method shows
the highest accuracy.

Table 2: Results of the performance metrics

Methods | Pixel accuracy | Mean accuracy | F1 | mIoU | Precision | Recall
---|---|---|---|---|---|---
The existing methods |  |  |  |  |  |
UNet++ | 96.7% | 61.7% | 56.0 | 56.4 | 60.9 | 51.9%
PSPNet | 95.5% | 64.6% | 56.0 | 56.5 | 59.5 | 52.8%
DeepLabv3+ | 96.9% | 70.9% | 62.8 | 58.5 | 66.7 | 59.3%
The proposed methods |  |  |  |  |  |
Baseline (HRNet+OCR) | 97.1% | 71.7% | 65.9 | 59.0 | 67.9 | 63.9%
Baseline+$Loss{}_{W}^{1}$ | 98.3% | 73.6% | 72.9 | 61.5 | 77.2 | 69.1%
Baseline+$Loss{}_{W}$+OAGC2 | 98.8% | 75.2% | 87.2 | 66.7 | 89.7 | 85.1%
Baseline+$Loss{}_{W}$+OAGC+WL3 | 99.0% | 81.6% | 91.8 | 71.7 | 93.5 | 90.2%

  * 1

"$Loss{}_{w}$" denotes weighted cross-entropy loss.

  * 2

"OAGC" denotes optimized adaptive gamma correction.

  * 3

"WL" denotes weak learning.

Further, four HRNet+OCR models were compared to test the effects of the
weighted cross-entropy loss, optimized adaptive gamma correction, and weak
learning on the performance. The results show that the accuracy is greatly
improved by applying the weighted cross-entropy loss, optimized adaptive gamma
correction, and weak learning. The pixel accuracy of the final model
(HRNet+OCR+$Loss{}_{W}$+OAGC+WL) is higher than 99%. The class-wise
intersection-over-union (IoU) and accuracy of the HRNet+OCR models are shown
in Figure 6. The IoU and accuracy of the highly under-represented classes 1L,
3L, and 6-10L are greatly increased by the weighted cross-entropy loss,
optimized adaptive gamma correction, and weak learning. The minimum IoU of the
final model is higher than 56%. Figure 7 shows representative images in the
test set and the corresponding prediction results.

(a)

(b)

Figure 6: The class-wise performance of the applied segmentation methods: (a)
Mean IoU (mIoU), (b) accuracy. $Loss{}_{W}$, OAGC, WL stand for the weighted
cross-entropy loss, optimized adaptive gamma correction, and weak learning,
respectively.

Microscope image |  Ground-truth label |  Predicted result |  Overlay
---|---|---|---

---

(a)

(b)

(c)

(d)

---

(e)

(f)

(g)

(h)

---

(i)

(j)

(k)

(l)

---

(m)

(n)

(o)

(p)

---

(q)

(r)

(s)

(t)

---

(u)

(v)

(w)

(x)

(y)

Figure 7: Representative optical microscope images and the predicted results
using the final model (HRNet+OCR+$Loss{}_{W}$+OAGC+WL).

##  4 Conclusions

For the first time, we developed a deep learning method for automatic
identification and classification of exfoliated graphene flakes from optical
microscope images with various microscope settings. Novel methods such as
optimized adaptive gamma correction are presented to improve the microscopy
images and address the challenges associated with the imbalanced classes,
limited dataset size, different background conditions, and different
resolutions of the microscopy images. The presented machine learning framework
based on object-contextual representations is promising for automatically
identifying 2D material flakes and classifying the thickness of the flakes
based on microscope images. The pixel accuracy of the trained model is higher
than 99%, outperforming the state-of-the-art deep convolutional neural
networks DeepLabv3+, PSPNet, and U-Net ++. The proposed method is robust to
the background colors of microscope images because the proposed machine
learning method can automatically identify the background color of the images
and categorize the images into appropriate groups. The proposed method is also
robust to the brightness and resolution of images because the optimized
adaptive gamma correction method can effectively improve the quality of
overexposed microscopy images. The presented machine learning framework is
capable of learning characterization from new images while preserving the
knowledge learned from old images. Weak learning method enables us to retrain
convolutional neural networks while preserving the learned knowledge by the
network. The method is useful especially when the training set is small but
there is a convolutional neural network trained on a similar task. The
presented method can handle imbalanced datasets. The IoU of different classes
is in the range of 57% and 99%, and the mIoU of all the classes is 59%. This
methodology shows that image augmentation, iterative stratification, and
weighted cross-entropy loss are able to significantly enhance the predictive
performance of convolutional neural networks on under-represented classes and
successfully identify the optical images with various backgrounds.

## Acknowledgement

This research was funded by National Science Foundation (grant number
CMMI-2046407).

[1] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, A. Kis, Single-
layer MoS2 transistors, Nature Nanotechnology, 6 (2011) 147-150.

[2] B. Huang, G. Clark, E. Navarro-Moratalla, D.R. Klein, R. Cheng, K.L.
Seyler, D. Zhong, E. Schmidgall, M.A. McGuire, D.H. Cobden, W. Yao, D. Xiao,
P. Jarillo-Herrero, X. Xu, Layer-dependent ferromagnetism in a van der Waals
crystal down to the monolayer limit, Nature, 546 (2017) 270-273.

[3] C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang,
Y. Wang, Z.Q. Qiu, R.J. Cava, S.G. Louie, J. Xia, X. Zhang, Discovery of
intrinsic ferromagnetism in two-dimensional van der Waals crystals, Nature,
546 (2017) 265-269.

[4] X. Xu, W. Yao, D. Xiao, T.F. Heinz, Spin and pseudospins in layered
transition metal dichalcogenides, Nature Physics, 10 (2014) 343-350.

[5] X. Xi, Z. Wang, W. Zhao, J.-H. Park, K.T. Law, H. Berger, L. Forró, J.
Shan, K.F. Mak, Ising pairing in superconducting NbSe2 atomic layers, Nature
Physics, 12 (2016) 139-143.

[6] X. Zhang, D. Sun, Y. Li, G.-H. Lee, X. Cui, D. Chenet, Y. You, T. F.
Heinz, and J. Hone, Measurement of lateral and interfacial thermal
conductivity of single- and bilayer MoS2 and MoSe2 using optothermal Raman
technique, ACS Applied Materials & Interfaces, 7 (2015) 25923-25929.

[7] Y. Li, F. Ye, J. Xu, W. Zhang, P. Feng, and X. Zhang, Gate-tuned
temperature in a hexagonal boron nitride-encapsulated 2D semiconductor
devices, IEEE Transactions on Electronic Devices, 65 (2018) 4068-4072.

[8] E. Easy, Y. Gao, Y. Wang, D. Yan, S.M. Goushehgir, E.-H. Yang, B. Xu, X.
Zhang, Experimental and computational investigation of layer-dependent thermal
conductivities and interfacial thermal conductance of one- to three-layer
WSe2, ACS Applied Materials & Interfaces, 13 (2021) 13063-13071.

[9] E. Easy, Y. Gao, Y. Wang, D. Yan, S.M. Goushehgir , E.-H. Yang , B. Xu, X.
Zhang , Thermal conductivities and interfacial thermal conductance of 1- to
3-layer WSe2, arXiv:2011.00730, (2020).

[10] F. Ye, Q. Liu, B. Xu, P.X.L. Feng, X Zhang, Very high interfacial thermal
conductance in fully hBN-encapsulated MoS2 van der Waals heterostructure,
ArXiv:2102.05239, (2021).

[11] Y. Li, A. Chernikov, X. Zhang, A. Rigosi, H. M. Hill, A. M. van der
Zande, D. A. Chenet, E.-M. Shih, J. Hone, and T. F. Heinz, Measurement of the
optical dielectric function of monolayer transition-metal dichalcogenides:
MoS2, MoSe2, WS2, WSe2, Physical Review B, 90 (2014) 205422.

[12] X. Zhang, Characterization of layer number of two-dimensional transition
metal diselenide semiconducting devices using Si-peak analysis, Advances in
Materials Science and Engineering, 2019 (2019) 7865698.

[13] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V.
Dubonos, I.V. Grigorieva, A.A. Firsov, Electric field effect in atomically
thin carbon films, Science, 306 (2004) 666.

[14] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevich, S.V.
Morozov, A.K. Geim, Two-dimensional atomic crystals, Proceedings of the
National Academy of Sciences of the United States of America, 102 (2005)
10451.

[15] H. Li, J. Wu, X. Huang, G. Lu, J. Yang, X. Lu, Q. Xiong, H. Zhang, Rapid
and reliable thickness identification of two-dimensional nanosheets using
optical microscopy, ACS Nano, 7 (2013) 10344-10353.

[16] X. Lin, Z. Si, W. Fu, J. Yang, S. Guo, Y. Cao, J. Zhang, X. Wang, P. Liu,
K. Jiang, W. Zhao, Intelligent identification of two-dimensional
nanostructures by machine-learning optical microscopy, Nano Research, 11
(2018) 6316-6324.

[17] S. Masubuchi, T. Machida, Classifying optical microscope images of
exfoliated graphene flakes by data-driven machine learning, npj 2D Materials
and Applications, 3 (2019) 4.

[18] Z.H. Ni, H.M. Wang, J. Kasim, H.M. Fan, T. Yu, Y.H. Wu, Y.P. Feng, Z.X.
Shen, Graphene thickness determination using reflection and contrast
spectroscopy, Nano Letters, 7 (2007) 2758-2763.

[19] C.M. Nolen, G. Denina, D. Teweldebrhan, B. Bhanu, A.A. Balandin, High-
throughput large-area automated identification and quality control of graphene
and few-layer graphene films, ACS Nano, 5 (2011) 914-922.

[20] P. Blake, E.W. Hill, A.H. Castro Neto, K.S. Novoselov, D. Jiang, R. Yang,
T.J. Booth, A.K. Geim, Making graphene visible, Applied Physics Letters, 91
(2007) 063124.

[21] E. Shelhamer, J. Long, T. Darrell, Fully convolutional networks for
semantic segmentation, IEEE Transactions on Pattern Analysis and Machine
Intelligence, 39 (2017) 640-651.

[22] L.-C. Chen, G. Papandreou, I. Kokkinos, K. Murphy, A.L. Yuille, DeepLab:
Semantic image segmentation with deep convolutional nets, atrous convolution,
and fully connected CRFs, IEEE Transactions on Pattern Analysis and Machine
Intelligence, 40 (2018) 834-848.

[23] O. Ronneberger, P. Fischer, T. Brox, U-Net: convolutional networks for
biomedical image segmentation, in: Medical Image Computing and Computer-
Assisted Intervention – MICCAI 2015, Springer International Publishing, Cham,
2015, pp. 234-241.

[24] V. Badrinarayanan, A. Kendall, R. Cipolla, SegNet: A deep convolutional
encoder-decoder architecture for image segmentation, IEEE Transactions on
Pattern Analysis and Machine Intelligence, 39 (2017) 2481-2495.

[25] S. Masubuchi, M. Morimoto, S. Morikawa, M. Onodera, Y. Asakawa, K.
Watanabe, T. Taniguchi, T. Machida, Autonomous robotic searching and assembly
of two-dimensional crystals to build van der Waals superlattices, Nature
Communications, 9 (2018) 1413.

[26] B. Han, Y. Lin, Y. Yang, N. Mao, W. Li, H. Wang, K. Yasuda, X. Wang, V.
Fatemi, L. Zhou, J.I.J. Wang, Q. Ma, Y. Cao, D. Rodan-Legrain, Y.-Q. Bie, E.
Navarro-Moratalla, D. Klein, D. MacNeill, S. Wu, H. Kitadai, X. Ling, P.
Jarillo-Herrero, J. Kong, J. Yin, T. Palacios, Deep-learning-enabled fast
optical identification and characterization of 2d materials, Advanced
Materials, 32 (2020) 2000953.

[27] Y. Saito, K. Shin, K. Terayama, S. Desai, M. Onga, Y. Nakagawa, Y.M.
Itahashi, Y. Iwasa, M. Yamada, K. Tsuda, Deep-learning-based quality filtering
of mechanically exfoliated 2D crystals, npj Computational Materials, 5 (2019)
124.

[28] E. Greplova, C. Gold, B. Kratochwil, T. Davatz, R. Pisoni, A. Kurzmann,
P. Rickhaus, M.H. Fischer, T. Ihn, S.D. Huber, Fully automated identification
of two-dimensional material samples, Physical Review Applied, 13 (2020)
064017.

[29] S. Masubuchi, E. Watanabe, Y. Seo, S. Okazaki, T. Sasagawa, K. Watanabe,
T. Taniguchi, T. Machida, Deep-learning-based image segmentation integrated
with optical microscopy for automatically searching for two-dimensional
materials, npj 2D Materials and Applications, 4 (2020) 3.

[30] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A.A.
Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska, D. Hassabis,
Overcoming catastrophic forgetting in neural networks, Proceedings of the
National Academy of Sciences, 114 (2017) 3521-3526.

[31] W. Chen, X. Mao, H. Ma, Low-contrast microscopic image enhancement based
on multi-technology fusion, 2010 IEEE International Conference on Intelligent
Computing and Intelligent Systems, IEEE, 2010, pp. 891-895.

[32] F.W. Leong, M. Brady, J.O.D. McGee, Correction of uneven illumination
(vignetting) in digital microscopy images, Journal of Clinical Pathology, 56
(2003) 619-621. [33] Z. Wu, Y. Gao, L. Li, J. Xue, Y. Li, Semantic
segmentation of high-resolution remote sensing images using fully
convolutional network with adaptive threshold, Connection Science, 31 (2019)
169-184.

[34] G. Cao, L. Huang, H. Tian, X. Huang, Y. Wang, R. Zhi, Contrast
enhancement of brightness-distorted images by improved adaptive gamma
correction, Computers & Electrical Engineering, 2018 (66) 569-582.

[35] K. Sechidis, G. Tsoumakas, I. Vlahavas, On the stratification of multi-
label data, Joint European Conference on Machine Learning and Knowledge
Discovery in Databases, Springer, 2011, pp. 145-158.

[37] R. Poli, J. Kennedy, T. Blackwell, Particle swarm optimization. Swarm
Intelligence, 2007, 1(1), pp.33-57.

[38] U. Shin, J. Park, G. Shim, F. Rameau, I.S. Kweon, Camera exposure control
for robust robot vision with noise-aware image quality assessment, 2019
IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS).

[39] D. Arthur, S. Vassilvitskii, k-means++: The advantages of careful
seeding, Stanford, 2006.

[40] Y. Yuan, X. Chen, J. Wang, Object-contextual representations for semantic
segmentation, in: A. Vedaldi, H. Bischof, T. Brox, J.-M. Frahm (Eds.) Computer
Vision – ECCV 2020, Springer International Publishing, Cham, 2020, pp.
173-190.

[41] J. Wang, K. Sun, T. Cheng, B. Jiang, C. Deng, Y. Zhao, D. Liu, Y. Mu, M.
Tan, X. Wang, W. Liu, Deep high-resolution representation learning for visual
recognition, IEEE Transactions on Pattern Analysis and Machine Intelligence,
2021, 43(10), 3349-3364.

[42] Z. Tian, L. Liu, B. Fei, Deep convolutional neural network for prostate
MR segmentation, Medical Imaging 2017: Image-Guided Procedures, Robotic
Interventions, and Modeling, International Society for Optics and Photonics,
2017, pp. 101351.

[43] X. Zhou, C. Yao, H. Wen, Y. Wang, S. Zhou, W. He, J. Liang, EAST: An
efficient and accurate scene text detector, 2017 IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2017, pp. 2642-2651.

[44] C. Shorten, T.M. Khoshgoftaar, A survey on image data augmentation for
deep learning, Journal of Big Data, 6 (2019) 60.

[45] Z. Zhou, M.M. Rahman Siddiquee, N. Tajbakhsh, J. Liang, UNet++: A nested
U-Net architecture for medical image segmentation, in: Deep Learning in
Medical Image Analysis and Multimodal Learning for Clinical Decision Support,
Springer International Publishing, Cham, 2018, pp. 3-11.

[46] H. Zhao, J. Shi, X. Qi, X. Wang, J. Jia, Pyramid scene parsing network,
Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition,
2017, pp. 2881-2890.

[47] L.-C. Chen, Y. Zhu, G. Papandreou, F. Schroff, H. Adam, Encoder-decoder
with atrous separable convolution for semantic image segmentation, Proceedings
of the European Conference on Computer Vision (ECCV), 2018, pp. 801-818.


††thanks: Ronit Chitre is with Aerospace Engineering, Indian Institute of
Technology, Bombay, Mumbai, India, and Arpita Sinha is with Systems and
Control Engineering, Indian Institute of Technology Bombay, Mumbai, India
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# An Enhanced RRT based Algorithm for Dynamic Path Planning and Energy
Management of a Mobile Robot

Ronit Chitre and Arpita Sinha

###### Abstract

Mobile robots often have limited battery life and need to recharge
periodically. This paper presents an RRT-based path-planning algorithm that
addresses battery power management. A path is generated continuously from the
robot’s current position to its recharging station. The robot decides if a
recharge is needed based on the energy required to travel on that path and the
robot’s current power. RRT* is used to generate the first path, and then
subsequent paths are made using information from previous trees. Finally, the
presented algorithm was compared with Extended Rate Random Tree (ERRT)
algorithm [4].

###### Index Terms:

RRT, path planning, robot energy management, autonomous systems

##  I Introduction

This paper addresses a path planning problem with a battery power management
system. Not all tasks are doable with one battery charge. We may need to
recharge the battery several times before completing the task. An example can
be a battery-powered autonomous robot harvesting a huge field. The robot must
know when to travel to a charging point in such cases. A conservative approach
will increase the total time to complete the job. The other extreme may lead
to the robot running out of battery before it reaches the recharge point. In
this paper, we develop an online algorithm that indicates the best time for
the robot to return for recharging.

We use RRT* to plan the path from the current robot location to the recharge
station to find the appropriate time to return to the base. However, the robot
follows a path to complete its task and the return-to-base path is taken only
when required. Since the return-to-base path planning needs to happen every
instant, RRT* generates a new tree every time. We propose to use the trees
built previously to reduce the computational time of RRT*. This concept is
similar to the RRT algorithms applied to a dynamic environment. However,
unlike the dynamic environment scenario where the robot moves on the RRT
algorithm’s path, here the robot follows a track to execute the task
independently of the RRT path generated.

We assume the robot knows the energy required to execute the task and to
travel back to the charging station. An application can be an onion harvester
robot, as shown in Fig. 1. The robot must harvest along the rows of onions
moving from one row to another. It also finds a path to the closest charging
station at each instant. Since the robot knows its current battery level and
can estimate the energy required for harvesting and traveling back to the
charging station, it can decide when to return.

Figure 1: Schematic of an onion field with onion harvester robot

We survey the literature in the next subsection, followed by the problem
formulation in Section II. An overview of RRT and RRT* is presented in Section
III. The proposed algorithm is explained in Section IV, and energy management
is addressed in Section IV-A. Simulation results are presented in Section V
followed by the concluding remarks in Section VI.

###  I-A Literature Survey

RRT and RRT* algorithms are studied extensively in the literature. Several
extensions are presented. These algorithms are sampling-based methods for
robotic path planning which are probabilistically complete. RRT generates a
path while RRT* gives an optimal path (in the limit number of nodes tend to
infinity) based on some cost function. Some of the advantages of RRT and its
variants include their applicability in complicated workspace space or
configuration space, its capability to include robot motion constraints, and
so on. An extensive survey on the extensions of RRT is available in [13] and
the references therein. We present the papers relevant to the work presented
in this paper.

In [1], the RRT* algorithm was used for replanning in a dynamic environment
with random, unpredictable moving obstacles. When an obstacle moves to a node
location that was included in the path, the path is replanned around the
obstacle using the node that is immediately after the node closest to the
obstacle. Authors in [2] and [9] also used a similar replanning approach to
planning a path around obstacles. Additionally, the work in [2] also limited
the number of nodes by removing childless nodes when the number of nodes
exceeded a limit which will reduce its complexity.

Paper [3] proposed the Dynamic Rate Random Forest or the DRRT. It, too,
discards nodes affected by an obstacle and reforms the tree. The Extended Rate
Random Tree algorithm (ERRT) mentioned in [4] proposed using waypoints that
are nodes from the tree generated in past iterations. When building the tree
either a random node is selected or a node from the waypoint array is selected
or the goal node is selected based on a certain probability distribution. Some
further improvements to this were made in [10] which used the waypoint cache
method and BG-RRT algorithm.

[5] developed a combination of RRT and $\text{A}^{*}$ in which multiple random
nodes and the corresponding nearest nodes are selected, but only the node that
minimizes a certain premade cost function is further examined. It also
experimented with using other types of norms apart from the Euclidean norm
while finding distance metrics. [12] combined the artificial potential field
method and RRT by using two trees that advance towards each other with respect
to an attractive or repulsive potential.

If $\text{RRT}^{*}$ takes too long to converge to a path, some other methods
can be used to remove redundant turns and optimize the path, like the ant
colony optimization algorithm that was used in [6]. [11] developed a new
approach to $\text{RRT}^{*}$, which involves first generating a path from an
initial point to a goal and then optimizing it by interconnecting directly
visible points and doing intelligent sampling. The robot used in this study
can not take sharp turns instantaneously and has limits on maximum and minimum
turning angles. These nonholonomic constraints need to be accounted for in the
path planning algorithm. Some ways to incorporate these constraints have been
discussed in [7] and [8].

_Contributions_ \- The problems addressed in the literature considered a
dynamic environment and different methods of replanning using RRT are
presented. We consider the case where the robot is following a pre-defined
path, and the RRT* is used to plan a path to the base for re-charging.
Therefore the starting point of the RRT* is changing at every instance. We
propose a fast algorithm called ”dynamic path RRT” that can find the cost to
reach the base so that the performance of the robot to carry out its assigned
task is not hampered. Hence, our problem setting is different from what exists
in the literature. We also compared the performance of our proposed algorithm
with some of the algorithms in the literature.

##  II Problem Formulation

We consider an autonomous robot executing some task in a known environment.
There are one or more charging stations. The robot follows a predefined path
until it requires recharging its batteries. The robot plans an online path to
the charging station using RRT*. We assume the power needed to move at a
constant speed is known, and the robot moves at a fixed desired speed. So, the
robot knows the energy it will need to return to the charging point. The robot
can also measure the energy remaining in the batteries.

To model the robot, we choose a classic car model as

$\displaystyle\frac{d\mathbf{x}}{dt}$ $\displaystyle=\mathbf{v}$ (1)
$\displaystyle\frac{d\psi}{dt}$ $\displaystyle=\omega=\frac{v}{L}\tan\delta$
(2)

where $\mathbf{x}$ represents position in the $xy$ plane, $\mathbf{v}$
represents the velocity in the $xy$ plane, $v=\|\mathbf{v}\|$ is its speed,
$\psi$ represents the heading and $\omega$ represents the angular velocity of
the robot, $L$ is the length of the robot and $\delta$ is its steering angle.
We assume $v$ is fixed, so the input to the robot is $\omega$.

$\alpha$$\mathbf{v}$Initial pointFinal point$R$$\theta$$d$ Figure 2:
Trajectory of the robot from one point to another

We use simple geometry to make the robot go from current to new positions.
Consider Fig. 2. Here, $\theta=2\alpha$ and $R=\frac{d}{2\sin\alpha}$. The
linear and angular velocities are related by

$\displaystyle\omega=\frac{v}{R}=\frac{2v\sin\alpha}{d}$ (3)

where both $\alpha$ and $d$ are measurable. Since there exist limits on the
steering angle $\delta\leq\delta_{\text{max}}$, $\omega$ will be restricted to

$\displaystyle\omega\leq\left|\frac{v}{L}\tan(\delta_{\text{max}})\right|$ (4)

So, if (3) demands an $\omega$ outside the above limits, the final point will
not be reachable. We associate a cost with the path from the current state to
the next state. The cost is equal to the energy required to travel the path.
For simplicity, we assume that the energy spent is proportional to the
distance covered. However, formulations can also be used. Therefore, the cost
($E$) of the path is

$\displaystyle E=kR\theta=2kR\alpha$ (5)

where $k$ is the gain relating distance to cost.

##  III Overview of RRT and $\text{RRT}^{*}$

RRT is a well-known probabilistic algorithm used for path planning in
robotics. $\text{RRT}^{*}$ is a modification of RRT that can generate an
optimal path. We give an overview of these algorithms to relate to our
proposed algorithm and for the completeness of the paper. A pseudocode for the
RRT algorithm is given below.

Algorithm 1 RRT algorithm

procedure RRT($q_{i},q_{g}$)

$T\leftarrow$ initialize_tree($q_{i}$, $q_{g})$

while goal_not_reached do

$q_{\text{rand}}\leftarrow$random_sample()

neighbourhood,
$q_{\text{nearest}}\leftarrow$get_nearest_node($q_{\text{rand}}$)

$q_{\text{new}}\leftarrow$ steer($q_{\text{nearest}}$, $q_{\text{random}}$)

if is_obstacle_free($q_{\text{nearest}}$, $q_{\text{new}}$) then

cost($q_{\text{new}}$) = cost($q_{\text{nearest}}$) +
energy($q_{\text{nearest}}$, $q_{\text{new}}$)

T$\leftarrow$ insert_node($q_{\text{nearest}}$, $q_{\text{new}}$)

if $q_{\text{new}}=q_{g}$ then

goal_not_reached = False

end if

end if

end while

end procedure

First, the initial and final nodes are initialized i.e. $q_{i}$ and $q_{g}$. A
random node $q_{\text{rand}}$ is selected from the workspace. Then by checking
all the nodes in the tree, the node closest to the random node is found:
$q_{\text{nearest}}$. Also, the random node is moved towards the nearest node
up to a certain step size. Then in the steer() function angular velocity to go
from the nearest node to the random node is computed, and if it lies outside
the bounds set by turning rate constraints, it is set to
$\text{sign}(\omega)\omega_{\text{max}}$. The coordinates of the new node can
then be fixed by moving from the nearest node with the constant $\omega$ that
is calculated. The minimum distance is the RRT step size and the maximum
distance is the radius of the ball used for finding the neighborhood. The
pseudocode for this is given in algorithm 1. If the path from the nearest node
to the new node intersects an obstacle then a different random node is picked.
In this work we consider obstacles to be line segments thus, it is easy to
check if the path is intersecting the obstacle. Before adding the new node to
the tree its cost is computed by adding the cost of its parent and the energy
required to go from the nearest node to the new node. The new node now becomes
a child node of $q_{\text{nearest}}$. If this new node happens to be the goal
then RRT is ended and a path is generated.

There is a very low probability of the goal node being selected while drawing
random samples. Thus instead, the algorithm is stopped when a node becomes
sufficiently close to the goal node. Another alternative [4] is to modify
random_sample() such that a completely random node is drawn with a probability
$1-p_{\text{goal}}$ and the goal node is drawn with probability
$p_{\text{goal}}$

$P(q)=\begin{cases}1-p_{\text{goal}},&\text{if }q\text{ is a random node}\\\
p_{\text{goal}},&\text{if }q\text{ is the goal node}\end{cases}$ (6)

RRT algorithm will converge on a path from the initial point to the goal.
However, the path may or may not be the optimal path. $\text{RRT}^{*}$
modifies this code to find the optimal path with two changes. First, in the
steer function, after the random node is brought within a certain step size of
the nearest node, all the nodes lying in a neighborhood around the random node
of a size greater than the step size are selected and put into an array. Thus
the nearest node will always be in this neighborhood. Instead of taking the
nearest node as a parent, the node with the least cost in that neighborhood is
chosen as a parent. Secondly, the nodes are ”rewired” that is after a new node
has been added it is checked if it can be connected with a neighboring node to
reduce the cost of that node.

##  IV Dynamic Replanning

We assume that the map of the environment is available to the robot. For
example, the onion-harvesting robot can map the field while harvesting since
the path to the base will be only through the harvested region. A suboptimal
path to the goal is generated as the robot moves, and the energy required to
transverse this path is computed. This energy is compared with the current
battery level of the robot. If sufficient energy is available, the robot
continues to harvest.Else the plucking is stopped, and an optimal path is
computed using $\text{RRT}^{*}$. The energy required for the new path is used
to make the decision on return to base.

We propose a new algorithm that relies on trees built in past iterations.
Consider a tree already generated by the algorithm, and a path is found. We
call this the old path and the tree being built in the current iteration is
called new tree. Now the initial node is placed in the new position of the
robot. As new nodes are added to the new tree, it checks if there are any
nodes from the old path nearby. If so, the new tree is built further by
replicating nodes from the old path.

Algorithm 2 Dynamic Path RRT Algorithm

procedure Dynamic_Path_RRT($q\\_i,q\\_g$)

$T\leftarrow$ initialize_tree($q_{i}$, $q_{g}$)

$F\leftarrow$ initialize_forest($T$)

while goal_not_reached do

$q_{\text{rand}}\leftarrow$random_sample()

$q_{\text{nearest}}\leftarrow$get_nearest_node($q_{\text{rand}}$)

$q_{\text{new}}\leftarrow$ steer($q_{\text{nearest}}$, $q_{\text{random}}$)

if is_obstacle_free($q_{\text{nearest}}$, $q_{\text{new}}$) then

cost($q_{\text{new}}$) = cost($q_{\text{nearest}}$) +
energy($q_{\text{nearest}}$, $q_{\text{new}}$)

T$\leftarrow$ insert_node($q_{\text{nearest}}$, $q_{\text{new}}$)

p $\leftarrow$ UniformRandom(0, 1)

if p $<$ $\text{p}_{\text{scan}}$ then

old_nodes $\leftarrow$ scan_forest(F, $q_{\text{new}}$)

check_connection(old_nodes, $q_{\text{new}}$)

end if

if goal_reached(T) then

goal_not_reached = False

end if

end if

end while

end procedure

Here it is important to introduce a new data structure - forest. It is the
collection of the paths generated in past iterations. Not all paths need to be
remembered since this might lead to excessive memory requirements. We propose
to discard the oldest path from the forest when a new path is added. The
scan_forest function checks if any of the nodes from the older iterations are
present in a neighbourhood around the new node.

The check connection function iterates through all the nodes gathered from
scan forest function and finds the node that can be attached to the newly
added node with minimum cost. It then calls the path building function. The
pseudocode for this is in algorithm 3.

Algorithm 3 Check Connection and Build Path

procedure check_connection(old_nodes, F, $q_{\text{new}}$)

candidate_nodes[] $\leftarrow$ initialize_array

for $q_{\text{old path}}$ in old_nodes do

$q_{\text{old path}}\leftarrow$ steer($q_{\text{new}}$, $q_{\text{old path}}$)

cost($q_{\text{old path}}$) = cost($q_{\text{new}}$) +
energy($q_{\text{new}}$, $q_{\text{old path}}$)

candidate_nodes.add($q_{\text{old path}}$)

end for

$q_{\text{old path}}$ = $\min_{\text{cost}}$candidate_nodes

path_building($q_{\text{new}}$, $q_{\text{old path}}$)

end procedure

The path_building function takes in $q_{\text{new}}$ i.e. the latest node
added to the new tree and $q_{\text{old path}}$ that is the minimum cost node
returned by ’check connection’. Path building attaches $q_{\text{old path}}$
to the new tree and then finds its child node which is also part of the old
path. It then attaches that child node to the new tree and this goes on in a
recursive process. The pseudocode for this is given in algorithm 4.

Algorithm 4 Path Building from Older Paths

procedure path_building($q_{\text{new tree}}$, $q_{\text{old path}}$)

new_tree $\leftarrow$ insert_node($q_{\text{new tree}}$, $q_{\text{old
path}}$)

for child_node in $q_{\text{old path}}$.children do

if child_node.part_of_old_path is True then

new_child_node $\leftarrow$ steer($q_{\text{old path}}$, child_node)

end if

path_building($q_{\text{old tree}}$, new_child_node)

end for

end procedure

###  IV-A Robot Energy Management

Now that the dynamic path RRT algorithm is ready the procedure to determine
robot’s decision to go to the recharge point needs to be fixed.

Algorithm 5 Energy management algorithm

procedure decide_return($\mathbf{x}$, F)

latest_path $\leftarrow$ dynamic_path_RRT($\mathbf{x}$, $q_{g}$)

if latest_path.cost $>$ safety_factor $\times$ $\mathbf{x}.\text{power}$ then

latest_path
$\leftarrow\text{RRT}^{*}(\mathbf{x},q_{g},\mathbf{x}.\text{power})$

if latest_path.cost $>$ safety_factor $\times$ $\mathbf{x}.\text{power}$ then

latest_path$\leftarrow$execute($\mathbf{x}$)

end if

end if

$\mathbf{x}$ += robot_velocity

end procedure

The pseudocode for this is given in algorithm 5. Here $\mathbf{x}$ denotes the
state of the robot i.e. its position, velocity, angle, and charge. Dynamic
path RRT is first used to generate a suboptimal path. Then it is checked if
the current battery energy is enough to execute that path with a pre-
determined safety factor ($<1$). If yes, then the robot moves to its new
position. If not it runs RRT* until a path with low enough cost is found or
until a maximum node limit is reached. If RRT* was unsuccessful the robot
returns to the charging station.

##  V Simulation Results

We simulated the algorithm in Python on a computer with an Intel CORE i7
processor and Ubuntu 22.04 OS. Two different environments were considered. We
compared the execution time of our algorithm with that of ERRT[4] and regular
RRT. In the simulations, we assumed the forward speed of the robot to be 1
unit and a maximum steering angle of $40^{\circ}$. The probability of scanning
for old nodes was 0.7 and the probability for the sampling goal was 0.2. We
used a step size of 0.05 unit for RRT and a waypoint sampling probability of
0.7 for ERRT. Rewiring was done in the first iteration when an initial tree
was generated but not done later.

Figure 3: First tree formed using pure $\text{RRT}^{*}$ Figure 4: Decrease in
number of nodes after using the dynamic path planner Figure 5: Average time
taken for different y coordinates

In the first scenario, the robot moves parallel to the $y$-coordinate. The
snapshot of the path generated at the initial time and some intermediate time
is shown in Figs. 3-4. The dynamic path RRT algorithm performed much better
than ERRT and regular RRT. Please refer to Fig. 5. It is worth mentioning that
in the time analysis of dynamic path RRT, the very first dot in the plot at
$y=0$ represents pure RRT*. That is why the time taken for the first iteration
is much higher than all others. Planning a trajectory by using information
from older iterations is far more efficient than using pure RRT in each case.

Figure 6: Average cost of paths vs y coordinate

Figure 6 shows how the cost of each path varies as the position of the robot
changes. This was done for different values of scan forest probability. This
shows initially, the cost of RRT* is similar to the cost of the dynamic path
planner. The costs of the path generated by varying $p_{\text{scan}}$ in the
algorithm are also similar. However, as the robot moves farther away, the
algorithm gives paths with higher costs. Thus as the robot moves further away
from its starting point, running an RRT* intermittently will help.

The energy management algorithm was also tested on a field with dimensions
$6\times 6$ and no obstacles. The robot goes from $(3,0)$ to $(3,3)$ and then
takes a 90-degree turn and starts moving towards $(0,3)$. There are two
recharging stations at $(1,5)$ and $(5,1)$ respectively. The snapshots at
different instances are shown in Fig. 7.

Figure 7: Top-Left: Initial path made by $\text{RRT}^{*}$; Top-Right: The path
made by the robot can not be used with its current battery thus it makes a new
path using $\text{RRT}^{*}$; Bottom-Left: New path made using the dynamic path
RRT; Bottom-Right: The same situation happens here too but this time it finds
a path to a different recharging station. It now considers this as its goal
point.

##  VI Conclusion

This paper presents an innovative approach that addressed the critical
challenges of path planning and energy management in autonomous mobile robots.
The dynamic path RRT algorithm is proposed that generates paths from the
robot’s current position to the recharge station efficiently by using
information from old iterations. It is then checked if the robot’s current
battery is sufficient to safely execute this path and a decision is taken
based on this condition. We compared the proposed algorithm with ERRT as well
as regular RRT and found that the proposed strategy is significantly faster.
However, the cost of the path is usually more than RRT*. We also analyzed the
proposed algorithm by varying the probability of building the connection with
the older paths. It is observed that the new paths are found faster as the
probability is increased, but the path cost also increases. A user can decide
on the probability based on the trade-off between time and cost.

## References

  * [1] Devin Connell and Hung Manh La, ”Dynamic Path Planning and Replanning for Mobile Robots using RRT*”, 2017 IEEE International Conference on Systems, Man, and Cybernetics (SMC), Banff, AB, Canada, 2017, pp. 1429-1434, doi: 10.1109/SMC.2017.8122814.
  * [2] Adiyatov and H. A. Varol, ”A novel RRT*-based algorithm for motion planning in Dynamic environments,” 2017 IEEE International Conference on Mechatronics and Automation (ICMA), Takamatsu, Japan, 2017, pp. 1416-1421, doi: 10.1109/ICMA.2017.8016024.
  * [3] D. Ferguson, N. Kalra and A. Stentz, ”Replanning with RRTs,” Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006., Orlando, FL, USA, 2006, pp. 1243-1248, doi: 10.1109/ROBOT.2006.1641879.
  * [4] Bruce, James and Manuela M. Veloso. “Real-time randomized path planning for robot navigation.” IEEE/RSJ International Conference on Intelligent Robots and Systems 3 (2002): 2383-2388 vol.3.
  * [5] J. Li, S. Liu, B. Zhang and X. Zhao, ”RRT-A* Motion planning algorithm for non-holonomic mobile robot,” 2014 Proceedings of the SICE Annual Conference (SICE), Sapporo, Japan, 2014, pp. 1833-1838, doi: 10.1109/SICE.2014.6935304.
  * [6] J. Qi, H. Yang and H. Sun, ”MOD-RRT*: A Sampling-Based Algorithm for Robot Path Planning in Dynamic Environment,” in IEEE Transactions on Industrial Electronics, vol. 68, no. 8, pp. 7244-7251, Aug. 2021, doi: 10.1109/TIE.2020.2998740.
  * [7] Khanal, Abhishek. “RRT and RRT* Using Vehicle Dynamics.” ArXiv abs/2206.10533 (2022): n. pag.
  * [8] Heß, Robin. (2013). Trajectory Planning for Car-Like Robots Using Rapidly Exploring Random Trees. IFAC Proceedings Volumes. 46. 44-49. 10.3182/20131111-3-KR-2043.00018.
  * [9] E. Shan, B. Dai, J. Song and Z. Sun, ”A Dynamic RRT Path Planning Algorithm Based on B-Spline,” 2009 Second International Symposium on Computational Intelligence and Design, Changsha, China, 2009, pp. 25-29, doi: 10.1109/ISCID.2009.155.
  * [10] Chengren Yuan, Guifeng Liu, Wenqun Zhang, Xinglong Pan, An efficient RRT cache method in dynamic environments for path planning, Robotics and Autonomous Systems, Volume 131, 2020, 103595, ISSN 0921-8890,
  * [11] F. Islam, J. Nasir, U. Malik, Y. Ayaz and O. Hasan, ”RRT$\ast$-Smart: Rapid convergence implementation of RRT$\ast$ towards optimal solution,” 2012 IEEE International Conference on Mechatronics and Automation, Chengdu, China, 2012, pp. 1651-1656, doi: 10.1109/ICMA.2012.6284384.
  * [12] W. Xinyu, L. Xiaojuan, G. Yong, S. Jiadong and W. Rui, ”Bidirectional Potential Guided RRT* for Motion Planning,” in IEEE Access, vol. 7, pp. 95046-95057, 2019, doi: 10.1109/ACCESS.2019.2928846.
  * [13] Iram Noreen, Amna Khan and Zulfiqar Habib, “Optimal Path Planning using RRT* based Approaches: A Survey and Future Directions” International Journal of Advanced Computer Science and Applications(IJACSA), 7(11), 2016. http://dx.doi.org/10.14569/IJACSA.2016.071114