TecLines : A MATLAB-Based Toolbox for Tectonic Lineament Analysis from Satellite Images and DEMs , Part 1 : Line Segment Detection and Extraction

Geological structures, such as faults and fractures, appear as image discontinuities or lineaments in remote sensing data. Geologic lineament mapping is a very important issue in geo-engineering, especially for construction site selection, seismic, and risk assessment, mineral exploration and hydrogeological research. Classical methods of lineaments extraction are based on semi-automated (or visual) interpretation of optical data and digital elevation models. We developed a freely available Matlab based toolbox TecLines (Tectonic Lineament Analysis) for locating and quantifying lineament patterns using satellite data and digital elevation models. TecLines consists of a set of functions including frequency filtering, spatial filtering, tensor voting, Hough transformation, and polynomial fitting. Due to differences in the mathematical background of the edge detection and edge linking procedure as well as the breadth of the methods, we introduce the approach in two-parts. In this first study, we present the steps that lead to edge detection. We introduce the data pre-processing using selected filters in spatial and frequency domains. We then describe the application of the tensor-voting framework to improve position and length accuracies of the detected lineaments. We demonstrate the robustness of the approach in a complex area in the northeast of Afghanistan using a panchromatic QUICKBIRD-2 image with 1-meter resolution. Finally, we compare the results of TecLines with manual lineament extraction, and other lineament extraction algorithms, as well as a published fault map of the study area. OPEN ACCESS Remote Sens. 2014, 6 5939

