1 Introduction

River channel extraction is an important step in remote sensing river image processing and has great significance for investigating and monitoring of water resources, climate monitoring, map matching and ship navigation [1]. Traditional measurement consumes lots of human and material resources, and it is inappropriate for timely monitoring. In recent decades, with the rapid development of remote sensing technology, we can achieve remote sensing images in the real time [2], and moreover, remote sensing imaging cycle is short, has high resolution and of low cost. Specially, SAR imaging technology has advantages of all-weather all-day, and strong penetration. Therefore, extracting river information from SAR image has become an important approach of river channel extraction.

Recently, active contour model has been studied and used for target extraction widely. The model divides the image by minimizing the gray-level difference of the pixel in the regions [3]; Euler–Lagrange equation guides the evolution of level set function. The advantages of active contour model are that the zero level set function splits or merges, while the topology of evolving curve changes [4, 5].

The Chan–Vese (CV) model [4] is one of the most popular region-based active contour model and segments the image with two regions which have different grayscale information [6,7,8,9,10,11,12,13], but it is not applicable to the image which presents intensity inhomogeneity. Recently, numbers of improved models emerged [14, 15] to solve these shortcomings. Niu et al. [16] proposed a region-based model for the segmentation of objects or structures in images by introducing local similarity factors, which relies on the local spatial distance within a local window and intensity difference, and the model has the disadvantages of the high time complexity and the sensitive initialization of the contour. With the purpose of coping with noisy images, numbers of region-based models have been presented, such as local binary fitting model [17], local image fitting model [18], region-scalable fitting model [19], local and global intensity fitting model [20], but they are sensitive to the initialization of active contour and easy to fall into the local optimum.

The paper proposes an improved active contour model for SAR river channel extraction. The regional energy function is defined by the logarithm kernel distance which is robust and non-Euclidean. In addition, the local median information of image is considered rather than global image statistics, and the fitting center is calculated by the mean of local median information. Furthermore, fixed energy weight affects the curve evolution, and the interior and exterior energies are calculated by the variation which integrates every pixel with fitting center in the paper. Finally, distance regularized term is introduced into the objective function to accelerate the convergence of the model.

The rest of the paper is organized as follows: Sect. 2 presents the region-based active contour model, which is driven by the logarithm kernel-based energy function; some experiments applying to the proposed model and other state-of-the-art models in SAR river channel extraction are presented in Sect. 3; Sect. 4 concludes the paper.

2 The region-based active contour model

2.1 Active contour model driven by the logarithm kernel-based energy function

Let \( \varOmega \in R^{2} \) be a bounded image domain, suppose that the pixel I(x,y) can be segmented into two classes which are \( \varOmega_{1} \) and \( \varOmega_{2} \), and then, define the energy function as:

$$ E = \sum\limits_{i = 1}^{M} {\int_{{\varOmega_{i} }} {\left\| {\varphi \left( {I\left( {x.y} \right)} \right) - \varphi \left( {c_{i} } \right)} \right\|} }^{2} {\text{d}}x{\text{d}}y $$
(1)

where \( \varphi \left( \cdot \right) \) is kernel function and \( c_{i} \) represents the fitting center of \( \varOmega_{i} \).\( \left\| {\varphi \left( {I\left( {x,y} \right)} \right) - \varphi \left( {c_{i} } \right)} \right\|^{2} \) represents the distance between I(x,y) and \( c_{i} \). The kernel distance has the following advantages: It has nonlinear transformation, transformed space can be linearly classified, and the distance is more robust and can improve the accuracy of distance measure. Here, kernel distance can be represented as:

$$ \left\| {\varphi \left( {I\left( {x,y} \right)} \right) - \varphi \left( {c_{i} } \right)} \right\|^{2} = \left( {\left\langle {\varphi \left( {I\left( {x,y} \right)} \right),\varphi \left( {I\left( {x,y} \right)} \right)} \right\rangle - 2\left\langle {\varphi \left( {I\left( {x,y} \right)} \right),\varphi \left( {c_{i} } \right)} \right\rangle + \left\langle {\varphi \left( {c_{i} } \right),\varphi \left( {c_{i} } \right)} \right\rangle } \right) $$
(2)

