CIRCULAR-SHAPED OBJECT DETECTION IN LOW RESOLUTION SATELLITE IMAGES

Circle and ellipse detection have been a widely discussed topic among the computer vision community. Applying circle detection methods to satellite images can bring valuable information on urban or industrial areas. A limitation of such methods lies in the resolution of the object to detect. The smaller they are in the image, the harder it becomes to detect them accurately. In this paper, we explore several circle detection methods adapted to low resolution satellite images. An algorithm based on level-lines detection and classification will be presented in this paper along with other well known algorithms. The methods are evaluated in the context of oil tank detection in Sentinel-2 images where those circular objects can have an observed diameter as small as 2 pixels.


INTRODUCTION
Recognizing circular objects in images can bring insightful information on a scene. This is all the more interesting on satellite images with a wide coverage of the Earth's surface. Several sophisticated algorithms proposing to tackle this problem have also been proposed such as the Ellipse and Line Segment Detector presented in (Pȃtrȃucean et al., 2017) and the Hough Circle Transform (Illingworth, Kittler, 1987).
In this paper we propose to compare these methods to another one based on Closed Level-Line Extraction inspired by the Fast Level-Line Transform (Monasse, Guichard, 2000a) on which the isoperimetric ratio is computed to detect circles. Indeed, circles are characterized among all other curves as maximizing the isoperimetric ratio. Extracting the edge map of an image and observing which closed shape has a high enough isoperimetric ratio can therefore be a simple way to detect circles.
The specific context of oil tank detection has been addressed in several works. Detecting oil depots in satellite images allows civilian and military actors to assess a region's economic asset. Several papers have proposed methods to tackle this problem such as the one described in (Han, Xu, 2012, Cai et al., 2014, Ok, Baseski, 2015, Zhang et al., 2015, Soundrapandiyan, 2017 and more recently in (Zalpour et al., 2019) and (Tadros et al., 2020).
These papers use two properties of oil tanks. The first one is that they are cylindrical buildings that appear as circular or ellipsoid objects in satellite images. The second one is that they usually are grouped in clusters. State-of-the-art methods take advantage of these geometrical characteristics. Han et al. (Han, Xu, 2012) use geometrical properties of circles to detect them from a saliency map. Then, they filter their detections and extract oil depots by using a graph-based clustering method. Cai et al. (Cai et al., 2014) enhance the saliency map method from (Han, Xu, 2012) and use the Hough Transform (HT) for circle detection and a Support Vector Machine Classifier (SVM) to keep only oil tanks. * Corresponding author -atadros@ens-paris-saclay.fr Several machine-learning based oil tank detectors have also been developed in recent years. The authors of (Zhang et al., 2015), (Soundrapandiyan, 2017) and (Zalpour et al., 2019) use feature extractors (SURF, HOG or pretrained neural networks) on top of which a classifier is trained on an annotated dataset. On the other hand, (Tadros et al., 2020) proposes a calibrationfree clustering method controlling the number of false detections.
In the work presented here we will focus on the use of circle detectors to recognize oil tanks in optical satellite images.

OBSERVATIONS
Tanks can appear in various ways in the Sentinel-2 images. They are commonly white and therefore appear as lighter disks. However, there is a high variability in the diameters of those tanks, ranging from 1 to 10 pixels. See figures 1a-1d. It is particularly difficult to single out tanks that have a diameter less than 2 pixels. When the diameter gets larger than 4 pixels, it is possible to perceive textures ( Figure 1a) or curvatures (Figure 1b) on the roof.
Those tanks also cast shadows whose length varies between 1 and 5 pixels depending on their height and on the season. This has the advantage of further highlighting white tanks and could be used as a complementary information for detecting tanks.
Shadows can also appear inside the so-called floating roof tanks ( Figure 1e). The roof of floating roof tanks are not fixed but floats on the surface of the stored liquid. The external walls are therefore higher than the roof and cast a shadow on it. See Figure 2. Tanks are, however, not always light and can be at the same grey level or even darker than the ground. See Figure 1f. These tanks are noticeably more difficult to detect: inside and outside shadows allow here to differentiate these types of tanks from other structures.
Furthermore, differentiating tanks from other structures can sometimes prove difficult as shown in Figure 3. Small square structures -such as houses -are for instance similar to small tanks. Roundabouts are very clear examples of round structures that are not tanks. Vegetation can also sometimes create small and darker circular shapes that could be mistaken for darker tanks.

METHODS
Our circle detector applied to oil tank detections is a two-step process. Firstly, we detect salient objects in the image, one notable technique being edge detection. Secondly, we check the circularity of each of those salient objects and filter out those that are not circular. These two steps are solved independently: we can use various combinations of methods of the first and second step.

