Exploring the Opencv Library for Image Processing in Long-Pulse Thermography

Abstract

The objective of the present paper is to propose a framework, utilizing the OpenCV library, for post-processing infrared images obtained using Long-Pulse Thermography (LPT), with the goal of segmenting the images into defective and sound areas. A series of thermograms of a carbon fiber/epoxy specimen, containing precision milled flat-bottom holes, was acquired using an LPT system comprised of an uncooled microbolometer imager and halogen optical sources. Flaw detectability and planar size estimation were used to evaluate the results obtained with the proposed post-processing framework, in comparison to raw images, and images subjected only to pre-processing algorithms.

Keywords:
Thermography; Nondestructive Testing; Composite Materials; Image Processing; OpenCV

1. Introduction

Infrared Thermography (IRT) can be considered one of the most promising Nondestructive Testing (NDT) methods for the inspection of composite structures, mainly due to its non-contact nature, high inspection rate and easy data interpretation^1,2. IRT for NDT applications, or simply Thermal NDT (TNDT)³, can be classified as passive or active, depending on the available control of the thermal radiation source. In passive mode, the object under analysis is naturally at a temperature higher or lower than the environment, creating a measurable thermal contrast on its surface. On the other hand, in active mode, the object is thermally excited by an external source and its thermal response is analyzed. Stimulation sources include optical radiation, ultrasonic wave propagation, eddy current and others¹. Particular attention will be given here to optical TNDT.

The most common optically stimulated thermography methods applied to aerospace structures are Pulsed Thermography (PT) and Lock-in Thermography (LT)². PT is carried out by providing a heat pulse to the object under evaluation and monitoring its surface temperature under the transient cooling or heating phases. Hence, the technique is also referred to as Transient Thermography (TT). Heating is usually done with flash lamps but can also be carried out with halogen lamps. In the latter case, the heating requires a longer time, and the technique is called Long Pulse Thermography (LPT) or Step Pulse Thermography (SPT). LT, in the other hand, consists of heating the object with a modulated heating source, usually halogen lamps. An interesting parallel can be traced between these techniques and a vibration test. PT would be equivalent to a hammer impact test, and the verification of the transient behavior of the inspected component, while LT would be equivalent to the periodic excitation of a shaker, with evaluation of the steady-state vibration response of the component. Obviously, in this comparison, the mechanical waves are replaced by heat waves.

Thermograms obtained using optical stimulation are often contaminated with different noise sources, such as reflected infrared energy from the surroundings, non-homogeneous heating, and variations in surface properties (emissivity) of the specimen¹. These noise sources influence flaw detectability, since they produce abnormal thermal patterns in the infrared images. In order to enhance the contrast between sound and defective areas, a number of signal processing algorithms have been developed^1,3. The most widely used algorithms are grouped in Figure 1.

Figure 1
Most widely used image processing algorithms in TNDT³.

Data from a typical TT analysis consists of a sequence of infrared images (thermograms) that display the surface temperature distribution and evolution in time. This allows for distinct data processing strategies, namely: one-dimensional (temporal), two-dimensional (spatial) or a combination of both. 1D algorithms are applied on a pixel-base and analyze the pixel temperature evolution (image sequence). On the other hand, 2D algorithms are applied to a single image and operate on a pixel-to-pixel basis (e.g., contrast adjustment, filtering, etc.). Most TNDT processing algorithms rely on the pixel temporal evaluation, since more information about defect parameters, such as defect depth, can be obtained³. However, 2D algorithms are fundamental for image segmentation and feature extraction, which are important steps in the implementation of automated image processing algorithms employing machine learning strategies⁴.

Bearing this in mind, the objective of the present paper is to explore the Open Computer Vision (OpenCV) library for spatial processing of thermograms obtained using LPT. The general image processing framework consists of: a) pre-processing the thermograms utilizing 1D algorithms: Thermal Signal Reconstruction (TSR) and Pulsed Phased Thermography (PPT); and b) post-processing using OpenCV: spatial filtering, contrast enhancement and segmentation. The goal of the post-processing step is to segment the image into defective and sound areas, removing image background and easing defect sizing.

2. Image Processing

This section describes the pre-processing and post-processing algorithms applied to the thermograms.

2.1. Pre-Processing Algorithms