where define K(x,y) as the inner product of x and y, and then, Eq. (9) can be transformed as:

$$ \left\| {\varphi \left( {I\left( {x,y} \right)} \right) - \varphi \left( {c_{i} } \right)} \right\|^{2} = {\kern 1pt} \left( {K\left( {I\left( {x,y} \right),I\left( {x,y} \right)} \right) - 2K\left( {I\left( {x,y} \right),c_{i} } \right) + K\left( {c_{i} ,c_{i} } \right)} \right) $$
(3)

\( K\left( {m,n} \right) = \left\langle {\varphi \left( m \right),\varphi \left( n \right)} \right\rangle \). In this paper, define \( K\left( {m,n} \right) \) by logarithm kernel function, where

$$ K\left( {m,n} \right) = - \;\log \left( {1 + \left( {m - n} \right)^{d} } \right)\quad {\kern 1pt} d \ne 0 $$
(4)

As \( K\left( {I\left( {x,y} \right),I\left( {x,y} \right)} \right) = K\left( {c_{i} ,c_{i} } \right) = 0 \), Eq. (3) can be transformed as:

$$ \left\| {\varphi \left( {I\left( {x,y} \right)} \right) - \varphi \left( {c_{i} } \right)} \right\|^{2} = - \;2K\left( {I\left( {x,y} \right),c_{i} } \right) $$
(5)

Then, define the energy function as:

$$ E = \sum\limits_{i = 1}^{M} {\int_{{\varOmega_{i} }} {\left\| {\varphi \left( {I\left( {x,y} \right)} \right) - \varphi \left( {c_{i} } \right)} \right\|} }^{2} {\text{d}}x{\text{d}}y = {\kern 1pt} \sum\limits_{i = 1}^{M} {\int_{{\varOmega_{i} }} {\left( { - 2K\left( {I\left( {x,y} \right),c_{i} } \right)} \right)} } {\kern 1pt} ,\quad M = 1,2 $$
(6)

In order to minimize the energy function, the level set function \( \phi \left( {x,y} \right) \) is introduced [5], and the total energy function which depended on logarithm kernel function takes the following form:

$$ F(c_{1} ,c_{2} ,C) = \mu \int\limits_{\varOmega } H^{{\prime }} (\phi )\left| {\nabla \phi } \right|{\text{d}}x{\text{d}}y + v\int\limits_{\varOmega } {H(\phi ){\text{d}}x{\text{d}}y} - \;2\lambda_{1} \int\limits_{{\varOmega_{1} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - c_{1} } \right)^{d} } \right)H(\phi ){\text{d}}x{\text{d}}y} - \;2\lambda_{2} \int\limits_{{\varOmega_{2} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - c_{2} } \right)^{d} } \right)(1 - H(\phi )){\text{d}}x{\text{d}}y} $$
(7)

The logarithm kernel distance, nonlinear transformation, maps from low-dimensional space to high-dimensional space. In addition, compared with traditional CV model, if the pixel I(x,y) (such as noise point) has great differences with \( c{}_{i} \), the external and internal energy would have significant value variation. However, because of the logarithmic computation, while I(x,y) has great differences with \( c{}_{i} \) in Eq. (7), the external and internal energy \( \int\limits_{{\varOmega_{2} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - c_{2} } \right)^{d} } \right)(1 - H(\phi )){\text{d}}x{\text{d}}y} \) and \( \int\limits_{{\varOmega_{1} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - c_{1} } \right)^{d} } \right)H(\phi ){\text{d}}x{\text{d}}y} \) would have gentle variation. Therefore, the active contour model suppresses noise well; moreover, logarithm kernel distance would have influences on the curve evolution and improves the accuracy of image segmentation.

2.2 Active contour model based on the local median information

In the traditional CV model, \( c_{1} \) and \( c_{2} \) are the average intensities inside and outside the curve, while the image includes some noise; \( c_{1} \) and \( c_{2} \) ignore local information of the image, and then, they are not accurate. Therefore, the calculation of \( c_{1} \) and \( c_{2} \) based on local median information of the whole image is proposed in this section.