Salient Object Extraction
3.1.1 Minimum of all directional top-hat transforms A fast local extrema extractor can be obtained by top-hat transforms. As most tanks are represented by lighter or darker small circles compared to the ground, one can extract them using the top-hat and bottom-hat transforms (Serra, 1983). The top-hat of a grey level image I extracts thin lighter elements in the image I and is obtained as follows: where S is a structuring element and opening(I, S) is the morphological grey opening transform of the image I with the structuring element S.
The bottom-hat of an image I extracts thin dark elements in the image I and is obtained as follows: where S is again a structuring element and closing(I, S) is the morphological grey closing transform of the image I with the structuring element S.
The top-hat or bottom-hat transform with a round structuring element of diameter d will highlight respectively lighter or darker elements that are thinner than d in at least one direction. It will therefore highlight small elements but also lines and rectangles thinner than d.
As we only want to extract small round elements and not lines or rectangles, we compute the minimum of directional top-hatŝ Tw or bottom-hatsT b : where S d is a line structuring element of size d. This ensures that we highlight elements that are thinner than d in all directions (contrary to any direction for round structuring elements).
In multi-channel images, the top-hats and bottom-hats transforms are applied independently on each channel. Figure 4 shows a Sentinel-2 image on which we have applied the top-hat transform with a round structuring element and the minimum of directional top-hats. The minimum of directional bottom-hat can detect darker tanks, but in most cases detect shadows as shown in Figure 5. These modified top-hat and bottom-hat transforms could therefore be used in complementary processes to enhance detection.
These transforms can either be used as a preprocessing method before another method, or directly used to detect tanks by applying a simple threshold on them. See Figure 6.

Canny-Devernay's Sub-Pixellic edge detector
Edge detection provides a useful image representation. Most of the current edge detection algorithms are inspired by Canny's method presented in (Canny, 1986). Let I be an image. To detect the edges of I its gradient is first computed, where h is a Gaussian kernel of standard-deviation σ. An edge point strength is defined as the norm of the gradient at this point ||g(·)||. Canny's method keeps points whose gradient magnitudes are local maxima along the direction of the gradient as proposed by Haralick (Haralick, 1982) .
After this non-maxima suppression step, Canny proposed to clean the edge map from spurious low contrasted structure by applying a high threshold H and a low threshold L on the edge points strength. If the strength of the gradient at a pixel is higher than H, the edge point is validated; if it is between H and L the edge is validated only if it is connected to another validated edge point.
We are primarily interested in detecting oil tanks larger than 20 meters. For medium resolution images acquired with satellites such as Sentinel-2, this means that the smallest tanks can have a diameter as small 2 pixels. Therefore, we need to detect edges at a sub-pixellic precision to be able to detect circles in the next step. To achieve this, we use the method presented by Devernay (Devernay, 1995) which adapts Canny's method by using a quadratic interpolation of the gradient norm to approximate the sub-pixel positions of edges.
More details on the algorithm can be found in (Devernay, 1995). Some improvements were proposed in (Grompone von Gioi, Randall, 2017). All in all, we end up with three parameters to set for the edge extraction step: σ, H, L.
Automatic parameters tuning In (Fang et al., 2009), the authors proposed to use Otsu's method as a way to automatically decide the value of H and L.
Otsu's method segments an image I by finding a threshold on the pixel's histogram that minimises the intra-class variance.
Let τO be the optimal Otsu's thresholding, (Fang et al., 2009) defines Canny's thresholds as For the experiments we will set L and H using a grid search approach and the Otsu-based method, that we will refer to as the "Otsu-Devenay" method in the rest of the paper.

Level-Line Extraction
The set of all image levellines constitutes a complete image representation, also known as topographic map. Each level line is actually a local holistic geometric feature that often corresponds to the boundary of a shape present in the image. The Fast Level Line Transform (FLLT) (Monasse, Guichard, 2000a) is a fast algorithm organizing all level lines in an inclusion tree and encoding each level line as a piecewise affine polygon. The Image Curvature Microscope (Ciomaga et al., 2017) (ICM) filters by affine curvature motion (Sapiro, Tannenbaum, 1993, Alvarez et al., 1993) the level lines obtained by this subpixel extraction. This procedure is more accurate than standard methods using a finite difference method. To extract the initial subpixel level-lines, a bilinear interpolation of the image is used (Caselles, Monasse, 2010). Then, an independent smoothing is performed on the levellines by applying the affine scale space (Ciomaga et al., 2010), through a geometric scheme (Moisan, 1998) that smooths the level lines and in particular removes pixelization artifacts. At the end, the Image Curvature Microscope delivers level lines endowed with an accurate curvature, which are much more informative on the topology of the objects in the image than simple edges. For the circle detection, we only extract level-lines with the highest local contrast in the tree of level lines. In other terms, the average image gradient across these level lines is higher than on the next and former level lines in the inclusion tree (Monasse, Guichard, 2000b). This method can replace the edge map to give a good representation of an image content.

