COASTLINE DETECTION USING FUSION OF OVER SEGMENTATION AND DISTANCE REGULARIZATION LEVEL SET EVOLUTION

Coastline detection is a very challenging task in optical remote sensing. However the majority of commonly used methods have been developed for low to medium resolution without specification of the key indicator that is used. In this paper, we propose a new approach for very high resolution images using a specific indicator. First, a pre-processing step is carried out to convert images into the optimal colour space (HSV). Then, wavelet decomposition is used to extract different colour and texture features. These colour and texture features are then used for Fusion of Over Segmentation (FOOS) based clustering to have the distinctive natural classes of the littoral. Among these classes are waves, dry sand, wet sand, sea and land. We choose the mean level of high tide water, the interface between dry sand and wet sand, as a coastline indicator. To find this limit, we use a Distance Regularization Level Set Evolution (DRLSE), which automatically evolves towards the desired sea-land border. The result obtained is then compared with a ground truth. Experimental results prove that the proposed method is an efficient coastline detection process in terms of quantitative and visual performances.


INTODUCTION
The identification of a coastline involves two steps.The first needs the definition of an indicator (also known as a key) that will represent the land-sea boundary.The second step includes detection of the chosen indicator within the available data source (Boak, E. H. and Turner, I. L. , 2005).The definition of an apparently simple notion such as the coastline, which is assumed to represent a linear boundary between land and sea, is very challenging because of the wide range of indicators that can be defined.For this reason, Robin (Robin, n.d.) found more than a dozen of indicators, while Boak and Turner (Boak, E. H. and Turner, I. L. , 2005) identified 19 generic lines from 45 different indicators (Faye et al., 2011).This diversity of indicators leads to the development of numerous methods for coastline detection, which can be achieved in various steps like pre-processing, segmentation and edge detection (Zhang et al., 2013).Pre-processing methods can be grouped into two categories, which are noise reduction and image correction.For multispectral imagery, the normalized difference water index (NDWI) proposed in (Mcfeeters, 1996) and the superfine water index (SWI) recently introduced in (Sharma et al., 2015) are important metrics for classification.For panchromatic images, the K-means and ISODATA methods are the most popular clustering methods (Liu et al., 2011) (García-Rubio et al., 2015).Nevertheless the initial positions of the cluster centroids guide the final K-means clustering result Ursani et al. (Ursani et al., 2009) proposed an approach, referred to as FOOS, which is a mean of finding the K optimal centroids initialisation.We implemented a FOOS based clustering by integrating colour and textures features obtained from wavelet decomposition.
For the edge detection step, the literature lets us to appreciate different approaches.The most famous of these approaches is the Canny edge detector (Canny, 1986).The snakes (Kass et al., 1988) and the level set (Pro, 2008) (Liu et al., 2016) have also been used, but, in general, for SAR imaging.Compared with snakes' method, in level set algorithm, the initial contour is controlled to move automatically to the edges without rigorous limits.However, level set evolution is slow, particularly in high-resolution images due to lots of variables and complicated iterative methods.This leads to the development of a new type of level set evolution method called distance regularized level set evolution (DRLSE) (Li and Xu, 2010).The distance regularization effect eliminates the need for re-initialisation, thus avoiding numerical errors.DRLSE also accepts the use of more wide-ranging and efficient initialisation parameters of the level set function.This leads to reduce the number of iterations, while ensuring sufficient numerical accuracy in order to improve the coastline detection as shown in the results of this study.The rest of the paper is as follow.Section 2 presents the proposed method.Section 3 shows a validation step by comparing the results of the proposed method with the ground truth.The main conclusion of this work is recapped in the last section.

