Detection and Delineation of Sorted Stone Circles in Antarctica

: Sorted stone circles are natural surface patterns formed in periglacial environments. Their relation to permafrost conditions make them very helpful for better understanding the past climates where they were formed and have evolved and also for monitoring current underlying processes in case circles are active. These metric scale patterns that occur in clusters of tens to thousands of circular elements, can be more comprehensively characterized if automated methods are used. This paper addresses their identiﬁcation and delineation through the development and testing of a set of automated approaches, namely, template matching, sliding band ﬁlter, and dynamic programming. All of these methods take advantage of the 3D shape of the structures conveyed by digital elevation models (DEM), built from ultra-high resolution imagery captured by unmanned aerial vehicles (UAV) surveys developed in Barton Peninsula, King George Island, Antarctica (62 o S). The best detection results achieve scores above 85%, while the delineations are performed with errors as low as 7%.


Introduction
Sorted stone circles are a type of patterned ground formed in periglacial regions of the Earth [1,2], and also on Mars [3,4], easily attracting the attention due to its remarkable circular geometry ( Figure 1). These metric circular shapes, occurring in clusters of up to thousands of individuals, are constituted by an inner fine domain of soil surrounded by stone rings ascending few centimetres above its interior. This radial segregation of the geological material according to size (named sorting) is performed at the surface after the upward swelling of the soil in the active layer of permafrost in seasonal freezing and thawing cycles [1,5,6].
Their characterization and monitoring can provide useful information about the underlying processes if circles are currently active [7,8] and also about the past climatic conditions that allowed their formation/evolution [9][10][11][12][13]. In addition, besides their evident geomorphic importance and their use as paleoclimatic indicators, the periodic burial and exhumation of geological materials may also play an important role in the soil carbon cycle in connection to permafrost conditions [6]. Although the processes underlying their formation and evolution are relatively well understood [14][15][16][17][18], their geometric characteristics, dating, and seasonal dynamics are mainly supported by sampling and monitoring few circles in the field [6,11], preventing obtaining robust statistics and an extended spatial analysis of large fields. Ultra high-resolution imagery (mm to cm) captured by unmanned aerial vehicles (UAV) can greatly contribute to a much more complete geomorphic description of this type of patterned ground in large areas, which cannot be perceived in current sub-metric satellite imagery, giving also context to data built at ground-level [11,19]. The use of UAV in polar regions, and particularly in Antarctica, is becoming a very important tool for many studies due to the possibility of gathering data on these less accessible regions with unprecedented detail. Although field operations can be easily constrained by weather conditions, the diversity of applications is widening fast, encompassing for instance, glacial retreat quantification [20], landforms mapping [21], vegetation mapping [22] and health state assessment [23], or wildlife inventorying and monitoring [24], just to name some of the most recent ones. In what concerns their use for studying stone circles, it is a topic we have started addressing recently [25] and which, as far as we know, has not been studied before using remote sensing.
In addition, the large amount of circles to detect requires the use of automated methods to make their segmentation, characterization and also the multitemporal monitoring to account for the dynamics of the processes involved. This paper deals now with the initial phase, that is, with the identification and delineation of the stone circles. However, the detection of circular patterns is a classic pattern recognition problem using, for instance, popular methods such as template matching and the Hough transform [26,27]. These approaches, requiring the pre-definition of templates to compute the associated voting scheme, are not the most adequate to deal with the stone circles, as they are not able to take into account the natural deformation these circular patterns can exhibit. The use of other methods that can deal with more than slight deformations is therefore the right choice. Thus, we use a pair of such methods, based on the sliding band filter [28] and on dynamic programming [29][30][31], which we compare to a traditional template matching method.
Thus, this paper aims at detecting sorted stone circles and delineating its boundaries in elevation data derived from imagery acquired by a UAV in Antarctica. It extends the preliminary work presented in [25] by considering: • the delineation of the stone circle contour; • an additional method based on dynamic programming; • the evaluation of detection and delineation problems in an extended dataset.

Study Area
The sites selected for developing the surveys are located in the ice-free areas of Barton Peninsula (62 • 14'S; 58 • 46'W) in King George Island, the largest one of the South Shetlands in Antarctica (Figure 2 left). This peninsula has an area of about 10 km 2 that started to be exposed about 8000 years ago after glacial retreat [32]. It shows a rugged topography, with a wide and gentle slope in the central area, ranging from about 90 to 190 m a.s.l., where the majority of stone circles fields are found [33]. It is also one of the regions in Antarctica where sorted stone circles are more ubiquitous and therefore an adequate location to study such type of patterned ground.