Shape Classification
The salient object extraction methods presented above, are efficient ways to simplify an image representation for our circle detection task. We can then apply methods adapted to this new representation to answer the problem.
3.2.1 Isoperimetric Ratio Thresholding Once the edges or the level-lines are extracted from the image, a geometrical property of circles can be used in order to detect them. First, only closed curve are kept. If a chain is closed, it defines a polygon. Let S and P be respectively the surface and the perimeter of the geometrical shape defined by this closed edge chain. To classify a closed edge chain as a circle, a threshold is set on the isoperimetric ratio ρ of the object where ρ = 4πS P 2 . The closer an object's shape to a circle, the closer the isoperimetric ratio ρ is to 1, following the isoperimetric inequality: The shape of oil tanks is never perfectly circular in satellite images as shown in section 2, especially when the resolution is relatively low. To remove wrong tank candidates, we reject objects with isoperimetric ratio lower than ρ th . Most buildings and objects appear in satellite images as either circles or polygons, and in particular rectangles. Among the latter, squares have the largest ρ value, equal to π 4 ≈ 0.79. To be sure to reject rectangular shapes, we set ρ th in a range between 0.8 and 1.

Hough Circle Transform
The Hough Circle Transform is a well known method for extracting features to detect circles in images.
In a 2D space, a circle is described by x0 and y0, the coordinates of its center, and its radius r; the points (x, y) of the circle satisfy the following equation: The Hough Space in which we will look for the best parameters is then a 3D space. We first extract the edges from the image and we discretize the 3D Hough Space, which can then be represented by a 3D accumulator matrix, to test various triplets of parameters. For each triplet, we count the number of edge points that fall on the circle that it describes. Circles are then selected by finding the local maxima of this accumulator matrix.

Ellipse and Line Segment Detector -ELSDc
The last detector we shall evaluate is the ELSDc algorithm (Pȃtrȃucean et al., 2017), which jointly detects line segments, elliptic arcs and circular arcs while controlling the number of false detections using the a contrario detection theory (Desolneux et al., 2008). This method works in three consecutive steps: 1) candidate generation; 2) candidate validation; and 3) model selection. What follows is a brief summary of the method; we refer to the original publication for more details.
Candidate generation To avoid testing exhaustively all the possible geometrical features in the image, a heuristic step is used to propose a reduced set of candidate. Starting from a seed pixel, a region growing procedure recursively merges pixels with similar gradient orientation and a rectangle is computed that encloses the obtained set of pixels. Then, new seed pixels are selected at the rectangle end-points and the procedure is repeated, resulting in a chain of rectangles. Finally, these chains are cut into parts according to curvature and three candidate geometrical features are fit to each part: a polygonal line p d , a circular arc c, and an elliptic arc t, see Figure 8.
Candidate validation A validation step is applied to candidates, based in the a contrario detection theory (Desolneux et al., 2008) which implies keeping only the structures which are too regular to be the result of an accidental arrangements of independent parts. This theory requires a background model H0 where NT is the number of tests performed in principle and the right-hand term is the probability of obtaining in the background model H0 an error KH 0 (e) smaller or equal to the observed one k(e). The smaller the NFA, the more unlikely the event e is to be observed by chance in the background model H0; thus, the more meaningful. Candidates with NFA ≤ ε, for a predefined ε value are kept as valid detections. It can be shown (Desolneux et al., 2008) that under H0, the expected number of tests with NFA ≤ ε is bounded by ε. As a result, ε corresponds to the mean number of false detections under H0.
In ELSDc, the error of a geometrical event e is measured by the image gradient orientation, which should be orthogonal to the geometrical feature. More precisely, for an event e supported by the set of pixels S(e), where, ∇I(i) is the gradient of the image at pixel i, and dir ⊥e (i) corresponds to the normal direction to the geometrical event e at position i. In this case, H0 corresponds to pixels in which the gradient orientations are independent random variables following a uniform distribution on [−π, π). It can be shown that in these conditions, KH 0 (e) follows the Irwin-Hall distribution and the probability term can be upper-bounded by where l(e) is the number of pixels supporting the event e.
The formulation is completed by the Number of Test, which is an estimation of the theoretical number of features of each kind. In ELSDc and for an n × m image, the number of tests is approximated by 2 d+1 (mn) (d+1)+ 1 2 for a polygonal line of d sides, and 18(mn) 4 for an elliptical event. The final NFAs are for a polygon p d of d sides and for an elliptical arc t. (Circular arcs share the same NFA as elliptical arcs.) Events with NFA ≤ ε are considered as detection. For this problem, ε = 1 (Pȃtrȃucean et al., 2017).
Model selection For some chains of rectangles, it could happen that more than one candidate (p d , c, or t) obtains NFA ≤ ε and therefore are meaningful. ELSDc simply keeps, among those three, the interpretation with smallest NFA.