The basis of the TSR method is the reconstruction, or regression, of thermal responses as a function of time by low-order polynomials. The technique was proposed by Shepard⁶ and quickly became widely used, for its simplicity and for bringing important advantages over the analysis of raw data, such as noise reduction and contrast enhancement. The cooling of an intact area is modeled considering the one-dimensional solution of the conduction equation:

$\ln \left(T-T_0\right)=\ln \left(\frac{Q}{e}\right)-\frac{1}{2}\ln \left(\pi t\right)$ (1)

Where $T$ is the variable temperature,$T_0$ is the initial temperature, $Q$ is the applied heat energy, $e$ is the thermal effusivity of the material and $t$ is the time. From Equation 1, it is possible to model the temperature evolution of a sound area as a fixed, straight line with slope equal to -0.5. This behavior is independent of material properties and applied heat flux. A defective area will present a temperature evolution that diverges from this linear response. The polynomial function used to fit the temperature evolution data can be written as:

$\ln \left(T-T_0\right)=a_0+a_1\ln \left(t\right)+a_2{\left[\ln \left(t\right)\right]}^2+\dots +a_m{\left[\ln \left(t\right)\right]}^m$ (2)

The number of coefficients used to fit the thermal data plays an important role in quality of the processed images⁵. A fourth order polynomial was used in the present study⁶.

One of the main advantages of the technique is the ability to generate images of the derivatives of the thermal signal (first and second derivatives). The derivative images significantly reduce non-uniform heating effect and provide earlier indications of thermal contrast^3,7.

The PPT technique is based on the decomposition of the temporal thermal response in the frequency domain by means of a Fourier series expansion. Originally proposed by Maldague and Marinetti⁸, the technique represents the parallel between pulsed thermography and lock-in thermography⁹. As previously discussed, pulsed thermography makes use of square thermal pulses, with different durations depending on the type of lamp used and the purpose of the analysis. Mathematically, it is possible to decompose such a pulse as the combination of sine waves. Therefore, it can be said that pulsed thermography consists of the simultaneous analysis of thermal signals with varying amplitude and frequencies during the transient phase (heating or cooling), while lock-in thermography consists of the individual analysis of a single thermal signal, with given frequency and amplitude, during steady state (thermal equilibrium)².

According to the superposition principle, a time-domain response can be decomposed into a frequency-domain response using a Fourier expansion⁷:

$T\left(x,t\right)=\sum a\left({\omega }_n\right)T\left(z,{\omega }_n,t\right)$ (3)

Where $T\left(z,{\omega }_n,t\right)$ is a plane thermal wave of angular frequency ${\omega }_n$ propagating in the $z$ direction and $a\left({\omega }_n\right)$ is a measure of the strength of this component at the given frequency. The transformation between the time and frequency domain can be performed using the one-dimensional discrete Fourier transform (DDT):

$F_n=\Delta t{\displaystyle \sum }_{k=0}^{N-1}\left(k\Delta t\right)e^{\frac{-j2\pi k}{N}}=R{e}_n+I{m}_n$ (4)

The amplitude and phase of the transform can be estimated from Equation 4:

$A_n=\sqrt{R{e}_n{}^2+I{m}_n{}^2}$ (5)

${\varphi }_n={\tan }^{-1}\left(\frac{I{m}_n}{R{e}_n}\right)$ (6)

The components given by Equation 5 and Equation 6 are used to reconstruct the thermal image sequence. The amplitude image is proportional to local optical and infrared surface features, while the phase image, or phasegram, is related to the heat propagation time delay. Phasegrams are independent of surface and optical features, being particularly interesting for subsurface defect characterization⁸.

2.2. OpenCV

The OpenCV library was originally developed by Intel aimed at real-time computer vision applications. The library has over 2500 algorithms implemented in C++ language, with wrappers and bindings for many other languages, such as Python, MATLAB and Java. Many of the implemented algorithms have been extensively optimized since the first release of the library in 1999 and support models from deep learning frameworks such as Torch and TensorFlow. This makes the library extremely powerful for the development of automated image processing software. Another appealing feature of the library is that it is cross-platform and free for use under the open-source BSD license¹⁰.