UAV Surveys
The surveys were developed during February 2018 in the frame of a field campaign with logistics support from the Portuguese Polar Program (PROPOLAR) and the Korean Polar Research Institute (KOPRI). The aerial flights were conducted with a quadcopter DJI Phantom 3 (SZ DJI Technology Co., Ltd., Shenzhen, Guangdong, China) equipped with an RGB camera (12 Mpixel, lens 35 mm equivalent, f/2.8) in the locations of Barton Peninsula where sorted stone circles have the most extended occurrence (Figure 2, right). They were developed at a low altitude (10 m above ground) to allow perceiving stones down to 5 mm in size. Each flight was planned in a double-grid mode with high lateral and frontal overlaps between adjacent images (80%) to provide an increased number of views for the same point. The area surveyed in each individual flight was conditioned by battery duration reaching up to 80 m × 80 m and acquiring up to 500 images.

Mosaic and DEM Construction
The set of images captured in each flight were processed with Agisoft Photoscan software (version 1.4, Agisoft LLC, St. Petersburg, Russia) to build an orthorrectified image mosaic and a digital elevation model (DEM), using ground control points measured with a D-GPS to provide accurate georreferencing. The main procedures underlying the construction of the two products are based on SIFT (Scale-Invariant Feature Transform) to detect the same feature in adjacent images, RANSAC (Random Sample Consensus) method to filter erroneous matches and SfM (Structure from Motion) technique to estimate the 3D geometry of the scene [34]. An example of such product is presented in Figure 3.
There are small variations in the spatial resolutions of the final mosaic and DEM from site to site due to the natural fluctuations in position of the UAV along the development of the aerial surveys. The spatial resolution for all mosaics is in the range of 2.6-3.2 mm/pixels, while, for all DEMs, is in 2.1-2.4 cm/pixel. These variations in resolution are small and thus can be considered negligible for the generic application of detection methods. However, only the DEM were used in this work, as the geometric information provided is the most relevant to detect and delineate the circles. The use of RGB color information can be more difficult to interpret since color is not directly related to geometry. However, RGB information may be helpful to solve ambiguous situations, when the transitions between the geological materials of the circles and background have a poor expression. On the contrary, if vegetation also occurs with the patterned ground, namely mosses surrounding the circles or lichens on top of the their stones (which does not happen in all field sites), some additional information (edges) can be used and therefore improve the performances. This approach deserves being exploited in the future.
The image mosaics will be also used to analyze other features like, for instance, the size of the stones or the type of associated vegetation [35].

Methods
This paper addresses two main problems: the detection of stone circles in elevation data and the delineation of the contour of each circle. The first problem will be tackled by three methods: template matching (TM), sliding band filter (SBF), and dynamic programming (DP). The second problem (contour delineation) is discussed later (Section 3.4), but, as we will see in the following, the last two methods (SBF and DP) implicitly solve the delineation problem as an intermediate step.

Template Matching (TM)
The patterns to be detected mostly have a circular shape with a size approximately constant within each site (Figure 4, left). Their dimension, which may change a bit in the field from site to site, depends on some local features, namely the depth of the active layer of permafrost and the characteristics of the soils.
We assume that the elevation image contain 3D patterns of known shape T(x, y), which is denoted as template. The fit between the elevation image I(x, y) and the template T(x, y) that slides along the whole image can be computed through the correlation which reaches a local maximum when the template and the circle are aligned with each other. The template that better describes the stone circles is a half-torus (Figure 4, right) described by where R is the radius, R is the torus width ( = 0.2), and x + = max(0, x) is the rectifier function. The detection of the stone circles is based on the analysis of the correlation matrix R(x, y) by detecting all the local maxima (non-maximum suppression), followed by a comparison with a threshold Th ( Figure 5).
Since the elevation image normally has some natural slope, the image I is first pre-processed by a high-pass filter to level it.

Sliding Band Filter (SBF)
The gradient field of natural circles exhibits radial vector patterns ( Figure 6). To detect such patterns, we adopt a modified version of the sliding band filter (SBF) proposed in [28] for the detection of cells in microscopy images. The problem we are addressing now is slightly different. The gradient of cell images points inwards while the gradient of stone circles points inwards and outwards, depending on the distance to the centre.
Let us assume that we know the centre of the circle c = [x, y] T and let us define N radial directions starting from c, with angles θ i = 2πi/N, i = 0, . . . , N − 1. Along each direction, we will measure the alignment between the gradient direction and the radial orientation θ i , given by where φ i (r) denotes the direction of the gradient field at a point c + r[cos(θ i ), sin(θ i )] T . The gradient orientation changes when the point is located at the circle boundary. To detect the boundary, we will perform template matching along each direction and choose the radius r 0 that corresponds to the highest peak. Therefore, the best alignment in direction θ i is where is the template (u(r) that denotes the unit step signal). Since there are N radial directions, the global alignment is the sum of all the directional alignments These equations provide the alignment score associated with the centre hypothesis c = [x, y] T . However, since, in practice, we do not know c, we compute the alignment score for each image point c and detect the circle centers through non-maximum suppression and thresholding.
The gradient vector is computed with Sobel masks and its magnitude is compared to a threshold: gradients with small magnitude are discarded.