Traditionally, lineament mapping is based on a visual or semi-automatic interpretation (photo-geology).Lineaments are often extracted manually by digitizing, which is subjective, time consuming, expensive and requires expertise, training and adequate scientific skills.In addition, it cannot produce results for large scale areas [28].Lineament extraction could be more advantageous if the results were reproducible.Automatic lineament extraction is therefore needed.
Increasing spatial and radiometric resolution in satellite images favors the development of automatic, or criteria-based, lineament extraction algorithms [29,30].However, algorithms are usually slow and commonly generate false features that are related to artificial linear features such as roads and power lines [31].Linear feature enhancement and detection mostly have been done in the spatial domain (based on direct manipulation of pixels in an image) by using Sobel [32], Prewitt [32,33], and LOG filters [12,[34][35][36], as well as morphological filters [37] and frequency domain (based on modifying the spectral transform of an image).The binary edge maps that have been derived in the edge detection step have been used as inputs for extraction algorithms, such as edge following (graph searching [29]), edge linking operators (standard Hough Transform [27,35,[38][39][40], and edge tracing algorithms (STA, START, and ALERT algorithm [41]).
The success of automatic lineament extraction procedures depends on the reliability and accuracy of edge detection mechanism [28,42,43].Until now it has been impossible to recognize or measure linear features of interest that have to be detected and to remove all remaining artifacts [43].Selection of image preprocessing procedures and appropriate edge detection methods are very important in tectonic linear feature extraction, because they have significant effects on the accuracy of the final results.
The main goal of this study is to develop a new MATLAB based toolbox (TecLines) for automatic lineaments mapping from satellite images and digital elevation models (DEM).TecLines contains a set of functions for detecting and extracting potential edges with integration between image frequency and spatial filtering, Tensor voting framework, Hough transformation and polynomial fitting methodology.
The extraction of lineaments using the TecLines toolbox is performed in two main steps: edge detection and edge linking.Edge detection methods should ideally generate sets of pixels lying only on edges.In practice, these pixels seldom completely characterize edges because of unwanted noise and breaks.The general shape of the objects may be initially unknown, but in many cases they can be approximated by piecewise linear segments.It is not easily feasible to fit linear segments to all the edges in an image and discard false segments.These problems are addressed in the edge-linking step.The goal of edge linking is to describe an edge as a linear segment of specified shape and estimate the missing edge pixels from the assumed equation of the curvilinear segment.Due to differences in the mathematical background of the edge detection and edge linking procedure as well as the breadth of topics, we present the approach in a two-part paper.We chose MATLAB because it provides a multipurpose environment for mathematical processing based on a high-level programming language.Additionally, the functions provided in MATLAB are easy to modify and open for improvements.
The specific objective of the first part is to describe and assess the comprehensive edge detection procedure by integration of frequency domain filters, edge detection methods, and tensor voting framework.We introduce this part, with a focus on tectonic lineament extraction (binary edge maps with efficient thickness, length and pixel connectivity, and increased degree of accuracy of the edge detection).Finally, we validate this functionality using a synthetic image with known discontinuities and a high-resolution satellite image (QUICKBIRD-2) from an active tectonic area: the Andarab fault zone in NE Afghanistan.

Data
In this paper, we demonstrate the performance of the TecLines for edge detection, where validation has been performed on a synthetic and a real dataset.

Synthetic Dataset
In recent decades, several approaches for performance evaluation of edge detection methods have been proposed according to the presence or absence of ground truth data [44][45][46][47][48][49][50].These approaches are based on the characteristic of the images (i.e., real images, synthetic images) [44,[49][50][51].Most of the edge detection methods that rely on ground truth use simple synthetic images [52][53][54][55], because it is easy to specify the ground truth edge locations [44,45].In these cases, the edge detection can be quantitatively evaluated based on the known ideal detection considered to be the ground truth [55].The synthetic DEM (Figure 1) used here is the result of landscape evolution algorithm created using set river incision and different uplift rates across tectonic faults.The drainage system adapts to the evolving surface conditions.

Study Area and Data
We evaluated the performance of TecLines on satellite images of the active Andarab fault zone, northeastern part of Afghanistan (Figure 2a).The selected area has a long history of damaging earthquakes [56].Furthermore, it encompasses regions interpreted as deformed by primarily transtensional forces in the Trans-Himalayan orogenic belt [57,58].The Andarab fault is dextral and coincides with an approximately 150 km long, east-west-trending valley north of the intersection between the Paghman and Hari-Rud faults.In the high valley of Darya, many evidences of recent tectonic movements were observed along the fault trace [57].In this study, we used panchromatic band of the Quickbird-2 (1 m spatial resolution) for 2 March 2006 (Figure 2b).This data is in UTM coordinate system, datum -WGS84‖ and zone 42 N.

Figure 1.
The synthetic DEM that is the result of landscape evolution algorithm created using set river incision and different uplift rates across tectonic faults.The drainage system adapts to the evolving surface conditions.

Methodology
TecLines is a new MATLAB based framework that contains various functions for automatic detection and extraction of tectonic lineaments from satellite images and digital elevation models (DEM).Besides import and export functions that support the raster and vectors in standard file formats, TecLines provides functions for image filtering in the frequency and spatial domains to produce primary binary edge maps.Final binary edge maps in TecLines are produced by performing the computation of the Tensor voting framework.In addition, TecLines extract line segments from final binary edge maps by employing standard Hough transformation functions.A set of functions serves the grouping and merging line segments, which will be resulted in final lineaments maps.
Comparing results with published/non-published lineament maps and lineament analysis are also possible in TecLines.An overview of TecLines divisions can be found in Figure 3.

Frequency Domain Filtering
Image enhancement methods can contribute significantly to the automated extraction of linear features by using noise reduction and edges enhancement in the input images.There are two main categories for image enhancement methods: (1) spatial domain methods, such as Mean, Median, and Mode filters, and (2) frequency domain methods, such as Gaussian and Butterworth high, low, and band pass frequency filters.The main difference between both categories is that the spatial domain methods are usually applied locally, while the frequency domain methods are usually done in the global context.Unfortunately, there is no common principle for determining what is -good‖ image enhancement.However, when image enhancement methods are used as a pre-processing phase for other image processing methods, then quantitative measures can determine which methods are most appropriate.
In this study, we used Butterworth band-pass filter in the frequency domain Equation (1) in order to enhance edges by suppressing low frequencies and filter noise by attenuating high frequencies [20].The Butterworth band-pass method has several advantages with respect to other methods such as the smooth transfer function without any discontinuity or clear cut-off frequency to reduce the sharp truncation effect and the easy control of the range of frequencies passes (bandwidth) dependently on the order of the filter.In addition, the Butterworth band-pass method can improve the performance of the Ideal band-pass method while reducing the ring and blurring effect.
where   and   are the low and high cut frequencies respectively, f is the distance from the origin, and n is the order of the filter.The Butterworth filter is available in MATLAB but it shows some limitations.According to MATLAB documentation, numerical problems can happen for frequency filter orders as low as 6 when the transfer function coefficient form ([b, a] form) is used.These problems are due to round-off errors.The reason is that higher order filters may need extremely precise tap coefficients to get the desired performance.Therefore, numerical error in the calculation of the tap coefficients can destroy the performance.In fact, the computed filter is unstable.Accordingly, we implemented a new MATLAB based code to apply the band-pass filter on the satellite images.This code allows edge enhancement by frequency domain filtering at the specified band-pass parameters such as the order of the filter and the upper and lower frequency cutoffs of the filter for a given bandwidth that satisfies the rate of transmission at the cutoffs while maintaining the desired performance.

Spatial Domain Filtering
The principal physical edges correspond to significant variation in the reflectance, illumination, orientation and depth of scene surface [59].Conceptually, edge detection refers to the process of identifying and locating sharp discontinuities in an image in three steps: differentiation, smoothing and labeling [36,[60][61][62].Several types of methods are available for edge detection in the spatial domain.These methods are classified in: gradient based method (first order derivative), laplacian based methods (second order derivative) and optimal edge detection methods [59,63].We successfully tested and implemented in TecLines three common methods: Sobel and Prewitt [32,33] (gradient based), LOG [36,64,65] (Laplacian based), and Canny [60] (optimum) edge detection methods.In the gradient-based edge detection methods, the magnitude of the gradient vector reflects the strength of the edge, or edge response, at any given point.The effect of the noise in the signal will appear on the small local maxima in magnitude (edge strength), thus the resulting map of local maxima is thresholding to distinguish convincing edges [66][67][68].
In the Laplacian based methods such as LOG method, the gaussian smoothing filter is used for decreasing the sensitivity to noise in the input image, in other words to slightly blur the image.Then, Laplacian is applied to detect regions of rapid intensity change.The disadvantages of these edge detectors are their sensitivity to noise also producing two-pixel thick edges [69].
The Canny edge detection method is an optimum method and provides a multi-stage algorithm to detect a wide-range of edges in images.Those should be at a minimum distance to the actual edge in the real image.In addition, the detected edges should have minimal response.In other words, the discovered edges should have only one response to a single edge and completely eliminate the possibility of multiple responses to an edge by using adaptive thresholds with hysteresis.

Morphological Image Processing
Binary edge maps produced by Sobel, LOG, and Canny edge detection methods can have small-scale edges and isolated (or island) edge pixels.A common approach for improving these results is to use mathematical morphology methods to eliminate extra edge pixels.In this study, we used opening morphological filter using bwareaopen command in MATLAB.Morphological opening is an erosion followed by a dilation.Erosion eliminates small-scale details by removing outlying pixels and isolated pixels.Dilation restores all remained edges to their original size.The opening method presents many advantages because it relies only on the relative ordering of the pixel value, not on their numerical values.The opening is anti-extensive.In addition, opening is idempotent operation because it can be applied multiple times without changing the result beyond the initial application.

Tensor Voting Framework
Tensor voting is a non-iterative methodology to the inference of statistically salient features from possibly sparse and noisy data.In the tensor voting framework [70], the input data is encoded as elementary second order, symmetric, non-negative definite tensors (position/orientation pairs), then support information (including proximity and smoothness of continuity) is propagated by vote casting According to these principles for each possible dimension, each input token, where the input token is an unstructured point cloud with no priori orientation, encodes a local potential orientation [71].Therefore, every token is a location where an orientation is defined.The representation is in the form of a second order, symmetric, non-negative definite tensor which essentially indicates the saliency of each type of perceptual structure (curve, junction or region in 2-D) the token may belong to its preferred normal and tangent orientations [71,72]: T + λ 2 (e 1 e 1 T + e 2 e 2 T ) The tensor can specify its preferred tangent, normal orientation, as well as saliency corresponding to its perceptual structures.In Figure 4, the major axis e 1 is the preferred normal orientation of a potential curve segment.The magnitude of the stick component λ 1 − λ 2 indicates curve saliency.Tensor voting was implemented using tensor fields to pre-compute and store the votes from both stick and ball voters in receivers at various distances and angles (Figure 5) [70 -72].Vote orientation corresponds to the smoothest local curve continuation from voter to the recipient, while vote strength decays exponentially with distance and curvature.The second order vote also is a stick tensor and has a normal lying along the circular arc from voter O to recipient P (Figure 5).According to the Gestalt principles, vote orientation corresponds to the smoothest local curve continuation between two points, O as a voter and P as a recipient, while vote strength decays exponentially with distance and curvature.The magnitude of the vote is described by a function of confidence in spherical coordinates that the voter O and the receiver P indeed belong to the same perceptual structure.This magnitude can be calculated according to [70]: where S is the distance between voter (O) and receiver (P), K is the curvature, σ is the scale of voting, which determines the effective neighborhood size and C as a function of the scale controls the degree of decay with curvature (Figure 5).In the feature extraction process, a saliency map is produced for each feature type, which is assigning a second order symmetric tensor that estimates the structure of the feature type and the associated saliency.The characteristics of the elementary tensors in 2-D are shown in Table 1.
Figure 5. Second order vote casting by a stick tensor located at the origin [71].Where S is distance between voter (O) and receiver (P), K is the curvature, σ is the scale of voting, which determines the effective neighborhood size and C as a function of the scale controls the degree of decay with curvature.In TecLines, we wrote and implemented a set of MATLAB based functions for tensor voting to compute the gradient vector and the tensor matrix at the edge pixels of the binary image, which were extracted by edge detection methods.Then, voting tensor field generated by casting votes on all the edge pixels of the image were used to construct the stick saliency map.Finally, we extracted the distinguished edge pixels based on the stick saliency map.

Accuracy Measurements
There are specific statistical measures that we can use to assess erroneous measures [48,73,74].The difference between edge pixels obtained by edge detection methods and a reference map (ground truth) can be inferred [75,76].True positive is the number of correctly detected edge pixels, false positive is the number of pixels erroneously classified as edge pixels, and false negative is the amount of pixels that were not classified as edge pixels [77,78].
where TP is the total number of true positive edge.FN and FP are total number of false negative and false positive edges, respectively.     are the total number of edge pixels in the image and reference dataset, respectively.P is the subset of positions in the image at which edges occur.Accuracy requires that edges should be detected as close as possible to their correct positions.In a given image, the edge positions and numbers can be vary according to resolution and procedures [75].In this study, by comparing the edges detected using TecLines with a reference map, the accuracy was computed as follows.
× 100 (7) where Ac stands for accuracy.The value of Ac for accurate edge detection methods should be close to 100%.In order to evaluate the performance of the proposed edge detection procedure a reference dataset is required.The references dataset for both synthetic and QuickBird 2 images are determined based on manual extraction (Figure 6), because it is easy to distinguish geological and non-geological lineaments by visual interpretation [35,41,51,79,80].

Performance Evaluation of the Edge Detection Methods on a Synthetic DEM
Accuracy assessment is based on 21 known lineaments, which have been shown by 1180 edge pixels in the synthetic DEM (Figure 6b).The range of parameters to be used in edge detection procedure should be large enough to cover a wide range of detection results.We implemented 27 sets of parameter combinations of sigma, low and high thresholds for Canny edge detection method, where sigma ∈ (300, 500, 700), low ∈ (0.01, 0.05, 0.1), and high ∈ (1, 3, 5) In the opening morphological filter, minimum number of connected neighboring for each pixel was 4, 8, and 24 pixels.In the first step, we extract edges from original dataset without anterior Butterworth band-pass filtering (Figure 7a).The accuracy assessment shows that 30% of the known lineaments are detected as true positive edge pixels and 70% are detected as false negative edge pixels.The percentage of pixels that are false positives is about 45%.The tensor voting framework is applied to improve the result by considering the adjustment between increasing the information and decreasing the noise in the detection procedure.As seen in Figure 7b, the results after tensor voting framework demonstrate an increase in accuracy.40% of the known lineaments are correctly detected.The percentage of pixels that are detected as false negative edge pixels is around 60%.In the next step, we used a synthetic DEM preprocessed using a Butterworth band-pass filter in frequency domain.Due to local intensity changes in the image, edges are better detected in the west of the image (Figure 7c).The accuracy assessment of the Canny method shows that 78% of the known lineaments are detected as true positive edge pixels and 22% are detected as false negative edge pixel.Drainage patterns witness the critical relationship between streams, tectonic, and erosion changes [26,81,82].Tectonic activity disrupts drainage networks.Studying the nature of this disruption can give clues about the magnitude and orientation of the original tectonic activity [83][84][85].Thus, they are potential instruments for tectonic-geomorphology analysis [26,86].The result of the visual analysis shows that nearly 60% of the streams are also detected as edge pixels.Due to noise, the results also include many fragmented edge pixels that lead to the extremely low accuracy.
As seen in Figure 7d, 93% of the known lineaments are correctly detected after performing tensor voting framework.The percentage of false negatives is less than 7%.However, these results consist of many discontinuous edges, instead of edge contours.Results show that employing a Butterworth band-pass filtering prior to Canny edge detector improves the its performance.Additionally, using the tensor-voting framework both improved binary edge extraction and merged neighboring edges with similar direction.The overall accuracy with the proposed approach is about 52% and is higher than the overall accuracy achieved by the other methods (Table 2).In particular, a comparison with true positive edge pixels directly validates the superiority of the implemented edge detection procedure based on frequency band-pass filtering, edge detection, and tensor-voting framework implemented in TecLines.

Performance Evaluation of the Edge Detection Methods on a Satellite Image
In the first step, we applied a Butterworth band-pass filter in the frequency domain.We used 0.2 as an upper-lower cutoffs frequency and 6 for order numbers of filter such that these preserved as many significant edges (tectonic linear features) present in the truth as possible (Figure 8).Frequencies above and below these cutoffs are in the stop-band.Therefore, TecLines cannot find it in the edge detection step.Clearly, using a wide bandwidth and high order number of filter can decrease this assumption and preserve small edges but also increase the number of the wrong edges resulting from noises.After filtering in frequency domain, we used Sobel, LOG and Canny methods to detect linear features (edges).We used a 3 × 3 pixel masking window size in order to ensure the detection of sufficiently small feature and preserve image details.We classified image pixels in two groups: edge points (marked as 0) and non-edge points (marked as 255).We used a fixed threshold of 60 for Sobel filter.We selected a sigma (standard deviation of the noise) of 1.2 and a threshold of 20 for LOG filter.We set the hysteresis thresholding to 0.25 and sigma to 0.08 for Canny filter.The number of edge points produced in thresholding procedure can be modified on a case-by-case basis.The preliminary binary edge maps produced by Sobel, LOG and Canny filters are shown in Figure 9.
Table 2. Quantitative measures and overall accuracy obtained by TecLines for synthetic digital elevation models (DEM).True positive (TP) is the number of correctly detected edge pixels.False positive (FP) is the number of pixels erroneously classified as edge pixels.False negative (FN) is the amount of pixels that were not classified as edge pixels.We used opening morphological operation in order to thin the edges.This operation removes all pixels that are below a threshold of connected pixels.The minimum number of connected neighbors for each pixel was set to 4 pixels.The final binary edge maps are produced after removing all residual pixels.We used the tensor-voting framework to improve edge detected map accuracy and merge neighboring edges with similar direction.This provides better results than other methods for detecting common edges.It should be noted that using tensor-voting framework could lead to loss of small edges surrounded by bigger ones.After the voting process, a saliency map is produced by assigning a second order symmetric tensor that estimates the structure of the feature type and the associated saliency (Figure 10).Then, edges are extracted based on that certain orientations of the features coexist at a given location (Figure 10).If we compare the extracted edges with and without the use of tensor voting framework it appears that the extracted edges become more apparent and continuous after tensor voting.This technique allows the production of longer and smoother edge lines than the common edge detection methods.The accuracy was computed by comparing the edges detected using TecLines to the set of edges in the reference dataset (ground truth), The reference dataset is based on visual image interpretation and manual extraction.We analyzed the results before and after tensor voting framework.The results are shown in Table 3.This table shows that the tensor-voting framework significantly improves the true positive percentage in all of three edge detection methods.The number of false positive pixels is also drastically reduced.The highest overall accuracy was achieved by applying a combination of the Canny method and tensor voting framework and reaches about 74.5%.The accuracy of sobel and LOG methods were also improved after performing tensor-voting framework.

Conclusions
The aim of this work is to develop a MATLAB based toolbox (TecLines) in order to extract discontinuities from satellite images and digital elevation models.The high edge detection accuracy achieved by using TecLines shows that the proposed edge detection procedure can be used for geological purposes using high-resolution satellite images and DEM.The combined process of Butterworth band-pass filter in the frequency domain, edge detection methods in spatial domain and tensor voting framework slightly improves the edge detection accuracy.Consequently using the tensor-voting framework, the accuracy of edge detection is improved by 20% to 30%.We show that the combination of Sobel or LOG methods and tensor voting framework have similar overall accuracy for both the QuickBird and the synthetic images.However, the combination of the Canny and tensor-voting framework is considered more effective for edge detection.There are several important points that should be taken into account for further investigation.The percentage of achieved accuracy of edge detection methods largely depends on the gradient thresholds.The values of thresholds are still based on trial and error.For low thresholds, more edges will be detected but the resulting image will be increasingly susceptible to noise and trigger the detection of irrelevant features.A high threshold may miss subtle edges, or result in fragmented edges.Therefore, it is recommended to investigate further interactive thresholding in future edge detection methods.The outcome of this study is based on the panchromatic band of the QuickBird 2, which has a very high spatial resolution (1 m resolution).The proposed edge detection procedure is also adequate for medium resolution satellite images such as Landsat.The Landsat images have some advantages over QuickBird 2, such as free to availability, more coverage area, as well as shorter repeat intervals.While these results for edge detection procedures from satellite images are promising, the further mathematical approach presented in part 2 is still needed to extract linear segments from edges.In the second part of TecLines paper series, we describe the procedure of edge linking using the Hough transform and polynomial curve fitting using B-spline and Tavares methods.The quantitative representation of the binary edge maps can be exported as shape files, Geotiff and ENVI headers.Users can write their own MATLAB functions in order to expand the TecLines capabilities.The proposed toolbox is available from the TecLines website [87].This provides new opportunities for education and research in tectonic lineament analysis.

Figure 2 .
Figure 2. (a) location of the study area in Northeast Afghanistan.(b) panchromatic band of the Quickbird-2 (1-m spatial resolution) for 2 March 2006, of the study area.

Figure 3 .
Figure 3. Overview of the essential components of lineament mapping using TecLines.

Figure 6 .
Figure 6.(a) The reference lineament map for real dataset that is based on manual extraction from panchromatic band of QuickBird-2; (b) The reference map of the synthetic DEM consists in the digitized traces of the modeled faults (black line).

Figure 7 .
Figure 7. (a) and (b) preliminary binary edge map produced by Canny edge detection method and binary edge map obtained by tensor voting framework; (c) preliminary binary edge map produced by Canny edge detection method with prior Butterworth band-pass filter; (d) Final resulting binary edge map was obtained by tensor voting framework.

Figure 8 .
Figure 8.(a) Quickbird-2 (1 m spatial resolution) for 2 March 2006 of the study area shown in pseudo color; (b) Filtered image after Butterworth band-pass filter.Colors represent the image brightness digital number (DN) values.

Figure 9 .
Figure 9. Preliminary binary edge maps resulting from Sobel (a), LOG (b) and Canny (c) methods and morphological filtering.We used a fixed threshold of 60 for Sobel, a sigma (standard deviation of the noise) of 1.2 and a threshold of 20 for LOG, and the hysteresis thresholding to 0.25 and sigma to 0.08 for Canny.

Figure 10 .
Figure 10.(Left Pannel) Saliency density maps produced individually for (a) Sobel, (b) Log, and (c) Canny edge detection methods, respectively.The saliency density maps were constructed by using voting tensor fields that are generated from casting votes on all of the edge pixels.(Right Pannel) final extracted edge map by using tensor-voting framework after (a) Sobel, (b) Log, and (c) Canny filtering.

Table 1 .
Elementary [71,72] in 2-D, where n and t represent the normal and tangent vector respectively and the geometric features extracted after 2-D tensor voting[71,72].

Table 3 .
Quantitative measures obtained by TecLines for panchromatic band of QuickBird-2.True positive (TP) is the number of correctly detected edge pixels.False positive (FP) is the number of pixels erroneously classified as edge pixels.False negative (FN) is the amount of pixels that were not classified as edge pixels.