# Pyramidal Watershed Segmentation Algorithm for High-Resolution Remote Sensing Images Using Discrete Wavelet Transforms

Authors: K. Parvathi; B. S. Prakasa Rao; M. Mariya Das; T. V. Rao
Journal: Discrete Dynamics in Nature and Society
 (2009)

---

## Abstract

The watershed transformation is a useful morphological segmentation tool for a variety of grey-scale images. However, over segmentation and under segmentation have become the key problems for the conventional algorithm. In this paper, an efficient segmentation method for high-resolution remote sensing image analysis is presented. Wavelet analysis is one of the most popular techniques that can be used to detect local intensity variation and hence the wavelet transformation is used to analyze the image. Wavelet transform is applied to the image, producing detail (horizontal, vertical, and diagonal) and Approximation coefficients. The image gradient with selective regional minima is estimated with the grey-scale morphology for the Approximation image at a suitable resolution, and then the watershed is applied to the gradient image to avoid over segmentation. The segmented image is projected up to high resolutions using the inverse wavelet transform. The watershed segmentation is applied to small subset size image, demanding less computational time. We have applied our new approach to analyze remote sensing images. The algorithm was implemented in MATLAB. Experimental results demonstrated the method to be effective.


---

## Body

## 1. Introduction

Image  segmentation is object oriented and hence useful in high-resolution image analysis. It  provides a partitioning of the image into isolated regions, each one representing a different image. In particular, multiscale methods are based on image transformations that reduce image resolution, and they are able to segment objects of different sizes depending on the chosen resolution. From high-resolution images, small details can be detected, while from low-resolution images, larger structures are detected. The multi-scale behavior of image features has been analyzed in different ways. Tracking of intensity extrema along scales defined by Lifshitz [1] and Lindeberg  [2] was used for image segmentation. Other approach for image segmentation is based on the watershed transform. This transform can be applied to the gradient magnitude image as defined by Meyer and Beucher [3], Vincent and Soille [4], and Bieniek and Moga [5] to obtain the segmented regions. Small fluctuations in the grey levels produce spurious gradients, which cause over segmentation. To overcome this problem many techniques based on watersheds have been proposed. Meyer [6] introduced the leveling approach, which applies morphological filters to reduce the small details in the image. Jung and Scharcanski [7], Jung  [8] and Haris et al. [9] proposed edge preserving statistical noise reduction approach as a preprocessing for the watershed transform, and a hierarchical merging process as a postprocessing stage. Jackway [10] and Gauch [11] used morphological scale spaces for image segmentation. Weickert [12] proposed partial differential equations for image denoising and edge enhancement, combined with watershed segmentation and region merging. J. B. Kim and H. J. Kim [13, 14] proposed a wavelet-based watershed segmentation technique, by projecting the segmented image into higher resolutions.  Nguyen et al. [15] proposed a combination of energy-based segmentation and watersheds, called water snakes. Recently Liu et al.  [16] modified the wavelet-based watershed segmentation algorithm, introduced the morphological filters to smooth image while filtering out noise as preprocessing of medical image segmentation. Almost all the proposed techniques are applied on the medical and nonmedical images.


## 2. Methodology

In this paper, we proposed a new segmentation algorithm for high-resolution remote sensing images, which can also be applied to medical and nonmedical images. We used a bi-orthogonal (bior) wavelet decomposition to describe a remote sensing image in multiple resolutions. A suitable resolution is chosen. The gradient image is estimated (or computed) by the simple grey scale morphology. To avoid over segmentation, we have imposed the selective minima (regional minima of the image) on the gradient image. The watershed transform is applied and the segmentation result is projected to a higher resolution, using the inverse wavelet transform until the full resolution of segmented image is obtained. A general outline of the proposed method   is shown in Figure1.Figure 1
Flow chart of algorithm.


### 2.1. The Wavelet Transform