The algorithms explored in the present study are part of the Image Processing, Histogram and Feature Detection modules of OpenCV. The proposed framework is depicted in Figure 2. It is noteworthy to mention that the proposed framework will be employed to a Region of Interest (ROI) covering the complete laminate. An alternative strategy would be to delimit smaller ROIs at different areas of the laminate in order to account for local temperature variations and noise, possibly yielding better segmentation results. The choice a broader ROI, covering the complete structure, aims to eliminate the necessity of the user to subjectively specify localized ROIs of the image for analysis.

Figure 2
Image post-processing framework.

Linear filters are the most common type of neighborhood operators⁴. The basic idea of a linear filter is to compute a new output value of a pixel based on the weighted combination of pixel values in its neighborhood. The entries in the weight kernel (mask) are often called filter coefficients. In the present work, smoothing filters (e.g., Gaussian blur), or low-pass kernels, were used to further reduce high-frequency noise from the thermograms. This can aid contrast enhancement, avoiding noise propagation.

The proposed framework explores the use of Contrast Limited Adaptive Histogram Equalization (CLAHE). CLAHE is a variant of ordinary Adaptive Histogram Equalization (AHE). AHE is used to improve the contrast of images through histogram equalization of several distinct sections of the image. This is achieved by transforming each pixel with a transformation function computed from a neighborhood region, as shown in Figure 3. In comparison to ordinary histogram equalization, AHE has the capability of improving local contrast and enhancing definitions of edges even if the image presents nonuniform lighting. However, conventional AHE has the tendency to amplify noise in relatively homogeneous regions of an image. CLAHE prevents such problem by limiting the contrast amplification factor based on the slope of the transformation function. In the implementation of CLAHE, two important parameters must be set: clip limit and grid size. The clip limit is the threshold for contrast limiting and the grid size is the size of the area used for histogram equalization.

Figure 3
Schematic of CLAHE applied to an image illustrating the different sections (grids) used to compute the transformation functions.

After contrast enhancement, the thermograms were segmented using two approaches: adaptive threshold and watershed. Among the adaptive threshold algorithms implemented in OpenCV, Otsu’s method¹¹ was chosen. It is a nonparametric automatic threshold selection algorithm for image segmentation, extensively used in document image binarization¹². The algorithm returns a single intensity threshold used to segment the image in two classes: foreground and background. This approach eliminates, to a certain degree, the subjectivity of a threshold determination by the user. This is achieved by minimizing the intra-class variance, defined as a weighted sum of variances of the two classes¹¹:

${\sigma }_w^2\left(k\right)={\omega }_0\left(k\right){\sigma }_0^2\left(k\right)+{\omega }_1\left(k\right){\sigma }_1^2\left(k\right)$ (7)

Where ${\omega }_0$ and ${\omega }_1$ are the probabilities of the two classes separated by a threshold $k$ and ${\sigma }_0^2$ and ${\sigma }_1^2$ are the variances of the two classes.

In addition to the Otsu’s method, the thermograms were also segmented using the watershed algorithm. The method has its origins in mathematical morphology, with its general concept introduced by Digabel and Lantuejoul¹³. The idea behind the watershed segmentation is the representation of a gray scale image as a topographic landscape, with ridges and valleys. The elevation values are given by the pixel intensity or gradient magnitude. By placing a water source in each regional minimum (valley), the topography is flooded until different sources meet (See Figure 4). At these contact regions, watersheds are placed, forming delimitations. OpenCV employs a marker-based algorithm, which reduces problems related to over segmentation by specifying locations of know background (sure background) and foreground (sure foreground)¹⁴. These locations are determined by performing morphological erosion and dilation, followed by distance transform.

Figure 4
One dimensional representation of the watershed transform, based on the flooding from sources at local minima.

The contour localization step was performed using a border following algorithm proposed by Suzuki and Abe¹⁵, also implemented in OpenCV. After computing the contours, parameters such as area and equivalent diameter can be directly calculated.

The different image segmentations strategies were evaluated by comparing the number of identified flat-bottom holes, as well as their planar size computed by each algorithm, in comparison to measurements obtained directly from the pre-processed images, manually measured using the ImageJ software¹⁶. Since the measurement is performed manually, an average of 3 repetitions was considered. Additionally, the effect of the histogram equalization is analyzed by means of absolute contrast comparison and segmentation results (i.e., detectability, defect geometry and planar size estimation). The absolute contrast is calculated using Equation 8 ²:

$C=\left|S_d-S_i\right|$ (8)

Where $S_d$ is the average pixel intensity at the defect, $S_i$ the average pixel intensity at a sound reference area in the vicinity of the defect. Pixel intensity is given by a single channel, 8-bit value from 0 to 255. Defective and sound pixels were selected delimiting a circular area in the center of the defect and an outer ring, with thickness of 5 pixels, surrounding the defect, respectively.

3. Materials and Methods

Thermograms of a carbon fiber/epoxy specimen, containing precision milled flat-bottom holes, were acquired using a LPT system comprised of an uncooled microbolometer imager operating in the Long Wave Infrared (LWIR) spectral band, with resolution 640x480 pixels and sensitivity of 70 mK, and halogen optical sources with a total power of 2kW (See Figure 5). The specimen was manufactured following recommendations of ASTM 2582-19¹⁷, with a total of 25 flat-bottom holes of varying diameter (20, 15, 10, 5 and 2 mm) and depth (0.5, 1.0, 1.25, 1.5 and 2 mm), resulting in simulated flaws with Aspect Ratios (AR) – ratio between diameter and depth – ranging from 40 to 4. The total thickness of the laminate is 2.5 mm.

Figure 5
a) Hardware and test setup b) Back surface of the specimen with flat-bottom holes with diameters varying from 2-20 mm and depths of 0.5-2 mm.

The specimen was placed at 50 cm from the LPT system and heated for a total of 10 seconds, followed by monitoring for 10 seconds during free convection cooling. In trial runs, this configuration presented the best contrast results. An in-house software was developed in C++ for controlling the excitation sources and retrieving the thermal data from the imager. Pre-processing and post-processing were also carried out using the developed software package.

4. Results

In order to analyze the effect of the pre-processing algorithms, thermograms in the early stages of cooling were processed and compared to the raw data (See Figure 6). It is possible to see how the pre-processing algorithms are capable of increasing defect contrast. Since the earlier thermograms represent lower frequencies, the PPT is capable of revealing deeper defects when compared to the TSR. However, the TSR first derivative yields higher contrast between defective and sound area when the maximum contrast frame is analyzed. Table 1 presents the value of maximum absolute contrast for the 20 mm flat-bottom holes. The results show that the TSR algorithm outperforms the PPT, with significant improvement when compared with the raw thermograms.

Figure 6
Comparison of raw and pre-processed thermograms, for the maximum absolute contrast frame.

Table 1
Absolute contrast values for the 20 mm flat-bottom holes, for the raw and pre-processed thermograms.

	AR = 40	AR = 20	AR = 16	AR = 13.33	AR = 10
Raw	127.80	88.30	52.75	11.18	4.76
TSR 1st Der.	182.89	132.13	77.97	55.13	13.95
PPT Amp.	172.86	104.10	47.89	35.77	16.22

Figure 7 shows the evolution of the absolute contrast for the 20 mm flat-bottom holes, for the raw and TSR pre-processed thermograms. It is possible to notice a significant improvement in contrast for all the simulated flaws. Another interesting finding is that for the pre-processed thermograms, the maximum absolute contrast of all the flaws occurs roughly at the same instant of time, while in the raw thermograms the maximum contrast time varies with depth. Shallower defect present maximum contrast earlier, while deeper defects later, as expected from classical heat transfer theory².

Figure 7
Evolution of the absolute contrast of the 20 mm flat-bottom holes during the transient cooling phase for the raw thermograms and thermograms pre-processed using the TSR algorithm (1st derivative).

Figure 8 shows the result of the TSR first derivate at distinct times during cooling. It is possible to notice that the shallower defects appear with greater contrast in the early thermograms, while deeper defects appear later. However, as depicted in Figure 7, the maximum contrast of all defects occurs at roughly the same instant of time, approximately 7.75 seconds. The images also reveal that the defects with diameter of 2 mm cannot be detected, even at 0.5 mm depth (aspect ratio of 4). Defects at 2 mm depth also present very weak contrast, even for larger diameters. These results may be improved by adjusting the testing parameters (e.g., heating time) and performing inspections at closer distances (higher resolution). Nevertheless, it is not the goal of the present work to determine the limitations of LPT, but rather investigate how post-processing algorithms can improve data analysis.