PROPOSED METHOD
The block diagram in Figure 1 represents a general scheme of the proposed method.The proposed method is implemented as a 3-step algorithm: (1) pre-processing, (2) segmentation, and (3) edge detection.The pre-processing step is used to achieve image transformation into HSV colour space and to extract colour texture features using the coefficients of the image decomposition by wavelet transform.The segmentation step aims to transform the pre-processed image in order to discriminate the different regions of this image.This allows isolating the dry sand from the wet sand.At the edge detection step we use the DRLSE method to extract the coastline.The data used in this study are downloaded from the Google Earth Pro application; see an example on figure 2. They are provided by DigitalGlobe and derived from pan-sharpening process.Google Earth imagery has been used in land-sea segmentation for coastline detection (Cheng et al., 2016), and in mapping water depth and land cover (Collin et al., 2014).To detect the coast-Figure 2. Image of test line indicator, we use both colour and texture features.Texture analysis was initially performed on grayscale images, thus discarding colour information.Nonetheless, many new works have been made to show the importance of colour in texture analysis (Ilea andWhelan, 2011) (González-Rufino et al., 2013).In this sense, to increase the accuracy of a good segmentation and a perfect detection of the coastline, we carried out an experimental study on different colour spaces: RGB, XYZ, Lab, HSV, YCbCr.By comparing the different results of this study, we have seen that the HSV space offers the best clustering performance while using colour and texture features; which is perfectly coherent with the results presented in the literature.In fact, Paschos (Paschos, 2001) experimentally analysed the impact of colour space (RGB, Lab and HSV) on colour and texture classification using a dataset of 50 colour images and as texture descriptors the Gabor filters.He concluded that HSV is the best colour space for texture classification.The HSV colour model defines colours in terms of hue, saturation, and value (or brightness).The colour and texture features considered in this study are generated using the wavelet transform, which makes it possible to perform a multi-scale analysis of the local structures (Mallat, 1989).The wavelet transform decomposes the signal into a family of translated and scaled wavelets.As signified in (1), a wavelet is a function of zero average.
For 1-dimension, the continuous wavelet transform (CWT) is written as shown in (2).
where ϕ (a,b) is a wavelet computed from the mother wavelet Ψ by translation and scaling as shown by equation ( 3).
where a and b are positive integers and refers respectively to the scale and the translation of the wavelet.The Discrete Wavelet Transform (DWT) is an implementation using a discrete set of scales and wavelet translations following certain rules.It can be considered as an all-pass filter, in which a two-band Quadrature Mirror Filter (QMF) bank uses orthogonal analysis filters in order to decompose data into low-pass and high pass bands.
The high-pass and the low-pass filters provide respectively the details (texture) and an approximation of the image.
In this study, we use the Haar wavelet function, which is given by (4).
Where, r and c refer respectively to the row and column number in the image f(x,y).j0 represents the arbitrary starting scale.Wϕ(j0, r, c) defines an approximation of the image f (x,y) at scale j0.W i Ψ (j, r, c) defines diagonal, vertical and horizontal details for scales j ≥ j0 .i= H, V, D .Such that Ψ H indicates variations measured in horizontal orientations, Ψ V means variations in vertical orientations and Ψ D shows the variations along the diagonal orientations.The criterion of the conservation of the spatial variability of the texture intensity is used to define the level of composition.This allowed us to conclude that a one-level decomposition is sufficient because beyond that, texture becomes too homogeneous (too smoothed).The coefficients obtained with the approximation and the different details of the image are then used in order to carry out a segmentation.