I’(x,y) can be defined as the median point in a neighbor template of 5*5 size, and the weights of the neighbor pixels are set to 1. The pixels in the template are arranged in descending order, and \( I^{'} \left( {x,y} \right) \) is in the middle of the order. Then, the calculation of I’(x,y) is expressed as:

$$ I^{{\prime }} \left( {x,y} \right) = {\text{med}}\left\{ {I\left( {x \pm k,y \pm k} \right),\left( {k \in W} \right)} \right\} $$
(8)

The fitting center \( e_{1} \) and \( e_{2} \) can be defined as the mean of I’(x,y), which can be given by:

$$ \left\{ {\begin{array}{*{20}l} {e_{1} = \frac{{\int_{\varOmega } {I^{'} \left( {x,y} \right)H\left( \phi \right){\text{d}}x{\text{d}}y} }}{{\int_{\varOmega } {H\left( \phi \right){\text{d}}x{\text{d}}y} }}} \hfill \\ {e_{2} = \frac{{\int_{\varOmega } {I^{'} \left( {x,y} \right)\left( {1 - H\left( \phi \right)} \right){\text{d}}x{\text{d}}y} }}{{\int_{\varOmega } {\left( {1 - H\left( \phi \right)} \right){\text{d}}x{\text{d}}y} }}} \hfill \\ \end{array} } \right. $$
(9)

The fitting centers merge local median information and global mean information of the whole image, which are more accurate than \( c_{1} \) and \( c_{2} \). Additionally, the fitting centers reduce the influences of intensity inhomogeneity which is presented in each region and improves the accuracy of image segmentation.

2.3 Active contour model based on adaptive energy weight

\( \lambda_{1} \) and \( \lambda_{2} \) are fixed constant parameters in the traditional CV model; however, the interior and exterior energies are varying in the curve evolution. Constant parameters not only have influences on the result of image segmentation, but also affect the speed of curve evolving. The interior and exterior energies are weighted by the variation which integrates every pixel with fitting center in this paper.

Here, variables \( D_{1} \) and \( D_{2} \) are given by:

$$ \left\{ {\begin{array}{*{20}l} {D_{1} = I\left( {x,y} \right) \cdot \left| {I\left( {x,y} \right) - e_{2} } \right|} \hfill \\ {D_{2} = I\left( {x,y} \right) \cdot \left| {I\left( {x,y} \right) - e{}_{1}} \right|} \hfill \\ \end{array} } \right. $$
(10)

\( D_{1} \) is related to I(x,y) which is in the interior area of the evolving curve, and \( e_{2} \) is the fitting center of the curve exterior area, and it is same to \( D_{2} \). Thus, \( D_{1} \) and \( D_{2} \) are more robust than \( \lambda_{1} \) and \( \lambda_{2} \) which are expressed in the CV model. Here, the total energy function takes the following form:

$$ F(e_{1} ,e_{2} ,C) = \mu \int\limits_{\varOmega } H^{{\prime }} (\phi )\left| {\nabla \phi } \right|{\text{d}}x{\text{d}}y + v\int\limits_{\varOmega } {H(\phi ){\text{d}}x{\text{d}}y} - 2D_{1} \int\limits_{{\varOmega_{1} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - e_{1} } \right)^{d} } \right)H(\phi ){\text{d}}x{\text{d}}y} - \;2D_{2} \int\limits_{{\varOmega_{2} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - e_{2} } \right)^{d} } \right)(1 - H(\phi )){\text{d}}x{\text{d}}y} $$
(11)

In addition, the variables (\( D_{1} \) and \( D_{2} \)) accelerate the convergence of the model. While the grayscale differences of the internal area are greater than the external area, internal energy is larger than external energy, for controlling the curve evolution, and we should increase the weight of internal area and reduce the weight of external area. In Eq. (11), the weights of external area and internal area are \( - \;2D_{1} \) and \( - \;2D_{2} \) and then \( - \;2D_{1} > - \;2D_{2} \), and another situation can be equal. Adaptive energy weight of \( D_{1} \) and \( D_{2} \), which integrates every pixel with fitting center, controls the curve evolution and improves the accuracy of image segmentation and furthermore accelerates the convergence of the proposed model.

2.4 Active contour model based on distance regularized term