The wavelet transform is a mathematical tool that can be used to describe images in multiple resolutions. The wavelet transform can be represented with an equation like Fourier transform.(2.1)F(w)=∫∞∞f(t)ejwtdt,C(scale,position)=∫∞∞s(t)ψ(scale,position,t)dt.
The  results of wavelet transform are many wavelet coefficients C, which are functions of scale and position,s(t) is the image function and ψ(scale,position,t) is the wavelet function. If scales and positions based on powers of two, the so-called dyadic scales and positions then our analysis will be more efficient and just as accurate. We can obtain such analysis from the Discrete Wavelet Transform (DWT). For many images, the low frequency content is the most important part. This is the one which gives the image identity. The high frequency content imparts flavor or nuance.According to Mallat’s pyramidal algorithm [17], the original image is convolved with low-pass and high-pass filters associated with a mother wavelet, and down sampled afterwards as shown in Figure 2. In order to explain the process of study, a satellite image of Visakhapatnam city coastal harbor area has been chosen.  It is PAN sensor with 5.8 m spatial resolution of IRS-1C satellite and it is as shown in Figure 3. Four images (each with half the size of the original image) are produced. The first subimage corresponding to low frequencies in both the directions (LL), is called the Approximation image (cA1) of the original image as shown in Figure 4(a). Second one corresponds to low frequencies in horizontal and high frequencies in the vertical direction (LH), is called the horizontal coefficients (cH1) of original image (Figure 4(b)). The high frequencies in horizontal and low frequencies in the vertical direction (HL), are called vertical coefficients (cV1) as shown in Figure 4(c) where as high frequencies in both the directions (HH), are called diagonal coefficients (cD1) as shown in Figure 4(d). The process is repeated on the LL subband (cA1) to generate the next level of the decomposition adopting the same process. Second level of decomposed images are shown in Figure 5. The Approximation coefficient image (cA2) is shown in Figure 5(a). The horizontal, vertical and diagonal coefficients of second level are also shown in Figures 5(b), 5(c), and 5(d) respectively. The original image, S is represented in terms of its wavelet coefficients as given byFigure 2
Discrete wavelet transform of an image (single level). L—Low Pass Filter; H—High Pass Filter. LL—Approximation coefficients (High scale, low frequency)—cA1. LH, HL and HH– Detail coefficients (low scale, high frequency)—cH1, cV1  and  cD1.Figure 3
Original image size (264×414×3).First level decomposition of Image: (a) approximation coefficients, (b) horizontal coefficients, (c) vertical coefficients, (d) diagonal coefficients.
(a)
cA1 (132×207)(b)
cH1(132×207)(c)
cV1 (132×207)(d)
cD1 (132×207)Second level decomposition of Image (cA1) : (a) approximation coefficients, (b) horizontal cients, (c) vertical coefficients, (d) diagonal coefficients.
(a)
cA2(66×104)(b)
cH2(66×104)(c)
cV2(66×104)(d)
cD2  (66×104)(2.2)S=cA1+{cH1+cV1+cD1}.
Iterating the process to J levels.(2.3)S=cAJ+{cHi+cVi+cDi}i=1J.
Since the analysis process is iterative in theory, it can be continued indefinitely. In reality the decomposition can be continued only until the individual details consist of a single sample or pixel.The Approximations are related to one another by(2.4)cAJ-1=cAJ+(cHJ+cVJ+cDJ).


### 2.2. Pyramidal Watershed Segmentation

The original image is decomposed into different resolutions using wavelet transform explained as above. In this paper resolution of (66×104 pixels) is chosen (to avoid under segmentation due to loss of image details in further decomposition process.). The image is computed up to the high-resolution (264×144). The Approximation image cA2 (66×104) is used as a starting point as shown in Figure 5(a). The Sobel operator (Sobel horizontal edge emphasizing filter) is used to estimate the image gradient. The image differences in the horizontal and vertical directions are obtained using multidimensional image filters applied separately in both the directions. cH2 (66×104) and cV2 (66×104) denote these two images as shown in Figures 5(b) and 5(c).  Then the gradient magnitude can be computed by(2.5)G2=(cH22+cV22)The direct application of watershed transforms causes over segmentation; the over segmented images may require further merging of some regions. The decision on such regions to be merged depends upon homogeneity and similarity criteria based on the extended minima transform defined by Soille [18], the general syntax is(2.6)BW=Imextendedmin[G,H].
Equation (2.6) computes the extended minima transform, which is the regional minima of the h-minima transform. G is the gradient image and H is a selective minima value (as mentioned in Table 1) depends on the required number of segmented regions and sizes of the desired objects. We detect all intensity valleys below a particular threshold with the imextendedmin function. The output of the imextendedmin function is a binary image. The location rather than the size of the regions in the imextendedmin image is important at each pyramid level. The obtained binary image was imposed on original image to contain only those valleys found by the imextendedmin function. All regions containing imposed minima will be detected by the watershed transform. The standard watershed transform by Vincent and Soille [19] is applied to get well-segmented image at low-resolution as shown in the Figure 6(b).Table 1
The results of image segmentation on each resolution.
Test imageScale levelsSelective minima value, HNumber of segmentsComputational time (sec)Satellite image 1cA0951880.55cA11000300.23cA2300070.15Satellite image 2cA0602150.75cA11200350.25cA2390040.10Brain imagecA0120620.45cA1500220.3cA2350020.10Flowers imagecA0130690.40cA180080.25cA2300030.10Peppers imagecA060380.35cA1120050.20cA2270030.15Final segmented Images (a) original image and (b) to (d) are the watershed imposed original images from level 2 to level 0 resolution respectively.
(a)(b)(c)(d)


