Pyramidal Watershed Segmentation Algorithm for High-Resolution Remote Sensing Images Using Discrete Wavelet Transforms
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.
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  and Lindeberg  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 , Vincent and Soille , and Bieniek and Moga  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  introduced the leveling approach, which applies morphological filters to reduce the small details in the image. Jung and Scharcanski , Jung  and Haris et al.  proposed edge preserving statistical noise reduction approach as a preprocessing for the watershed transform, and a hierarchical merging process as a postprocessing stage. Jackway  and Gauch  used morphological scale spaces for image segmentation. Weickert  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.  proposed a combination of energy-based segmentation and watersheds, called water snakes. Recently Liu et al.  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.
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 Figure 1.
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.
The results of wavelet transform are many wavelet coefficients , which are functions of scale and position, is the image function and 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 , 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, is represented in terms of its wavelet coefficients as given by
|(a) cA1 ()|
|(b) cH1 ()|
|(c) cV1 ()|
|(d) cD1 ()|
|(a) cA2 ()|
|(b) cH2 ()|
|(c) cV2 ()|
|(d) cD2 ()|
Iterating the process to levels.
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.2. Pyramidal Watershed Segmentation
The original image is decomposed into different resolutions using wavelet transform explained as above. In this paper resolution of ( 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 (). The Approximation image cA2 () 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 () and cV2 () denote these two images as shown in Figures 5(b) and 5(c). Then the gradient magnitude can be computed by
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 , the general syntax is
Equation (2.6) computes the extended minima transform, which is the regional minima of the h-minima transform. is the gradient image and 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  is applied to get well-segmented image at low-resolution as shown in the Figure 6(b).
2.3. Image Projection
Once the segmented image is generated at the layer, it must be projected into next layer 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.
where is the Original Image, and represent Scale,Position respectively. Coefficients function and 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 image . 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 Figure 6. 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. Figure 7 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.
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 Table 1. 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.
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.
The authors thank the All India Council of Technical Education for utilizing the data obtained in the research project for this paper. The authors also thank Andhra Pradesh State Remote Sensing Application Center for providing the high-resolution satellite data. The medical and nonmedical images are downloaded from Google.com/images.
L. M. Lifshitz, Image segmentation using global knowledge and a priori information, Ph.D. thesis, University of North Carolina, Chapel Hill, NC, USA, 1987.
T. Lindeberg, Discrete scale space theory and the scale space primal sketch, Ph.D. thesis, Royal Institute of Technology, Stockholm, Sweden, May 1991.
F. Meyer and S. Beucher, “Morphological segmentation,” Journal of Visual Communication and Image Representation, vol. 1, no. 1, pp. 21–46, 1990.View at: Google Scholar
F. Meyer, “Leveling and morphological segmentation,” in Proceedings of International Symposium on Computer Graphics, Image Processing, and Vision (SIBAGRAPI '98), pp. 28–35, Rio de Janeiro, Brazil, October 1998.View at: Google Scholar
C. R. Jung and J. Scharcanski, “Robust Watershed segmentation using the wavelet transforms,” in Proceedings of the XV Brazilian Symposium on Computer Graphics and Image Processing (SIBGRAPI '02), pp. 1530–1834, February 2000.View at: Google Scholar
C. R. Jung, “Multiscale image segmentation using wavelets and watersheds,” in Proceedings of the XVI Brazilian Symposium on Computer Graphics and Image Processing (SIBGRAPI '02), pp. 1530–1834, March 2003.View at: Google Scholar
K. Haris, S. N. Efstratiadis, N. Maglaveras, and A. K. Katsaggelos, “Hybrid image segmentation using watersheds and fast region merging,” IEEE Transactions on Image Processing, vol. 7, no. 12, pp. 1684–1699, 1998.View at: Google Scholar
P. Jackway, “Morphological multiscale gradient watershed image analysis,” in Proceedings of the 9th Scandivania Conference on Image Analysis, pp. 87–94, 1995.View at: Google Scholar
J. M. Gauch, “Image segmentation and analysis via multiscale gradient watershed hierarchies,” IEEE Transactions on Image Processing, vol. 8, no. 1, pp. 69–79, 1999.View at: Google Scholar
J. B. Kim and H. J. Kim, “A wavelet-based watershed image segmentation for VOP generation,” in Proceedings of International Conference on Pattern Recognition, vol. 16, no. 3, pp. 505–508, Queboc city, Canada, August 2002.View at: Google Scholar
H. Liu, Z. Chen, X. Chen, and Y. Chen, “Multiresolution medical image segmentation based on wavelet transform,” in Proceedings of the 27th Annual International Conference of the Engineering in Medicine and Biology Society (EMBS '05), vol. 7, pp. 3418–3421, Shanghai, China, September 2005.View at: Google Scholar
P. Soille, Morphological Image Analysis Principles & Applications, Springer, Berlin, Germany, 1999.View at: MathSciNet