AUTOMATED DETECTION OF OIL DEPOTS FROM HIGH RESOLUTION IMAGES: A NEW PERSPECTIVE

This paper presents an original approach to identify oil depots from single high resolution aerial/satellite images in an automated manner. The new approach considers the symmetric nature of circular oil depots, and it computes the radial symmetry in a unique way. An automated thresholding method to focus on circular regions and a new measure to verify circles are proposed. Experiments are performed on six GeoEye-1 test images. Besides, we perform tests on 16 Google Earth images of an industrial test site acquired in a time series manner (between the years 1995 and 2012). The results reveal that our approach is capable of detecting circle objects in very different/difficult images. We computed an overall performance of 95.8% for the GeoEye-1 dataset. The time series investigation reveals that our approach is robust enough to locate oil depots in industrial environments under varying illumination and environmental conditions. The overall performance is computed as 89.4% for the Google Earth dataset, and this result secures the success of our approach compared to a state-of-the-art approach.


INTRODUCTION
Industrial facilities to store oil and/or petrochemical products are described as oil depots.The primary structure of an oil depot is the tankage, either above ground or underground, wherein such valuable products are stored.One of the key imperatives is the safety of industrial facilities.The recent experiences expose that these regions are highly vulnerable, especially to natural disasters.The latest devastating instance occurred in Sendai, Japan caused massive damages to the region once the oil depots of the largest refinery in Japan were set ablaze by the earthquake.Therefore, risk evaluation of such regions prior to natural disasters is crucial, and could be performed with the help of remotely sensed images.Aerial/satellite images could be useful for locating individual oil tanks, providing information about their content (to some extent), virtual modelling the environment nearby, etc.This may eventually help services responsible for the emergency planning, rescuing operations, and protecting individuals in the nearby residential areas.
This paper is devoted to the automated detection of individual oil depots visible in a single near-nadir high resolution aerial/satellite image.The motivation to use near-nadir image is that the projection of the boundary of the top of a vertical cylindrical shape formed by an oil depot is close to a circle in the image.Therefore, we focus on detecting circles in images acquired from industrial facilities containing oil depots.
Many previous studies in this context assume that oil depots are bright features, whose foreground is clearly separable from the neighbouring background, which eventually simplifies the problem.However, these studies do not evaluate their approaches for difficult conditions (e.g.complex foreground and/or background, occlusion, varying seasonal effects such as illumination, shadow, smoke, and snow cover).However, such cases are quite common for aerial/satellite images; thus, reduces the applicability and generalization of the previously developed approaches.
In this paper, we propose a new approach to automatically detect circular objects from high resolution satellite images.The approach considers the symmetric nature of the circular oil depots and allows us to detect them even in difficult conditions.The proposed approach alleviates the weaknesses of the fast radial symmetry transform (Loy and Zelinski, 2003) with a new framework, and involves an automated thresholding method to focus on circular regions as well as a new measure to validate circles.To show the feasibility of our approach, we use six GeoEye-1 (0.5 m GSD) test images.Besides, our dataset consist of 16 Google Earth images (0.5 m GSD) of an industrial test site acquired in a time series manner (between the years 1995 and 2012).These Google Earth images are taken from both aerial and spaceborne platforms and provide an excellent test domain for circle detection in different/difficult environmental and illumination settings (Fig. 1).
The remainder of this paper is organized as follows.The previous studies are summarized in Section 2. The fast radial symmetry transform is summarized in Section 3, and the details of the proposed approach are presented in Section 4. Our test dataset and evaluation strategy are given in Section 5.The results are reported and discussed in Section 6.The concluding remarks and future directions are provided in Section 7.