The calculation of logarithm function and local median information would increase the burden of computation. Distance regularized term avoids re-initialization of the level set function and accelerates the convergence of the proposed energy function. Here, we introduce distance regularized term into the objective function, and the final objective function can be expressed as:

$$ F^{r} \left( {e_{1} ,e_{2,} C} \right) = F\left( {e_{1} ,e_{2} ,C} \right) + \mu \int\limits_{\varOmega } {p\left( {\nabla \phi } \right)} \text{d}x $$
(12)

where \( p\left( s \right):\left[ {0, + \infty } \right) \to R \), a potential function with two minimum values at \( s = 0 \) and \( s = 1 \), and the potential function takes the following form:

$$ p(s) = \left\{ {\begin{array}{*{20}l} {\frac{1}{{(2\pi )^{2} }}\left( {1 - \cos 2\pi s} \right),\;{\text{if}}\; \le 1} \hfill \\ {\frac{1}{2}\left( {s - 1} \right)^{2} ,\;{\text{if}}\;s > 1} \hfill \\ \end{array} } \right., $$
(13)

While the potential function is used, it is possible to verify the boundedness of the diffusion rate. Then, the Gateaux derivative of \( F_{p} \left( \phi \right) = \mu \int\limits_{\varOmega } {p\left( {\nabla \phi } \right)} {\text{d}}x \) is given as:

$$ \frac{{\partial F_{p} \left( \phi \right)}}{\partial \phi } = - \;\mu {\text{div}}\left( {d_{p} \left( {\left| {\nabla \phi } \right|} \right)\nabla \phi } \right) $$
(14)

The key differences between the diffusion with potential functions [23] and \( p\left( s \right) \) are the boundedness of the corresponding diffusion rates and the diffusion behavior for the case \( |\nabla \phi | < 1/2 \). In summary, the total energy function is shown as:

$$ F^{r} (e_{1} ,e_{2} ,C) = v\int\limits_{\varOmega } {H(\phi ){\text{d}}x{\text{d}}y} + \mu \left[ \begin{aligned} \frac{1}{{\left( {2\pi } \right)^{2} }}\int\limits_{\varOmega } {\left[ {1 - \cos 2\pi \left| {\nabla \phi \left( x \right)} \right|} \right]} {\text{d}}xH\left( {1 - \left| {\nabla \phi \left( x \right)} \right|} \right) \hfill \\ + \frac{1}{2}\int\limits_{\varOmega } {\left( {\left| {\nabla \phi \left( x \right)} \right| - 1} \right)^{2} {\text{d}}xH\left( {\left| {\nabla \phi \left( x \right)} \right| - 1} \right)} \hfill \\ \end{aligned} \right] + 2D_{1} \int\limits_{{\varOmega_{1} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - e_{1} } \right)^{d} } \right)H(\phi ){\text{d}}x{\text{d}}y - } 2D_{2} \int\limits_{{\varOmega_{2} }} {\log \left( {1 + \left( {I\left( {x,y} \right) - e_{2} } \right)^{d} } \right)(1 - H(\phi )){\text{d}}x{\text{d}}y} $$
(15)

Minimizing the energy function in Eq. (15), the final partial differential equation is expressed as:

$$ \frac{\partial \phi }{\partial t} = H^{{\prime }} \left( \phi \right)\left\{ {\mu \left\{ {\begin{array}{*{20}l} {\text{div}\left\{ {\frac{{\text{sin}\left[ {2\uppi\left| {\nabla \phi \left( x \right)} \right|} \right]}}{{2\uppi\left| {\nabla \phi \left( x \right)} \right|}}\nabla \phi \left( x \right)} \right\}H\left[ {1 - \left| {\nabla \phi \left( x \right)} \right|} \right] + } \hfill \\ {\text{div}\left[ {\nabla \phi \left( x \right) - \frac{\nabla \phi \left( x \right)}{{\left| {\nabla \phi \left( x \right)} \right|}}} \right]H\left( {\left| {\nabla \phi \left( x \right)} \right| - 1} \right)} \hfill \\ \end{array} } \right\} - v{ + }2D_{1} \text{log}\left( {1 + \left( {I\left( {x,y} \right) - e_{1} } \right)^{d} } \right) - \;2D_{2} \text{log}\left( {1 + I\left( {x,y} \right) - e_{2} } \right)^{d} } \right\} $$
(16)

