Semi-automatic multi-level approach for extraction of tidal channels from Aerial Photographs and hyperspectral data Bharat Lohani and Hemendra Singh Bist Department of Civil Engineering, IIT Kanpur Kanpur 208016 (India) E-mail:blohani@iitk.ac.in Tel: +91-512-2597623 David C. Mason Environmental Systems Science Centre The University of Reading Reading RG6 6AL (UK)
Abstract
Tidal channels play a critical role in morphodynamics of tidal environment. In order to incorporate tidal channels for understanding this dynamics and for management of tidal zones it is essential to map these channels accurately at high speed. Mapping of tidal channels using land surveying is a cumbersome, time-consuming, and error-prone task. The remotely sensed data provide a suitable alternative for channel mapping. In a number of studies, employing remotely sensed data, the channels are mapped using manual interpretation and tracing. However, the manual approach suffers with several limitations including subjectivity in the outcome and inaccuracy. Furthermore, the manual approach is not suitable if the channel mapping is to be carried out for a large area or involving several images spread in time. In view of the above, development of an automatic procedure for tidal channel extraction from remotely sensed data is much warranted. The aim of this paper is to generate a method which will work for aerial photographs, MSS and hyperspectral data. The results in this paper, however, are shown using an aerial photograph. The technique developed is a semi-automatic multi-level approach, which seeks to classify the image in a feature constrained space. The low level processing employs edge detection followed by hysteresis thresholding to identify edge locations. The centrelines of all features (bounded by edges) are identified using distance with destination transform. The user is prompted to identify a few channel pixels, using which channel spectral clusters are formed by ISODATA classification. The Mahalanobis distance is used to classify entire image using these clusters. At this stage, beginning from the seed pixels the channel centreline segments are located and connected to other similar components. The final step includes fattening the located channel centrelines to generate full channels. This paper presents the results obtained till this stage and will outline the further work proposed to make the procedure complete. 1. Introduction The study of tidal channels plays a critical role in understanding the tidal morphodynamics of intertidal zone. The tidal channel morphology is still not well understood and differing schools of thought exit to explain their development (French and Stoddart, 1992; Pye, 1992; Steel, 1996; Pethick, 1992; Rodriguez-Iturbe and Rinaldo, 1997; Rinaldo et al., 1999). It is important to understand tidal channels and their evolution for various applications involving coastal management, salt marsh reclamation, artificial salt marsh development, and for managing lagoon cities like Venice (Fagherazzi et al., 1999). Physical measurement of tidal channels, i.e. mapping of their planimetric extent and cross-section, is essential for their study. Conventionally, these measurements are obtained using land surveying. Despite availability of advanced equipments, land surveying fails to meet the requirement when accurate and dense measurements are needed within a short time. The satellite and airborne data have potential to be used for this purpose. These data are also capable of providing the historical measurements for channels. In most of the geomorphological studies employing these data the channel networks are extracted using manual tracing. This approach is accurate and appropriate if the area of study includes only a few images spread in time and space. However, with the availability of a large number of these images, as are being produced by various air- and space-borne sensors, and also the need to include these in morphological studies spanning wider in space and time, the manual methods of channel mapping become too cumbersome. Furthermore, as at every step of manual tracing the input of an operator is involved, the mapped channels also suffer from subjectivity. Considering the above, there have been attempts to develop computer based methods to extract channel information from remotely sensed and other spatial data. Lohani (1999) provides a detailed review of algorithms for extracting channels (mostly terrestrial) from low resolution Digital Elevation Models (DEMs). The tidal channels, in view of their different morphology than their terrestrial counterparts, require different treatment and the algorithms designed for terrestrial channels fail to extract tidal channels (Lohani, 1999). Furthermore, in view of small sizes of tidal channels (up to a width of 30 cm) there is a need to use high resolution data. The recently emerged technique of airborne altimetric LiDAR provides high resolution DEM and have been employed by Lohani and Mason (2001) and Mason et al. (2003) for tidal channel network extraction. With the geometric information contained in them, LiDAR data should be better for channel extraction (Mason et al, 2003). However, there are circumstances in which only spectral data are available, e.g. mapping channels using historical images and cases when only spectral sensor is flown. Also, it is common now that with most of the LiDAR flights a spectral sensor is also flown for strengthening the data captured which increases the accuracy of information extracted. In view of the significant role that spectral data can play and their availability this papers attempts to use these for extraction of tidal channels. 2. Data used The data used in this paper is an aerial photograph (Figure 1 (a)) of the tidal zone of Wash in the UK. The image consists of channels which can be interpreted manually owing to their 1) spectral signatures, 2) typical shape, and 3) connectivity. The aerial photograph is separated in red, green and blue bands. It can be observed from the image that despite the manual interpretation being easy, computer processing for extraction of channels is not a straightforward process. This is basically due to the high variability in characteristics (spectral signature, shape, connectivity etc.) of channels which makes it difficult to extract channels using a direct low level image processing technique. ![]() Figure 1: (a) Original image; (b) Non-maxima suppressed edge image; (c) Distance transform image 3. Methodology The methodology consists of several steps run in parallel and sequence. These steps are discussed below and shown in the flow diagram (Figure 2). These steps are coded in C under Sun Solaris environment. 3.1 Edge detection and thresholding For the results presented in this paper, Sobel edge operator was employed on the Red band of aerial photograph. However, investigation is on to use the colour edge detector, which would make use of channel information stored in all bands to detect maximum edges from a scene. An edge threshold THRESH is calculated above which 20% of edge strength lies. At this stage to reduce the thickness of edges, those edges which are non-maxima in the direction perpendicular to the edge direction are suppressed. Hysteresis thresholding (Canny, 1986) is then used with the upper threshold being THRESH and the lower threshold being THRESH/2. The resulting binary image is shown in Figure 1 (b). While performing Sobel operator an edge direction image is also generated, which will be used in latter steps. ![]() Figure 2: Flow chart of methodology adopted 3.2 Distance-with-destination transform Using the thresholded image of step 3.1 a distance-with-destination transform is generated, which is a modified form of Euclidean distance transform (Mason et al. 2003). In the Euclidean distance transform (Castleman, 1996) each TURE (edge) pixel assumes a value of zero while for all other FALSE (non-edge) pixels the values given are the distance of the FALSE pixel from the nearest TRUE pixel. In distance-with-destination transform, in addition to above, each pixel also has information of the direction where from the minimum distance was propagated. This information will also be used in latter stages. The distance transform image is shown in Figure 1(c). 3.3 Centreline image generation In distance transform image the saddle points mark the centrelines of channels and as well as the centrelines of non-channel areas which are bounded by TRUE pixels. Suppression of non-maxima distances (in this case along the direction from where minimum distance is propagated) yields an image where the saddle points are retained as TRUE while the rest of the image is turned FALSE. This image is further processed by skeletonization algorithm given by Zang and Suen (1984) which erodes the contour pixels (i.e. those which are having at least one FALSE neighbour) without breaking the connectivity, and without eroding the endpoints. This centreline image is shown in Figure 3(a). 3.4 Cluster generation using ISODATA At the beginning of programme the user is shown an image of the study site and prompted to click at channel points. The results presented in this paper are generated using 25 such points clicked arbitrarily. As seen in Figure 1(a), the channels are characterized by varying spectral responses. A channel appears different in an image from its surroundings due to its different sediment type, water level, moisture content, presence of shadow in particular if the channel is empty and aligned perpendicular to the azimuth of Sun direction, and the vegetation growth on its banks. Furthermore, at any cross section the channel may have one or more of these characteristics present. The assumption at this stage is that the user clicks at all types of spectral responses forming the channel. ![]() Figure 3: (a) Centreline image; (b) Mahalanobis classified image; (c) Channel centreline segments connected with sample points The pixel at clicked location is further grown to a small neighbourhood provided the digital numbers (DNs) do not change beyond ±2 DN. The sample data are collected for all image bands employed (three in present case). The samples thus collected are operated by ISODATA clustering technique (Rees, 2001). The software begins with a specified number of clusters, four in the present case (should be equal to the possible spectral responses that the channels exhibit). At the end of the ISODATA clustering four clusters along with their mean vector and variance-covariance matrix are generated. 3.5 Mahalanobis distance classification of image For each cluster formed in the aforesaid step a distribution function is generated. For this, the Mahalanobis distance of each pixel in cluster is computed from the cluster mean vector. All pixels in the clusters (N) are sorted as per their Mahalanobis distance and the pixels are given an Estimated P Value (EPV) as: EPV1 = 1 i.e. for pixel with least Mahalanobis distance. EPVi = 1-(i-1)/N where i = 2..N. To classify an unknown pixel of image the Mahalanobis distance between this and all the clusters are computed. It is checked that to which cluster the unknown pixel is nearest to. Furthermore, an EPV value is found for the Mahalanobis distance computed comparing it with the distribution generated in above step. The pixel is classified as channel pixel if its EPV is more than the least EPV of the cluster. All pixels of image undergo this process and get classified in either a channel pixel or a non-channel pixel. The classified image is shown in Figure 3(b). 3.6 Channel segment identification from sample points This step is the beginning of feature constrained classification which attempts to identify only the channel areas using the centreline generated, classified image and the sample points. The algorithm uses the coordinates of sample points clicked and locates the nearest centreline pixel which is also classified as channel. Search for this pixel from sample pixel progresses using the direction image generated in step 3.2, which ensures that the centreline pixel inside the channel is located and the search does not go out of channel. Having located this pixel all connected centreline pixels which are also TRUE in classified image are switched on. The end result of this process is shown in Figure 3 (c) where each segment is definitely a channel segment because of the control exercised in their generation. 3.7 Joining of channel segments The channel segments generated in above step are unconnected and incomplete. The channel joining mechanism attempts to join these segments amongst themselves and also with other channel segments present in the image. The pixels in channel segment image (Figure 3(c)) for which only one neighbour is TRUE are marked as endpoints. Besides the segments shown in Figure 3(c), there are many other centreline pixels which are also classified as channel. All these pixels are flagged as candidate pixels. In the process of joining, an endpoint may be joined with other endpoint (the second endpoint becomes candidate pixel in this case) or other candidate pixels. The direction of centreline at an endpoint is found as the relative direction with previous pixel in the segment. An isosceles triangle of 6 pixel height with the vertex at endpoint and triangle axis in the direction of the computed direction is employed to facilitate joining. To narrow down the search in the direction of channel the base-height ratio of this triangle is kept small (0.8 for vertical and horizontal directions and 0.6 for triangles in diagonal directions). Search is conducted for a candidate pixel scanning parallel to the base of triangle beginning from the vertex. The endpoint is joined with the first candidate pixel encountered using the shortest distance. The result of this processing is shown in Figure 4(a). ![]() Figure 4: (a) Result of joining endpoints with candidate pixels; (b) Channel segments of Figure 4(a) expanded till the edges along direction through which minimum distance in distance transform propagated ; (c) Channel image after dilation and erosion. 3.8 Fattening of channel segments In order to generate full channel extent the channel centreline segments traced in step 3.7 are expanded till their edges. Beginning from all pixels of the segments the algorithm tracks in the direction of edges along the route through which minimum distance was arrived at the segment pixel (using the direction generated in step 3.2). All pixels in this track are switched on. The tracking stops at the edges (Figure 4(b)). Furthermore, the resulting image is dilated twice and eroded twice to fill the gaps in the channels. Final output is shown in Figure 4(c). 4. Discussion A comparison of original image (Figure 1 (a)) and final channel image (Figure 4(c)) shows that the algorithm developed has been able to identify many of the channels along with their areal extent. However, there are some minor channel segments which have not been identified. The success of this algorithm is basically because of the feature constrained classification approach where the control is exercised by the sample points selected by the user. Furthermore, the Mahalanobis classified image and centreline image help in growing the sample point further so all connected components can be located. This approach should be able to identify all connected channel components. Failure of present method to locate some minor channels is primarily because of less sophistication in segment joining mechanism at present. This mechanism will stop at the first appropriate pixel found in the joining triangle, though the same pixel may not be part of the main channel centreline but an isolated or a very small segment. In addition, as the height of joining triangle is restricted to 6 pixels, the endpoints and candidate pixels apart by more than 6 pixels will not be joined and further channel tracing stops. This is evident in the lower part of image where despite the presence of channels (i.e. TURE pixels in both centreline and Mahalanobis classified image) no tracing has been done as the gap between two segments is wider. Attempt is being made to improve channel joining mechanism which will find the channel direction not on the basis of only two pixels but full segment and proceed in that direction. Furthermore, before flagging any candidate pixel as channel it will be assessed if this pixel continues further or ends there only. The pixels which continue further will be given higher weightage for joining with endpoint. 5. Conclusion The algorithm presented has been successful in identifying many channels in the image. More number of sample points will improve the performance of method, as these will produce a better classification and also the control on channel tracing will be more rigorous. This method of feature constrained classification can also be employed for other similar features which possess curvilinear shape and exhibit unique but spatially and temporally varying spectral characteristics. The algorithm presented needs several modifications and work is going on at present in this direction. 6. Acknowledgement: This project is supported by the grants from the University of Padua, Italy under CORILA. 7. References:
| ||
| © GISdevelopment.net. All rights reserved. |