Figure 8
TSR first derivative images at 3.00, 5.00 and 7.75 seconds after heating.

Considering that the TSR first derivative image at 7.75 seconds presented the best absolute contrast when compared to the raw and PPT thermograms, it was used to evaluate the proposed post-processing framework. First, the thermogram was filtered using a median kernel (cv::medianBlur) with size of 11x11 pixels. The filter reduces high frequency noises, mitigating amplification during histogram equalization. After filtering, the image contrast was enhanced using CLAHE (cv::createCLAHE) with clip limit of 1.3 and grid size of 15x15 pixels. These parameters yielded the best results, with an average increase in absolute contrast of 37.4% (See Figure 9).

Figure 9
Effects of median filter and CLAHE.

Lastly, image segmentation was evaluated using Otsu’s method (cv::threshold) and the watershed approach (cv::watershed). From Figure 10, it is possible to see that the segmented images retained an acceptable number of defects, despite using a global automatic strategy. During this step, the ROI was reduced to the laminate only, in order to avoid identification of other contours, such as the stand support and contour of the laminate itself. The Otsu’s binarization output a threshold value of 167, which resulted in a good segmentation, but induced a considerable level of noise, especially at the laminate edges, which emphasizes nonhomogeneous heating effects. The right edge displayed higher levels of noise when compared to the left edge. This noise can be removed by further processing, such as erosion and dilation, for example. Nevertheless, the result allowed the successful application of the border following algorithm for automatic defect sizing. The watershed method was applied using a global threshold of 185, followed by morphological opening (cv::morphologyEx), for sure background estimation. Sure foreground was computed using the distance transform (cv::distanceTransform) of the sure background image. The watershed markers were then calculated using the subtraction of sure foreground and sure background.

Figure 10
Measurement procedure performed using the ImageJ software on the TSR first derivative image after histogram equalization (CLAHE), in comparison with the Otsu’s method and the watershed segmented images.

Flaw size estimation was performed using a border following algorithm (cv::findContours), which allows the computation of the contour area by estimating the equivalent diameter of the smallest circle that contains the complete contour. Segmentation is essential for this step, since the border following algorithm cannot be applied to non-binary images. Table 2 summarizes the planar size estimation of the flat-bottom holes for the segmented images in comparison to the sizing retrieved manually using ImageJ. Out of the 25 simulated flaws, a total of 12 flaws were automatically identified and quantified using the proposed framework. The flat-bottom holes with diameter of 20, 15 and 10 mm at 2.00 mm depth can be detected with weak contrast, as well as the flaws with 5 mm in diameter at 0.50, 1.0, 1.25 and 1.50 mm in depth. The watershed segmentation presented more consistent sizing results when compared to the Otsu’s binarization, both retrieving the same defects.

Table 2
Planar size estimation retrieved with the border following algorithm.

Depth [mm]	Diameter [mm]	Otsu’s Method [mm]	Watershed [mm]	ImageJ [mm]
0.50	20.00	22.89	21.01	20.81
0.50	15.00	17.35	14.98	16.13
0.50	10.00	11.10	10.75	9.66
1.00	20.00	23.64	20.58	21.49
1.00	15.00	16.56	15.52	15.06
1.00	10.00	11.38	11.05	10.35
1.25	20.00	21.06	20.28	19.15
1.25	15.00	16.51	15.33	15.01
1.25	10.00	11.39	10.49	10.35
1.50	20.00	21.82	20.83	18.98
1.50	15.00	15.60	15.55	13.46
1.50	10.00	11.09	9.74	11.39

5. Conclusions

The present work proposed a simple framework, based on the OpenCV library, for post-processing thermograms acquired during LPT inspections. The framework consists of single image techniques used to filter the thermograms and segment the image into defective and sound areas (background). The thermograms were pre-processed using the TSR and PPT algorithms. Results revealed that the framework can aid defect detection by automatically segmenting the image and can also improve flaw size estimation. The following conclusions may be drawn from this study:

• The TSR algorithm outperformed the PPT algorithm in terms of absolute contrast.

• Post-processing by means of adaptive histogram equalization (CLAHE) can significantly increase defect contrast.

• Image segmentation for automatic defect detection can be successfully achieved using Otsu’s method and the marker-based watershed algorithm.