Dynamic Programming (DP)
The sliding band filter finds the best alignment in each direction θ i , according to Equation (4). This is a sub-optimal approach and often leads to abrupt transitions of the radius at consecutive directions. To overcome this difficulty, we will perform a joint estimation of the contour points at all the directions, using dynamic programming (DP).
A contour is represented as a sequence of radius variables r 0 , . . . , r N−1 . Since the contour is closed, we will extend the sequence with an additional variable r = (r 0 , . . . , r N ) and assume that r 0 = r N .
In addition, two types of penalties will be considered: the alignment cost at direction i and radius r, inspired in the SBF alignment, and the transition cost which becomes infinite if the radius difference exceeds a threshold: |r − r | > m. Therefore, the total cost E associated with a contour starting at r 0 = k and ending at r N = k is given by the quantity The minimization of E(r), under the constraint r 0 = r N = k can be efficiently solved by dynamic programming [36,37]. This procedure involves two steps which are known as forward recursion and backward recursion.
The forward recursion computes the optimal cost to go from the direction 0 and radius k to the direction t and radius j E t (j) = min The optimal costs can be recursively computed in an efficient way through Since we also want to retrieve the optimal path, we need to store the value of r that minimizes [E t−1 (r) + c(r, j)]. This is done by defining a set of pointers φ t (j) = arg min r [E t−1 (r) + c(r, j)] , j = 0, . . . , M − 1.
After computing E t (j), φ t (j), t = 1, . . . , N, the optimal contour can be recovered by backtracking. The backward recursion is initialized as r * N = k. Then, we recursively compute The algorithm with the constraints r 0 = r N = k is summarized in Table 1. It finds the optimal path under the constraint r 0 = k. Since we do not know the best initial radius k, all the k values should be tested and the best (the one that minimizes E) should be chosen. However, this approach is computationally heavy, therefore we have adopted a faster algorithm based on two repetitions only, as proposed in [29]. Table 1. Dynamic Programming algorithm with constraints r 0 = r N = k.

Backward Recursion
The minimum value of E is a score that measures how well the image patch considered can be fitted by a smooth contour. It is, therefore, a score that can be used for the detection of the stone circles. In a typical application, the patch slides along the image. Every image pixel is considered to be the centre of a stone circle being the detection obtained by non-minimum suppression followed by thresholding.

Contour Delineation
The SBF and DP methods provide an estimate of the stone circle contour in polar coordinates. Therefore, nothing else has to be done to delineate the object contour.
In SBF, each contour point is independently obtained by the maximization of the filtered alignment along each radial direction Equation (4). In the DP method, the contour is the output of a joint minimization of all the contour points performed by dynamic programming. The extraction of the contour corresponds to the backtracking step of the DP algorithm.

Datasets
The datasets used for detecting and delineating the stone circles were obtained by the authors in Antarctica (Section 2.2). These images contain the natural variability found in the field related to several circle features, namely, dimension, circularity, ring thickness, and height and degree of completeness.
The detection dataset is constituted by 20 elevation images of 1000 × 1000 pixels extracted from the larger DEMs built in each site, whereas the delineation dataset is constituted by 24 individual circles extracted from the previous dataset.

Performance Evaluation
The evaluation of the performances of detection and delineation algorithms is performed separately with each of their own metrics.
For the detection algorithms, ground truth (GT) information was built by detecting visually all the circles in each DEM image and annotating manually their centers using an image analysis software (MATLAB, version 9.6 (2019a), MathWorks, Natick, Massachusetts, USA). The algorithms also identify a set of centers, so we need to establish a correspondence between the detected and ground truth centers. This was achieved through a simple matching algorithm defined as follows. First, we build a matrix of distances from all the GT centers to all the detected centers. The minimum element of the matrix is found (excluding the diagonal elements) and the corresponding points are matched, provided that their minimum distance is smaller than a given threshold ∆. If two points are matched, they are considered a true positive (TP) detection and the corresponding line and column are removed from the distance matrix. The process is repeated until there are no more candidates for matching. The ground truth points not matched are considered false negatives (FN), and the detected points that were not matched are considered false positives (FP).
The performance of the detection algorithm is obtained through the well known measures The Precision and Recall combined provide the single measure F-score For the delineation algorithms, the evaluation is performed through the comparison of the delineated contour with a ground-truth contour that is, a manually drawn contour. For computing the distortion between each pair of contours, we adopt the approach used in a similar delineation problem [38], assuming that the detected contour should be within a band of a given width ∆ around the ground-truth contour. This way, points located outside this band are considered errors, with their percentage used as a measure of distortion between both contours.
Let A denote the set of boundary points detected by the delineation algorithm, and let B denote the ground-truth boundary defined by an expert. We define the distance from a point p ∈ A to the ground-truth contour B as the distance from p to the nearest point q in B The amount of errors is given by the following equation: where #A is the number of points of set A; δ(x) = 1 if x = true, and δ(x) = 0 if otherwise. We have chosen for this study ∆ = 0.15R, where R is the stone circle radius. If d(p, B) = 0, the point p is assumed to be correct, and there is no error (cp). If 0 ≤ d(p, B) ≤ ∆, the error in p is considered a small error (se), and if d(p, B) ≥ ∆, p is considered a gross error (ge).