Using the following iterative calculation [4] to discretize Eq. (16), we will achieve the results of image segmentation:

$$ \phi^{{k^{'} + 1}} \left( x \right) = \phi^{{k^{'} }} \left( x \right) + \Delta t\frac{{\partial \phi^{{k^{'} }} \left( x \right)}}{\partial t} $$
(17)

where \( k^{'} \) and \( \Delta t \) are number of iterations and time. The implementation steps of the proposed model are:

  • Step 1 Initialize the level set function and contour.

  • Step 2 According to Eqs. (8) and (9), calculate \( e_{1} \) and \( e_{2} \).

  • Step 3 According to Eq. (10), calculate \( D_{1} \) and \( D_{2} \).

  • Step 4 Minimize Eq. (15), and according to Eq. (16), find the zero level set which is the evolution curve.

  • Step 5 If \( k \) is greater than the number of iterations, return the zero level set (contour); else, return to Step 2.

3 Experimental results and analysis

In order to verify the effectiveness of the proposed model, large numbers of SAR river images have been done in this section. The SAR river images are collected from RADARSAT-2 (C band) and TerraSAR-X (X band) satellites by random cropping, about 500 images. The example images of SAR database are illustrated in Fig. 1. What is more, comparing to the CV model [4], the local binary fitting model (LBF) [17], the local image fitting model (LIF) [18], the region-scalable fitting energy model (RSF) [19], the geodesic active contour (GAC) [8], the bias field correction level set model (BFC) [10], the distance regularized level set evolution (DRLSE) [21], active contour algorithm inspired from the cross-entropy (ACCE) [22]. The experiments presented in this section were all performed on Intel Core i5-4210U CPU 2.70Ghz/4.00 GB RAM. The algorithm was implemented in Matlab R2012a.

Fig. 1
figure 1

Example images of SAR river dataset

Full size image

We test the proposed model on the SAR dataset (ours), and about 500 images are used. Limited by the layout of the article, we will introduce five SAR river images to demonstrate the experimental results in this section. The size of these images are \( 147 \times 147 \), \( 166 \times 173 \), \( 137 \times 207 \), \( 178 \times 238 \), \( 215 \times 216 \), \( 214 \times 218 \). In the paper, the parameters are set to \( \lambda_{1} = \lambda_{2} = 10 \), \( \mu = 0.1 \times 255 \times 255 \), \( \nu = 0 \), \( \Delta t = 0.2 \), \( \varepsilon = 1 \). In the DRLSE energy model, set to \( \alpha = - 4.0 \). In the proposed model, set \( p = 5 \), and p is the index of logarithm kernel function. Then, Figs. 2, 3, 4, 5 and 6 show the segmentation results.

Fig. 2
figure 2

Image segmentation result of river image 1. a is the original image, b–j are the segmentation results of CV model, LBF model, LIF model, RSF model, GAC model, BFC model, DRLSE model, ACCE model and the proposed model

Full size image
Fig. 3
figure 3

Image segmentation result of river image 2. a is the original image, b–j are the segmentation results of CV model, LBF model, LIF model, RSF model, GAC model, BFC model, DRLSE model, ACCE model and the proposed model

Full size image
Fig. 4
figure 4

Image segmentation result of river image 3. a is the original image, b–j are the segmentation results of CV model, LBF model, LIF model, RSF model, GAC model, BFC model, DRLSE model, ACCE model and the proposed model

Full size image
Fig. 5
figure 5

Image segmentation result of river image 4. a is the original image, b–j are the segmentation results of CV model, LBF model, LIF model, RSF model, GAC model, BFC model, DRLSE model, ACCE model and the proposed model

Full size image
Fig. 6
figure 6

Image segmentation result of river image 5. a is the original image, b–j are the segmentation results of CV model, LBF model, LIF model, RSF model, GAC model, BFC model, DRLSE model, ACCE model and the proposed model

Full size image

With the purpose of evaluating the segmentation results quantitatively, DSC (Dice Similarity Coefficient) standard, an area-based method, is used for the analysis.

$$ {\text{DSC}}\left( {S_{S} ,S_{R} } \right) = \frac{{2 \cdot {\text{Area}}\left( {S_{S} \cap S_{R} } \right)}}{{{\text{Area}}\left( {S_{S} } \right) + {\text{Area}}\left( {S_{R} } \right)}} $$

where \( S_{S} \) represents the segmentation results of every model and \( S_{R} \) represents the standard segmentation result which is made by manual work (ground truth). While the value of DSC is closer to 1, the segmentation result is better. The ground truth of SAR river image is illustrated in Fig. 7, and the DSC of every segmentation result is illustrated in Table 1.

Fig. 7
figure 7

Ground truth of SAR river image

Full size image
Table 1 DSC of every segmentation result
Full size table

From the experimental segmentation results, observe that the traditional CV model cannot segment river image accurately, and it is that CV model is based on global gray value and ignores local information of the whole image and the grayscale difference; LBF, LIF and RSF models include non-convex level set function, which leads to be sensitive to the initialization of level set function, and sometimes achieves inaccurate results. We select the initial contour three times randomly and choose the best segmentation result to be illustrated in this paper; GAC model is edge-based active contour model, and the model utilizes boundary stop function which is based on gradient information to guide the contour to move toward the boundary, and SAR river image presents intensity inhomogeneity. Therefore, Figs. 2, 3, 4, 5 and 6f are worse than Figs. 2, 3, 4, 5 and 6b. BFC model is a novel level set method in the presence of intensity inhomogeneity, comparing Fig. 2g with Fig. 2b–f, and the BFC model segments better. However, while the target or background has large gray scale fluctuations, such as Figs. 3g and 4g, it fails. The distance regularization term eliminates the need for re-initialization and then avoids inducing numerical errors, and consequently, segmentation results are similar to traditional CV model. ACCE model aims to achieve the segmentation results where the cross-entropy between the original image and the segmented image is minimum, and therefore, segmentation results are better than traditional CV model.

Segmentation results shown in Figs. 2, 3, 4, 5 and 6j achieve the best performance, and the proposed model, which is defined as a logarithm kernel function distance, suppresses the noise which would have influences on the curve evolution; Furthermore, the fitting center (different from other fitting centers) is based on local median information of the whole image, and consequently, it is more robust than the CV model. In particular, for SAR river image, the proposed model achieves competitive performance superior to other state-of-the-art models. Moreover, Fig. 8 gives the segmentation results of synthetic image and real-world river image to validate the absolute accuracy of the proposed model. The proposed model, regional energy function that depends on the logarithm kernel distance, the local median information of image, the interior and exterior energies integrate each pixel with fitting center, and all of them have good performances for the curve evolution and the characterization of the energy inside and outside the curve. Therefore, the proposed model achieves good segmentation results on synthetic and real-world river image.

Fig. 8
figure 8

Segmentation result of synthetic image and real image, a synthetic images, c real-world river image, b and d segmentation results

Full size image

The running time of nine methods are listed in Table 2. In this section, the iteration is set to 500 to make sure all the images are segmented successfully. In Table 2, the running time of RSF LBF and LIF models is long, and this is because these models are based on regional or gradient information. The calculation of bias field is introduced to cause that the running time of BFC model is longer than CV model. DRLSE and ACCE models avoid re-initialization of the level set function, and thus, the running time is short. The proposed model is slower than the CV model, and it is that the proposed model includes the computation of logarithm kernel distance, local median information and adaptive weight; however, on the other hand, distance regularized term avoids re-initialization of the level set function, and moreover, adaptive weight accelerates the convergence of the model.

Table 2 Running time of 9 methods (unit: s)
Full size table

4 Conclusions and future work

A region-based active contour model driven by logarithm kernel-based energy function is proposed in this paper. Different from other level set models, the transformed space of the kernel distance can be linearly classified; the fitting center which merges local median information and global mean information reduces the influence of the homogeneous gray inside the region. Adaptive energy weight integrates every pixel with fitting center controls the curve evolution and improves the performance of segmentation results; finally, the model introduces distance regularized term to avoid the re-initialization of level set function.

There are some drawbacks in the proposed model which need to be improved. The convergent condition cannot be set adaptive in the proposed model, and the curve converges slowly in some SAR river image. We should try to solve these problems in the future work.