### 2.3. Image Projection

Once the segmented image is generated at theJ layer, it must be projected into next layer (J-1) in order to obtain the full resolution image segmentation. Direct projection of the segmented image offers very poor results. To overcome this problem, we have used an Inverse Discrete Wavelet Transform (IDWT) to implement projection from low to high-resolution layer. The IDWT is the summation of the products of the decomposed coefficients and wavelet function at respective scales.(2.7)S(t)=∑j∈J∑k∈JC(j,k)ψj,k(t),
where S(t) is the Original Image,  j and k  represent Scale,Position respectively.  C(j,k) Coefficients function and ψj,k(t) wavelet functions as used in DWT.To prevent noise from being introduced back into the unsampled image during Inverse transform, we have used only some detail coefficients (those related to edges). This means that all detail images (horizontal, vertical, diagonal) are set to zero, except those whose position correspond to the watershed lines of imageS. We have applied the Inverse wavelet transform on segmented Approximation image Figure 6(b), to get next layer image. The procedure of imposing the selective minima on the gradient image followed by application of watershed transform and image projection to the next layer has been repeated till we obtained the full resolution segmented original image at high-resolution (level 0) as explained in the flow chart (Figure 1). The final segmented images from low-resolution (level 2) to high-resolution (level 0) are shown in Figure 6.


## 2.1. The Wavelet Transform

The wavelet transform is a mathematical tool that can be used to describe images in multiple resolutions. The wavelet transform can be represented with an equation like Fourier transform.(2.1)F(w)=∫∞∞f(t)ejwtdt,C(scale,position)=∫∞∞s(t)ψ(scale,position,t)dt.
The  results of wavelet transform are many wavelet coefficients C, which are functions of scale and position,s(t) is the image function and ψ(scale,position,t) is the wavelet function. If scales and positions based on powers of two, the so-called dyadic scales and positions then our analysis will be more efficient and just as accurate. We can obtain such analysis from the Discrete Wavelet Transform (DWT). For many images, the low frequency content is the most important part. This is the one which gives the image identity. The high frequency content imparts flavor or nuance.According to Mallat’s pyramidal algorithm [17], the original image is convolved with low-pass and high-pass filters associated with a mother wavelet, and down sampled afterwards as shown in Figure 2. In order to explain the process of study, a satellite image of Visakhapatnam city coastal harbor area has been chosen.  It is PAN sensor with 5.8 m spatial resolution of IRS-1C satellite and it is as shown in Figure 3. Four images (each with half the size of the original image) are produced. The first subimage corresponding to low frequencies in both the directions (LL), is called the Approximation image (cA1) of the original image as shown in Figure 4(a). Second one corresponds to low frequencies in horizontal and high frequencies in the vertical direction (LH), is called the horizontal coefficients (cH1) of original image (Figure 4(b)). The high frequencies in horizontal and low frequencies in the vertical direction (HL), are called vertical coefficients (cV1) as shown in Figure 4(c) where as high frequencies in both the directions (HH), are called diagonal coefficients (cD1) as shown in Figure 4(d). The process is repeated on the LL subband (cA1) to generate the next level of the decomposition adopting the same process. Second level of decomposed images are shown in Figure 5. The Approximation coefficient image (cA2) is shown in Figure 5(a). The horizontal, vertical and diagonal coefficients of second level are also shown in Figures 5(b), 5(c), and 5(d) respectively. The original image, S is represented in terms of its wavelet coefficients as given byFigure 2
Discrete wavelet transform of an image (single level). L—Low Pass Filter; H—High Pass Filter. LL—Approximation coefficients (High scale, low frequency)—cA1. LH, HL and HH– Detail coefficients (low scale, high frequency)—cH1, cV1  and  cD1.Figure 3
Original image size (264×414×3).First level decomposition of Image: (a) approximation coefficients, (b) horizontal coefficients, (c) vertical coefficients, (d) diagonal coefficients.
(a)
cA1 (132×207)(b)
cH1(132×207)(c)
cV1 (132×207)(d)
cD1 (132×207)Second level decomposition of Image (cA1) : (a) approximation coefficients, (b) horizontal cients, (c) vertical coefficients, (d) diagonal coefficients.
(a)
cA2(66×104)(b)
cH2(66×104)(c)
cV2(66×104)(d)
cD2  (66×104)(2.2)S=cA1+{cH1+cV1+cD1}.
Iterating the process to J levels.(2.3)S=cAJ+{cHi+cVi+cDi}i=1J.
Since the analysis process is iterative in theory, it can be continued indefinitely. In reality the decomposition can be continued only until the individual details consist of a single sample or pixel.The Approximations are related to one another by(2.4)cAJ-1=cAJ+(cHJ+cVJ+cDJ).