• The watershed segmentation approach yielded more accurate defect size estimations when compared to Otsu’s method.

The OpenCV library offers a wide variety of powerful, built-in algorithms that can be further used to enhance image processing of thermograms for TNDT. Moreover, the library also offers tools for the development of machine learning algorithms, such as Convolution Neural Networks (CNNs) and Support Vector Machines (SVMs), that can be used to develop automatic defect classification strategies.

6. Acknowledgements

The authors would like to acknowledge the Brazilian funding agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and the São Paulo Research Foundation (FAPESP) for the grants CNPq 306193/2017-5 and 163101/2017-5 and FAPESP 2019/00917-8 and 2018/26122-9, that supported this research.

7. References

• ¹
  Ciampa F, Mahmoodi P, Pinto F, Meo M. Recent advances in active infrared thermography for non-destructive testing of aerospace components. Sensors. 2018;18(2):609.
• ²
  Meola C, Boccardi S, Carlomagno GM. Infrared thermography in the evaluation of aerospace composite materials [Internet]. Duxford: Elsevier; 2017 [cited 2020 Mar 1]. Available from: https://linkinghub.elsevier.com/retrieve/pii/C20130182268
  » https://linkinghub.elsevier.com/retrieve/pii/C20130182268
• ³
  Vavilov VP, Burleigh DD. Review of pulsed thermal NDT: physical principles, theory and data processing. NDT Int. 2015;73:28-52.
• ⁴
  Szeliski R. Computer vision: algorithms and applications. London: Springer; 2010. 979 p.
• ⁵
  Shepard SM. Active thermography for NDT. São Paulo: Conaendi; 2018.
• ⁶
  Shepard SM. Advances in pulsed thermography. In: Rozlosnik AE, Dinwiddie RB, editors. SPIE Proceedings: Thermosense XXIII; 2001; Orlando, FL. Bellingham: SPIE; 2001. p. 511-5.
• ⁷
  Rodríguez F JL. Detection and characterization of subsurface defects by infrared pulsed thermography [thesis]. Florianópolis: Universidade Federal de Santa Catarina; 2013.
• ⁸
  Maldague X, Marinetti S. Pulse phase infrared thermography. J Appl Phys. 1996;79(5):6.
• ⁹
  Ibarra-Castanedo C, Genest M, Piau JM, Guibert S, Bendada A, Maldague X. Active thermography techniques for the nondestructive testing of materials. In: Chen CH, editor. Ultrasonic and advanced methods for nondestructive tesing and material characterization. Singapore: World Scientific Publishing; 2007. Chapter 4; p. xix.
• ¹⁰
  Open CV Team. OpenCV [Internet]. 2020 [cited 2020 Mar 1]. Available from: www.opencv.org
  » www.opencv.org
• ¹¹
  Otsu N. A threshold selection method from gray-level histograms. IEEE Trans Syst Man Cybern. 1979;9(1):62-6.
• ¹²
  Farrahi Moghaddam R, Cheriet M. AdOtsu: an adaptive and parameterless generalization of otsu’s method for document image binarization. Pattern Recognit. 2012;45(6):2419-31.
• ¹³
  Digabel H, Lantuejoul C. Iteractive algorithms. In: Proceedings of the 2nd European Symposium on Quantitative Analysis of Microstructures in Material Science, Biology and Medicine; 1978; Caen, France. Stuttgart: Riederer; 1978. p. 85-99.
• ¹⁴
  Preim B, Botha C. Visual computing for medicine: theory, algorithms, and applications. 2nd ed. Amsterdam: Elsevier; 2014. 812 p.
• ¹⁵
  Suzuki S, Abe K. Topological structural analysis of digitized binary images by border following. Comput Vis Graph Image Process. 1985;30(1):32-46.
• ¹⁶
  Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012;9(7):671-5.
• ¹⁷
  American Society for Testing and Materials – ASTM. ASTM E2582-19: standard practice for infrared flash thermography of composite panels and repair patches used in aerospace applications [Internet]. West Conshochen: ASTM International; 2019 [cited 2020 Mar 1]. Available from: www.astm.org

Publication Dates

• Publication in this collection
  27 Sept 2021
• Date of issue
  2021

History

• Received
  26 July 2021
• Reviewed
  21 Aug 2021
• Accepted
  25 Aug 2021