PREVIOUS STUDIES
Automated extraction of circular objects from images is an open research area of computer vision and remote sensing (Ok, 2014).Pioneering method in computer vision for detecting circular shapes in images is the Standard Circular Hough Transform (SCHT) (Duda and Hart, 1972).As stated by their method, strong edges in an image space votes for the accumulation in parameter space.This accumulation process generates peaks in the parameter space, where each peak is expected to indicate a circle location.A major disadvantage of SCHT is the 3D parameter space (two parameters for the centre of the circle and one parameter for the radius) used during the accumulation.This is because 3D accumulation requires significant storage to handle large images, and also reduces the efficiency of the processing.Thus, a number of previous studies aimed at reducing the storage/time requirements of the SCHT (e.g.Chiu and Liaw, 2005;Guo et al., 2006;Chung et al., 2012;Huang et al., 2012).Convolution filters for SCHT and the invariance kernels were also proposed as alternative approaches for the search performed in the parameter space (e.g.Atherton and Kerbyson, 1999;Zelniker and Clarkson, 2006;Rhodes and Bai, 2011).Circle detection is not limited to Hough transformation.For instance, (Loy and Zelinski, 2003) incorporated local radial symmetry to detect circles, whereas (Qiao and Ong, 2004) extracted circles by searching meaningful arcs in image space.In a different work, (Chauris et al., 2011) proposed a circlet transform which is formulated in the Fourier domain with specific filters.Neural networks (e.g.D' Orazio et al., 2004), genetic algorithm (e.g.Ayala-Ramirez et al., 2006) and several other optimization routines (e.g.Cuevas et al., 2012) were also tested for circle detection.
There are only a small number of studies dealing with oil depot detection from high resolution aerial/satellite images.(Weisheng et al., 2005) detected storage tanks from SPOT-5 pansharped images using an improved Hough Transform, and a correlation-based template matching.(Li, 2006) tested an approach based on segmentation and feature-based classification, while (Chen, 2009) employed a circle detection algorithm based on shape parameters and a region-growingbased clustering.(Han et al., 2011) utilized Hough transform to detect circles in QuickBird images and a graph-based search is developed to eliminate false detections.Soon after, (Han and Fu, 2012) proposed a saliency model for the detection of circular storage tanks with well-defined boundaries.(Zhu et al., 2012) developed a coarse-to-fine strategy in which the coarse level aimed at finding the image patches with oil tanks, and the fine level detected the circular regions.(Kushwaha et al., 2013) proposed a supervised strategy to detect bright oil tanks from satellite images, where the oil tanks were detected based on the analysis of their statistical and textural information.
A major problem of those above approaches is that the test images involving the oil depots are mostly selected from wellcontrast areas where the roofs of the tanks are consistently bright and the background around the tanks is significantly darker.However, this might not be case for many industrial facilities, especially if the images are not taken under suitable imaging and environmental conditions (Figs. 1b and 1c).In very recent works, two studies have taken these issues into consideration.In the first study, (Ok, 2014) proposed an original approach for the detection of oil depots.In his work, shadow information was the key input to detect circular structures, and the representative boundaries of the shadow areas are used to infer circular regions.However, because the approach completely relies on shadow evidence, the proposed approach may lose its efficiency when this evidence is not visible or complete in image space.In the second study, (Zerman et al., 2014) proposed an algorithm to detect oil depots from high resolution satellite images.In their work, the circular targets were detected using the slightly modified version of the fast radial symmetry transform (FRST) (Loy and Zelinski, 2003), where near-infrared information was used to eliminate the false alarms.However, their approach requires many ad-hoc thresholds, and Loy and Zelinski's FRST may fail to detect circular structures in aerial/satellite images because of a number of reasons (cf.Section 3).

FAST RADIAL SYMMETRY TRANSFORM
For given one or more radii  ∈ , where R is the set of radii, (Loy and Zelinski, 2003) proposed a transform to detect radially symmetric features.In the first step, the transformation computes the gradient of an image (g), generates an accumulated orientation projection image O r (p) and a magnitude projection image M r (p) for each pixel p.The method uses two labels (Fig. 2), positively-affected pixel  + , and negatively-affected pixel  − , to compute the orientation and magnitude images.Thereafter, the orientation and magnitude images are combined into a single projection image  r () using the following formula: where In Eq. 1, α is the radial strictness parameter which determines the level of radial symmetry desired, n is a factor to normalize M r (p) and O r (p) for different input radii.The radial symmetry (S r ) at radius r is computed as a convolution, and the final transform (S) is the average of the radial symmetry considering the set of radii R.
The transform performs well for close range scenes and quite efficient to compute.However, the direct implementation of the approach for aerial/satellite images considering our task (that is the detection of individual oil depots) may not produce satisfactory results because of the following reasons:


The boundary of the top of an oil depot may be composed of both positively-affected pixels and negatively-affected pixels, depending on the illumination and shadow region around the depot.Thus, if that happens, the affected pixels may cancel each other or at least reduce the impact of orientation and magnitude projection values. It is not easy to set the parameters of the transform for different images.The maximum observable radius is large (≈ 50 m) for oil depots in high resolution images and Figure 2. Positively-affected pixel  + , and negatively-affected pixel  − , affected by the gradient element g(p) for a radius of r = 2.The circle shows all the pixels that can be affected by the gradient at p for a radius r (Loy and Zelinski, 2003).
therefore, setting a small subset of radii might not represent the whole possible set R of oil depots.The normalization factor n must be learned from very large training sets, which still might not yield an appropriate single value to represent the whole image set. It is not easy to achieve a valid threshold for S to locate the centre locations of circles, since an average of radial symmetries is computed from a set of S r images considering the set of input radii R.  It is not directly possible to define circle radius from S. This is because each radial symmetry (S r ) is computed using a convolution kernel depending on the radius r, and the final symmetry map is produced after combining all S r .Thus, different radii information is also accumulated in S.


It is not possible to evaluate the number of pixels at certain locations and specific radii which contribute to the projection matrices.This is because of the accumulation, and such information, if evaluated, could be helpful to eliminate false alarms, e.g. a false alarm object with clear boundary (large gradient) that contains an arc in a small part of its boundary could be identified.
Because of above-stated reasons, the direct implementation of the fast radial transform might not produce the desired results for aerial/spaceborne datasets.Our approach intends to alleviate the weaknesses of FRST defined in this section.In this respect, we propose a new strategy (cf.Algorithm 1).

PROPOSED APPROACH
Our approach starts with the computation of the gradient image.
We utilize the 5-tap coefficients introduced by (Farid and Simoncelli, 2004) for better estimation of the magnitude and orientation components of the gradient.Next, we normalize the orientation component to achieve unit direction vectors,  ⃗ .We define the range of minimum and maximum radii (r min ,r max ) after taking into account all values of radii of a single oil depot can retain, R = {r min , r min +1, … , r max }.At each radius r, first, we compute positively-affected pixel  + , and negativelyaffected pixel  − : This will allow us to compute positively-and negativelyaffected orientation and magnitude images: In Eq. 4,  r + and  r − denote the positively-and negativelyaffected orientation images, respectively, and  r + and  r − denote the positively-and negatively-affected magnitude images, respectively.Thus, the orientation O r (p) and magnitude M r (p) images (Fig. 3) can be computed as Note that, subtractions in Eq. 5 intensify the orientation and  where  r  and  r  denote the maximum values in the magnitude and orientation projection images, respectively.Finally we compute  r (Fig. 4) using a convolution (*) where A is a Gaussian kernel with a fixed σ.Once we compute the radial symmetry image S r at a certain radius r, we expect to achieve higher pixel values around certain locations which probably indicate circle centres (Fig. 4b and 4c).Thus, to focus on these high valued regions, we propose to use an automated thresholding of S r with the following rule: where    is the computed threshold value, (. ) denotes the median operator, and α ≥ 2. In this way, we compute the interest regions after thresholding the radial symmetry image (S r >    ) and finding the connected regions (b j ) in the resulting binary image using a connectivity of 8 neighbouring pixels (Fig. 4d  and 4e).
The pixels belonging to a connected component ( j ) involve a notable accumulation of symmetry which might correspond to centre location (  j ) of a candidate circular object in image space.The maximum affected pixel of each connected component must denote the corresponding centre location of the circle.Thus, we search for the maximum of  r ( j ) and label that pixel as centre of a candidate circle (  j :    ( j )).To understand whether the accumulation arises from a certain amount of pixels around the radius r of the circle centre   j , we compute the reverse of positively-and negatively-affected orientation images.To do that, we check and label all pixels which indicate pixel   j and its 8 neighbouring pixels to form the orientation image.Thus, we are able to trace all pixels supporting the candidate circle.As expected, these supporting pixels (  j ) must be around the boundary of the candidate circle; thus, a new measure, support ratio (ρ j ), can be computed.First we find the ideal boundary pixels (  j ) of the candidate circle estimated from [  j , ], and include its 8-neighbours to form a search region (  j ) .Thus, the support ratio (ρ j ) can be computed as where |(  j )| and |(  j )| denote the number of supporting pixels in the search region and cardinality of the search region, respectively.If the candidate circle belongs to a circular object in image, we expect to have relatively large support ratio (ρ j ≥ ρ min ).Otherwise, the candidate circle is not verified.
As a final step, a very important issue is to check whether a circle belonging to an oil depot is previously detected in any part of the pixels of the candidate circle.This test is not trivial, because some roofs of the oil depots may contain concentric evidences, and a number of radii may point to a single centre location.For that purpose, we collect all previous candidates ( k ) (if exist) and check whether the current candidate ( j ) has the highest projection value  r () (Eq.6).Note that this check must be performed without any normalization because normalized projection image (  r () ) takes into account the maximum values ( r  and  r  ) considering the projection images; thus, it does not allow us a direct comparison between  j and  k .Therefore, we use the projection image  r () and apply the test    ( j ) >    ( k ).If this test is passed, we verify  j as a candidate circle in image space (  j ∈  ) and remove the other candidate circles ( k ∉ ) belonging to that region, where C indicates the set of circles in image space.The processing loop continues until the last radius (r max ) is processed in R. Thus, after that step, all circle candidates remaining in set C are accepted as circular objects (Fig. 4f).

DATASET, EVALUATION, AND PARAMETERS
We selected 6 GeoEye-1 test images (0.5 m GSD) to assess the performance of our approach (Fig. 5).Besides, our dataset further involve 16 Google Earth images (RGB) of a part of an industrial facility (Latitude: 40.632811°,Longitude: -74.228924°) collected in a time series manner (between the years 1995 and 2012) (Fig. 6).These Google Earth images are taken from both aerial and spaceborne platforms with varying ground sampling distance (GSD).In this study, GSD of all 16 images is compiled as 0.5 m for numerical assessment, regardless of their original GSD.Besides, Google Earth images compiled with different GSD were also evaluated to investigate the robustness of our approach for scale variations.
The reference images of each image were prepared manually by a qualified human operator.In our case, a pixel-based assessment is not representative, because correct detection of relatively large oil depots in a scene always leads to superior pixel-based measures.Therefore, we prefer an object-based assessment which consists of two stages: (i) finding the one-toone matches between the reference circle objects and the circle objects in the output, and (ii) evaluating the performance of the  (Atherton and Kerbyson, 1999), and (third column) results of the proposed approach.Green, red and blue colours represent TP, FP and FN pixels, respectively.approach using the matches established.We assign one-to-one matches in a way that the difference between the centre positions of the objects in the reference and the output is minimized.Once all one-to-one matches are established, we follow the well-known three measures (Precision, Recall, and  1 -score) to evaluate the object-based performance of the proposed approach (Ok, 2014): In Eqs. 10 -12, TP denotes an output object that takes a part in the final matching list established through one-to-one mapping, FP indicates an object found by the proposed approach but does not participate in any matches, FN implies a reference object that does not involve in the matches, and the operator | .| is the set cardinality.In Eq. 12, the F 1 -score can be used to evaluate the overall object-based performance.
One-to-one matches between the objects in the reference and the output can also be used to assess the quality of the detected circles in terms of their parameters.To do that, we compute the Euclidean distance between the centre pixels of each final match together with the Euclidean differences of their radii, and that information is used to compute the root mean square (RMS) distances between the reference and output circle objects.
The implementation and processing was performed in MATLAB.All experiments were performed on a computer with an Intel i7 processor with 2.40 GHz and 16 GB RAM.Our approach contains four parameters (Table 1).The experiments show that a single parameter set is sufficient to properly handle very different datasets, and therefore, we fix each parameter to a constant.We set the minimum radius to a value where the accumulation starts producing reliable results for circle detection (that is r min = 3 pixels), and the maximum radius to the largest observable radius for oil depots (that is r max = 100 pixels).We fix the radial strictness parameter α to 4 because we observed that setting provided the best balance between precision and recall.We preferred isotropic filtering presented by (Geusebroek et al., 2003) to efficiently compute the Gaussian smoothing, where we set σ = 4. Finally, we fixed the support ratio, ρ min = 0.4, by which we seek at least 40% supporting evidence around the circle boundary.

RESULTS AND DISCUSSION
We visualize the detection results of GeoEye-1 test images in Fig. 7.These results demonstrate that our approach can provide remarkable results for the detection of oil depots.Our results are quite good and representative, and almost all oil depots are very well detected despite their complex characteristics, e.g.roof colour, texture, and size.The numerical results in Table 2 favour these facts.We achieved overall precision and recall as 97.3% and 94.4%, respectively.The computed F 1 -score for these six test images is around 96%.We computed RMS of circle centre differences of reference and output images between 0.65 and 1.6 pixels with an overall centre RMS of 1.11 pixels (Table 2).Besides, the RMS of radius distances are almost better for each test image, leading to a total radius RMS value of ≈ 0.9 pixel.
The results of Google Earth time series investigation are shown in Fig. 6, and all related numerical results are presented in Table 3.The outcomes in Fig. 6 show how our approach can provide a useful output for the identification of oil depots even in very different/difficult illumination and acquisition conditions.The proposed approach has distinct capability to focus oil depots that are characterized by severe environmental conditions, e.g.snow cover (#3) and smoke (#14-15).In addition, note that our approach is not sensitive to the type of platform (aerial/spaceborne) used (as long as the images are collected from a near-nadir viewing angle).According to the numerical results in Table 3, we achieved overall precision and recall as 94.7% and 84.6%, respectively.The computed F 1 -score for these 16 test images is 89.4%.However, we also state that some of the F 1 quality measures are relatively low (≈ 80%, e.g.#1, #5, #7, #8, and #15).This is mainly because of small sized oil depots whose boundaries are missed by the proposed approach.
Overall RMS of centre differences of reference and output circle objects is nearly 5 pixels (Table 3).Such poor result, of course, mostly due to three images, #12, #14 and #15 whose centre RMS values are computed above 7 pixels.If we investigate the results of these images from Fig. 6, the problem is apparent: for a number of oil depots (one depot in #12, five depots in #14, and one depot in #15), the exact positions of rooftops could not be found.For images #14-15, this is related   Table 3. Performance results of Google Earth test images.
to smoke over the depots during image acquisition, which significantly weakens the gradient magnitude around the roof boundary of the depot.For image #12, we observe blurred boundaries (probably because of relatively low GSD in original image), where the gradient magnitude of a shadow boundary around the depot produced a stronger radial symmetry than the rooftop boundary.This, of course, is also one of the main reasons of false positive objects shown in several images (upper-left of images #2, #8, and #12).Despite these problems, we believe that centre locations of rooftops are achieved close to their reference positions in most of the test images.The RMS results are much better for radius parameter, with an overall performance of ≈ 1.7 pixels.As one can easily predict, test image #15 is the worst case (3.12 pixels) yet again because of heavy smoke visible over the oil depots.
Figs. 8 and 9 show the robustness of our approach to scale variations.As seen in those figures, our detection performance slightly drops when the images are compiled with relatively coarse GSD.This is because of two reasons: (i) it is not possible to detect small-sized depots with a coarser GSD (e.g. 2 m), and (ii) the orientation of the gradient computed from a coarse resolution image is not as accurate as the orientation component generated from a high resolution image.However, we believe that more tests on this topic are necessary to expose a strong conclusion.
Our approach successfully alleviates the weaknesses of FRST summarized in Section 3. In the proposed approach, positivelyaffected and negatively-affected pixels support each other to achieve improved orientation and magnitude images.Besides, all parameters of the approach can be fixed for different datasets.In addition, the new approach offers an automated way to focus centre locations and radii of circular structures.Finally, a new efficient test is employed to remove false alarms.
We compare our results with the approach presented in (Atherton and Kerbyson, 1999).Comparative assessments in Fig. 6 and Table 3 prove that the results provided by the proposed approach are more reliable and robust.The processing times required by the proposed approach are provided in Table 4.According to the computational times computed, it is possible to detect circular depots from images with sizes 1500x1500 pixels approximately in one minute with our approach.

CONCLUSIONS AND FUTURE WORK
In this paper, a new approach is presented to automatically detect circles belonging to oil depots in industrial environments from a single satellite image.Assessments performed on six GeoEye-1 test images and 16 Google Earth time series images reveal that our approach is capable of detecting circles in very different/difficult images.The time series investigation reveals that our approach is robust enough to monitor oil depots in industrial environments under varying illumination and environmental conditions.
In the future, we will focus more to develop an automated validation strategy for the detected circles to understand whether they belong to a depot rooftop or not.In this respect, it might be useful to integrate the shadow evidence (Ok, 2014) or the evidence of statistical and structural characteristics of oil depots to discover their compound structures (Akcay and Aksoy, 2011).Besides, a different interesting tasks is to separate oil depots based on their types (e.g.flat roof or sphere etc.) and try to understand (to some extent) whether the depot is full or not; therefore, we plan to expand our studies in these research directions.Finally, our approach can be extended to ellipses (Ni et al., 2012), and it might also be possible to process oblique images.
Figure 1.Three test images from Google Earth time series images of an industrial facility.(a) Panchromatic aerial image acquired on 3.29.1995,(b) and (c) satellite images (RGB) acquired on 2.1.2004and 11.3.2012,respectively.
(a) GeoEye-1 test image, (b) orientation O r and (c) magnitude M r images (r = 45 pixels).Note that the images are enhanced for better visualization.magnitude images.Thereafter, we compute the projection image  r () and its normalized form  r () using  r () =  r (). r ()  ( (a) GeoEye-1 test image,  r image for (b) r = 45 pixels and (c) r = 90 pixels.Interest regions of (b) and (c) are shown in (d) and (e), respectively and (f) final result of the approach.Note that the S r images are enhanced for better visualization.
Figure 6.(first column) Google Earth test images (0.5 m GSD,#1-8), (second column) results of the approach presented in(Atherton and Kerbyson, 1999), and (third column) results of the proposed approach.Green, red and blue colours represent TP, FP and FN pixels, respectively.
Effect of different image scales: Google Earth test image #9 compiled at (a) 0.5 m, (b) 1 m, and (c) 2 m GSD, and the results are shown in images (d), (e), and (f), respectively.Effect of different image scales: Google Earth test image #13 compiled at (a) 0.5 m, (b) 1 m, and (c) 2 m GSD, and the results are shown in images (d), (e), and (f), respectively.

Table 4 .
Computational time required by the proposed approach for GeoEye-1 test images.