Detection
The detection algorithms were tested, assuming different configurations for their parameters. To be specific, the three parameters for TM method are R (radius of template), (template width) and Th (threshold), which were tested with the configurations: R ∈ {25, 30 The outputs of each method with the optimal parameters are illustrated with three different images in Figure 7.
These results show the ability of automatic methods to detect stone circles. The TM method performs worse than SBF and DP, since it has more false positive detections. The best method is DP, but the difference with respect to SBF is small. It can also be observed that the performance of the automatic detectors is better in the first and third examples (left and right columns) than in the second (middle column) where the circles have stronger deformations and exhibit a less uniform spatial distribution pattern. Table 2 shows the detection scores (precision, recall and F-score). Not unexpectedly, DP and SBF methods achieve the highest average F-scores with 85.2% and 82.0%, respectively, clearly outperforming TM method with 66.7%. This difference derives from TM making a strong rigidity assumption about the shape of the circles while SBF and DP are flexible and make appropriate trade-offs between shape and deformation.  Nevertheless, the much larger average computation time per image of DP and SBF methods (23.9 and 14.9 min, respectively) when compared to TM (0.1 min) requires their evaluation under less restrictive conditions. Just decreasing the number of directions of analysis in each method from 128 to 16 turns the computation time to be only 15% and 12% of the initial duration for DP and SB, respectively, and still keeping the F-score performances at a high level: 84.3% for DP and 81.3% for SBF.

Delineation
The delineations performed by each method are illustrated in Figure 8. Both methods find the contour of the stone circles, but DP leads to smoother contours since it performs a global estimation of the contour and has less gross errors.
The delineation does not always provide an accurate contour. Two examples are shown in Figure 9 in which SBF and DP fail. In the first case (SBF), there are strong elevation artefacts that produce spurious detections in some directions. In the second case (DP), the shape of the contour can no longer be approximated by a circle and the shape model (circle) does not fit the image. Table 3 shows the performances achieved by each method in the whole dataset of circles. The evaluation is based on gross errors only, since we consider small errors to be acceptable. Although both methods perform well, DP clearly outperforms SBF, with gross errors of 7.1% and 18.3%, respectively. On the contrary, the average computation time for each individual circle with DP (50.2 s) is almost ten times slower than SBF (5.5 s).
This difficulty can be overcome by substantially decreasing the number of directions (from 360 to 64). This change has, similarly to the detection problem, a strong effect on the reduction of the average computation time of DP method with almost no effect on the quality of its delineations: it is more than 15 times faster (from 50.2 to 3.2 s), keeping gross errors practically unchanged (from 7.1% to 7.2%). Although less significant, there is also an average computation time reduction for SBF (from 5.5 to 4.4 s) but now with a slight decrease of gross errors (from 18.3% to 18.0%).

Conclusions
This paper provides new advances for the study of sorted stone circles through automated methods applied to elevation maps built from UAV surveys. Three methods were proposed for circle detection (template matching, sliding band filter, and dynamic programming) with the last two being also used for contour delineation.
Good results were obtained by the sliding band filter and dynamic programming methods, both in the detection and delineation problems. However, the dynamic programming method performs better in both problems. An F-Score of 85.2% was achieved in a database of 20 elevation maps. In the delineation problem, a percentage of 7.1% of gross errors was obtained, being all the other points correctly estimated.
The flexibility inherent to sliding band filter and dynamic programming detection methods makes them a good option to deal with the natural variability associated with the stone circles, since they can encompass variations in their circularity, thickness, and height of the stone ring and also in their percent completeness.
For future work, we intend to exploit the use of the RGB images to improve the performances and also to apply other methods (e.g., deep leaning) to the detection and delineation of the stone circles. These methods require the enlargement of the datasets, which will be achieved after the development of more UAV surveys in other ice-free regions of Antarctica.
Author Contributions: All authors contributed to the development and testing of the methodologies and to writing the manuscript. P.P. and S.H. developed the field surveys in Antarctica. All authors have read and agreed to the published version of the manuscript.