Segmantation stept
K-means clustering is widely used in remote sensing image processing for various tasks such as image segmentation, image classification, and image compression, due to its simplicity of implementation.It organizes the data of an image into K clusters.The algorithm returns a data partition, in which the objects within the same cluster are as close as possible to each other and as far as possible from the objects of the other clusters.Each cluster of the partition is defined by its objects and its centroid.The initial position of the centroids conditions the final result.So, to solve this problem various methods were proposed.These methods are based on centroids initialization or multiple restart of the algorithm.Ursani et al. (Ursani et al., 2009)proposed an approach, referred to as FOOS (Fusion Of Over Segmentation), which is a mean of finding the K initial cluster centres that, unlike McQueen initialization, lead to a near-optimal solution, even if the desired clusters do not lie in the spherical Gaussian distribution.We propose a FOOS-based land-sea segmentation by integrating colour and textures features obtained from wavelet decomposition.The FOOS method for clustering into K clusters consists of the following steps: • Fusion of the two clustering in order to obtain P×Q classes • Select the K greatest classes common to both clusterings.
• Take centroids of these K classes as the initial centroids of the clustering K classes and rerun the k-means.
Besides, we introduce median filtering before the fusion of the two over-segmentation results.The median filter is used in this case to remove the isolated points before merging the two oversegmentations.This makes it possible to reduce errors considerably and to find the optimal centre of gravity used as centroids for the final k-means clustering.The clustering result is represented by the figure 3 where, we can clearly discriminate the dry sand and the wet sand.In the following step, a mean of finding the boundary between these two classes will be presented.

Edge detection stept
For the detection of the aforementioned boundary, we use a level set formulation that gives an active contour model, which is called an implicit active contour or geometric active contour model.For the level set formulation of parametric active contour, the readers are referred to (On the relationship between parametric and geometric active contours, 2000).The idea of the level set method is to represent an edge as the zero level set of a higher dimensional ) and to formulate the motion of the contour as the evolution of this LSF.However, the LSF develops irregularities during its evolution.This destroys the stability of the level set evolution.To ensure the LSF stability, reinitialization, can be applied to periodically replace the degraded LSF with a signed distance function; see (Sethian, 1999) for more details.Reinitialization can also be performed by using the fast marching method as shown in (Osher and Fedkiw, n.d.).
Here, we use a variational level set formulation in which the regularity of the LSF is intrinsically maintained during its evolution.
The level set evolution is derived as the gradient flow that minimizes an energy functional with a distance regularization term and an external energy that drives the motion of the zero level set toward desired locations.This process is called distance regularized level set evolution (DRLSE) (Li and Xu, 2010).In this formulation, an energy functional of a LSF φ, defined in an area Ω, is expressed as shown in ( 7): where Rp(φ) is the regularization term.It is defined with a potential function p such that the derived level set evolution has a unique forward-and-backward (FAB) diffusion effect, which is used to maintain the desired shape of the level set function, particularly a signed distance profile near the zero level set.The distance regularization effect eliminates the need for reinitialization and thereby avoids its induced numerical errors.Lg(φ) and Ag(φ) refer to the external energy.The energy Lg(φ) computes the line integral of the function g along the zero level contour of φ.!It is minimized when the zero level contour of φ is located at the coastline.The energy functional Ag(φ) is introduced to speed up the motion of the zero level contour in the level set evolution process, which is necessary when the initial contour is placed far away from the coastline.µ ,λ and α are the weighting coefficients of the following functions: For more details to the definition of φ0(x) , readers are referred to (Li and Xu, 2010).For the implementation, we use the same parameter as in (Li and Xu, 2010).µ= 0.2, λ=5, α=-3 and ∆t =1, but the initial contour is manually set.R0 is the part of Ω were φ0(x) < 0 ; see the red line in figure 4. For the narrowband implementation of the DRLSE , (Li and Xu, 2010) denote by φi, j a LSF defined on a grid.A grid point (i,j) is called a zero crossing point, if either φi−1,j and φi+1,j are of opposite signs or φi,j−1 and φi,j+1 are of opposite signs.The set of all the zero crossing points of the LSF is denoted by Z.Then, the narrowband is constructed as shown in ( 14) where N (r) (i,j) is a (2r + 1)(2r + 1) is square block cantered at the point (i,j), r can be set to be the smallest value r = 1, in which case the narrowband Br is the union of the 3×3 neighbourhoods of the zero crossing points.The results of this process are represented in figure 5.The narrowband implementation consists of the following steps: Step 1) Initialize the LSF φ(x) to the function φ0(x) Then, create the initial narrowband B 0 r = (i,j)∈ Z 0 N (r) (i,j) , where Z 0 is the group of the zero crossing pixels of φ0.
Step 2) Update φ k+1 i,j = φ k i,j +τ L (φ k i,j ) on the narrowband B k r as suggested by the equation ( 20) of (Li and Xu, 2010) Step 3) Determine the set of all the zero crossing pixels of φ k+1 i,j on B k r , denoted by Z k+1 .Then, update the narrowband by setting Step 4) For every point (i,j) in B k+1 r , but not in B k r , set φ k+1 i,j to h if φ k i,j > 0 , or else set to φ k+1 i,j − h , where h is a constant, which can be set to r + 1 as a default value Step 5) Terminate the iteration if the zero crossing pixels stop chaging for ?m consecutive iterations, or k exceeds a predifined maximum number of iterations, otherwise, go to Step 2.