## 2.2. Pyramidal Watershed Segmentation

The original image is decomposed into different resolutions using wavelet transform explained as above. In this paper resolution of (66×104 pixels) is chosen (to avoid under segmentation due to loss of image details in further decomposition process.). The image is computed up to the high-resolution (264×144). The Approximation image cA2 (66×104) is used as a starting point as shown in Figure 5(a). The Sobel operator (Sobel horizontal edge emphasizing filter) is used to estimate the image gradient. The image differences in the horizontal and vertical directions are obtained using multidimensional image filters applied separately in both the directions. cH2 (66×104) and cV2 (66×104) denote these two images as shown in Figures 5(b) and 5(c).  Then the gradient magnitude can be computed by(2.5)G2=(cH22+cV22)The direct application of watershed transforms causes over segmentation; the over segmented images may require further merging of some regions. The decision on such regions to be merged depends upon homogeneity and similarity criteria based on the extended minima transform defined by Soille [18], the general syntax is(2.6)BW=Imextendedmin[G,H].
Equation (2.6) computes the extended minima transform, which is the regional minima of the h-minima transform. G is the gradient image and H is a selective minima value (as mentioned in Table 1) depends on the required number of segmented regions and sizes of the desired objects. We detect all intensity valleys below a particular threshold with the imextendedmin function. The output of the imextendedmin function is a binary image. The location rather than the size of the regions in the imextendedmin image is important at each pyramid level. The obtained binary image was imposed on original image to contain only those valleys found by the imextendedmin function. All regions containing imposed minima will be detected by the watershed transform. The standard watershed transform by Vincent and Soille [19] is applied to get well-segmented image at low-resolution as shown in the Figure 6(b).Table 1
The results of image segmentation on each resolution.
Test imageScale levelsSelective minima value, HNumber of segmentsComputational time (sec)Satellite image 1cA0951880.55cA11000300.23cA2300070.15Satellite image 2cA0602150.75cA11200350.25cA2390040.10Brain imagecA0120620.45cA1500220.3cA2350020.10Flowers imagecA0130690.40cA180080.25cA2300030.10Peppers imagecA060380.35cA1120050.20cA2270030.15Final segmented Images (a) original image and (b) to (d) are the watershed imposed original images from level 2 to level 0 resolution respectively.
(a)(b)(c)(d)


## 2.3. Image Projection

Once the segmented image is generated at theJ layer, it must be projected into next layer (J-1) in order to obtain the full resolution image segmentation. Direct projection of the segmented image offers very poor results. To overcome this problem, we have used an Inverse Discrete Wavelet Transform (IDWT) to implement projection from low to high-resolution layer. The IDWT is the summation of the products of the decomposed coefficients and wavelet function at respective scales.(2.7)S(t)=∑j∈J∑k∈JC(j,k)ψj,k(t),
where S(t) is the Original Image,  j and k  represent Scale,Position respectively.  C(j,k) Coefficients function and ψj,k(t) wavelet functions as used in DWT.To prevent noise from being introduced back into the unsampled image during Inverse transform, we have used only some detail coefficients (those related to edges). This means that all detail images (horizontal, vertical, diagonal) are set to zero, except those whose position correspond to the watershed lines of imageS. We have applied the Inverse wavelet transform on segmented Approximation image Figure 6(b), to get next layer image. The procedure of imposing the selective minima on the gradient image followed by application of watershed transform and image projection to the next layer has been repeated till we obtained the full resolution segmented original image at high-resolution (level 0) as explained in the flow chart (Figure 1). The final segmented images from low-resolution (level 2) to high-resolution (level 0) are shown in Figure 6.