EXPERIMENTAL VALIDATION
In this section, we will first quantitatively evaluate the methods we previously presented on a validation site where the ground truth is known. We will in a second part discuss the influence of the preprocessing. Another simple preprocessing step has been considered, consisting of 2× zoom by zero-padding in the Fourier domain.

Quantitative results
We have selected a validation site of a refinery containing all interesting cases: large and small tanks, dark tanks and floating roof tanks. See Figure 9a. We have annotated the center of all 231 visible tanks. See Figure 9b.
Given the detection results of an algorithm on this image, we evaluate its performance as follows. First, we take the center of each detection. We have then two lists: P det , the list of detected points, and Pgt, the list of points in the ground truth. We then compute the recall, precision and F1 precision metrics as described in Algorithm 1.
Algorithm 1: Recall, precision and F1 measurements Input: P det , Pgt, tol = 3 S ← all pairs of points Pi ∈ P det and Pj ∈ Pgt D ← distance of each pair in P I ← Argsort D from lowest to highest Remove all elements in I where if Pj has already been assigned then continue; As our methods are parametric, we evaluate their performance on a set of plausible parameters as shown in Table 1. In section 3 we stated that the isoperimetric threshold should be above 0.8 to avoid detecting squares, but it appears that most of the shapes detected by the level-lines extractor have isoperimetric ratios under this value as we can see in Figure 10. We extend then the range of tested values for ρ th to 0.2, see Table 1.
We observe that the level-line-based method has the best F1 score and it achieves the highest recall when it is preceded by the top-hat transform. The ELSDc method seems to attain a really high precision, but on the other hand it has the lowest recall, meaning that it has only a few accurate detections. We can also note that the simplest method based on Devernay's edge extraction and isoperimetric thresholding is the second best method despite its simplicity. The Otsu-Devernay method seems to present results close to the ones that we can observe with the Devernay method, while preventing us from hand-tuning the parameters H and L for the edge extraction. Detection result with the best configuration for each method are shown in Figure 11.
Overall, the dark tanks tend to be much more often missed than the lighter ones. It is not very surprising as the edges are much less visible than on lighter tanks. Floating roof tanks are often detected despite the fact that their roof doesn't appear as perfect disks but as ellipses. This is not a very important issue as ellipses tend to have an isoperimetric ratio closer to circles than rectangular shapes.
In Figure 12, we show several qualitative examples of tank detection on other sites using the level-line-based method that achieved the best f1-score. Overall, the observed result in the validation site seems to generalize well to other areas. Two preprocessing on the images have been used on before applying the methods presented in Section 3. The top-hat, which help enhancing contrasted areas, and the zero-padding zoom.

Preprocessing influence.
Top-hat Transform: Based on Table 2, we note that adding the top-hat processing always helps to improve the recall but can lead to a precision loss. For all the methods this can be interpreted quite easily. This transformation increases the gradient magnitude at the edges by increasing the contrast. This helps Canny-Devernay's extraction method which is used for both the isoperimetric threshold and the HCT. The ELSDc also uses the gradient information to build its region candidates. The level-line-based method can also benefit from this processing. However, the top-hat also bring to light small structures that would otherwise have been much less visible, creating false positive and therefore reducing the overall precision.

Generalisation to other images.
Those models have been calibrated for subpixel circle detection. We have also tested on PlanetScope images, which have a resolution of 3m per pixel, to see if the methods proposed are also suited for more informative images. As seen in Table 2, only the Dervernay and leve-line based methods seems promising. Nevertheless, the later method that seems to perform well most of the time, struggles when it comes to more resoluted images, while the Devernay method continue to give decent results. As shown in Figure 13, the level-line method detects a lot of false alarm, while Canny-Devernay's still seems accurate.

CONCLUSION
We have presented several circle detection methods adapted for tank detection on Sentinel-2 images. We have quantitatively compared them and we have shown how well they generalize to other sites.
Results are promising despite the low resolution of Sentinel-2 images. There are some areas of improvements, such as the detection of darker tanks. Applying these detectors on image series instead of a single image and using additional inputs such as Sentinel-1 images might further improve their performance.
Moreover, using time series might reduce the illumination problems as a missed detection at one date might be detected on another date when the illumination conditions are better.
Code 1 and dataset are available online.