EXPERIMENTAL RESULT AND DISCUSSION
An automated or semi-automated detection of a coastline ought to be validated using a ground truth.This ground truth can be achieved using terrestrial surveys process or a manually detected and validated coastline.In this study, the deficiency of a ground truth motivated us to generate our own one.This ground truth is obtained by averaging several manual detection results achieved by different persons; see figure 6.Since a manual detection of an objet depends on the visual capability of the operator, we collected different outlooks and then took the mean of them.By analysing the figure 7, one can clearly find that the red line is very close to the blue one.Visually, it is very difficult to make a difference between them.This shows a good correlation between the two coastlines.
By looking closer, we can see a slight offset between the two curves in figure 8.We notice that the semi-automated coastline vary along the length of the image, but over the width, it is concentrated between the pixels of rank between 260 and 275.This is quite conformed to the ground truth.
Figure 8.Our approach vs manually detected coastline The computed Root Mean Square Error (RMSE) is about 0.3716, corresponding to an offset of 0.3716× 0.6 = 0.2m.0.6 m is in fact the image spatial resolution.This value of RMSE let us approve the good correlation between the two coastlines.To approve the goodness of the proposal method, we found the coastline with different values of the initial contour.
In figure 9, the initial contour is set right to the coastline.The domain R0, previously define has moved.Despite this change of the initial contour,we found the same results with a good detection.Therefor, the initial contour ought to be along the coastline to ensure a good detection.

CONCLUSION
The current study provides a semi-automated coastline detection method using very high resolution images like those from Google Earth.It is done by image processing and computer vision techniques aiming at region segmentation and edge detection.We integrate both colour and texture information and perform a multiresolution analysis based on the wavelet transform.After the preprocessing step, which aims to find the best features for segmentation, we apply image segmentation in order to split the image into different regions.The process of DRLSE is adopted to improve the accuracy of the extracted coastline.The obtained results clearly show the effectiveness of this approach.Likewise, many works show coastline detection process, without specification of the key indicator that is used.The proposal brings a contributions in this sense.As for future work, we intend to conduct a diachronic study over several years to describe the coastline dynamics and thus to appreciate how it has moved over time

Figure 1 .
Figure 1.Block diagram of the proposed method

Ω
gH(−φ)dx (10) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-3/W4, 2018 GeoInformation For Disaster Management (Gi4DM), 18-21 March 2018, Istanbul, Turkey g is an edge indicator of the analysed image I. is the Dirac function and H(−φ) the Heaviside function.Gσ a Gaussian filter which is used to reduce the image noise.The equation of the curve evolution is thus obtained by deriving the total energy as shown in (12): ∂φ ∂t = µdiv(dp(|∇φ|)∇φ) + λδ (φ)div(g ∇φ |∇φ| ) + αgδ (φ) (12) We use LSF that take negative values inside the zero level contour and positive values outside.The initial LSF is defined by the function (13). φ0

Figure 9 .
Figure 9. initialization with a different value of φ0(x)