## 3. Experimental Results

All the computation involving wavelet transform, edge detection and watershed transform are implemented using MATLAB. Application of the wavelet transform takes very short time, this quick response is mainly due to the property of wavelet decomposition and reconstruction which have fast algorithms. The algorithms are based on convolutions with a bank of filters. We used two different satellite images. Satellite image 1 is a PAN image of IRC-1C (5.8 m resolution). It shows Visakapatnam sea coast with fishing jetties, harbor channels and dense city with concrete structures. Image 2 is a Cartosat-1 satellite PAN image (2.5 m resolution) of Hyderabad city. The final segmentation results are shown in Figure6. The Figure 6(a) are the original images and Figures 6(b) to 6(d) are the watershed imposed original images from level 2 to level 0, respectively. At the lower resolution (level 2), we selected the minima to detect minimum number of intensity valleys, which gives less number of segmented regions. The final high-resolution image is having more number of segmented regions compared to lower resolutions.   At level 2 (Figure 6(b)),  segmentation is seen on broad features like sea and the coastal features in the left image while dense city structures are separated form isolated buildings and field areas from the right side satellite image. At level1 all the regions with different gray level intensities of major structure and surroundings are clearly segmented (Figure 6(c)), while in the Hyderabad city area park with vegetation, buildings and building shadows and buildings with utility space are all segmented. At level 0 all the regions with different gray level intensities are segmented very clearly (Figure 6(d)) between the various structures, vegetation, parks and orchards in both the images.We also applied our algorithm on different medical and nonmedical images. Figure7 shows the results of image segmentation using the above approach. The first column shows the original image, second, third and fourth column show the watershed imposed images of the wavelet decomposition from level 2 to level 0 respectively. The results show not only the elimination of over segmentation but also reducing the number of segments. In MRI medical image, all the regions with different intensities are clearly segmented. We can easily extract desired parts like tumor in the brain image. In flowers image, the borders (edges) of all flowers are traced and petals segmented with respect to their intensity variation. Even in the peppers image, all the spices are clearly segmented.Final segmented images on wavelet decomposition from level 2 to level 0. (a) Original image and, (b) to (d) are the watershed imposed original images from low-resolution (64×64) to high-resolution (256×256) respectively.
(a)(b)(c)(d)


### 3.1. Performance Evaluation

To evaluate the performance of the presented method, simulation was carried out on the sequence of frames of different images. The pyramid image generated by the 2-scale Bior (Bi-orthogonal) wavelet transformation and selective minima imposed on gradient image for watershed segmentation. We evaluated the results of the segmentation of the present method by using three common objective measurements viz., the selective minima, number of segments and computation time at each scale. The full resolution image has more number of regions, longer computation time and lowest minima value where as the low-resolution image has less number of regions, shorter computational time and higher minima value. Each set of the results of image segmentation on each resolution is shown in Table1. In the present method the use of a low-resolution image is highly preferable as a starting point of pyramidal segmentation since it gives the best objective quality, when the number of regions and computational time are taken into consideration.


## 3.1. Performance Evaluation

To evaluate the performance of the presented method, simulation was carried out on the sequence of frames of different images. The pyramid image generated by the 2-scale Bior (Bi-orthogonal) wavelet transformation and selective minima imposed on gradient image for watershed segmentation. We evaluated the results of the segmentation of the present method by using three common objective measurements viz., the selective minima, number of segments and computation time at each scale. The full resolution image has more number of regions, longer computation time and lowest minima value where as the low-resolution image has less number of regions, shorter computational time and higher minima value. Each set of the results of image segmentation on each resolution is shown in Table1. In the present method the use of a low-resolution image is highly preferable as a starting point of pyramidal segmentation since it gives the best objective quality, when the number of regions and computational time are taken into consideration.


## 4. Conclusions

In this paper, we proposed a new technique to improve the watershed segmentation on remote sensing images based on Discrete Wavelet Transforms. The lower resolution must be chosen in accordance with the size and number of desired objects to avoid under segmentation. The experimental results indicate that the proposed technique performs well for the remote sensing images as well as for medical and low-resolution images too. The full resolution image has more number of regions, longer computation time and lowest minima value where as the low-resolution image has less number of regions, shorter computation time and higher minima value. We would like to extend this work using other soft computing techniques.



---

*Source: 601638-2009-09-23.xml*