Image structure analysis for seismic interpretation
Proefschrift
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema, voorzitter van het College voor Promoties, in het openbaar te verdedigen op dinsdag 4 juni 2002 om 13.30 uur door
Peter BAKKER
doctorandus in de natuurkunde geboren te Linz Oostenrijk
Dit proefschrift is goedgekeurd door de promotor: Prof. dr. ir. L.J. van Vliet Samenstelling promotiecommissie: Rector Magnificus, prof.dr.ir. L.J. van Vliet, dr. P.W. Verbeek, prof.dr.ir. A. Gisolf, prof.dr.ir. R.L. Lagendijk, prof.dr.ir. F.A. Gerritsen, dr. G.C. Fehmers, dr. W.J. Niessen, prof.dr. I.T. Young,
voorzitter Technische Universiteit Delft, promotor Technische Universiteit Delft, toegevoegd promotor Technische Universiteit Delft Technische Universiteit Delft Philips Medical Systems Shell International Exploration and Production B.V. University Hospital Utrecht Technische Universiteit Delft, reserve lid
This project was financially supported by the Netherlands Ministry of Economic affairs, within the framework of the Innovation Oriented Research Programme (IOP Beeldverwerking, project number IBV97005).
Advanced School for Computing and Imaging
This work was carried out in graduate school ASCI. ASCI dissertation series number 78.
ISBN: 90-75691-08-4 c 2002, Peter Bakker, all rights reserved.
Contents 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Traditional interpretation of 3-D seismic data . . . . 1.2 Improving the efficiency of the interpretation process 1.2.1 Structure enhancement for horizon tracking . 1.2.2 Seismic attributes for detection . . . . . . . . 1.3 Image processing . . . . . . . . . . . . . . . . . . . . 1.3.1 Adaptive filtering . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
1 2 4 4 5 6 9
2. Linear structures . . . . . . . . . . . . . . . . . . . . . . 2.1 The representation of linear structures . . . . . . . . 2.1.1 Orientation representation . . . . . . . . . . . 2.1.2 Gradient structure tensor . . . . . . . . . . . 2.2 Orientation adaptive filtering . . . . . . . . . . . . . 2.3 Edge preserving filtering . . . . . . . . . . . . . . . . 2.3.1 Generalized Kuwahara filtering . . . . . . . . 2.3.2 Improving orientation estimation near borders 2.3.3 Edge preserving orientation adaptive filtering 2.4 Application: Automatic fault detection . . . . . . . . 2.5 Application: Structure enhancement . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
11 12 13 14 18 21 22 26 29 32 37
3. Curvilinear structures . . . . . . . . . . . . . . . 3.1 The GST for 3D plane-like curvilinear structures 3.1.1 The quadratic surface approximation . . 3.1.2 The quadratic GST for 3D surfaces . . . 3.1.3 Experimental tests and results . . . . . . 3.2 The GST for 2D curvilinear structures . . . . . 3.3 Curvature adaptive filtering . . . . . . . . . . . 3.4 Appendix A: Symmetry . . . . . . . . . . . . . 3.5 Appendix B: Implementation . . . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
41 42 42 43 45 51 52 55 56
4. Line-like curvilinear structures . . . . . . . . . . 4.1 The GST for 3D line-like curvilinear structures . 4.1.1 The quadratic curve approximation . . . 4.1.2 The quadratic GST for space curves . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
61 61 61 63
iv
Contents 4.2
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
66 67 70 71 71 74 77 79
5. Structural analysis using a non-parametric description . . 5.1 The tracking of line-like curvilinear structures . . . . . . . . 5.1.1 Application to the tracking of growth rings . . . . . . 5.1.2 Application to the tracking of sedimentary structures 5.2 Non-parametric adaptive filtering . . . . . . . . . . . . . . . 5.3 Non-parametric confidence estimation . . . . . . . . . . . . . 5.3.1 Application to channel detection . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
81 82 83 84 88 92 93
. . . . . . . . . . . . . . . . . . . . search . . . .
. . . . . . .
97 98 100 101 103 104 108
4.3 4.4 4.5 4.6
Experimental tests and results . 4.2.1 Circle image . . . . . . . 4.2.2 Helix image . . . . . . . 4.2.3 Ellipse image . . . . . . Discussion . . . . . . . . . . . . Application: Channel detection Appendix A: Symmetry . . . . Appendix B: Bias . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
6. Coherency estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Coherency based on the eigenstructure of the covariance matrix . 6.2 Coherency estimation using the GST . . . . . . . . . . . . . . . . 6.3 The presence of structural dip . . . . . . . . . . . . . . . . . . . . 6.3.1 Dip estimation by sampling the dip dependency . . . . . . 6.3.2 Comparison between the dip estimates of the GST and dip 6.4 An experimental comparison for fault detection . . . . . . . . . .
7. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.1 Image processing approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.2 Application to seismic interpretation . . . . . . . . . . . . . . . . . . . . . 114 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Samenvatting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Dankwoord . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Curriculum vitae . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
1. Introduction The growing global population and living standards cause an increasing demand for energy. Despite the many efforts made to exploit ‘new’ energy sources such as solar energy and bio-mass, oil and gas continue to be the primary sources of energy. The total amount of oil and gas produced each year is still increasing, and is very likely to continue to do so for at least 30 years. The oil industry is searching for new reservoirs on land and offshore, in increasingly difficult environments. Furthermore, the trend in the oil industry has changed from producing ‘at any cost’ in the seventies to a cost efficient and environmentally aware production today. Since the first seismic surveys, in the 1920’s, the seismic reflection method has played an important role in the exploration of oil and gas. The seismic method is a powerful remote sensing technique that can image the subsurface over depths from several meters to several kilometers. The basic idea is to first generate an acoustic wave field by a localized source. This field travels down the subsurface and partly reflects at locations where the acoustic rock properties change. The reflected wave field is measured by an array of localized receivers on the surface. The seismic method can be divided into three parts. It starts with the acquisition that consists of collecting raw data directly from the receivers. Usually, several different ‘shots’ are recorded of the same location. These different shot-records are stacked, or averaged, to reduce the influence of noise. Next, all the acquired data is processed to isolate the signal that corresponds to the travel-time of the reflected wave field from the surface to the reflectors. Signals due to diffraction and multiples should be suppressed. Multiples are wave fields that have been reflected more than once. Migration techniques are used to get a sharp, ‘focused’ image of the subsurface. Finally, the identification of possible oil and gas reserves is done by interpretation of the processed image using geological models and information from well measurements such as logs and bore-hole images. The first 3D seismic survey was shot over a field in Texas in 1967. Since then, there has been an increasingly rapid expansion in the application of this technology. A 3D seismic image I(x, y, t) has two spatial coordinates (x, y) parallel to the surface and one perpendicular time coordinate (t). The time usually corresponds to the two-way traveltime, i.e. the time it takes for a wave field to travel from the surface to a reflector and back. The typical spatial sample spacing is 12.5 m and the time resolution varies from 25 m in the shallow part of the data to 100 m in deeper parts. Although many geological features are still below seismic resolution, seismic images can in favorable circumstances reveal the
2
Chapter 1 Introduction
internal structure of a reservoir. A 2D cross-section of a 3D seismic image is shown in figure 1.1. The 3D image is the shallow part of a processed 3D seismic survey, showing the sea-floor. The total time interval of the image is 0.75 s, assuming that the acoustic velocity is 2000 m/s, this interval corresponds to 750 m. The traditional 2D surveys gave only a few cross-sections of the subsurface, but the 3D surveys give a full representation of the subsurface. This offers great advantages for both the processing and the interpretation of the seismic data. The development of 3D migration techniques has significantly improved the quality of seismic images. Densely sampled horizontal or time slices made it possible for the interpreter to map the lateral changes in the subsurface with a higher accuracy. Although much progress has been made in the interpretation process since the first seismic surveys, the current seismic interpretation methods still do not exploit the data to its full potential. Subtle geological features are often difficult to recognize in a seismic image for the human expert. Furthermore, the total amount of data that has to be interpreted grows more rapidly than the total number of interpreters, demanding an increasing efficiency. Image processing techniques can be used to enhance certain features in the data, and they can provide the basis for the automation of certain interpretation tasks. In this thesis we will present image processing techniques that can contribute to the automation of the interpretation of seismic data. 21.3 km t
Sea 0.75s
x Figure 1.1: A vertical 2D cross-section of a 3D seismic image. Three horizons have been tracked and they are displayed as white lines and pointed to by white arrows.
1.1
Traditional interpretation of 3-D seismic data
The introduction of 3D seismic images caused a radical change in the interpretation. Instead of having a dozen of paper sections, interpreters now have densely sampled 3D images and powerful workstations at their disposal. One of the direct advantages is the ability to freely choose the view angle, and it was soon discovered that a horizontal view or a time slice is very useful for the interpretation of depositional structures. Furthermore, the 24-bit color displays of modern workstations provide intuitive color maps with a large dynamic range.
1.1 Traditional interpretation of 3-D seismic data
3
Prior to the description of the traditional interpretation, we will introduce some terminology with the help of a simple geological model of the subsurface, in which we limit ourselves to sedimentary rocks. The formation of sedimentary rocks begins with the physical deposition of sediments on the surface. The physical properties such as the grain size, and the chemical properties of the sediments depend on the environmental conditions. During the continuation of this physical deposition process, the older sediments are buried deeper into the subsurface. The increasing of the pressure and the temperature with depth causes the conversion of the sediments into sedimentary rock. The sedimentary layers can be tilted, bend of fractured by the plate tectonic forces in the earth’s crust. A sequence of sedimentary layers formed in one type of sedimentological environment is called a stratigraphic sequence. A change of the sedimentological environment, for example caused by a sea-level change, means the beginning of a new stratigraphic sequence. Seismic interpretation begins with mapping the large scale structure of the area. This structural interpretation mainly consists of creating horizons and fault planes. Horizons are surfaces that are created by the interpreter by selecting a reflector and following it over the volume. As an example three horizons are shown as white lines in figure 1.1. There are several possible reasons why a reflector is selected to be interpreted as a horizon. The simplest reason is that the reflector is outstandingly clear and strong, making it easy to track. Sequence boundaries are important horizons to distinguish between the different geological periods. Another example of an important horizon is the top of a reservoir. A fracture in the subsurface rock caused by tectonic forces is called a fault. Faults cause discontinuities in the layered structure that make the creation of horizons more difficult. To be able to continue a horizon over a fault, it is necessary to know the amount of vertical displacement between both sides of the fault. It may be possible that one reflector seamlessly continues over a fault into a different reflector. To avoid that these two reflectors are incorrectly interpreted as belonging to one horizon, the entire fault-surface should be known. A common next step in the interpretation process is to map the structure inside a sedimentary layer. The horizons created in the structural interpretation can be used to approximate the depositional surfaces within the time interval where the reflectors are approximately parallel to the horizons. With a depositional surface we mean the earth’s surface in a previous geological time where paleodeposits were made upon. The sedimentary structures manifest themselves as lateral changes in the seismic response1 along the horizon, and can be recognized by their morphology. A channel for example, a sedimentary structure formed by flowing water, has a very distinctive meandering morphology.
1
In essence the shape of the seismic wavelet. The amplitude and the phase of the wavelet are the most commonly used descriptors for this shape.
4
1.2
Chapter 1 Introduction
Improving the efficiency of the interpretation process
The volume of the seismic datasets has increased tremendously over the last three decades, the area covered by the surveys has grown and the sample density has become larger. The last few years a lot of research is done on time lapse seismic. The goal of this technique is to monitor during the production, the changes inside a reservoir over time. Currently, these datasets contain only a few time samples, but installing a permanent array of receivers on top of the reservoir could lead to full 4D seismic images. Another development that potentially adds an extra dimension is the sampling of the amplitude versus the offset (AVO). This offset is the lateral distance between the start point of the acoustic wave and the point of reflection. These developments have already contributed to an increase of at least a factor 1000 of the total amount of data per survey. The demand for more data points will continue in the future. Seismic dataset currently require a couple of Terabytes of computer storage space. It is clear that the analysis of this increasing amount of data, requires an increase in the efficiency of the interpretation process. One way forward is to improve the visualization of the data. Interpreting 3D data by the visual inspection of 2D cross-sections is not optimal. Improved 3D visualizations are provided by immersive virtual reality rooms. The possibility to literally walk through the data also presents a new form of human-computer interaction. Another possibility to improve the visualization is to make a 2D animation of the 3D data by stepping through the consecutive slices at a high frame rate. Although an improved 3D visualization can considerably speedup the manual interpretation, it falls outside the scope of this thesis. The approach chosen in this thesis is to develop image processing tools for the automation of interpretation tasks. In the remainder of this section we will indicate the application domain of the tools presented in this thesis.
1.2.1
Structure enhancement for horizon tracking
Three dimensional seismic data made it possible to delineate faults and sedimentary structures laterally. Time slices, however, can cut through different depositional surfaces, showing only parts of these structures. This makes the interpretation of time slices difficult. The traditional solution is to approximate the depositional surface by picking a reflector and follow it over the seismic image. This process is commonly referred to as horizon tracking. Since it is very time consuming to manually mark every individual point of the horizon, algorithms have been developed to automate this task. A poor data quality or complex faulting can hamper the tracking considerably, and the tracking software often fails at many points on the horizon. The success of the automatic horizon tracker could be improved by enhancing the structural
1.2 Improving the efficiency of the interpretation process
5
information in the seismic image. We will present an adaptive filter method that suppresses random fluctuations and small deviations. This filter should be flexible enough to avoid filtering across faults, thereby preventing the merging of different reflectors from both sides of the fault.
1.2.2
Seismic attributes for detection
Geological features in the data that are difficult to estimate by the human visual system, can be made explicit by the computation of seismic attributes. In the field of image processing and pattern recognition attributes are referred to as features. Examples of conventional seismic attributes are the complex seismic trace attributes: instantaneous amplitude, phase, frequency, and bandwidth [TKS79]. A seismic trace is a vertical line of data in a seismic image I(x = x0 , y = y0 , t). Sedimentary structures are usually found by computing a seismic attribute in a certain time window around a horizon. Computing the variance of the seismic response within a time window, for example, projects all lateral changes in a seismic sequence onto a two dimensional surface. This allows the interpreter to check a volume of data for the occurrence of sedimentary structures by inspecting only one 2D display. The computation of the 3D attributes dip and azimuth along a horizon are very useful for the structural interpretation [DGS+ 89]. Since a few years computers have become powerful enough to compute 3D seismic attributes over an entire data volume, creating attribute volumes. In 1995 Bahorich and Farmer introduced a 3D seismic attribute they called the coherency cube [BF95]. Their algorithm computes the optimally lagged cross-correlation between neighboring traces. This was the first time a coherency measure was presented as a seismic attribute in the literature. Three dimensional seismic attributes like coherency, dip-magnitude, and dip-azimuth, can be used to highlight lateral changes. A time slice of a 3D attribute volume is much easier to interpret than a time slice of a seismic volume. Because each point in the attribute volume is computed using a 3D neighborhood, the detection of geological features is less dependent on the exact location, and therefore less dependent on carefully picked horizons. The main use of seismic attributes is to create an intuitive display that allows an interpretation task to be performed both more efficiently and more effectively. In this thesis we want to go one step further. We want to develop seismic attributes or features that allow the automatic detection of certain geological features. We will focus on the automatic detection of faults and the automatic detection of channels.
6
Chapter 1 Introduction
1.3
Image processing
Now that we have sketched our application domain, it is time to introduce the image processing approach. Although seismic interpretation is both the motivation and the application of this work, our goal is to develop image processing techniques that are general enough to be applicable to other natural images as well. We will describe the required properties for this more general class of images, and we will give some practical examples. Nevertheless, seismic images will continue to form the thread that runs through this introduction. Seismic images often show patterns with a layered structure due to the depositional nature of the subsurface, see for example figures 1.1 and 1.2a. In image processing a pattern with a certain regularity or structure is called a texture. The description of the ‘layered’ textures in seismic images can be split up in two parts. One part is the geometrical description of the structure, the other part is the description of the signal perpendicular to the layered structure. Examples of geometrical properties are the orientation and the curvature of the layered structure. An example of a property of the perpendicular signal is its characteristic frequency. In the case of a seismic image, the perpendicular signal is determined by the change in the acoustic impedance of the subsurface rock, convolved with the seismic wavelet. This convolved signal is usually described by using a time-frequency representation [TKS79, Ste97, Des97]. The main subject of this thesis is the geometrical description of the structure of layered textures.
(a) seismic
(b) wood
(c) finger print
(d) interference
Figure 1.2: Examples of natural images containing ‘layered’ textures.
The geometrical description is insensitive to the properties of the perpendicular signal2 , and is therefore applicable to all images with a layered structure. Some practical examples of images belonging to that class are shown in figure 1.2. Besides a seismic image (a), a 2D cross-section of a CT image of a tree-trunk is shown (b). This image can be used for the counting of growth rings. Furthermore, a finger print image (c) and an image of an interference pattern of a vibrating plate (d) are shown. 2
The only demand on the perpendicular signal is that it is not constant. In other words, its spectrum should contain frequencies other the zero.
1.3 Image processing
7
The individual layers of a layered texture in an image can locally be approximated with isophotes. Isophotes are curves and surfaces with a constant intensity value. For the geometrical description of these isophotes, we do not analyze the image intensities directly. We analyze their spatial relationship to characterize the shape of the isophotes. The first task for the geometrical description is to define properties that can describe this shape. The relevant properties depend, however, on the scale of the analysis. At the smallest scale, the level of the individual points of the image, only the intensity of the image point itself can be measured. The definition of geometrical properties at this scale is therefore not possible. At the largest or global scale, the description of the isophotes in general becomes very complex. The geometrical properties should therefore be defined and estimated locally, which means at some intermediate scale. The scale of the geometric analysis is not directly related to the frequency content of the data, but merely determines the spatial extent of the data that has to be modeled. In general, the optimal scale for the estimation of some feature may vary over the image. The feature may also be present at several scales simultaneously. Both cases require the estimation of the feature at multiple scales.
Parametric description The parametric description of the local geometry is based on the heuristic that decreasing the analysis scale decreases the geometric complexity. In other words, at a small scale a layered structure can be described using a simple geometric model. This complexity can be compared with the number of polynomial terms necessary to approximate an arbitrary function around a certain point. The larger the area around the point that needs to be modeled correctly the more terms are needed. Since we want to describe the differential geometry, we first compute the gradient at each point, then we apply a parametric model to the gradient information. In chapter 2 we start the analysis with the simplest gradient model, namely a straight model. A straight structure is fully determined by its orientation, and the representation and estimation of orientation is therefore discussed in detail in chapter 2. Orientation representation and estimation is an intensively studied subject in image processing [Knu89, SF96, Mar97, GVV97, KHRV99, CST00]. An important property of an orientation representation is its modality. Imagine two straight layered textures3 at different angles as drawn in figure 1.3. This configuration can occur in seismic images at sequence boundaries. The textures inside the windows 1 and 2, can be described using a single orientation, thus using a uni-modal representation. The description of the texture inside window 3, requires two orientations, and therefore a multi-modal orientation representation. Uni-modal representations are in general more computationally efficient, and we will therefore avoid multi-modal orientation estimation if possible. Examining the images in figure 1.2, one can conclude that modeling these images as locally 3
Straight layered textures are more commonly referred to as oriented textures.
8
Chapter 1 Introduction
1
3
2
Figure 1.3: Schematic display of two touching layered textures. The textures are analyzed inside a circular window at three different locations indicated by circles.
straight is feasible. It is also clear from these images that the curvature of the structure becomes more important for larger scales. In chapters 3 and 4 we will extend the local straight models by incorporating curvature.
Non-parametric description A parametric description of the geometry at an even larger scale requires more parameters. The estimation of these higher order parameters increases the computational demand, and the applicability of these parameters is limited. In chapter 5 we therefore turn to a nonparametric description. We will use the orientation and curvature estimates for a piece-wise description of layered structures.
Detection Faults manifest themselves in seismic images as discontinuities in the layered structure, as can be seen in figure 1.2a. Modeling this seismic image as locally straight gives an accurate description of the largest part of the data, but fails on the faults. In other words, if we define a resemblance measure between a planar reflector and the local data, then we can use this resemblance measure for the detection of faults. A fault detection method based on this idea is presented at the end of chapter 2. Closely related to this approach is the estimation of seismic coherency. Different measures for the estimation of coherency are recently published in [MKF98, GM99, MK00], and some of the techniques presented in these papers will be discussed in chapter 6. These articles show that the seismic coherency is lower at both faults and channel banks. For the delineation of faults one should use a large temporal window and a small spatial window, while the delineation of channels and other depositional structures requires a large spatial window size and a small temporal window. This is due to the fact that faults can continue across several different geological time periods, but depositional structures are restricted to a single time period. The typical vertical extent of a fault is 100-1000m and for a depositional structure 10-100m. The image processing techniques developed in chapters 4 and 5 are applied to the detection of channels.
1.3 Image processing
1.3.1
9
Adaptive filtering
Once we have a description of the structure of a texture, we can use this description for the enhancement of this structure. We now define the description as the model of the data and any deviation of the data from this supposed model as the noise in the data. Using these definitions, we can approach structure enhancement as a noise reduction problem. Noise in images is traditionally suppressed by low-pass filtering. This assumes that the signal at some point changes slowly compared to the noise, i.e. the signal and the noise are spectrally separable. In the case of layered textures this is only partly true. The signal changes slowly inside the layers, and low-pass filtering across the layers suppresses both the signal and the noise. The suppression of the signal can be avoided by locally adapting the shape of the low-pass filter window in such a way that the window ‘fits’ inside of a single layer. The local shape of the layers is estimated using the parameter estimates of the local geometric model. Parameters σ
Input image
φ
κ
.......
Adaptive filter
Output image
Figure 1.4: Generic adaptive filtering scheme. Parameters (σ, φ, κ, . . . ) are locally estimated on the input image to locally control an adaptive filter.
The processing scheme of adaptive filtering is depicted in figure 1.4. First the control parameters are estimated. Candidate parameters are, in the case of structure enhancement, scale, orientation and curvature (σ, φ, κ). Next, the parameters are used to adapt the shape of the filter window. In the case of the straight model, this results in disc-shaped windows steered to the local orientation. This kind of adaptive filtering can also be used as a method for providing the analysis window that yields an optimal property estimate. As we have seen above, the optimal windows for the detection of faults and channels have an anisotropic shape. The fault detection window is elongated along one dimension that should be steered perpendicular to the horizons. The window for the detection of channels, on the other hand, is elongated in two dimensions and this disc-shaped window should be steered parallel to the horizon of interest.
2. Linear structures Classical image processing filters for smoothing and edge detection, are designed for very simple local neighborhoods: homogeneous regions separated by discontinuities. However, a local neighborhood can also contain a pattern. If this pattern has some form of ordered structure, then it is called a texture. Typical examples of textures are the patterns on different types of cloth, in wood, or brick walls. In this chapter we are going to study a specific type of texture, namely oriented textures. An oriented texture is a pattern with a linear structure. We start the analysis of oriented textures by giving a definition of linear structure.
(a)
(b)
(c)
Figure 2.1: Images that consist of multiple domains. In the first image (a) the domains are characterized by the intensity. The domains in images (b) and (c) contain oriented textures, which are characterized by the dominant orientation.
definition A linear structure in N-D space is shift invariant in at least one orientation, but not in all orientations. For the mathematical definition of N-D linear structures, we describe a neighborhood as a function f (x). A neighborhood f (x) is shift invariant along x1 , if f (x) ≡ f (x1 , x2 , · · · , xn ) = f (x1 + d, x2 , · · · , xn )
(2.1)
f (x) = f (x2 , x3 , · · · , xn ).
(2.2)
for all values of the scalar d. Since the neighborhood f does not change as a function of x1 , we can write if as a function with a reduced dimensionality Thus we also can define linear structures by the reduction of dimensionality. A neighborhood f has a linear structure if there exists a basis β, such that fls (a) ≡ fls (a1 , a2 , · · · , an ) = f (a1+i , a2+i , · · · , an ) , i ∈ {1, 2, · · · , n − 1}
(2.3)
12
Chapter 2 Linear structures
with a = [x]β , the coordinate vector in basis β. We call the axes of the coordinate system (a1 , a2 , · · · , an ), the principal axes of the linear structure. In the two-dimensional case, there are two linear independent orientations, allowing only one possible linear structure. Namely, with one shift invariant orientation, and the other orientations not shift invariant. Two examples of two-dimensional oriented textures are given in figure 2.1b and c. In n-dimensional data, there are n − 1 types of linear structures that differ by the number of independent shift invariant orientations (i). In this thesis we differentiate between two types of linear structures: line-like and plane-like. A linelike structure is shift invariant along only one orientation (i = 1), and a plane-like linear structure is shift invariant along two orientations (i = 2). The simplest examples of lineand plane-like linear structures are a single line and a single plane respectively. A line-like oriented texture can also be regarded as a collection of individual lines, and therefore be referred to as a line-bundle. Note that plane-like linear structures are only meaningful if the dimensionality n > 2. Granlund and Knutsson [GK95] defined the simple or one-dimensional neighborhood as the building block of their analysis. In a simple neighborhood the intensity fs only changes along one orientation and can be mathematically represented by fs (x) = f (x · n).
(2.4)
We do not use this structure definition directly, because it has no clear geometrical meaning. Simple neighborhoods are identical to line-like linear structures in 2D, and identical to plane-like linear structures in 3D. The geometrical interpretation of line-like and plane-like linear structures is the same for arbitrary dimensional data. The geometrical interpretation of a simple neighborhood on the other hand, is different for each dimensionality.
2.1
The representation of linear structures
In this section we address the problem of finding a representation for a neighborhood with a linear structure, that is suitable for further processing. A linear structure can be characterized by the orientation of its principal axes and the number of shift invariant orientations (i). A property of linear structures is that they are invariant under pointinversion fls (x) = fls (−x). (2.5) The representation should have this property as well. In two dimensions this property corresponds to the fact that 180 degrees rotation of a line amounts to no change at all. This structure should therefore be represented by orientation, which is defined modulo π (180◦ ) as opposed to direction which is defined modulo 2π (360◦ ). An orientation estimate is a valuable feature that allows the segmentation of images with several oriented textures. If we want to detect the borders between the regions in the
2.1 The representation of linear structures
13
images in figure 2.1, then we can use traditional edge detectors for the borders between the constant intensity domains. The detection of borders between different oriented textures requires the estimation of orientation. Later on in this thesis we will use orientation to steer an adaptive filter. Another example of using the output of an orientation estimator for further processing, is to measure the change in orientation, i.e. the curvature [GWVV99]. Stabilizing the orientation estimate by averaging is often required for these tasks. The representation should therefore be continuous.
2.1.1
Orientation representation
We start by studying the two dimensional case. The most obvious representation of orientation is a scalar quantity. For example the rotation angle φ, corresponding to the rotation that aligns the principal axes of the structure with the coordinate system. While it is easy to restrict the scalar quantity to the interval [0, π), it is impossible to make it continuous at the same time. For each arbitrary interval there is a ’jump’ between the two end-points of that interval, i.e. it cannot be made cyclic. Think of visualizing orientation with grey-values; at some angle there is always a jump from black to white, see figure 2.2.
Figure 2.2: Image of an oriented texture and its local orientation. Visualizing orientation with grey-values suffers from two unavoidable black to white jumps.
Another possibility is to represent orientation by a vector x, e.g. the gradient vector. Although this vector gives a continuous representation, it is defined mod 2π. The vectors x and −x do not map to the same value. Limiting the vectors to a two dimensional half plane makes the representation discontinuous. The correct procedure is to double the angle. r cos φ r cos 2φ → (2.6) r sin φ r sin 2φ This vector representation was introduced by Granlund [Gra78]. A representation of 3D orientation was given by Knutsson in[Knu85]. He addressed the problem by stating that a suitable representation should meet three basic requirements: 1. ’uniqueness’: the vectors x and −x should map to the same value, 2. ’uniform stretch’: preservation of the angle metric of the original space.
14
Chapter 2 Linear structures 3. ’polar separability’: while the angle may change, the magnitude of the mapped vector should only be a function of the magnitude of the original vector.
The mapping he found was a complicated mapping to a 5D vector. Later Knutsson [Knu89] found a mapping M that maps the vector x onto the tensor T defined by M: T≡
xxT . ||x||
(2.7)
This tensor mapping also meets all three criteria and is therefore suitable for further processing. The mapping consists of applying the dyadic product and normalizing the resulting tensor with the magnitude of the vector. To be able to compare the mapping M with the mapping in eq. 2.6, we show the two-dimensional case in polar coordinates. r 1 + cos 2φ cos φ sin 2φ r → (2.8) sin φ sin 2φ 1 − cos 2φ 2 The major advantage of the tensor mapping is that it also holds for higher dimensional images. Another interesting property of the tensor mapping is its inverse. The space spanned by the tensors T that are created by mapped vectors according to eq. 2.7, is only a sub-space of the entire tensor space. The least squares approximation Tls to a tensor outside of this sub-space T0 is given by Tls = λ1 e1 eT1 ,
(2.9)
where λ1 is the largest eigenvalue of T0 and e1 the corresponding eigenvector.
2.1.2
Gradient structure tensor
The structure tensor is defined as T≡
xxT , ||x||n
(2.10)
where ( ) indicates some weighted local average. The computation the structure tensor consists of two steps. First the orientation is estimated for each point in the image using orientation selective filters. Next, the filter outputs are mapped to the tensor representation and averaged. The normalization of the tensors with ||x||n , determines how the orientation averaging is weighted with the corresponding intensity contrast. The structure tensor can be used as a tool for the local analysis of linear structures. Several different implementations and applications of the structure tensor have been reported in the literature. The earliest publications are [KW87, RS91, Hag92], and the list of applications has steadily grown over the last decade. An efficient implementation of the structure tensor is the gradient structure tensor (GST). The first step is to estimate the gradient g = ∇I at scale σg . We compute the gradients by
2.1 The representation of linear structures
15
convolving the image with the first order derivatives of a Gaussian. For an N-dimensional image the components of the gradients are gi = I(x) ⊗
∂ G(x; σg ) , i ∈ {1, . . . , N } ∂xi
(2.11)
The second step consists of mapping the gradient to the tensor using the dyadic product and average the tensor components Tij at scale σT . The gradient structure tensor is defined by T ≡ ggT , (2.12)
and the mapping used is equivalent to eq. 2.10, with n = 0. This means that the elements of the GST can be interpreted as gradient energies. We compute the local average, or spatial integration, by convolving the tensor components with a Gaussian kernel. Tij = Tij ⊗ G(x; σT )
(2.13)
The tensor scale is usually chosen three to ten times the gradient scale: 3σg < σT < 10σg . Averaging the tensor has a number of advantages. Rapid changes in the orientation due to noise on the gradient vector are suppressed, yielding a smooth orientation estimate. Furthermore, the GST not only contains information about the gradient energy in the orientation of the maximum gradient, but also in all perpendicular orientations. This extra information allows for the differentiation between different types of local structures, as we will see below. As an example, we show the GST for a 3D local neighborhood f (x, y, z) using the derivative notation by indexes. 2 fx fx fy fx fz fx g = fy T = fx fy fy2 fy fz (2.14) 2 fz fx fz fy fz fz
The computation of the GST in 3D requires nine convolutions, three for the gradient and six for the tensor smoothing. Since the tensor is symmetric only N (N + 1)/2 components have to be processed, where N is the dimensionality of the image.
The relevant information contained by the GST is extracted by computing its eigenvalues and eigenvectors. Due to the way it is constructed the GST is symmetric and semi-positive definite. This means that the eigenvalues are real and positive. By diagonalizing the tensor, the eigenvalue analysis finds the the coordinate system that is aligned with the principal axes. This aligned coordinate system is given by the set of eigenvectors, and we call it the gauge-coordinate system. The gradient energies along the principal axes are given by the eigenvalues, and the largest eigenvalue corresponds to the dominant orientation. For high dimensional cases (N > 3), the eigenvalue problem must be solved numerically, but in the 2D and 3D cases the eigenvalues can be found analytically. The structure tensor can also be based on other orientation filters than the gradient filters we use. Knutsson [Knu89] used angular separable quadrature filters to estimate orientation.
16
Chapter 2 Linear structures
The advantage of quadrature filters over derivative filters is that they are phase invariant, i.e. give a response on both even and odd structures. Derivative filters only give a response on odd structures, e.g. on the edges and not on the ridge of a line. This limited support no longer constitutes a problem since the tensor averaging combines the gradient information from both slopes of a line, without the cancellation of opposite vectors. The minimum number of convolutions needed to estimate 3D orientation using the quadrature filters in 3D is 12, while the computation of the gradient only requires 3 convolutions. Furthermore, mapping of the gradient to the tensor is much less complicated than the mapping of the quadrature filters. The structure tensor analysis assumes a uni-modal orientation distribution. The dominant orientation corresponds to the weighted average orientation. In the case of multiple maxima in orientation histogram, the dominant orientation will, in general, not coincide with a maximum.
Interpretation of the eigenvalues The richness of the structure tensor analysis becomes apparent if we study the eigenvalues. Throughout the thesis we assume that the eigenvalues are ordered, i.e. λi > λi+1 . In two dimensions we can distinguish three different cases, corresponding to different types of local neighborhoods. They are given in the table below: λ1 0
λ2 0
>0 >0
0 >0
description Both eigenvalues are zero. No gradient energy, no contrast. Constant intensity with no measurable structure. One eigenvalue is zero. Linear structure. Both eigenvalues are greater than zero. The underlying structure deviates from the linear structure model, e.g. due to noise or curvature or multiple orientations. If λ1 = λ2 , then the structure is isotropic.
In practice, the eigenvalues should not be checked against zero, but against a threshold value that is determined by the noise and the signal energy in the image. Without a quantitative estimate of the noise level, it will not be possible to distinguish between noise and isotropic structures. The total gradient energy is given by the trace of the GST Eg = T r(T) =
X
λi .
(2.15)
i
Contrast independent measures can be constructed dividing the eigenvalues by the total energy Eg . In this thesis we use the following contrast independent confidence measure Can =
λ1 − λ 2 . λ1 + λ 2
(2.16)
2.1 The representation of linear structures
17
It takes values between 0 and 1 and indicates how much the local data resembles a linear structure. The next table sums up the behavior of Can . isotropic: Linear structure: (anisotropic)
λ 1 ≈ λ2 λ1 λ 2
Can ≈ 0 Can ≈ 1
The more isotropic a structure becomes (Can → 0), the more difficult it becomes to estimate the orientation of that structure. Therefore we use Can as the confidence measure of the orientation estimation. Since a linear structure can also be viewed of as an anisotropic structure, Can is also referred to as the anisotropy. The analysis of the eigenvalues of the structure tensor in 3D is similar to the analysis in 2D. The extra dimension introduces the plane-like linear structures and the different cases are given in the table below: λ1 0
λ2 0
λ3 0
description all eigenvalues are zero. No gradient energy, no contrast. Constant intensity with no measurable structure. >0 0 0 Two eigenvalues are zero. Plane-like linear structure. >0 >0 0 One eigenvalue is zero. Line-like linear structure. If λ1 = λ2 , then the cross-section of the line-like structure is isotropic, e.g. cylindrical. > 0 > 0 > 0 All eigenvalues are greater than zero. The underlying structure deviates from the linear structure model. If λ1 = λ2 = λ3 , then the structure is isotropic. Again we like to have contrast independent measures. With three eigenvalues we are able to construct two mutually independent confidence measures [KBV+ 99] Cplane =
λ1 − λ 2 λ2 − λ 3 , Cline = , λ1 + λ 2 λ2 + λ 3
(2.17)
and they both take values between 0 and 1. The confidence measure Cplane can be interpreted as how much the neighborhood resembles a plane-like structure, and Cline indicates the resemblance to line-like structures. These measures can differentiate between the following local structures. isotropic: Plane-like: Line-like:
λ 1 ≈ λ2 ≈ λ3 λ 1 λ2 ≈ λ3 λ 1 ≈ λ2 λ3
Cplane ≈ 0 Cplane ≈ 1 Cplane ≈ 0
Cline ≈ 0 Cline ≈ 0 Cline ≈ 1
Application of the GST to seismic data The seismic method is a remote sensing technique for the inspection of the outer layer of the earth’s crust with a maximum depth in the order of kilometers. This relatively thin
18
Chapter 2 Linear structures
outer layer mainly consists of sedimentary rock which has a stratified structure. If there is sufficient contrast in the acoustic impedance of the different layers that were deposited on top of each other, then the stratification will be visible in the seismic image. Depositional surfaces, or more specifically the corresponding reflectors and seismic sequences, can locally be modeled as a plane-like linear structure. The GST is an efficient image processing tool that incorporates this model. We have seen that the GST yields parameter estimates of the local linear model and confidence measures for the resemblance of the model to the actual data. Figure 2.3a contains a two dimensional seismic image of a vertical cross-section of an anticline. In a two dimensional vertical cross-section, the planar structure of the reflectors reduces to a line-like structure, and can therefore be analyzed with the 2D GST. We have processed this image with the 2D GST (σg = 1, σT = 4) and the resulting estimates for orientation, gradient energy Eg , and confidence Can are shown in figure 2.3b-d. The GST yields a smooth and continuous estimate of the slowly varying structural orientation due to the curvature of the anticline. The strong reflectors are clearly highlighted in the energy estimate and the confidence estimate gives lower values if the reflectors are discontinuous.
2.2
Orientation adaptive filtering
Random measurement noise, which is present in every real world image, hampers manual interpretation by human experts as well as automatic segmentation and analysis by computers. Therefore many image processing techniques are developed to reduce noise. These methods are either based on spatial correlation or a spectral analysis. Additive uncorrelated noise on signals that only contain low frequency components, can be reduced in a straight forward manner by low-pass filtering. If the spectrum of the noise-free image as well as the spectrum of the noise are known, then the best linear filter to separate the two signals is the Wiener filter [Wie49, Pra72]. These spectra are in general not available. In the absence of a priori knowledge, the noise is usually assumed to be uncorrelated. Edges or transitions, which are imported for the analysis of images, often form lines or surfaces that can be locally modeled as a linear structure. In this linear structure model we have local a priori knowledge of the signal, namely shift invariance along one or more orientations. A linear filter such as the Wiener filter is not able to exploit this local information. The reduction of noise in an oriented texture domain or along an individual contour needs an anisotropic smoothing operator that adapts to the local orientation. Candidate approaches are elongated steerable filters [Hag92] [Fre92], and anisotropic diffusion [Wei98]. For this thesis we used a very flexible adaptive filter model. The adaptive behavior is obtained by local transformation of the image data. The transformation is formalized as a coordinate transformation based on locally estimated parameters. The transformed data
2.2 Orientation adaptive filtering
19
(a)
(b)
(c)
(d)
Figure 2.3: Analyzing a seismic image of a vertical cross-section of an anticline (a), using the 2D GST. Resulting in estimates of local orientation (b), gradient energy (c), and confidence or anisotropy (d).
20
Chapter 2 Linear structures
is computed by interpolation of the image data. Many different control parameters can be incorporated by simply adjusting the coordinate transform. Since this adaptive mechanism does not depend on the filter type, it can be used for non-linear filters as well, e.g. median or variance filters. Steering an elongated filter window requires a robust and continuous representation of orientation. Haglund [Hag92] showed that the structure tensor is suitable for controlling an orientation adaptive filter. Our orientation adaptive filter consists of an elongated Gaussian filter steered by the eigenvectors of the GST.
2 σ1 2 σ2 (a)
(b)
(c)
Figure 2.4: (a) Creating an elongated filter by combining isotropic filters. (b) The parameters σ1 and σ2 of a 2D elongated filter. (c) Steering an elongated filter to match the local structure.
An elongated filter can be created efficiently by combing several isotropic Gaussian filter responses, see figure 2.4. The responses are sampled at regular intervals along a straight line. If the intervals are smaller than two times the scale of the isotropic filter (∆ < 2σ), the resulting mask is essentially flat. A practical implementation of this filter consists of two stages. First an isotropic Gaussian filter is applied by simple convolution, setting the filter widths σ1 . . . σn equal to σ. This filter step can be interpreted as selecting the appropriate scale. A steered filter completes the task by elongating the filter along the shift invariant orientations. The steered filter is one-dimensional (σ1 ) for line-like structures and two-dimensional (σ1 ,σ2 ) for plane-like structures. When relevant information is present at the finest scale, the first of stage of the filter may be omitted. This is only allowed if the image is band-limited and sampled according to the sampling theory [Sha49]. The sampling or Nyguist theory requires that the the sampling frequency is higher than two times the maximum frequency of the analog signal fsample > 2fmax . To demonstrate the capabilities of an orientation adaptive filter, we applied it to a 2D image of a circular oriented texture, shown in figure 2.5. The parameter settings we used for the GST are σg = 1, σT = 5, and for the steered filter σ1 = 4. For comparison we also applied isotropic Gaussian filters at four different scales σ = {1, 2, 3, 4}. Due to the circular shape of the oriented texture in the test image, we can immediately verify from the result in figure 2.5c that the adaptive filter is rotation invariant. Furthermore, a limitation of the filter comes to light in the high curvature area. The elongated filter is no longer able to match the shape of the local texture in this area, and the filter error increases with increasing curvature. An efficient way to implement steerable filters is proposed in [FA91, Per95]. The presented technique allows a convolution based filter to be steered to an arbitrary orientation by
2.3 Edge preserving filtering
(a)
21
(b)
(c)
Figure 2.5: A demonstration of the need for anisotropic smoothing filters to improve the SNR of a noisy image of an oriented texture (a). (b) Isotropic Gaussian filters applied to (a) at four different scales. (c) The result of orientation adaptive filtering. The sizes of the filter are indicated by white drawings.
writing it as a linear combination of a set of linear basis filters. The basic idea is to sample the angular response of the filter and to interpolate between these samples. A draw-back of this method is that a narrow elongated filter (σ1 σ2 ) needs a lot of basis filters, since it has a high angular resolution. Furthermore, this technique does not work for non-linear filters, and the extension of this technique to curvature adaptive filtering is not trivial.
2.3
Edge preserving filtering
The goal of filtering is to get a better description of the local image properties by replacing the values of the individual points by a function of their neighborhood. The influence of random fluctuations in the image, for example, can be reduced by integration over a neighborhood. This integration assumes that the image properties of the individual points inside the neighborhood, are approximately constant. The assumption of constant image properties is not met when the neighborhood overlaps a border between two regions with different properties. Estimation of the images properties near such a border will result in a weighted average of the properties from the two different regions. This should obviously be avoided. The traditional example of this problem is low-pass filtering of an image with constant gray-value regions. Although it will reduce the noise, it will also blur sharp edges. A filter that avoids overlapping a border, is therefore said to be edge preserving. There are two principal mechanisms for filter windows to avoid a border. One of them is to make the filter window smaller near a border, i.e. scale adaptive filtering. The border is interpreted as a fine scale image feature and the scale of the filter should be adapted to this smaller scale. The second mechanism avoids borders by allowing decentralized windows that are ”repelled” by the border. A schematic visualization of the two mechanisms avoiding an orientation border is given in figure 2.6. Both mechanisms require an estimation of
22
Chapter 2 Linear structures
(a)
(b)
(c)
Figure 2.6: Two principal mechanisms to avoid a border. The behavior of a normal filter (a). Avoiding border by adapting the scale of the windows (b), and by using decentralized windows (c).
the border locations. The central issue is to transform this border estimate into the control parameter scale or displacement. An example of a filter based on decentralized windows that solves this issue is the generalized Kuwahara filter [BVV99].
2.3.1
Generalized Kuwahara filtering
A traditional filter for edge preserving smoothing of images containing constant grey-value regions, is the Kuwahara filter [Kea76]. Kuwahara divided a square symmetric neighborhood into four slightly overlapping windows, each containing the central pixel, see fig.2.7a. In each window the mean mi and the variance si of the intensity is computed. The edge preserving estimate of the mean mep is obtained by replacing the central pixel by the mean value of that window that has the smallest intensity variance. mep = mr , r = arg{min si } , i ∈ 1, 2, 3, 4 i
(2.18)
The mean operation (uniform filter) reduces the noise and the variance is used to select the most homogeneous window. The proposition on which the filter is based is that the variance in a window that overlaps an edge s2edge , is greater than the variance in a window in a homogeneous area s2homogenous , s2edge > s2homogenous .
(2.19)
The window with the smallest variance, therefore avoids edges. This filter has been further developed in [NM79] by increasing the number of windows to eight and changing the shape of the windows to pentagons and hexagons. The Kuwahara filter is an implementation of an edge preserving smoothing filter that uses decentralized windows. The crucial part of this class of filters is finding the optimal displacement. The generalized Kuwahara filter [BVV99] finds the optimal displacement
2.3 Edge preserving filtering
23
1
2
3
4
(a)
(b)
Figure 2.7: a) traditional Kuwahara filter, b) a realization of the generalized Kuwahara filter with circular windows. The dashed lines bound the neighborhood and the solid lines the windows.
by locating the highest confidence value. Since the confidence of the parameter estimation is lower near borders, the filter avoids them. By evaluating all decentralized windows that still contain the current pixel, the filter is made rotation invariant and maintains connectivity. Each decentralized window has a parameter and a confidence estimate (p i , ci ), the parameter from the window with the highest confidence value yields the edge preserving estimate pep . pep = pr , r = arg{max ci } (2.20) i
Making a parameter estimate edge-preserved using the generalized Kuwahara filter, therefore, requires the estimation the parameter and a corresponding confidence value. A realization of the generalized Kuwahara filter with round windows is depicted in figure 2.7b. In principal, there are no restrictions to the shape or size of the windows, except that they should both be fixed. Note that the shape of the window determines the shape of the neighborhood. To demonstrate the generalized Kuwahara filter we applied it to an image with three grey-value domains. We added uncorrelated noise to a SNR of 3 dB. The parameter to be estimated is the Gaussian mean and the confidence estimate minus the gradient magnitude (−|∇I|) is the confidence value, both computed at scale σ = 5. As can be seen from the results in figure 2.8, the Gaussian filter blurs the edges and the generalized Kuwahara filter preserves sharp edges. The sharp edges created by the generalized Kuwahara filter do not per definition coincide with the original edges of the noise-free image, since the certainty measure, on which the new edge location is based, is estimated on the noisy image. A mostly undesired side effect of the Kuwahara filter is a blemished result in regions without clear edges and the creation of false contours. These artifacts are due to the fact that the Kuwahara filter always selects. In homogeneous regions, however, the differences between the confidence values are due to noise, and the selection should not take place. A solution to this problem is to control the selection process by a function of the confidence value f (c). We used a linear mixing of the normal parameter estimate and the edge preserved version pnew = p(1 − f (c)) + pep f (c), 0 < f (c) < 1. (2.21)
24
Chapter 2 Linear structures
(a) Input image
(b) −|∇I|
(c) Gauss
(d) Generalized Kuwahara
(e) Controlled selection
(f) PM Diffusion λ = 0.06
(g) PM Diffusion λ = 0.08
(h) PM Diffusion λ = 0.1
Figure 2.8: Reducing the noise in an image with grey-value domains at 3 dB. The Gaussian filter (c) blurs the edges, but the Generalized Kuwahara filter (d) preserves them. Controlled selection using the confidence values (b) is shown in (e). For comparison the result of Perona and Malik diffusion [PM90] is shown in (f,g,h). 100 iterations are used in all cases.
2.3 Edge preserving filtering
25
The task of the function f is to suppress false border detections due to noise, and enhance true border detections. Clipping the confidence between the 5 and the 95 global percentile gives visually good results, see figure 2.8e. The lower clipping bound is based on the assumption that in 5 percent of the data, the contrast of the noise is larger than the contrast of the signal. Related methods Another approach to finding the optimal displacement of decentralized windows, is presented in [FS99]. In this paper Fischl and Schwartz split up the computation of the displacement vector field v(x) in three parts v(x) = m(x)d(x)o(x),
(2.22)
where m is the magnitude, d is the direction or sign, and o is the orientation of the displacement vectors. They state that the orientation should be perpendicular to the border, and the direction should be such that the decentralized window moves away from the border without crossing it. The magnitude or amount of displacement should move the window until it is inside a homogeneous region or until the influence of another border becomes dominant. They worked out this algorithm for the filtering of constant grey-value regions. The use the orientation of the gradient and the direction is found by locating the position of the maximum gradient magnitude. This yields their initial estimate vi of the displacement. The magnitude is found extending the initial vector along its direction until the inner product between the initial vector at ”head” and ”tail” becomes non-positive. In the case of a straight symmetric border, the displacement found by this algorithm is identical to the displacement of the generalized Kuwahara filter that uses the gradient magnitude as the confidence value. The complexity of the generalized Kuwahara filter is lower, since it is a one step algorithm and it does not explicitly enforce a displacement perpendicular to the border. Fischl and Schwartz compared their displaced window method to Perona and Malik diffusion [PM90]. Treating the image intensity as a conserved quantity being allowed to diffuse over time, allows image enhancement to be performed by solving the diffusion equation. Perona and Malik diffusion is a nonlinear version of diffusion that is controlled by the geometry. The magnitude of intensity gradient determines the amount of diffusion that is allowed in a certain direction. The corresponding diffusion equation is given by ∂I = div(g(|∇I|2 )∇I) (2.23) ∂t with 1 (λ > 0) (2.24) g(s2 ) = 1 + s2 /λ2 For a visual comparison we have applied the Perona and Malik diffusion filter to the test image in figure 2.8, for three values of the parameter λ (0.06,0.08,0.1). The Fischl and
26
Chapter 2 Linear structures
Schwartz algorithm is slightly more complex than our generalized Kuwahara filter, but one to two orders of magnitude faster than 50 iterations of the Perona and Malik diffusion. Fischl and Schwartz have also shown that their algorithm gives results that are visually comparable to the results of Perona and Malik.
2.3.2
Improving orientation estimation near borders
The structure analysis using the GST is based on the assumption that there is locally one dominant orientation, i.e. a uni-modal representation. The orientation estimation of the GST fails when there is locally more than one orientation present, for example at the border of two oriented textures. The resulting orientation estimation near that border is a weighted average of the two dominant orientations from both sides. A powerful solution for allowing multiple orientations in one neighborhood in the image, is to add orientation as a new dimension [Wal87, CH89, GVV97]. This orientation dimension needs to be sampled by applying a bank of directionally selective filters. The number of filters needed depends on the angular resolution. The drawbacks of this approach are the higher computational complexity and the demand for more computer memory, especially in 3D. Orientation in 2D can be represented by one angle, thus adding orientation as a dimension f (x, y) → f (x, y, φ) adds one dimension. The representation of orientation in 3D requires two angles. Adding orientation to a 3D image f (x, y, z) → f (x, y, z, φ, θ) amounts to a total of five dimensions. An image with multiple oriented textures can be correctly filtered without a full multimodal representation of the data, as long as the textures do not overlap. A correct orientation estimation near border can be obtained by applying the generalized Kuwahara filter to the output of the GST. We will now demonstrate this combination by example. We apply the GST to the three domains image and the orientation and confidence estimates are shown resp. in figure 2.9c and b. It can be verified from this result that the confidence estimate is lower near the borders. Applying the generalized Kuwahara filter with the GST orientation as the parameter estimate and Can as the confidence value, causes the GST windows to avoid the borders. The resulting improved orientation estimate is shown in figure 2.9d. A comparison of this result with the true orientation, shows that at some locations there is a small shift in the position of the border. Limitations of the GST confidence as orientation border detector The success of the method described above, highly depends on the ability of GST confidence estimate to locate the borders. We will now analyze a neighborhood around an orientation border point. First, the neighborhood is split up in two parts A and B, in such a way that neither of the these parts contain the border, see figure 2.10. Then the structure tensor is computed on both parts, according to
2.3 Edge preserving filtering
27
(a)
(c)
(b)
(d)
(e)
Figure 2.9: The orientation (c) and confidence estimate (b) of the GST applied to the threedomains image (a) at 9 dB. An improved orientation estimation near borders using the generalized Kuwahara filter with (c) as parameter and (b) as confidence. The true orientation (e) is given for comparison.
28
Chapter 2 Linear structures
e1A A
B (a)
e1
φ e 1B
φa φ b
(b)
(c)
Figure 2.10: Orientation estimation at a border point indicated by the small circle in (a). The neighborhood around this point is divided into parts A and B (a). The combination of the structure tensors from these parts is shown in (b). The tensors are depicted by the ellipses spanned by their eigenvectors, and the dashed line denotes the combined tensor.
Ti =
Z
x∈i
w(x)g(x)g(x)T dx , i ∈ {A, B}.
(2.25)
We assume for the confidence estimates Ci of these tensors λ1i − λ2i ≈ 1 , i ∈ {A, B}, (2.26) Ci = λ1i + λ2i i.e. the patterns in sub-neighborhoods A and B have approximately a one-dimensional structure. This means that A and B can be characterized by a single dominant orientation λ1 e1 , and the orientation border can be described by one angle φ. e1A · e1B φ = arccos (2.27) ||e1A || ||e1B || The ability of the GST confidence to detect the orientation border is determined by the angular dependence of the confidence estimate C(φ), of the combined tensor T = TA +TB . An expression for C(φ) can be found by analyzing the eigenvalues of T. The tensor T is brought to its diagonal form D, by writing it in the basis β spanned by the eigenvectors of T. D = [T]β = [TA ]β + [TB ]β = Q−1 (φa )DA Q(φa ) + Q(φb )DB Q−1 (φb )
(2.28)
where Q is a rotation in the two-dimensional sub-space P spanned by {e1A , e1B }, and φ = φa + φb as shown in figure 2.10. If λ2i = λ3i = · · · = λN i , i ∈ {A, B}, we can chose e2A and e2B such that they are elements of sub-space P. This gives the following equations λ1 = λ1A cos2 φa + λ1B cos2 φb + λ2A sin2 φa + λ2B sin2 φb λ2 = λ1A sin2 φa + λ1B sin2 φb + λ2A cos2 φa + λ2B cos2 φb
(2.29)
2.3 Edge preserving filtering
29
We can use these eigenvalues to construct C(φa , φb ). λ1 − λ 2 (λ1A − λ2A ) cos 2φa + (λ1B − λ2B ) cos 2φb = λ1 + λ 2 λ1A + λ2A + λ1B + λ2B = CA cos(2φa ) + CB cos(2φb )
C(φA , φb ) =
(2.30)
The angles (φa , φb ) can be found by solving 0 = (λ1A − λ2A ) sin(2φa ) − (λ1B − λ2B ) sin(2φb ) φ = φa + φb .
(2.31)
In the case that λ1A = λ1B ≡ λ01 , λ2A = λ2B ≡ λ02 → φa = φb =
φ , 2
we get λ01 − λ02 C(φ) = 0 cos(φ). λ1 + λ02
(2.32)
The angle φ is difference between the dominant orientations of the oriented textures at both sides of the border. As a consequence, at a border with a small angular difference, the confidence estimate C(φ) will only decrease by a few percent. For example, if φ = π6 , then the confidence will only decrease 10%. In the presence of noise it will be difficult to detect such a low-angle border.
2.3.3
Edge preserving orientation adaptive filtering
Earlier we have shown that noise in oriented texture domains can be reduced by orientation adaptive filtering. This filter fails near borders between oriented textures. The previous section shows that it is possible to improve the orientation estimate near borders. Using this improved orientation estimate, we are able to correctly steer an elongated window over images with several touching oriented textures. Still, we have to make sure that the adaptive filter does not overlap borders, see figure 2.11a. This is in fact the same problem we encountered during the orientation estimation. The difference is that now we can only allow displacement along the shift invariant orientations, as depicted in figure 2.11b. A steered version of the generalized Kuwahara filter immediately suggests itself. The generalized Kuwahara filter should have the same dimensionality as the steered filter. Onedimensional for line-like structures and two-dimensional for plane-like structures. However, an efficient estimation of the parameter and confidence value in the steered displaced windows is not possible, since the orientation can be different between each point in the image. To keep the computation cost low, we are forced to reduce the number of displaced windows. By varying this number, we experimentally verified the trade-off between ”visual
30
Chapter 2 Linear structures
(a)
(b)
Figure 2.11: Reducing noise near a border between two oriented textures, a) orientation adaptive filter, b) edge preserving orientation adaptive filter.
quality” and computation cost. With visual quality we mean the judgment of the human visual system. A reasonable compromise is to use 5 (1D) or 9 (2D) displaced windows evenly spread over the neighborhood. We have applied the steered Kuwahara filter to reduce the noise in the two-dimensional three-domains test image. The signal-to-noise ratio in the image is 9 dB, and the results are shown in figure 2.12d. In each displaced 21*1 pixel window, we measured the Gaussian weighted mean with σ = 5 as the parameter estimate and the variance as the confidence value. The filter is steered using the improved orientation estimate of the GST with σg = 1.4 and σT = 5. The result after controlling the selection process according to eq. 2.21, is shown in figure 2.12e. For comparison we have also shown the result of the normal orientation adaptive filter of the same size, and steered with both the GST orientation and the improved orientation estimate. From figure 2.12 is it clear that the steered Kuwahara filter gives the best result near borders. The created false contours inside the domains are reduced by controlling the selection of the Kuwahara filter on the basis of the GST confidence. A significant speedup could be obtained by using the GST confidence as the confidence value and computing the optimal displacements using a steered version of generalized Kuwahara filter. Applying these displacements to the result of an orientation adaptive filter leads to edge preservation, see figure 2.12f. However, since the confidence value was not computed in the same window as the parameter estimate, a small error in the orientation estimate could cause large errors in filter result. Another remark we should make is that we composed the three-domains image using oriented textures of different scales, to show that the steered filters are not very sensitive to the scale selection. Since the scale separation between the domains in this test image is so clear, a multi-scale analysis would in this case certainly lead to better results. For comparison we have applied the coherence enhancing diffusion filter as described by Weickert in [Wei98], to the three-domains test image as well. This filter is based on the anisotropic diffusion equation ∂I = div(D · ∇I). (2.33) ∂t
2.3 Edge preserving filtering
31
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Figure 2.12: Noise reduction using an orientation adaptive filter steered by the GST (b), and the improved orientation estimate (c). The result of the steered Kuwahara filter (d), and a space-variant mix of (c) and (d) on the basis of GST confidence (e). The faster steered generalized Kuwahara based on the GST confidence and (c). Coherence enhancing diffusion with 5,10, and 20 iterations (g,h,i).
32
Chapter 2 Linear structures
The diffusion is made coherence enhancing by basing the diffusion tensor D on the structure tensor according to D = d1 v1 v1T + d2 v2 v2T , (2.34) with v1 , v2 the eigenvectors of the structure tensor and d1 = α ( d2 = α
d2 = (1 − α)e
(2.35) λ 1 = λ2 −C (λ1 −λ2 )2m
λ1 > λ 2
(2.36)
with λ1 , λ2 the eigenvalues of the structure tensor and α, C, m parameters of the filter. In figure 2.12g,h,i we have shown the results of this filter (α = 0.01, C = 0.001, m = 1) after 5,10 and 20 iterations. These result are less ‘edge preserving’ than the generalized Kuwahara results.
2.4
Application: Automatic fault detection
In this section we will apply the image processing techniques described in this chapter to the automatic detection of faults. The first step in each automation process is to study the way things are done manually. In the case of fault detection this means that we should summarize the methods used by the interpreter. The manual detection of faults consist of two steps. First the fault location is found using visual cues for reflector discontinuities. If there is enough support, the fault is inspected more closely to estimate the fault parameters. This is done by matching seismic sequences from both sides of the fault. We will limit ourself to the automatic localization of faults.
(a)
(b)
Figure 2.13: Two schematic cross-section of faults showing a different vertical displacement of a seismic sequence. A normal fault (a), and a reverse fault (b). The faults are depicted by thick solid lines. The solid arrows indicate the direction of the displacement along the fault, and the dashed arrows show the forces in the subsurface rocks.
Faults are fractures in the subsurface rock caused by tectonic forces. These forces often lead to a vertical displacement (throw) between the rock formations on both sides of
2.4 Application: Automatic fault detection
33
the fracture. In seismic images, faults are recognized by abrupt vertical displacement of seismic sequences along some plane, i.e. the fault plane. Two schematic examples of the appearance of faults in seismic images are depicted in figure 2.13. In practice, it is rare that the reflectivity is identical on either side of the fault. Therefore there can be reflector discontinuities visible in the data, even if there is no measurable throw. The most important observation from an image processing point of view, is that faults cause discontinuities in the reflectors. Flat continuous reflectors can be well described as plane-like linear structures and the value of Cplane from eq.(2.17) will be close to 1. Near a fault, there are also intensity (seismic response) changes perpendicular to the fault plane, causing a decrease in the value of Cplane . We have computed Cplane over a part (1283 voxels) of a 3D seismic image, that contains many faults, and the result is shown in figure 2.14b. From this result we see that the faults are only partly highlighted by Cplane . The fault is not highlighted at locations where there are no clear reflectors present nor at locations where a the reflector from one side of the fault is continued by another reflector at the other side. In three dimensional data it is possible distinguish between random structures and faults using the GST. Since a fault plane is two dimensional, there is still one orientation with only small intensity changes. Namely perpendicular to both the reflector and the fault plane. The value of Cline will therefore be close to 1 near faults. We have summarized the observations in the table below: Random/isotropic: Fault:
λ 1 ≈ λ2 ≈ λ3 λ 1 ≈ λ2 λ3
Cplane ≈ 0 Cline ≈ 0 Cplane ≈ 0 Cline ≈ 1
We have combined Cplane and Cline to create the fault confidence Cf ault , which given by Cf ault = (1 − Cplane )Cline =
2λ2 (λ2 − λ3 ) . (λ1 + λ2 )(λ2 + λ3 )
(2.37)
This new fault confidence measure takes values between 0 and 1, and should give fewer false positives on random structures. A very distinguishing feature of faults is that they have a large vertical extent. Until now we have computed Cplane using isotropic windows. Increasing the isotropic window size to capture the vertical extent, will not yield a better signal to noise ratio1 . Since a fault is a 2D signal inside a 3D window, the amount of noise inside the window grows more rapidly than the amount of signal. By elongating the tensor windows vertically along the faults, we increase the SNR and we exploit the distinguishing vertical extent. Because we only know the orientation of the reflector planes, we steer the elongated tensor windows perpendicular to these planes. This is in many cases a reasonable approximation of the orientation of the fault. We computed Cf ault on the 3D seismic image and the results are shown in figure 2.14c. This new method gives a more continuous fault highlighting with less false positives. 1
Here we define the faults as the signal and everything else as the noise.
34
Chapter 2 Linear structures
(a) xt-slice
(b) Cplane
(c) Cf ault
(d) fault detection
Figure 2.14: Detecting faults by measuring reflector discontinuity in a 3D seismic image. A xt-slice of this image is shown in (a). The result of the GST confidence measure Cplane , computed using isotropic windows, is shown in (b) and the measure C f ault computed using elongated tensor windows is shown in (c). The fault detection using based on (c) is shown in (d).
2.4 Application: Automatic fault detection
35
(a) xy-slice
(b) yt-slice
Figure 2.15: Results of the fault detection algorithm on a 3D seismic image. The fault detections are indicated by black pixels. (a) A time slice (xy), (b) a yt-slice.
36
Chapter 2 Linear structures
Figure 2.16: Result of the fault detection algorithm on a 3D seismic image. The images show a xt-slice of this image, with the fault detections indicated by black pixels.
The final task of our fault detection algorithm is to segment the seismic image into fault and non-fault. This segmentation is based on Cf ault and consists of two steps. First we apply a non-maximum suppression. The candidate fault points are limited to the local maxima of Cf ault , measured in one-dimensional windows steered perpendicular to the fault-planes. This perpendicular orientation is the maximum gradient orientation in the Cf ault image. Next, we take the Cf ault value of these local maxima into account using two threshold values. The reduced Cf ault image is segmented using a low threshold value Tl . The final result is obtained by accepting only those connected segments whose maximum Cf ault value is above a high threshold value Th . The results of our automatic fault detection algorithm are show in figure 2.14d, 2.15 and 2.16. For all the results we have used the threshold values Tl = 0.1 and Th = 0.4. The segmentation is visualized by making the fault points black in the seismic image. To allow a manual verification of the result, we have depicted the original images as well. Now that the fault locations have been found, the estimation of fault parameters could be performed as a subsequent processing step. One could for example perform a matching pursuit between the seismic sequences at both sides of the fault to estimate the throw. This information could then be used to determine the extent of the fault, e.g. by splitting or merging of detected fault planes.
2.5 Application: Structure enhancement
2.5
37
Application: Structure enhancement
A display of a subsurface layer with the same viewpoint as one has when watching the earth’s surface from an airplane, is called a map view. Map views have been used by seismic interpreters for many years, and have become an essential tool for the inspection of lateral changes in the geological structures. Horizons, the earth’s surface in previous geological times, form the natural bases for map views. The manual construction of horizons is a very time consuming task and is therefore automated by horizon-tracking algorithms. Ideally, a horizon-tracker is initiated by picking a single point on the target reflector, and it automatically finds all points in the image that belong to that same reflector. In todays practice, the interpreter has to provide the algorithm with many start points to get an accurate result. Small deviations from the layered structure are often the cause of tracking errors. These deviations can be the result of noise, pre-processing artifacts or subtle geological features that are not important for the structural analysis. Enhancing the structural information by suppressing these small deviations could greatly improve the performance of horizontrackers. Reflectors can be locally described as planar linear structures using the structure tensor. Choosing a large tensor scale σT , makes this description insensitive to small deviations. Reflector surfaces are deformed by many geological processes, e.g tectonic forces, intrusive salt, or erosion. The tensor scale should therefore be chosen small enough, not to notice the deviations from the linear model due to curvature in the global structure, think of synclines and anticlines. In practice, the smallest scale should be coupled to the lowest frequency of seismic wavelet. A rough estimation can be made assuming a low signal frequency of 10 Hz, corresponding to a time window of 100 ms to contain one cycle. If the the wave velocity is 2000 m/s, then the window has a length 200 m. A time sampling of 4 ms gives a window of 25 pixels corresponding to σT = 4. Having estimated the local orientation of the structure using the GST, we can remove the small deviations by orientation adaptive filtering. A steered disc shaped filter (σ1 = σ2 σ3 ) enhances the planar structure. We applied this intrinsically 2D orientation adaptive filter to a 3D seismic image and the result is shown in figure 2.17b. The scale of the filter was σ1 = σ2 = 2 and it was steered parallel to the reflectors using the first eigenvector of the GST (σg = 1, σT = 6). The adaptive filter has clearly reduced the deviations, and the reflectors are now much easier to follow. The filter has also smoothed the reflectors across the faults, possible merging different reflectors. This should avoided and the filter should be made edge preserving, or more accurately fault preserving. Faults can be preserved using the generalized Kuwahara filter combined with the fault confidence Cf ault from the previous section. We applied a one-dimensional version of the generalized Kuwahara filter steered perpendicular to the fault-plane using the second eigenvector of the GST. We used the disc shaped filter as the parameter estimate and Cf ault as the confidence value. The result is shown in figure 2.17c. We removed the false contour by controlling the selection process using eq.2.21. The final result is shown in figure 2.17d
38
Chapter 2 Linear structures
and figure 2.18 for orthogonal slices through the data. These slices are same as the ones shown in the fault detection application.
2.5 Application: Structure enhancement
39
(a) xt-slice
(b)
(c)
(d)
Figure 2.17: The result of the adaptive filter controlled by the improved orientation estimate applied to a seismic image containing a lot of faults (a) is shown in (b). The fault preserving version of (b) is shown in (c). (d) shows a spatially variant mix of (b) and (c) based on the GST confidence.
40
Chapter 2 Linear structures
(a) xy-slice
(b) yt-slice
Figure 2.18: Results of the structure enhancement on a 3D seismic image. (a) A time slice (xy), (b) a yt-slice.
3. Curvilinear structures The GST model presented in chapter 2, describes images locally as a linear structure and yields estimates for orientation and confidence. The previous chapter shows that the GST has a broad range of applications. We have shown that the orientation estimate is stable and suitable for controlling subsequent processing steps. The confidence estimate decreases at image locations where the neighborhood deviates from the linear model. An example of such a neighborhood is the area around a border between two oriented textures. The confidence estimate can therefore be used to these highlight borders. Bending a straight structure also causes a decrease of the confidence value, because the local structure is no longer shift invariant. We have limited the influence of curvature in the previous chapter, by limiting the size of the analysis window, i.e. the tensor scale σT . However, the noise level in the image determines the minimum scale needed to obtain a confident description of a neighborhood. In the case of a curved structure with a high curvature or a high noise level, the linear model does not suffice. Verbeek et al. presented in [VvVvdW98, WVVG01] a curvature corrected confidence estimate for 2D curvilinear structures by extending the structure tensor. A coordinate transform is applied to locally transform the image in such a way that the rotational invariance of the local structure, becomes a translational one. Next, the structure tensor is applied in the new coordinates and transformed back to the original coordinates. This yields expressions for the eigenvalues of the tensor with the curvature κ as a free parameter. Maximization of the confidence measure yields the transformation that ‘closest resembles’ the local structure. From the parameter of the optimal transformation κopt , an estimate of local curvature is obtained. We will extend this method to the curvature corrected confidence estimation of 3D curvilinear structures. In three dimensions there are two types of linear structures, namely plane-like and linelike. These structures are generalizations of an individual plane and an individual line. The curvilinear counterparts of planes and lines are surfaces and curves. Evidently, there are also two types of curvilinear structures in 3D. Plane-like curvilinear structures are the generalization of surfaces, and line-like curvilinear structures are the generalization of individual curves. A plane-like linear structure was defined in the previous chapter as shift invariance along a plane. A stack of straight isophote surfaces therefore has a plane-like linear structure. A stack of curved isophote surfaces, on the other hand, has a plane-like curvilinear structure. To be able to describe plane-like curvilinear structures using the GST, the gradient of the isophote surfaces should be continuous, since the GST analysis is
42
Chapter 3 Curvilinear structures
based on differential geometry. In this chapter we will extend the local straight model of the GST to a local quadratic surface model. The curvature corrected confidence estimate for 2D curvilinear structures will be derived as a special case of quadratic surfaces. The curvature corrected description of line-like curvilinear structures differs significantly and will therefore be presented separately in the next chapter.
3.1 3.1.1
The GST for 3D plane-like curvilinear structures The quadratic surface approximation
In the previous chapter we modeled a local neighborhood of an arbitrary surface S(x) = 0 as a plane. This modeling can be interpreted as a first order polynomial approximation of this surface around a point P on this surface S(x) ≈ b · x + c = 0,
(3.1)
with x the local neighborhood coordinate with P as its origin, and b the unit normal vector (kbk = 1) of the plane. Apart from translation, a plane can be described by two angles. In this chapter we will extend this description with two principal curvatures. This can be written as the addition of a quadratic term, making it a second order polynomial approximation S(x) ≈ xt · A · x + b · x + c = 0, (3.2) with A a symmetric matrix with two non-zero eigenvalues [FIS89]. In essence we limit ourselves to the quadratic surfaces that, after rotation and translation, can be written as S(x) ≈ λ1 s21 + λ2 s22 − s3 = 0,
(3.3)
with s1 , s2 , s3 coordinates in the frame spanned by the eigenvectors of A. Now, we want to determine the relation between the two eigenvalues of A and the two principal curvatures (κ1 , κ2 ) of the surface S at point P (s1 = s2 = s3 = 0). The principal curvature κi (i ∈ 1, 2) can be found by computing the curvature of S(si ) = λi s2i using the definition [BSMM00] df (x) f 00 0 , f = (3.4) κ= [1 + (f 0 )2 ]3/2 dx This yields κi = 2λi , and we can write 1 1 S(x) ≈ κ1 s21 + κ2 s22 − s3 = 0. 2 2
(3.5)
Thus, the quadratic surface is either an elliptic paraboloid in case the signs of κ1 , κ2 are equal, or a hyperbolic paraboloid (saddle), if κ1 , κ2 have opposite sign.
3.1 The GST for 3D plane-like curvilinear structures
43
We define a quadratic surface structure as shift invariant along a quadratic surface and not shift invariant along the normal on each point of this surface. A quadratic surface can be described by two curvatures (κ1 , κ2 ) and a local frame. This frame is oriented along the normal vector and the two vectors that correspond to principal curvatures, i.e. the frame for which A is diagonal. Finding the optimal quadratic surface model therefore requires the ‘simultaneous optimization of five parameters, namely κ1 , κ2 , and three angles to determine the local frame. A significant speedup of the optimization problem can be obtained by utilizing the local frame Q given by the eigenstructure of the GST. 2 2 fx fx fy fx fz fu 0 0 [T]uvw = QT · fx fy fy2 fy fz · Q = 0 fv2 0 (3.6) 2 2 fx fz fy fz fz 0 0 fw
The definition of Q is given in appendix B of this chapter. In eq.(3.6), we labeled the gauge-coordinates of the GST: ‘uvw’. The coordinate u corresponds to dominant gradient orientation (the normal) and v, w to the directions of respectively κ1 , κ2 , with κ1 > κ2 . This uvw-frame is schematically depicted in figure 3.1a. Without curvature, fv2 = λ2 and fw2 = λ3 are determined by noise. Increasing curvature (κ1 , κ2 ) causes an increase of the eigenvalues (resp. λ2 , λ3 ). Utilizing the gauge-coordinates of the GST, the optimization of the model reduces to optimization of two independent parameters κ1 and κ2 . u v
w
u’ v’
(a) uvw
w’
(b) u0 v 0 w0
Figure 3.1: Deforming quadratic surfaces to planes using a coordinate transform.
3.1.2
The quadratic GST for 3D surfaces
The first step toward a curvature corrected GST, is to define the coordinate transform uvw → u0 v 0 w 0 that deforms a curvilinear surface into a plane, see fig 3.1. Using the quadratic surface approximation, the transform and its inverse are given by u0 = u − 21 κ1 v 2 − 21 κ2 w 2 v0 = v w0 = w
u = u0 + 12 κ1 v 0 2 + 12 κ2 w 0 2 v = v0 w = w0
44
Chapter 3 Curvilinear structures
This approximation is valid for small values of u, v or small values of κ1 , κ2 . The test results will show the practical limitations of this approximation. The inverse-transformation is used to express the derivatives in the u0 v 0 w 0 -coordinates as a function of the derivatives in the uvw-coordinates. f u0 = u u0 f u + v u0 f v + w u0 f w = f u fv0 = uv0 fu + vv0 fv + wv0 fw = κ1 vfu + fv fw0 = uw0 fu + vw0 fv + ww0 fw = κ2 wfu + fw
(3.7)
We now consider curved surfaces in the uvw-space with as shown in figure 3.1a. By applying the traditional GST method to the u0 v 0 w 0 -space for arbitrary κ1 , κ2 , we get the gradient structure tensor for quadratic surfaces fu20 f u0 f v 0 f u0 f w 0 TQS = fu0 fv0 (3.8) fv20 fv 0 fw 0 2 f u0 f w 0 f v 0 f w 0 f w 0
Using eq.(3.7) we can express TQS in the uvw-coordinates. fu20 = fu2 fv20 = fv2 + κ21 a + 2κ1 b fw2 0 = fw2 + κ22 c + 2κ2 d
(3.9)
f u0 f w 0 = f v 0 f w 0 = f u0 f v 0 = 0 with the abbreviations a ≡ v 2 fu2 , b ≡ vfu fv
c ≡ w 2 fu2 , d ≡ wfu fw .
(3.10)
The off-diagonal elements of TQS are all equal to zero, due to the symmetries in the model. The derivation and the details of the symmetry considerations are given in appendix A. Since the matrix of [TQS ]uvw is diagonal, the eigenvalues of TQS are equal to the diagonal elements of this matrix λ1 = fu2 , λ2 (κ1 ) = fv2 + κ21 a + 2κ1 b , λ3 (κ2 ) = fw2 + κ22 c + 2κ2 d.
(3.11)
For a plane-like linear structure the eigenvalues relate as λ1 λ2 ≈ λ3 . The optimal resemblance to the quadratic surface, can therefore be found by minimizing λ2 and λ3 , with respect to κ1 , κ2 . The curvatures κ1,min , κ2,min that minimize λ2 ,λ3 yield an estimate of local curvature and the are given by κˆ1 = κ1,min =
−b −d , κˆ2 = κ2,min = . a c
(3.12)
3.1 The GST for 3D plane-like curvilinear structures
45
Substitution of these curvatures in eq.(3.11), yields the following expressions λ1 = fu2 , λ2 = fv2 −
b2 d2 , λ3 = fw2 − . a c
(3.13)
These curvature corrected eigenvalues can be combined to create a curvature corrected confidence measure. Since λ2 , λ3 of the GST are determined by the curvature and the noise level, it is not guaranteed that λ2 > λ3 , after removing the influence of curvature. This should be checked explicitly, to avoid an overestimation of the confidence value C=
λ1 − λ 2 . λ1 + λ 2
(3.14)
Note that for the limit κ1 , κ2 → 0, the equations for λ1 , λ2 reduce to those of the traditional GST. The computation of the structure tensor for quadratic surfaces TQS , consist of two steps. First, the GST is computed to obtain the gauge-coordinates uvw, and next the extra terms (a, b, c, d) of eq.(3.10), need to be computed. The direct computation of (a, b, c, d) results in a spatially variant operation, since the gauge-coordinates change over the image. For large filter sizes, this can lead to a high computational demand. Appendix B shows how these terms can be computed efficiently using linear combinations of convolutions.
3.1.3
Experimental tests and results
The structure tensor for quadratic surfaces TQS , described above models a local neighborhood as a linear structure after a quadratic transformation. The parameters of the optimal transformation κˆ1 , κˆ2 yield an estimate of the local curvature. The resemblance between the local model and the image data is given by the confidence estimate C. The value of the confidence estimate can be decreased by two factors; noise and a deviation of the local structure from the supposed model. In this section we will experimentally test the structure tensor for quadratic surfaces TQS , by means of several measurements on test images. Quadratic structure We will start with a test of the correctness of the theory and its implementation using the ‘ideal’ test image It1 given by 1 1 0 It1 (x, y, z; κ1 , κ2 ) = It1 (u, v, w; κ1, κ2 ) = u − κ1 v 2 − κ2 w 2 , 2 2
(3.15)
where the transformation of the xyz to the uvw-coordinates is given by an arbitrary rotation. The neighborhood around the origin (x, y, z) = (0, 0, 0) of the test image has a
46
Chapter 3 Curvilinear structures
quadratic structure. The tensor TQS should therefore give an exact description of the image structure at the origin. We will experimentally verify this proposition in two steps. First we verify that the first eigenvector of the structure tensor is normal to the shift invariant surfaces and that the second eigenvector is aligned with the orientation of maximal curvature in the tangent plane. At the origin of the test image It1 this means that the first eigenvector should be parallel to the u-axis. The second eigenvector should be parallel to the v-axis if κ1 > κ2 and parallel to the w-axis if κ2 > κ1 . We found that the eigenvectors of the GST are exactly parallel to the uvw-axes for all rotations and curvatures 1 . There is however a value of the curvature κi > κmax , where first eigenvector of the GST swaps with the second eigenvector. The value of κmax depends on the tensor scale σT , and we experimentally verified that 1 κmax (σT ) = . (3.16) σT For example, we found the maximum curvature κmax = 0.2 for σT = 5 and κmax = 0.5 for σT = 2. Next we applied TQS to the test image It1 to verify the exactness of the description at the origin. We found for the curvature and confidence estimates that C=1, κ ˆ i = κi , i ∈ {1, 2},
(3.17)
within the precision of the measurement; 10−7 of 32-bit floating point numbers. Therefore we have shown that theory and the implementation of TQS are correct. Spherical structure To test the robustness with respect to noise and with respect to deviations of the structure from the local model, we use a second test image of noisy concentric spherical surfaces defined by p 2π It2 (x, y, z; p, σn ) = It2 (r; p, σn ) = cos r + N (0, σn2 ) , r = x2 + y 2 + z 2 , (3.18) p where σn is the standard deviation of the normally distributed noise N . First we applied TQS with (σg = 1, σT = 5) to the test image It2 (r; 8, 0). This noise-free version allows us to determine the estimation errors due to the deviation of the spherical structure from the quadratic model. We measured the resulting estimates for curvature and confidence according to eqs.(3.12) and (3.14) as a function of r, and the results are shown in figure 3.3a and b. In figure 3.3a we compared the confidence estimate (QGST) with the confidence estimate of the structure tensor for linear structures (GST). We have also plotted in this figure the improvement in the confidence estimate
1
∆C = (C)QGST − (C)GST ,
All curvatures for which the corresponding surfaces can be properly sampled.
(3.19)
3.1 The GST for 3D plane-like curvilinear structures
47
due to the curvature correction. For the evaluation of the curvature estimate we have plotted the ‘true κ’ in figure 3.3a. This true curvature is the curvature of the surface of a sphere along some arbitrary tangent, and is given by 1 κ= , r
(3.20)
with r the radius. The estimation error is therefore the difference between κ and the estimated value. When we analyze the estimation error in figure 3.3a, we can distinguish two parts. For the first part (r > 15), the error increases slowly for decreasing r from 10−5 up to 10−3 when r approaches 15. In the second part (r < 15), the estimation error increases rapidly for decreasing r and becomes very larger. The confidence estimate shows approximately the same behavior; close to 1 for r > 15, and quickly decreasing to 0 for r < 15. We have checked this behavior for several values of σT , and we found that this rapid decrease in confidence occurs for r < 3σT .
(a) r = 0
(b) r = 15
(c) r = 31
Figure 3.2: Three Gaussian windowed (σ=5, 31 3 pixels) neighborhoods of the noise-free test image around points at r=0,15,31. The corresponding confidence estimates using the quadratic model are C=0.00, 0.96, 1.00.
To give an impression of the underlying data, we selected three neighborhoods of 313 pixels from the noise-free test image It2 (r; 8, 0) around points at r=0,15,31. The neighborhoods are multiplied with a Gaussian window (σ = 5) to show the input-data for the local analysis with the quadratic structure tensor TQS . A two-dimensional cross-section of each of these neighborhoods, is shown in figure 3.2. At r=0 the neighborhood consists of concentric spheres, that constitute an isotropic structure and the confidence value is C=0. At r=15 the center of the spheres is at the edge of the neighborhood and the structure consists of half spheres. The quadratic model is able to describe this neighborhood with a confidence C=0.96. The value of the confidence estimate of the third neighborhood is C=1.00. The deviations of the local structures in It2 from the quadratic model, consist of two parts. The exact description of a spherical surface requires a higher order polynomial. Second, the quadratic model assumes a constant curvature along the u-coordinate. The u-axis is spanned by the first eigenvector of the GST. If we describe the spherical test image with spherical coordinates It2 (r, φ, θ), than the normal vector at each point in the image is parallel to the r-axis. The u-axis is parallel to the r-axis, since the first eigenvector of the GST is parallel to the normal vector Using eq.(3.20) we can write the curvature as a
48
Chapter 3 Curvilinear structures
0.12
1
0.1
0.8
True κ κ1
0.6
Curvature
Confidence
0.08
QGST GST QGST−GST
0.4
0.04
0.2
0
0.06
0.02
0
20
40
60
80
0
100
0
20
40
Radius
80
100
(b) Curvature, noise-free
(a) Confidence, noise-free
20 dB
1
0.9
Confidence
60 Radius
10 dB
0.8
0.7
3 dB QGST GST
0.6
0.5
15
20
25
30
35
40
45
50
Radius (c) Confidence, noise dependence
Figure 3.3: Experimental results of estimation of curvature and confidence (σ T = 5) on a noisefree spherical structure using the quadratic model (a,b). Experimental results of the test of robustness with respect to noise (3,10,20 dB) of the confidence estimate (c). The measurements in (c) are averaged over 20 different noise realizations, and the error-bars denote the corresponding standard deviation.
3.1 The GST for 3D plane-like curvilinear structures
49
0.07 True κ 20 dB 10 dB 3dB
0.06
κ1
0.05
0.04
0.03
0.02
0.01
15
20
25
30
35
40
45
50
Radius
(a) Curvature κˆ1 0.07 True κ 20 dB 10 dB 3dB
0.06
κ2
0.05
0.04
0.03
0.02
0.01
15
20
25
30
35
40
45
50
Radius
(b) Curvature κˆ2
Figure 3.4: Experimental results of the test of the robustness with respect to noise of the curvature estimates κˆ1 (a) and κˆ2 (b). The measurements are performed at scale (σ T = 5) for 3 noise levels, 3,10,20 dB. The measurements are averaged over 20 different noise realizations and the error-bars denote the corresponding standard deviation.
50
Chapter 3 Curvilinear structures
function of u
1 . (3.21) |u| The curvature is clearly not constant along the u-coordinate. Since the GST is a quadratic form, the sign of the eigenvectors is not defined by the image structure, and therefore we have written the curvature as a function of the absolute value of u. For use of the eigenvectors in subsequent processing steps, their signs are chosen in a consistent manner. κ(|u|) =
To explain why the curvature estimate approaches zero when the analysis window approaches the center of the test image, we consider the neighborhood (r = 0) in figure 3.2a. One half of this neighborhood (u > 0) contributes to a positive curvature and the other half (u < 0) to a negative curvature, resulting in a curvature estimate κ ˆ 1,2 = 0. For radii r smaller than the effective radius of the analysis window, approximately 3σ for a Gaussian window, the structure of spherical surfaces can no longer be straightened by a coordinate transform of the form u0 = u − F (v, w) , v 0 = v , w 0 = w, (3.22) with F (v, w) and arbitrary function. Even the optimally transformed structure becomes isotropic as r → 0.
Next, we tested the robustness with respect to noise, by applying the tensor TQS at a constant scale (σg = 1, σT = 5) to the test image It2 (r; 8, σn ) at three different noise levels σn2 = {0.04, 0.4, 2}, corresponding to the SNR values 20,10,3 dB. Again we measured the curvature and confidence estimate as a function of the radius (15 < r < 50). The measurements are averaged over 20 different noise realizations, and the results are shown in figure 3.3b and 3.4. The error-bars in these figures denote the corresponding standard deviation. From figure 3.3b it is clear that the confidence value decreases with an increasing noise level, for both the quadratic (QGST) and the linear (GST) model. If we consider a neighborhood where f = p + n, with p the pattern and n the noise, then we get the following noise dependency of the GST confidence (C(n))GST =
fu2 − fv2 p2u − p2v = , fu2 + fv2 p2u + p2v + 2n02
(3.23)
with n2u = n2v = n02 and pu nu = pv nv = 0. It is clear that as the standard deviation of the noise becomes larger the confidence measure becomes smaller. The noise dependency of confidence estimate of the QGST can be found using eq.(3.13). (C(n))QGST = with
2
p2u − p2v + K(n) p2u + p2v + 2n02 − K(n)
(3.24)
2
vfu fv v(pu pv + pu nv + pv nu + nu nv ) vpu pv 2 = . K(n) = = v 2 fv2 v 2 p2v + v 2 n2v v 2 (p2v + 2pv nv + n2v )
(3.25)
3.2 The GST for 2D curvilinear structures
51
Equation (3.25) shows that K(n) decreases for an increasing noise level. Therefore it is clear from eq.(3.24) that (C(n))QGST decreases if n2 increases. Furthermore the improvement of the confidence estimate due to the curvature correction ∆C, decreases as well. This is consistent with the results shown in figure 3.3b. The experimental results in figure 3.4 show that the magnitude of the curvature estimate becomes smaller with an increasing noise level. This under-estimation of the curvature can be explained in a similar fashion as the decrease in the confidence estimate due to noise. If we substitute f = p + n in eq.(3.12), we obtain κ1 (n) =
vfu fv vpu pv wfu fv wpu pv = , κ2 (n) = = , v 2 fv2 v 2 p2v + v 2 n2v w 2 fv2 w 2 p2v + w 2 n2v
(3.26)
and we see that both κ1 (n) and κ2 (n) decrease if n2 increases.
Discussion Summarizing the results, we can say that the estimators for curvature and confidence give error-free results on the ‘ideal’ test image, and we can therefore conclude that the theory and its implementation are correct. By applying the quadratic structure tensor to a second test image that has a structure of spherical surfaces, we showed that deviations from the quadratic model only result in significant errors, if the radius of curvature is smaller than the effective radius of the analysis window. We experimentally determined the noise dependency of the estimators, and found that both the confidence estimate and the improvement in the confidence estimate due to the curvature correction, decrease for an increasing noise level. For signal-to-noise ratios smaller than 10 dB, the curvature estimates becomes significantly biased, i.e. they give an under-estimation.
3.2
The GST for 2D curvilinear structures
The two-dimensional curvature corrected GST can be derived from the structure tensor for quadratic surfaces TQS , in a straight forward manner by considering a 2D sub-section of the uvw-space. Convenient sub-sections for the derivation of the 2D estimators are the uv-plane or the uw-plane. If we consider the uv-plane, then the transformation that straightens the curved structure in this plane is the coordinate transform of quadratic surface model with κ2 = 0, w = 0. The inverse transformation can be derived similarly, since κ1 and κ2 of the quadratic surface model are independent. u0 = u − 21 κ1 v 2 v0 = v
u = u0 + 21 κ1 v 0 2 v = v0
52
Chapter 3 Curvilinear structures
Applying the traditional GST to the u0 v 0 -space, yields the GST for 2D curvilinear structures " # 2 f u0 f u0 f v 0 Tcurv2D = . (3.27) f u0 f v 0 fv20 Expressions for the tensor element in the uv-space are found using the inverse transformation and they are given in eq.(3.9). The optimal resemblance between the model and the local image structure is found by maximizing the confidence measure C=
λ1 − λ 2 , λ1 + λ 2
(3.28)
and thus by minimizing λ2 (κ1 ) with respect to κ1 . This minimization yields the following results b2 −b λ1 = fu2 , λ2 = fv2 − , κˆ1 = κ1,min = (3.29) a a with the abbreviations (3.30) a ≡ v 2 fu2 , b ≡ vfu fv Applying the GST to the u0 w 0 -plane gives the same results, replacing κ1 ↔ κ2 and λ2 ↔ λ3 . The direct derivation of the curvature corrected structure tensor for 2D curvilinear structure including test results and an application to finger-print recognition is presented in [VvVvdW98]. The implementation shown in appendix B is also applicable to the two dimensional terms (a, b). In [WVVG01], Weijer et al. derive two-dimensional curvature corrected confidence measures by applying normalized curvilinear models to the gradient vector field. This allows for the quantitative comparison of different models, i.e. coordinate transformations. They compared the quadratic model (second order polynomial) with the circular model, and show that the latter gives better results on circular oriented textures for small radii of curvature. They also present an application to interference patterns of a vibrating plate and to images of growth rings of trees.
3.3
Curvature adaptive filtering
In the previous chapter, we showed that noise in oriented texture domains can be reduced by anisotropic low-pass filtering along the shift invariant orientations, i.e. by orientation adaptive filtering (OAF). The motivation for using adaptive filter windows is that increasing the window size isotropically does not yield a higher SNR, but increasing the size of an adaptive window that matches the local structure does. We have already encountered a limitation of elongated OAF windows. They are not able to match a curved oriented texture, causing an increasing filter error with increasing curvature. Using the parameter estimates of the curvature corrected structure tensor, we can make the filter windows both
3.3 Curvature adaptive filtering
53
curvature and orientation adaptive (COAF). In this section we present a two-dimensional curvature adaptive filter, and its application to noise reduction in curved oriented textures. Although we present the filtering method in 2D, we have applied it in 3D as well. The COA filter windows can be created in the same way as we created the OA filter windows in the previous chapter. The appropriate scale σ is selected by applying an isotropic filter, and the filter is elongated by sampling the isotropic responses at regular intervals along a curved line, see figure 3.5. This can be implemented by locally transforming the data using the quadratic coordinate transform of Tcurv2D from the xy to the u0 v 0 -coordinates, and sampling of the filter along the u0 -axis. The weighted average of the samples yields the adaptive filter result. Below we will use a Gaussian weighting with a width σ1 .
(a)
(b)
Figure 3.5: (a) Creating a curved elongated filter by combining isotropic filters. (b) Adapting the curved elongated filter to match the local structure.
To demonstrate the COA filter, we will use it to reduce the noise in a test image It of a curved oriented texture. p 2πr It (x, y; p, σn ) = It (r) = cos + N(0, σn2 ) , r = x2 + y 2 , (3.31) p where N is normally distributed noise. A realization of It (x, y; p, σn ) of 128*128 pixels with (p = 8, σn2 = 1/2) is shown in figure 3.6a. The first step in reducing the noise in curved oriented textures is to apply the tensor Tcurv2D . The resulting parameter estimates on It at scale (σg = 1, σT = 5), are shown in figures 3.6d-f. Since the tensor is based on a quadratic gradient model, the model error becomes large on circular structures with a radius R ≤ 3σT . This is reflected by the lower confidence values in the center of figure 3.6f. The sign of the curvature estimate in figure 3.6e is defined relative to the gaugecoordinates uv of the structure tensor, and is therefore only meaningful in combination with the orientation estimate. Note that the sign of the curvature estimate changes at the same locations as the ‘jumps’ in the orientation estimate occur. The improvement in the confidence estimate ∆C given by ∆C = (C)curv2D − (C)GST ,
(3.32)
is shown in figure 3.6g. From this result it is clear that the improvement has the highest values for small radii. The maximum improvement in confidence for the noise-free case, occurs at radius R = 8.
54
Chapter 3 Curvilinear structures
(a) Test image
(b) COAF
(c) OAF
(d) Orientation
(e) Curvature
(f) Confidence
(g) ∆C
Figure 3.6: Using adaptive filters (b),(c) to improve the SNR of curved oriented textures (a). The parameter estimates derived from T curv2D are shown in (d)-(f), and the improvement in the confidence estimate by curvature correction ∆C is shown in (g).
3.4 Appendix A: Symmetry
55
Using the curvature and the orientation estimate of the tensor Tcurv2D , we can control a COA filter. The result of a Gaussian COA filter with σ1 = 5 is shown in figure 3.6b. For visual comparison, we have shown the result of a Gaussian OA filter (σ1 = 5) in figure 3.6c. A quantitative comparison of the filter results can be made by measuring the signal-to-noise ratios of the filtered test images. The noise in the test image is white and additive, and therefore we estimate the standard deviation of the noise σn of a filtered image f ilter(It ) by σ ˆn = STDDEV[f ilter(It (x, y; p, σn )) − It (x, y; p, 0)]. (3.33) Since the intensities of the noise-free image are bounded between -1 and 1, we use the following definition for the signal-to-noise ratio max(I) − min(I) SN R = 20 log10 dB, (3.34) σn with I the noise-free image. The resulting global SNR’s of the test image after reducing the noise using different filter types are given in table 3.1. Filter Ideal COAF OAF Isotropic None
σn SNR (dB) 0.147 22.7 0.162 21.8 0.238 18.5 0.559 11.1 0.710 9.0
Table 3.1: A quantitative comparison of adaptive filters for noise reduction in a circular oriented texture.
All the filters of table 5.1 use a Gaussian window. The SNR of √ the test image without filtering (None), is 9 dB. The isotropic filter is applied at scale σ = 5. This scale is chosen such that the area of the isotropic window is equal to the area of the windows of the adaptive filters. we applied the ‘ideal’ filter to get an estimation of the maximal attainable SNR using a local filter. The shape ideal filter window matches the shape local structure, which is in this case circular. As the ideal filter we applied a circular shaped Gaussian filter at the same scale as the adaptive filters, controlled by the true local parameter values of the noise-free test image. The SNR’s of COAF and OAF are computed from the results shown in figure 3.6. The visible improvement of the COAF over the OAF is supported by the quantitative results in table 5.1.
3.4
Appendix A: Symmetry
The averaging of a tensor component Tij in 3D, can be written as a triple integral of the component times a symmetric window function w Z Z Z Tij (u, v, w) = w(u, v, w)Tij (u, v, w)dwdvdu. (3.35)
56
Chapter 3 Curvilinear structures
The integration of an anti-symmetric function fasym (u, v, w) ≡ −fasym (u, −v, w), equals zero Z Z Z fasym (u, v, w)dwdvdu = 0. (3.36)
Writing an anti-symmetric function as
fasym (x) = sign(x)fsym (x),
(3.37)
it is easy to verify that the multiplication of (anti-)symmetric functions leads to the following rules 0 0 0 fasym fsym = fasym , fsym fsym = fsym , fasym fasym = fsym .
(3.38)
The gradients in the quadratic surface model have the following symmetries with respect to the uvw-coordinates fu (u, v, w) = fu (u, −v, w) , fu (u, v, w) = fu (u, v, −w) fv (u, v, w) = −fv (u, −v, w) , fv (u, v, w) = fv (u, v, −w) fw (u, v, w) = fw (u, −v, w) , fw (u, v, w) = −fw (u, v, −w)
(3.39)
The combination of eqs.(3.39), (3.36), and (3.38) yields the following results fu0 fv0 = κ1 vfu2 fv + fu fv = 0 fu0 fw0 = κ2 wfu2 fw + fu fw = 0
(3.40)
fv0 fw0 = κ1 κ2 vwfu2 fv fw + κ1 vfu fv fw + κ2 wfu fv fw + fv fw = 0
3.5
Appendix B: Implementation
The elements of the three-dimensional quadratic structure tensor abfc fd , afc fd , a, b, c, d ∈ {u, v, w}
(3.41)
can be implemented using convolutions, by transforming them from uvw-coordinates back to xyz-coordinates before averaging. The terms in eq.(3.41) consist of two components, a (quadratic) coordinate and a structure tensor component. For an arbitrary coordinate system (x1 , x2 , · · · , xn ), we can consider the quadratic coordinates xi xj as elements of a tensor of order 2. The definition of tensors and their elementary operations are given in [BSMM00]. The elements of the second order tensor X, used to describe the quadratic coordinates are given by Xij = xi xj . (3.42) Equation (3.42) is the definition of the dyadic product, a special case of the tensor product that maps two tensors of order 1 to a tensor or order 2. Therefore we can define the quadratic coordinate X using the tensor product by X = xx.
(3.43)
3.5 Appendix B: Implementation
57
Since we will only use the tensors under a coordinate transformation, there is a practical identity between a vector and a tensor of order 1, and between a matrix and a tensor of order 2, see [Gol80]. The coordinate tensors x, X and the gradient tensor T can therefore be transformed to a different coordinate system [FIS89] using v = QT · x
V = QT · X · Q T
[T]uvw = Q · T · Q, where the dot denotes the standard matrix product and with the definitions x u T ≡ [T]xyz , X = xx , V = vv , x = y , v = v . z w The orthonormal transformation matrix Q is given by xu xv xw Q = [eu ev ew ] , eu = yu , ev = yv , ew = yw , zu zv zw
(3.44) (3.45) (3.46)
(3.47)
(3.48)
with (eu , ev , ew ) the eigenvectors of the traditional GST.
Now we are ready to transform the quadratic structure tensor elements to xyz-coordinates. We start with the terms that contain quadratic coordinates abfc fd = (ea · X · eb )(ec · T · ed )
(3.49)
= ea · (X(ec · T · ed )) · eb
(3.50)
= ea · (XT ◦ ec ed ) · eb , a, b, c, d ∈ {u, v, w}
(3.52)
= ea · (X(T ◦ ec ed )) · eb
(3.51)
This results in a linear combination of the 36 terms ijfk fl , i, j, k, l ∈ {x, y, z}.
(3.53)
In eq.(3.51) we introduced the operator (◦) that performs the sum over two indices after an element-wise multiplication, thereby reducing the order n of the tensor by two (n − 2). It maps the two second order tensors T, ec ed to a scalar. If Tij are the elements of T and Eij are the elements of ec ed then the scalar s = T ◦ ec ed is given by the sum of the element-wise multiplication X s= Tij Eij (3.54) i,j
For clarity we give a two-dimensional example: a11 a12 b11 b12 ◦ = a11 b11 + a12 b12 + a21 b21 + a22 b22 . a21 a22 b21 b22
(3.55)
58
Chapter 3 Curvilinear structures
In eq.(3.52) we create the fourth order tensor XT, that contains all the elements of eq. (3.53). A fourth order tensor in 3D has 34 = 81 elements, but since X, T are both symmetrical only 36 are unique. If we define the elements of this tensor Mijkl = Tij Xkl , then the elements Lkl of the second order tensor XT ◦ ec ed are given by X Lkl = Mijkl Eij (3.56) i,j
The transformation of the linear coordinate terms can be performed similarly afc fd = (x · ea )(ec · T · ed ) = ea · (x(ec · T · ed ))
= ea · (x(T ◦ ec ed ))
(3.57)
= ea · (xT ◦ ec ed ) , a, c, d ∈ {u, v, w},
and results in the linear combination of the 18 unique terms ifk fl , i, k, l ∈ {x, y, z},
(3.58)
of the third order tensor xT. Writing the transformation of the quadratic structure tensor elements in tensor form, has a number of advantages over the straight forward expansion. It is a compact representation. The coefficients for the linear combination of the terms in eq.(3.53) and (3.58) do not have to be computed individually, but can be computed by applying the tensor product and the matrix product to the eigenvectors of the GST. The most important advantage, however, is that eq.(3.52) and (3.57) hold for the twodimensional case as well. The only difference is that a, b, c, d ∈ {u, v}. The averaging of the terms in eqs.(3.53) and (3.58), is linear and shift invariant, since the xyz-coordinates are the same for each point in the image. The computation of these terms can now be written as the following convolutions Tkl ⊗ K(x) = Tkl ⊗ (ijG(x; σT )) , i, j, k, l ∈ {x, y, z},
(3.59)
where G is a Gaussian at scale σT . Since the convolution kernel K is xyz-separable, the convolution can be perform efficiently in the spatial domain using FIR-filters. The convolution can also be performed in the Fourier domain using a fast Fourier-transform, for example the FFTW described in [FJ97] or see http://www.fftw.org. For a comparison of the computational speed of convolutions performed with FIR-filters and convolutions using a FFT-algorithm see [YGV98]. We finish with a rough estimation of the difference in computation demand in 3D, between the implementation using spatially variant operations and the method described above. For the computation of the four extra terms in eq.(3.10) of the structure tensor for quadratic surfaces, four spatially variant operations are needed, or 36+18=52 spatially invariant convolutions. We limit our estimation to the FIR implementation of the convolutions.
3.5 Appendix B: Implementation
59
Furthermore, we assume that the computational demands of the coordinate transformation is equal in both methods. A xyz-separable convolution kernel with a width of N pixels, requires approximately 3N multiplications and additions. Where as the non-separable version requires approximately N 3 of these operations. If we combine these complexities with the number of convolutions needed (4N 3 = 156N ), then the filter width that gives an equal number of operations for both methods is N ≈ 6. Thus for a filter width N > 6, approx. σ > 1 for a Gaussian window, the spatially variant implementation uses approximately 4N 3 − 156N more operations than the method presented in this appendix.
4. Line-like curvilinear structures The Gradient Structure Tensor (GST) models an image locally as a linear structure. In three dimensions there are two types of linear structures: plane-like and line-like. Curvature in the local structure causes a decrease of the confidence value of the GST. In the previous chapter we have shown that a curvature corrected confidence estimate can be obtained, by extending the local linear model to a local quadratic model. The structure tensor for quadratic surfaces was presented and tested. In this chapter we will extend the structure tensor to a curvature corrected descriptor of line-like curvilinear structures. In chapter 2, we defined a line-like linear structure as shift invariant along a single orientation. The natural generalization to line-like curvilinear structures is shift invariance along a space curve. The isophotes or level-sets of the local image form parallel space curves. Examples of images that contain these structures are 3D medical images of blood-vessels, such as CT or MRI. Another example is time series of 2D images. Following an object in the 2D image through time yields a space curve, the curvature of which is related to the acceleration of the object. Channels in seismic images can also manifest themselves as space curves. For the local analysis of these structures, we will locally approximate a space curve with a quadratic curve. This curve will be straightened with a coordinate transform, and the GST will be applied to this straightened structure. Transforming the elements of the structure tensor back to original coordinates and maximizing the confidence value, yields an estimate for the curvature of the space curve and a curvature corrected confidence estimate. We will extensively test the new estimators by applying them to images of various space curves. We will conclude this chapter by applying the curvature corrected confidence estimate to the detection of channels. This chapter is based on the previously published article from the author [BVV01].
4.1 4.1.1
The GST for 3D line-like curvilinear structures The quadratic curve approximation
In chapter 2, we implicitly modeled a local neighborhood of an arbitrary space curve c, by a straight line. This modeling can be interpreted as a first order polynomial approximation of this curve around point p on this curve. The approximation of curve c can be expressed
62
Chapter 4 Line-like curvilinear structures
in parametric form by c(t) ≈ bt + p,
(4.1)
with b the tangent vector and t the parameter of c. In this chapter we will extend this approximation with a quadratic term, making it a second order polynomial approximation. This is equivalent to the local modeling of c with a quadratic curve, and can be written in parametric form as c(t) ≈ at2 + bt + p, (4.2) where a contains the quadratic coefficients. There exists a coordinate system for each quadratic curve such that only one element of a is non-zero. Writing the curve in this system gives c(t) ≈ e1 a1 t2 + e2 t, (4.3) with ei the ith basis vector of the coordinate system. This local coordinate system in point p corresponds to the Frenet-frame [Spi79] of the quadratic curve, where e2 is the tangent and e1 the normal. Using the arc length parameterized definitions we obtain
dc(s) ds dc(t) dt
dc = , e1 = c0 (s) ≡ = dt dt ds dt ds c00 (s) , e3 = e 1 × e 2 e2 = 00 kc (s)k
(4.4) (4.5)
The vector e3 is used for the 3D analysis and is called the bi-normal. This direction is needed for the computation of torsion τ , and is therefore not required for the quadratic curve approximation.
de de dt de
1 1 3 κ ≡ kc00 (s)k = = −τ e2 (4.6)
=
, ds dt ds ds The relation between a1 and the curvature κ in point p on the curve is found using the definition of curvature for parameterized curves. k2 =
[(2a1 t)2 + 1](2a1 )2 − [(2a1 )2 t]2 t=0 −−→ k = 2a1 [(2a1 t)2 + 1]3
(4.7)
Using this relation we can write 1 c(t) ≈ e1 κt2 + e2 t. 2
(4.8)
A quadratic curve can be determined by one curvature κ and the Frenet-frame. The optimization of the quadratic curve model therefore consists of the simultaneous optimization of four parameters. If the structure tensor is applied to a noise-free image of a straight line, then the eigenvalues of the tensor for a point on the line relate as λ1 = λ2 > λ3 = 0, and the third eigenvector is aligned along the tangent. Due to the degeneracy, there is still one degree of freedom in
4.1 The GST for 3D line-like curvilinear structures
63
the first two eigenvectors. Bending the line increases the third eigenvalue. Assuming that the total gradient energy is constant we have λ1 > λ02 > λ03 > 0.
(4.9)
The second eigenvector is now aligned along the normal, and the first eigenvector along the bi-normal. However, if the cross-section of the line is not circular symmetric, then we have λ1 > λ2 , even without bending. In this case it is not guaranteed that the second eigenvector is parallel with the normal, and therefore we can only use the third eigenvector as an estimate of the tangent. Still, if we use the third eigenvalue of the GST as the tangent, we reduce the optimization of the quadratic curve model to the simultaneous optimization of two parameters. Namely, the curvature and an angle to determine the normal. In the next section we will use an alternative description of the quadratic curve. We use the gauge-coordinates of the GST ‘uvw’ defined by the local frame of the GST Q, 2 2 fx fx fy fx fz fu 0 0 (4.10) [T]uvw = QT · fx fy fy2 fy fz · Q = 0 fv2 0 , 2 2 0 0 fw fx fz fy fz fz
where the w-axis is parallel to the tangent. The definition of Q is given in appendix B of the previous chapter. A schematic example of a quadratic curve and its local uvw-frame is shown in figure 4.1b. Instead of estimating the curvature κ and the angle between the normal and the v-axis, we will estimate the curvatures κ1 and κ2 . The curvature κ1 is defined as the curvature of the projection of the quadratic curve c on the uw-plane, and κ2 as the curvature of the projection of c on the vw-plane, see 4.1a. The relation between the projected curvatures and the total curvature is given by q (4.11) κ = κ21 + κ22 .
4.1.2
The quadratic GST for space curves
The first step towards a curvature corrected GST, is to define the coordinate transform uvw → u0 v 0 w 0 that deforms an arbitrary curve into a straight line, see fig 4.1b. Using the quadratic curve approximation and the κ1 , κ2 description, the coordinate transform and its inverse are given by u0 = u − 21 κ1 w 2 v 0 = v − 12 κ2 w 2 w0 = w
u = u0 + 21 κ1 w 0 2 v = v 0 + 21 κ2 w 0 2 w = w0
(4.12)
This approximation is valid for small values of κ1 and κ2 . The test results will show the practical limitations of this approximation. The inverse-transformation is used to
64
Chapter 4 Line-like curvilinear structures
w
w’
w κ1
κ2 v
u
v u
v’ u’
c
(a) projections
(b) transformation
Figure 4.1: A schematic representation of the projections of the curve c on the uw and the vw-plane and the corresponding curvatures κ 1 , κ2 (a). Deforming a quadratic curve to a straight line using a coordinate transform (b).
express the derivatives in the u0 v 0 w 0 -coordinates as a function of the derivatives in the uvw-coordinates. f u0 = u u0 f u + v u0 f v + w u0 f w = f u fv 0 = uv 0 fu + v v 0 fv + w v 0 fw = f v fw 0 = uw 0 fu + v w 0 fv + w w 0 fw = κ1 wfu + κ2 wfv + fw
(4.13)
We will now consider a curved line through the origin of the uvw-space with its tangent in the origin parallel to the w-axis, see figure 4.1b. By applying the traditional GST to the u0 v 0 w 0 -space for arbitrary κ1 and κ2 , we obtain the gradient structure tensor for quadratic curves fu20 f u0 f v 0 f u0 f w 0 (4.14) TQC = fu0 fv0 fv20 fv 0 fw 0 2 f u0 f w 0 f v 0 f w 0 f w 0
Using eq.(4.13) we can express the elements of TQC in the uvw-coordinates. fu20 = fu2 fv20 = fv2 fw2 0 = fw2 + κ21 a + 2κ1 κ2 b + κ22 c + 2κ1 d + 2κ2 e f u0 f w 0 = f v 0 f w 0 = f u0 f v 0 = 0
(4.15)
4.1 The GST for 3D line-like curvilinear structures
65
with the abbreviations a ≡ w 2 fu2 , b ≡ w 2 fu fv , c ≡ w 2 fv2
d ≡ wfu fw , e ≡ wfv fw .
(4.16)
Note that the terms a, b, c, d differ from those presented in the previous chapter. The off-diagonal elements of TQC are all equal to zero, due to the symmetry f (u, v, w) = f (u, v, −w) in the model. The derivation and the details of the symmetry considerations are given in appendix A. Since the matrix of TQC is diagonal, the eigenvalues are equal to the diagonal elements. λ1 = fu2 , λ2 = fv2 , λ3 (κ1 , κ2 ) = fw2 + κ21 a + 2κ1 κ2 b + κ22 c + 2κ1 d + 2κ2 e
(4.17)
For a line-like linear structure the eigenvalues relate as λ1 ≈ λ2 λ3 . The optimal resemblance to a quadratic curve structure can therefore be found by minimizing λ3 with respect to κ1 and κ2 . In contrary to the quadratic surface model, κ1 and κ2 are coupled and the minimization should be performed two-dimensional. The curvatures κ1,min and κ2,min that minimize λ3 yield an estimate of local curvature and the are given by κ ˆ 1 = κ1,min =
be − cd bd − ae , κ ˆ 2 = κ2,min = . 2 ac − b ac − b2
(4.18)
Substituting these curvatures in eq.(4.17), we get the following expressions λ1 = fu2 , λ2 = fv2 , λ3 = fw2 − K.
(4.19)
K ≡ aˆ κ21 + bˆ κ1 κ ˆ 2 + cˆ κ22
(4.20)
with the abbreviation The curvature corrected third eigenvalue can be used to create the curvature corrected line-confidence measure λ2 − λ 3 Cline = . (4.21) λ2 + λ 3 Note that for the limit κ1 , κ2 → 0, the equation for λ3 reduces to the third eigenvalue of the traditional GST. An estimate of the total curvature κ can be constructed from the estimates of the projected curvatures using eq.(4.11) q ˆ 21 + κ ˆ 22 . (4.22) κ ˆ= κ The computation of the structure tensor for quadratic curves TQC , consist of two parts. First, the GST is computed to obtain the gauge-coordinates uvw, and next the extra terms (a, b, c, d, e) of eq. (4.16) need to be computed. The direct computation of (a, b, c, d, e) results in a spatially variant convolution, since the gauge-coordinates change over the image. For large filter sizes, this can lead to a high computational demand. The extra
66
Chapter 4 Line-like curvilinear structures
terms have the same form as the correction terms of the structure tensor for quadratic surfaces, namely abfc fd , afc fd , a, b, c, d ∈ {u, v, w}. (4.23) The terms (a, b, c, d, e) can therefore be computed as linear combinations of convolutions, as shown in Appendix B of chapter 3.
4.2
Experimental tests and results
The structure tensor for three-dimensional space curves TQC described above, locally models an arbitrary space curve as a line-like linear structure after a quadratic coordinate transformation. The resemblance between the model and the image data is given by the confidence value Cline . The parameters of the model that maximize the confidence yield an estimate of the curvature κ ˆ . The confidence of the model fit can be decreased by two factors: noise and a deviation of the image structure from the model. In this section, we will analyze the estimates for curvature and confidence of the tensor TQC , by performing several measurements on test images. For these measurements we use images of several types of space curves. The shape of an arbitrary space curve in 3D is determined by the curvature and the torsion along the curve. Since our quadratic model only corrects for curvature we will first examine images of circles with various radii. Next we will measure the influence of torsion by applying the tensor TQC to images of a helix. The quadratic curve model assumes a curve structure that is symmetric along the w-axis. As a final test we will examine the influence of local asymmetry using images of an ellipse. Creating an image of a space curve means sampling the Euclidean space in the neighborhood of the curve. To be able to sample the curve it should have a finite thickness1 . One should make sure that the spectrum of space curve is band-limited, to be able to sample it according to the sample theorem of Shannon [Sha49]. An example of a practically band-limited image of a curve, is a curve with a Gaussian intensity cross-section Icross (u, v) = e−r
2 /2σ 2 c
, r 2 = u2 + v 2 ,
(4.24)
where r can be obtained by computing for each point in the image I(x, y, z), the distance d = r to the space curve c. A straightforward way to compute the distance to the space curve c is to find a parameterization for the curve c(t), and sample it for equal intervals ti . The distance d for each point in the image can now be computed by d = minkc(ti ) − xk , x = (x, y, z) i
(4.25)
As a 2D example, a 128*128 pixel image of a circle with a radius of 32 pixels is shown in figure 4.2. The computed distance to the circle, shown in figure 4.2a, is used as the input 1
Strictly speaking it is no longer a space curve since a mathematical curve is infinitely thin.
4.2 Experimental tests and results
67
(r) to eq. (4.24) with σc = 2 to create to image shown in figure 4.2b. Although we describe the curves using a parameterization, the estimation of curvature and confidence does not make use of this parameterization in any way. Furthermore, an arbitrary translation and rotation of the curves does not affect the estimates for curvature and confidence.
(a)
(b)
Figure 4.2: The creation of a band-limited image of a space curve. For each point in the image, the distance to the curve is computed (a), to create a space curve with a Gaussian profile (b).
We chose the width of the intensity cross-section σc = 2 voxel for all images. For both the gradient regularization and the local averaging of the structure tensor a Gaussian window is used. We computed the gradient at scale σg = 1 and the tensor smoothing at scale σT = 4 for all measurements. The scale σT is chosen such that an optimal SNR over all measurements is achieved. We will elaborate on the scale selection of the tensor analysis window, in the discussion section below. To test the robustness of the estimators with respect to noise, uncorrelated normally distributed noise N (0, σn2 ) is added to the test images. The noise dependency is measured by repeating the measurements for three different noise levels σn2 = {0.01, 0.1, 0.5}, corresponding to the signal-to-noise ratios 20,10,3 dB. Since the intensity of the test images are bounded, 0 ≤ I(x, y, z) ≤ 1, we use the following definition for the signal-to-noise ratio max(I) − min(I) dB = −20 log10 (σn )dB. (4.26) SN R = 20 log10 σn All measurements are averaged over 40 different realizations of the noise, and the corresponding standard deviations are indicated by error-bars. To show the improvement in the confidence estimation, we will also compute the confidence estimate of the traditional GST.
4.2.1
Circle image
The first type of test images Icircle (x, y, z; r, σc , σn ) are images of circle with radius r, created using equations (4.25) and (4.24), with normally distributed noise (σn ) added. The parameterization used for the circle in 3D is given by ccircle(t; r) = (r cos t, r sin t, 0).
(4.27)
68
Chapter 4 Line-like curvilinear structures
The shape of a circle is fully determined by its constant curvature. The curvature of the circle can be derived from the parameterization using equation (4.6), and is given by the simple relation 1 (4.28) κ= , r where r is the radius of the circle. We will measure both the curvature and the confidence estimate as a function of r. To give a visual impression of the circle images, figure 4.3 shows for each noise level the xy-slice (z = 0) of the test image Icircle(x, y, z; 20, 2, σn ).
(a) ∞ dB
(b) 20 dB
(c) 10 dB
(d) 3 dB
Figure 4.3: Normally distributed noise is added at three different levels for the measurement of the noise dependency.
The test image Icircle could also be interpreted as an image of a solid torus, with a Gaussian shaped intensity profile. A surface rendering of this torus is shown in figure 4.4.
2σ
R
Figure 4.4: A surface rendering of one half of a solid torus with radius R and Gaussian intensity cross-section of width σ.
The results of the curvature and confidence estimation on the circle images are depicted in figure 4.5. From the results in 4.5b we see that the curvature correction of the confidence estimate yields a significant improvement. The estimators show approximately the same noise dependency as the estimators of the quadratic surface model. The confidence estimate decreases for an increasing noise level, and the curvature estimate becomes biased for low signal-to-noise ratios, i.e. the curvature is under-estimated. The increase of the bias in the curvature estimation for decreasing SNR, is explained in appendix B. There is another bias visible in the curvature estimation in figure 4.5a. For small radii the curvature is over-estimated. This is due to the model error we make by locally modeling the circle with a quadratic curve. Since we keep the size of the analysis window constant, the quadratic curve must describe a larger part of the circles with smaller radii. The
4.2 Experimental tests and results
69
0.2
True κ 20 dB 10 dB 3 dB
κ
0.15
0.1
0.05
0
5
10
15
20
25
30
35
40
radius 1
20 dB 0.8
10 dB
Cline
0.6
0.4
0.2
0
3 dB
GST PGST 5
10
15
20
25
30
35
40
radius Figure 4.5: Curvature κ and confidence C line estimation as a function of the radius on the circleimages (σT = 4). The measurements are performed at 3,10,20 dB, and averaged over 40 different noise realizations. The error-bars denote the corresponding standard deviation. The small horizontal shift between the different labels is artificially introduced to improve their visibility.
70
Chapter 4 Line-like curvilinear structures
fourth order term of the polynomial approximation of a circle becomes more important and causes the over-estimation. Note that (Cline )GST ≤ (Cline )QGST for all radii. This can be explained by the fact that adding an extra parameter to the model, always results in a better fit.
4.2.2
Helix image
To test the influence of torsion on the estimators, we will apply the structure tensor for quadratic curves to images of a helix Ihelix (x, y, z; r, h, σc , σn ). The parameterization we used for the creation of the helix images is given by chelix (t) = (r cos t, r sin t, ht).
(4.29)
the parameter h is called the pitch and its geometrical meaning is indicated in figure 4.6. The curvature κ and the torsion τ of the helix are given by κ=
h r , τ= 2 . 2 +h r + h2
r2
(4.30)
The curvature of the helix approaches the curvature of a circle for h → 0. for an increasing pitch h, the curvature decreases and the torsion increases. We will measure both the curvature and the confidence estimate as a function of the radius r for a fixed pitch h = 10. The curvature κ(r) as function of the radius is maximal for r = h.
z h r
x
y
Figure 4.6: A Graph of a helix with its two constant parameters, radius r and pitch h, indicated.
The results of the curvature and confidence estimation on the helix-images are depicted in figure 4.8. In figure figure 4.8a we can see that curvature estimate at 20 dB is not significantly biased. The curvature values in the helix-images are κ ≤ 0.05. The curvature estimates on the circle-images for this curvature range are unbiased too. We have checked the curvature estimate on noise-free helix-images for a broad range of torsion values and we obtained an unbiased result in all cases. The improvement in the confidence estimate is comparable to the improvements measured on the circle-images, as is the noise dependency of the estimators.
4.3 Discussion
4.2.3
71
Ellipse image
For the derivation of the estimators from the structure tensor for quadratic curves, it was assumed that the curve is locally symmetric along the w-axis: f (u, v, w) = f (u, v, −w). This assumption was valid for the circle and the helix images. To test the influence of local asymmetry on the estimators, we will apply them to an image of a ellipse Iellipse(x, y, z; a, b, σc , σn ). b a Figure 4.7: A Graph of an ellipse with its radii a and b indicated.
The parameterization used for the three-dimensional ellipse-curve is given by cellipse(t) = (a cos t, b sin t, 0),
(4.31)
where a, b are the radii of the ellipse, see figure 4.7. The ellipse-curve is only symmetrical at the extrema of the ellipse, i.e. for t ∈ {0, 21 π, π, 32 π}. The non-constant curvature κ along the ellipse-curve is given by κ(t) =
(a2
ab . sin t + b2 cos2 t)3/2 2
(4.32)
The curvature κ(t) of the ellipse approaches the curvature of a circle for |a−b| → 0. We will measure both the curvature and the confidence estimate as a function of the parameter t for fixed radii a = 60, b = 30. Due to the symmetry of the ellipse, we can limit the measurements to one quadrant of the ellipse 0 < t < 21 π. If we examine the the results in figure 4.9, we see that at 20 dB the curvature estimate is unbiased for the extrema of the ellipse (t = 0, 21 π). For some of the estimated values in between, we see that the local asymmetry causes a small over-estimation of the curvature. This bias is partly due to a bias in the tangent estimate. For a neighborhood around an asymmetric curve the third eigenvector of the GST is not parallel to the tangent. The presence of noise causes an under-estimation of the curvature. The noise dependency of the estimators is comparable to the noise dependency measured on the circle-images, and so is the improvement of the confidence estimate.
4.3
Discussion
For the description of line-like curvilinear structures we have presented a curvature corrected gradient structure tensor TQC . A line-like curvilinear structure was defined as shift invariant along a space curve. The curvature correction is based on a local approximation
72
Chapter 4 Line-like curvilinear structures
0.06
True κ 20 dB 10 dB 3 dB
0.05
κ
0.04
0.03
0.02
0.01
0
5
10
15
20
25
30
35
40
radius 1
20 dB 0.8
10 dB
Cline
0.6
0.4
0.2
0
3 dB
GST PGST 5
10
15
20
25
30
35
40
radius Figure 4.8: Curvature κ and confidence C line estimation on the helix-image as a function of the radius r, with a constant pitch h = 10 (σ T = 4). The measurements are performed at 3,10,20 dB, and averaged over 40 different noise realizations. The error-bars denote the corresponding standard deviation. The small horizontal shift between the different labels is artificially introduced to improve their visibility.
4.3 Discussion
73
0.08
True κ 20 dB 10 dB 3 dB
κ
0.06
0.04
0.02
0 −0.1
0.3
0.7
1.1
1.5
t 1
20 dB 0.8
10 dB
Cline
0.6
0.4
3 dB GST PGST
0.2
0 −0.1
0.3
0.7
1.1
1.5
t Figure 4.9: Curvature κ and confidence C line estimation on the ellipse-image as a function of the parameter t, with constant radii a = 60, b = 30 (σ T = 4). The measurements are performed at 3,10,20 dB, and averaged over 40 different noise realizations. The error-bars denote the corresponding standard deviation. The small horizontal shift between the different labels is artificially introduced to improve their visibility.
74
Chapter 4 Line-like curvilinear structures
of a space curve with a quadratic curve. The complete representation of the shape of a three-dimensional space curve requires two parameters for each point on the curve, namely curvature κ and torsion τ . Our model, on the other hand, only incorporates curvature. Furthermore, mirror symmetry with respect to the plane normal to the tangent vector of the curve is assumed. To test the robustness of the estimators for curvature and confidence derived from TQC with respect to deviations of the image structure from the local model, we applied them to images of a circle, a helix, and an ellipse. Furthermore, we added noise to the space curve images to measure the robustness with respect to noise. For all measurements we kept the thickness of the space curves and the scale of the tensor window constant. Increasing the tensor scale σT suppresses the influence of noise, i.e. the standard deviation of the measurements decreases. However, a larger analysis window requires the modeling of a larger part of the space curve. The model errors become more apparent and the curvature estimate becomes biased. Another draw-back of increasing σT is a decrease in signalto-noise ratio. Consider an individual line in an isotropic three-dimensional window with radius R, with noise. Increasing the window size causes the total amount of signal energy to increase linearly with the radius of the window. The total noise energy, however, increases with R3/2 , assuming uncorrelated normally distributed noise. The estimators show approximately the same noise dependency as the estimators of the quadratic surface model. The confidence estimate and its improvement with respect to the traditional GST confidence, decrease for an increasing noise level. The curvature estimate becomes biased for low signal-to-noise ratios, i.e. it gives an under-estimation. The measurements on the helix-images show that a constant torsion has no measurable effect on the estimators. The measurements on the ellipse-images show that local asymmetry does influence the estimation of curvature, i.e. it causes an over-estimation of the curvature.
4.4
Application: Channel detection
The occurrence of sedimentary structures, such as channels, is an important cue for the geological model of a subsurface region. The manual localization of channels in often vast amounts of 3D seismic data, is very time consuming. One way of automating this process is to find a suitable mathematical structure model for channels. The resemblance between the model and the structure that is being analyzed should provide a measure of how much this structure resembles a channel. A characteristic feature of channels is their meandering shape. Therefore, a parametric model based on shape parameters such as curvature does suggest itself. In this section we will study the applicability of the structure tensor to the detection of channels. The image data we are going to analyze is a region (256*128*64 voxels) of a 3D seismic image I(x, y, t) around a channel. The intensity of the image corresponds to the amplitude
4.4 Application: Channel detection
75
of the reflected seismic wavelet. The xy-axes are spatial and the t-axis indicates the travel time of the seismic wavelet. The t-axis can be inverted to depth (spatial) if the velocity of the acoustic wave is known at each depth. To give an impression of how this channel manifests itself in this data, we have visualized it by 2-D cross-sections in figure 4.10. The (time) xy-slice shows the horizontal extend of the channel, that has the characteristic meandering shape. The yt-slice shows a cross-section of the channel. Figure 4.10d depicts the t-axis along the centre of the channel showing the depth (time) variation of the channel. As a preprocessing step, we create an attribute volume by computing the standard deviation within a window of (5*5*9) voxels, for each point in the seismic image. A cross-section of this volume, corresponding to figure 4.10a, is depicted in figure 4.10c. The advantage of this attribute volume is that the channel manifests itself as a space curve. The meandering nature of channels suggests that the curvature corrected model gives a more accurate description of a channel than the traditional straight model. Therefore we expect that the curvature corrected confidence estimate has a significantly higher value than the confidence estimate of the traditional GST. To test this hypothesis we have applied both the GST and structure tensor for quadratic curves (QGST) to the seismic image. In both cases we used a gradient smoothing σg = 1 and a tensor smoothing σT = 3.5. The choice of the tensor scale is a trade-off as pointed out in the discussion in section 4.3. For the analysis of the space curve images we used a tensor scale that was twice the width of the curves (σT = 2σc ). Here we use the same rule, and we roughly estimated the width of the channel. The results of the QGST confidence estimate and the improvement in the confidence estimate are shown in figure 4.10e and f. Due to the contract-stretch it is not possible to determine the relative improvement from these image. The purpose of (f) is to show where the improvement is achieved. The relative improvement in the bright spots in (f) is approximately 50%. For the evaluation of the result, we defined three regions in the seismic image. The first region contains the channel and is indicated in figure 4.11a with a thick white contour. The second region contains a straight sedimentary structure and is depicted by a black contour in figure 4.11a. The last region consist of points from an area without sedimentary structure, and represents the noise in the data. This area is indicated by a white rectangle in figure 4.11b. For all three regions we computed the cumulative frequency distribution of the confidence values in that region of both the GST and the QGST confidence estimate. These cumulative distributions are displayed in figure 4.11c. We will use the cumulative distributions to verify our hypothesis that the curvature correction yield a significant better description of channels. For notational convenience we denote the improvement in the confidence value as ∆C = (Cline )QGST − (Cline )GST .
(4.33)
Inspection of the confidence distributions of the channel region in figure 4.11c, shows that the improvement in this region ∆C is larger than zero. However, this does not necessarily mean that the improvement is significant. As we have already seen in the test results
76
Chapter 4 Line-like curvilinear structures
A
A
A’
A’ (a) xy-slice (t = 33)
(c) attribute volume (t = 33)
(e) (Cline )QGST (t = 33)
(b) yt-slice (x = 133)
(d) t along channel
(f) ∆C (t = 33)
Figure 4.10: Seismic data volume used to demonstrate channel detection (a,b,d), and a slice of the attribute volume (c). The resulting confidence estimate and the confidence improvement are shown in (e) and (f).
4.5 Appendix A: Symmetry
77
that the QGST confidence estimate always yields a higher value than the GST confidence estimate. Since ∆C > 0 in the noise region, the significance of ∆C along the channel should be proven by correcting it for the improvement due to a better description of the noise ∆Cn . For the following reasoning we will assume that the noise has the same distribution at each position in the volume. The improvement ∆C ≈ 0.03 in the noise region (SN R = −∞) gives an estimation of the improvement due to noise ∆Cn . However it should be taken into account that the value of ∆Cn is not constant for all Cline values. Imagine for example a noise-free neighborhood around a straight line (SN R = ∞). In that case the line confidence Cline = 1 and improvement ∆C = 0 and thus ∆Cn = 0. To get an impression of how ∆Cn changes for higher values of Cline we have included the ”straight structure” region. In this region the curvature correction improves the fit by giving a better description of the noise. Therefore we have ∆Cn ≈ 0.01 at Cline = 0.55. We can now conclude that the noise contribution ∆Cn to the improvement ∆C ≈ 0.05 in the channel region is smaller than 0.01. The curvature correction therefore yields a significant better description of the meandering structure of the channel. Summarizing the results we conclude that the detection of channels, based on the estimators of structure tensor alone is not feasible. This is most likely due to the local nature of the tensor analysis. Increasing the size of the analysis window isotropically does not result into a better description, for the same reasons we mentioned in the discussion of the analysis of space curve images. For a more accurate analysis of channels, windows that adapt to the shape of the channel are needed. We have shown that the curvature corrected structure tensor for space curves TQC , gives a better description of channels than the GST. Therefore TQC should be used to control a subsequent processing step with adaptive windows. Windows that are able to match the shape of a channel will be presented in the next chapter.
4.5
Appendix A: Symmetry
The symmetry in the quadratic curve model f (u, v, w) = f (u, v, −w) with respect to the uvw-coordinates, gives rise to the following symmetries in the gradients fu (u, v, w) = fu (u, v, −w) fv (u, v, w) = fv (u, v, −w) fw (u, v, w) = −fw (u, v, −w)
(4.34)
We will use the fact that the integration of an anti-symmetric function equals zero, to show that the off-diagonal elements of the quadratic structure tensor for curves, equal zero. The combination of eqs.(4.34) and the equations (3.36), (3.38) for (anti-)symmetric functions
78
Chapter 4 Line-like curvilinear structures
(a) xy-slice (t = 33)
(b) xy-slice (t = 47)
1
0.8
0.6
Noise
Straight structure
Channel
0.4
GST QGST
0.2
0
0
0.2
0.4
0.6
0.8
1
Cline (c) cumulative distributions
Figure 4.11: The cumulative frequency distributions of the confidence values (c) in selected regions of the seismic volume (a)+(b). The noise region is indicated by a white rectangle in (b). The region with the straight structure is indicated with a black contour and the channel region with a thick white contour in (a).
4.6 Appendix B: Bias
79
from appendix A of chapter 3, yields the following results fu0 fw0 = fu fw + κ1 wfu2 + κ2 wfu fv = 0 fv0 fw0 = fv fw + κ1 wfu fv + κ2 wfv2 = 0
(4.35)
f u0 f v 0 = f u f v = 0 The last equation requires mirror symmetry with respect to an arbitrary plane perpendicular to the uv-plane.
4.6
Appendix B: Bias
From figure 4.5a it is clear that noise introduces a bias in the curvature estimation. The expectation value of the curvature estimator E[ˆ κ] is smaller than the true curvature. This can be explained by examining κ ˆ 1 and κ ˆ2 . κ ˆ1 =
κ ˆ2 =
w 2 fu fv wfv fw − w 2 fv2 wfu fw w 2 fu2 w 2 fv2 − (w 2 fu fv )2 w 2 fu fv wfu fw − w 2 fu2 wfv fw w 2 fu2 w 2 fv2 − (w 2 fu fv )2 q κ ˆ= κ ˆ 21 + κ ˆ 22
(4.36)
(4.37) (4.38)
The terms w 2 fi2 , i ∈ {u, v}, increase when noise is added, while the other terms do not change. Since the increasing terms appear quadratic in the denominator and linear in the nominator, the curvature estimate κ ˆ becomes smaller for an increasing noise level. The sign of κ ˆ 1 and κ ˆ 2 is lost in eq.(4.38).
5. Structural analysis using a non-parametric description The local structural analysis of images is based on the idea that complex image structures can be described at a small scale using a simple geometric model. We started the structural analysis of images by assuming local shift invariance, and we called an image structure with shift invariance in one or more orientations a linear structure. The structure tensor was presented as a versatile tool for the analysis of these structures. Increasing the analysis scale leads to the description of a larger part of the complex structure. The curvature of the structure becomes gradually more important, and therefore we introduced the curvature corrected structure tensor based on a quadratic model. The application of the confidence value of this model to the detection of channels showed that the local structure analysis is not sufficient for this detection task. An important feature of the channel structure is its spatial continuity. At the scale needed to capture this continuity, the quadratic model is not flexible enough to describe the structure. Increasing the flexibility of a parametric model means increasing the number of parameters, and thereby increasing the computational demand. Instead of introducing a more complex parametric model with more parameters, we will turn to a non-parametric model. The local parameter estimates obtained from the simple geometric models can be used for a piece-wise description of curvilinear image structures. The shape of the structure can be found by the tracking of this piece-wise description. The shape of two-dimensional and line-like curvilinear structures can be described using a single space curve. This shape can therefore be found by the tracking of an individual curve. The shape of a plane-like curvilinear structure is defined by a surface. The tracking of a surface is significantly more complex than the tracking of a curve, and falls outside the scope of this thesis. The nonparametric description of plane-like curvilinear structures is therefore not discussed in this thesis. This chapter begins with the description of the tracking method. The method will be demonstrated and tested by applying it to the tracking of growth rings and channels. Next, the tracking method is used for the creation of a non-parametric adaptive filter. The results of this filter applied to the noise reduction in an oriented texture, will be compared to the results of the adaptive filters presented in chapter 3. We will conclude this chapter with the estimation of a non-parametric confidence measure, and its application to channel detection.
82
5.1
Chapter 5 Structural analysis using a non-parametric description
The tracking of line-like curvilinear structures
In chapter 4 we generalized line-like linear structures to line-like curvilinear structures. An image neighborhood has a line-like curvilinear structure if the isophotes or the level-sets in that neighborhood can be described as an individual space curve or a bundle of parallel of space curves. In this section we will present a method for the tracking of these structures. With the word ‘tracking’, we refer to the iterative process of following the spatial extent of a selected object. In the case of curvilinear structures, the objects are the isophotes. Since the isophotes form space curves, we start with the tracking of space curves. Consider a space curve c(t) and two points t=P and t=Q on this curve, as depicted in figure 5.1a. We select point P as the start position. For the tracking of the curve we need to find the displacement vector ∆c = c(Q) − c(P ),
(5.1)
to make the first step to the next point on the curve, Q. Since the parameterization of the curves formed by isophotes is in general not known, we need to estimate the displacement vector ∆c using a local geometric model. From differential geometry [Spi79], we know that the shape of an arbitrary space curve in 3D is fully determined by the curvature and the torsion along the curve. Therefore, if we have an estimation of the shape parameters at point c(P ), then we can use these parameters to estimate the displacement vector ∆c.
t4
c
P
n
c t
∆c
Q (a)
t3 t1
( φ 1, κ 1)
t2
t5
( φ 3, κ 3) (b)
Figure 5.1: Tracking a space curve c by following the displacement vectors ∆c. The displacement vectors are computed using the locally estimated parameters orientation and curvature (φi , κi ).
We start by applying this idea to the tracking of a space curve c in two dimensions. If we assume that the tangent vector t and the normal vector n are known in each point, then the displacement vector can be approximated by 1 (5.2) ∆c ≈ n κ(∆t)2 + t(∆t), 2 with ∆t the step size. Suppose that we have a sampled version of the space curve c(ti ), then we can use eq.(5.2) to estimate the displacement vectors between the sample points t i and ti+1 . Following these vectors ‘head-to-tail’ as shown in figure 5.1b, we track the curve.
5.1 The tracking of line-like curvilinear structures
83
Translating this method to the tracking of a two-dimensional curvilinear image structure is not difficult. An estimate for the tangent and the normal vector for each point of the structure, is given by the eigenvectors of the structure tensor {e1 , e2 }. The curvature estimate κ ˆ of the curvature corrected structure tensor can be used as the curvature in eq.(5.2). The frame of the GST in 2D {e1 , e2 } can be fully determined by one angle between the first eigenvector and the x-axis of the image. Therefore we can track the image structure by following the path that is piece-wise defined by the local orientation and curvature estimates (φi , κi ), as depicted in 5.1b. The tracking of a space curve in three dimensions requires more parameter estimates. A full local description of the curve is given by the Frenet-frame and the shape parameters curvature and torsion [Spi79]. We can approximate the curve using the quadratic model presented in chapter 4. The quadratic curve was parameterized using the local frame of the structure tensor {e1 , e2 , e3 } and two curvatures κ1 , κ2 . The third eigenvector e3 is an estimate for the tangent, and the curvatures κ1 and κ2 are defined along resp. e1 and e2 . In chapter 4, we have experimentally shown that the quadratic curve model is not significantly hampered by torsion. If we start at an arbitrary point ti on the space curve c then we can track the curve by following the displacement vectors 1 1 ∆c(ti ) ≈ e1 κ1 (∆ti )2 + e2 κ2 (∆ti )2 + e3 (∆ti ), 2 2
(5.3)
where ∆ti is the step size that can be chosen freely. Increasing the step size however, can increase the model error and thereby the tracking error. So far we have assumed that the start point from which we start the tracking is known or manually selected. If the starting position is not known and manual selection is not desirable, the curve needs to be detected first. The automatic detection of a space curve constitutes a more difficult problem than the tracking of the curve. Under favorable circumstances the confidence estimate of the local models can be used for detection. We will investigate the detection capabilities of the confidence estimates later on in this chapter.
5.1.1
Application to the tracking of growth rings
To demonstrate the tracking method described above, we applied it to the tracking of the growth rings of a tree. A 2D cross-section of a CT image of a tree-trunk that shows the growth rings, is depicted in figure 5.2a. We started the analysis of this 2D cross-section by applying the curvature corrected structure tensor at scale (σg = 1, σT = 3). The resulting estimates for orientation, curvature and confidence are depicted in figures 5.2e-g. We selected three starting points for the tracking algorithm, and these points are indicated by white arrows in figure 5.2d. The tracking results based on the orientation and the curvature estimate are shown in figure 5.2b as black curves super-imposed on the original image. For comparison we have shown the tracking result using only the orientation estimate in figure
84
Chapter 5 Structural analysis using a non-parametric description
5.2c. The difference between the tracked paths of the two methods is shown in 5.2d. A step size of one pixel is used in both cases, and we allowed sub-pixel positions. Evaluation of the tracking of the two inner rings (1,2), shows that following the estimated orientation and curvature according the eq. (5.3), results in a closed path that matches the corresponding growth rings. If the curvature estimate is not used, then the tracked path is no longer closed but has a spiral shape. Evaluation of the tracking of the third outer ring (3), shows that the tensor analysis is based on the differential geometry of the intensity values, i.e. it is a gradient model. The parameter estimators are only sensitive to intensity changes not to the intensity value itself. The dominant gradients in the neighborhood of the third ring are those belonging to the change in the background, not those on the slopes of the local intensity minimum that constitutes the growth ring. The resulting path therefore matches the shape of the background transition. Again we see that the path based on orientation alone is not able to keep up with the curvature of the structure. Eventually it follows a more outward laying growth ring uninfluenced by the background transition. The purpose of this section was to show that the parameter estimates of the curvature corrected structure tensor can be used for the piece-wise definition of paths that follow the local structure. A more accurate tracking of the growth rings could be obtained by utilizing more features. To give a simple example: one could check the intensity value along the path. If the intensity value of a new point differs too much from some average value along the path, then the path has probably taken a wrong turn.
5.1.2
Application to the tracking of sedimentary structures
In the previous chapter we showed that the confidence value of the quadratic structure tensor applied to a 3D seismic variance attribute volume, can detect channels. The curvature correction improved the detection of the strongly bent parts of the channel. However, this confidence based detector gives a lot of false positives, i.e. there are still a lot of non-channel structures that give rise to a high confidence value. The analysis of the channels at a relatively small scale, approximately two times the width of the channel, is not sufficient to distinguish between channel and non-channel. Taking the spatial continuity of the channel into account could improve the channel detection. One way to incorporate the continuity of channels into the detection, is to extend the analysis window along the structure. To accomplish this, a technique for the tracking of channels is needed. In this section we will test the ability of the tracking method described above, to track channels and related sedimentary structures. We will apply the tracking method to a horizontal 2D cross-section (512*256 pixels) of a 3D seismic variance attribute image. The planar cross-section shown in figure 5.3a, is chosen such that it approximates the depositional surface, and it thereby makes the sedimentary structures visible. We selected three clearly visible channels and show them in figure 5.3b. The structure tensor analysis can be interpreted as imposing a local parametric model on
5.1 The tracking of line-like curvilinear structures
(a) Wood
85
(b) Tracking O+C
(c) Tracking O
(e) Orientation
(f) Curvature
3
1 2
(d) Difference
(g) Confidence
Figure 5.2: Tracking the growth rings in a 2D cross-section of a CT image of a tree-trunk. The results are shown using only orientation (c) and using orientation and curvature (b). (d) shows the difference between the two methods and the starting points indicated by white arrows.
86
Chapter 5 Structural analysis using a non-parametric description
the gradients. To give an impression of the gradient information, we computed the gradient magnitude at scale (σg = 2) and depict it in figure 5.3c. Instead of manually selecting a starting point for the tracking algorithm for each structure we want to track, we will use all ‘edge points’ as start positions. The edge points are all points in the image that, according to an edge detector, belong to an edge. The edge points are usually found by locating the points where the gradient magnitude has a maximum [Can86]. The traditional way to find this maximum is to locate the zero-crossings of the second derivative in the gradient direction. We applied a one-dimensional window (7*1 pixel) steered along the gradient direction to each point in the gradient magnitude image. If the intensity value of the current pixel is equal to the maximum value inside this window, then it belongs to an edge. The result of this edge detection is shown in figure 5.3d. Next, we apply the curvature corrected structure tensor at scale (σg = 2, σT = 4) to the 2D attribute image shown in figure 5.3a, to estimate the orientation and curvature (φ, κ) at each point in the image. The angle φ is the angle between the x-axis of the image and the second eigenvector of the GST, and gives an estimate of the tangent direction. Consider an arbitrary edge point as starting position P0 and the corresponding orientation estimate φ0 , see figure 5.4. Since the structure continues in both the φ0 and −φ0 direction, we initiate the tracking algorithm in both directions. The only stop criterion for the tracking is a maximum number of steps n. The result of this tracking experiment using a fixed step size of one pixel and n = 100, is shown in figure 5.3e. The tracked paths are visualized by increasing the intensity value of the output image at all the points of all the paths with a constant value 1. Thus, if forty paths cross the point (x, y), then the intensity value at that point I(x, y) = 40. If we examine the tracking result near the channel on the left and the channel in the middle, we see that the tracked paths are parallel to the channel banks. This means that the tracking method is able to follow these channels. The tracking of the bottom-left channel is more difficult and fails near the ‘hairpin turns’. At the locations where the channel banks constitute salient edges in the attribute image, the tracking results in lines with a relatively high intensity value. This intensity value is not related to the gradient magnitude value, but indicates that all the paths started on that edge, follow the same route. To enhance the edges where most of the tracked paths coincide, we added two extra stop criteria based on heuristics to the tracking algorithm. If during the tracking the gradient magnitude value of a new point differs more than 40% from the average value along the path, the tracking is stopped. Assuming that an edge ends if the gradient magnitude value changes rapidly, this should prevent the path from continuing when it is no longer on the edge. If the angle between two subsequent displacement vectors is larger than 40 degrees, the tracking is stopped as well to eliminate incoherent paths. The result of this conditional tracking algorithm is shown in figure 5.3f. The edge where the paths are consistently defined are clearly visible in this image. The higher the intensity value the more consistent the tracking. Again we see that the tracking algorithm is able to follow the banks of the
5.1 The tracking of line-like curvilinear structures
87
(a) Attribute image
(b) Three channels
(c) Gradient magnitude
(d) Edge detect
(e) Tracking
(f) Conditional tracking
Figure 5.3: Tracking of the structures in a 2D cross-section (a) of a seismic attribute image. Three channels in this cross-section are manually highlighted in (b). The tracking result using the edge points (d) as the starting positions is shown in (e). The consistently tracked edges are enhanced by applying extra stop criteria (f). In both cases we use a maximum number of step n = 100, and a step-size of one pixel.
88
Chapter 5 Structural analysis using a non-parametric description
P1 P−n
P−2
P2
P0 P−1 ( φ , κ ) 0 0
Pn
Figure 5.4: Tracking a structure by following the displacement vectors. The structure is tracked from the starting point P0 in both directions φ0 and −φ0 . The tracking is stopped after n steps.
channel on the left and the channel in the middle, but fails on the bottom-right channel.
5.2
Non-parametric adaptive filtering
In the previous chapters we have shown how the parameter estimates of the local structure models can be used to control adaptive filters. In chapter 3 we applied four different filter types to the noise reduction in curved oriented textures. These filter types were in order of increasing complexity: isotropic, orientation adaptive (OAF), curvature and orientation adaptive (COAF), and ideal. The ‘ideal’ filter is the filter that exactly matches the underlying structure of the data. To obtain the most accurate filter result, the adaptive filter should approximate the ideal filter as much as possible. On the other hand, the computational burden of the filter and the required parameter estimates should remain practical. In figure 5.5a we have drawn a space curve c, and three isotropic analysis windows at three different scales (1-3). As mentioned before, the goal of the adaptation of the filter, is to approximate the shape of the structure a well as possible. Figuratively speaking we want to ‘deflate’ these isotropic windows so they wrap around the curve tightly. To be able to do so, we need a mathematical description of the curve inside the isotropic window. Describing the space curve at scale 1 using a linear model will probably result in a high confidence value. At a larger scale (2) the curvature of the curve becomes more important and a quadratic model is required to obtain a confident description. Making the parametric model flexible enough to describe the part of the space curve inside the third window, requires more parameters. The estimation of these extra parameters increases the computational burden, and the applicability of the higher order parameter estimates is limited. Instead of using a more flexible parametric description, we choose to use a non-parametric description. The parameter estimates of the quadratic local model at scale 2 are used to describe the curve at larger scales (e.g. 3) as piece-wise quadratic. The actual shape of the
5.2 Non-parametric adaptive filtering
c
2 1
89
3 (a) Window scales
(b) Non-parametric window
Figure 5.5: Analyzing a space curve at three different scales (a). The largest scale requires a non-parametric window shape (b).
structure is found by tracking. The tracking is performed as described in section 5.1.2 and is schematically depicted in figure 5.4. If P0 is the starting point and t0 is the estimated tangent at this point, then the tracking algorithm is initiated in both the t0 direction and in the −t0 direction. The tracking is stopped after a fixed number of steps n and no additional stop criteria are used. The non-parametric adaptive filter is created in the same way as the adaptive filters presented in chapter 2 and 3. First the appropriate filter scale is selected by applying an isotropic filter to the image. Next, a non-parametric filter is created at each point in the image by sampling and a weighted averaging of the isotropic outputs along a tracked path that is piece-wise defined by the parameter estimates of the (curvature corrected) structure tensor. The resulting non-parametric filter shape for the curve c of figure 5.5a is shown as an example in figure 5.5b. For the visual and quantitative comparison of the non-parametric filter with the other adaptive filters, we will use the 2D test image It of a curved oriented texture given by It (x, y; σn ) = It (r) = cos
p 2π r + N(0, σn2 ) , r = x2 + y 2 , 8
(5.4)
where N is normally distributed noise. A realization of It (x, y; σn ) with σn2 = 0.5, in a 128*128 pixel image, is shown in figure 5.6a. The orientation and curvature estimates used to control the filters are computed by applying the tensor Tcurv2D at scale (σg = 1, σT = 5), and they are depicted in figure 5.6e and f. The result of the non-parametric adaptive filter (NPAF) applied to the test image, is shown in figure 5.6b. For the visual comparison of the different adaptive filters we have shown the results of the COAF and the OAF in figure 5.6c and d. In all three filter results we have indicated one instance of the filter with a white drawing. From these drawings one can see that in the high curvature area the non-parametric filter is able to match the local structure, but the COAF is not. As a result the NPAF yields visibly better result in this area. We performed a quantitative comparison of the filters by computing the signal-to-noise ratios (SNR) of the filtered images. Since the noise-free version of the test image I =
90
Chapter 5 Structural analysis using a non-parametric description
(a) Test image
(b) NPAF
(c) COAF
(d) OAF
(e) Orientation
(f) Curvature
Figure 5.6: Using adaptive filters (b-d) to improve the SNR of a curved oriented texture (a). The orientation and curvature estimate used to control the filters are shown in (e,f). One instance of each is filter is indicated with a white drawing superimposed on the filter result.
5.2 Non-parametric adaptive filtering
91
It (x, y; 0) is known and (max(I) = 1, min(I) = −1), we compute the SNR using 2 SN R = 20 log10 dB. σn
(5.5)
The standard deviation σn of the remaining noise in the filtered images f ilter(It ), is estimated by σ ˆn = STDDEV[f ilter(It ) − I]. (5.6)
The resulting global signal-to-noise ratios and standard deviations estimates of the test image after noise reduction are given in table 5.1. In the right column we have shown the signal-to-noise ratios estimated in the high curvature area, i.e. for r < 20. Filter Ideal NPAF∗ NPAF COAF OAF Isotropic None
σn SNR (dB) 0.147 22.7 0.147 22.7 0.152 22.4 0.162 21.8 0.238 18.5 0.559 11.1 0.710 9.0
SNR1 (dB) 21.7 21.7 21.2 18.5 14.5 11.1 9.1
Table 5.1: A quantitative comparison of noise reduction filters applied to a circular oriented texture. The quantities σn and SNR are computed globally, SNR1 is computed for the high curvature area (r < 20).
All the filters of table 5.1 use a Gaussian window, and the unfiltered result ‘None’ is added to indicate the SNR of the original test image. The isotropic filter is applied at scale √ σ = 5 and all the other filters at scale σ1 = 51 . The scale of the isotropic filter is chosen such that the area of the corresponding window is equal to the area of the window of the adaptive filters. The ideal filter result indicates the highest attainable SNR at this filter scale. The ideal filter shape for this test image is a circular line piece, and the ideal filter is controlled by the true parameter values of the noise-free image. From the SNR values in table 5.1 it is clear that the NPAF is the adaptive filter that yields the highest SNR, and that this SNR is close that of the ideal filter. In figure 5.6 we saw that the non-parametric filter is able to match the circular structure. To show that this filter can mimic the circular filter, we have applied it using the true parameter values of the noise-free image as control parameters (NPAF∗ ). The experimental results show that there is no measurable difference between the SNR of the NPAF∗ and the ideal filter. One would expect that the differences between the filters are the biggest in the high curvature area. Therefore we added an extra column to table 5.1 with the SNR1 estimated in the high curvature area, r < 20. It is clear that the SNR of the OAF changes the most, and this filter performs even worse than the isotropic filter in this area. The NPAF still yields the best results, close to the performance of the ideal filter. 1
The suffix 1 in σ1 indicates that the filter is one dimensional as described in chapter 2
92
5.3
Chapter 5 Structural analysis using a non-parametric description
Non-parametric confidence estimation
In the previous section we computed a Gaussian weighted average of isotropic filter responses along tracked paths to create an adaptive filter. These paths could be used in the same way, to combine the analysis windows of the structure tensor. The advantage of a larger tensor window is that the continuity of the structure that is being analyzed, has more influence on the confidence value. This allows for the segmentation between continuous structures and fragmented structures. Straightforward application of the non-parametric adaptive filters described above to the averaging of the structure tensor elements, results in an unnatural description of the tracked structure. It describes the structure inside the non-parametric window with a single parametric model, see figure 5.7a. The natural extension of the structure tensor would be a piece-wise description as depicted in figure 5.7b. Both the presented linear and quadratic model can be used as the basis of this piece-wise model. For the remainder of this section, we will use the quadratic model.
e1 e1
e2 (a) Single model
e1 e2
e1
e2 e2
(b) piece-wise model
Figure 5.7: The combination of structure tensor analysis windows along a tracked path. Averaging of the elements of the GST using a non-parametric window assumes a linear structure inside this window (a). A piece-wise model (B) gives a better description of a tracked structure.
Since the piece-wise model gives a non-parametric description of the data, it does not yield new (geometric) parameter estimates. However, the confidence value of the piecewise model can be computed by averaging the eigenvalues of the parametric models it is composed of. First, the curvature corrected structure tensor is applied to the image using isotropic windows. This yields the eigenvalues of the tensor and an estimate of the orientation and curvature of the local structures. The orientation and curvature estimates define, for each point in the image, a path along the local structure that can be found by tracking. The tracking is performed similar to the tracking for the non-parametric filter, and is schematically depicted in figure 5.4. If P0 is the starting point and t0 is the estimated tangent at this point, then the tracking algorithm is initiated in both the t0 direction and in the −t0 direction. The tracking is stopped after a fixed number of steps n and no additional stop criteria are used. The eigenvalues of the non-parametric model (λi )N P are found by averaging the eigenvalues (λi )m of the N parametric models m that
5.3 Non-parametric confidence estimation
93
correspond to the points of the tracked path (λi )N P
N 1 X = (λi )m , N = 2n + 1. N m=0
(5.7)
These averaged eigenvalues can be used to compute the confidence value of the nonparametric model. In the two-dimensional case one can compute the confidence value using (λ1 )N P − (λ2 )N P . (5.8) (C)N P = (λ1 )N P + (λ2 )N P This confidence value indicates how well a structure can be modeled as piece-wise quadratic.
5.3.1
Application to channel detection
One of the practical goals of this thesis is to develop tools for the automatic detection of channels and related sedimentary structures. In chapter 4 we used the confidence value (C)QGST of the structure tensor based on the quadratic curve model, for the detection of channels in a seismic volume. The channels give rise to a higher confidence value than the background, but there are still a lot of other structures, in the data with a confidence value comparable to that of the channels. The visually most distinguishing feature between the channels and the other structures is the larger spatial continuity of the channels. This continuity is not exploited by the local analysis of the structure tensor. The quadratic parametric model did however give a confident local description on most parts of the channels. Earlier in this chapter, we utilized the corresponding parameter estimates for the tracking of channels. This tracking method was able to follow the channels, except for the location where a channel makes a radical turn. The non-parametric model described above uses these tracked paths to increase the analysis scale. Therefore we expect that the confidence value of this model (C)N P , is more suitable for the detection of channels than the confidence value (C)QGST . To experimentally verify this hypothesis we use both the confidence values (C)QGST and (C)N P for the detection of the channels in a 2D cross-section (512*256 pixels) of a seismic variance attribute volume. This cross-section shown in figure 5.8a and is identical to the 2D image in figure 5.3a, that was used for the tracking of channels. The manual channel detection is shown in figure 5.8b, and will be used for the evaluation of the automatic detection. We applied the curvature corrected structure tensor at scale (σg = 2, σT = 4) to the 2D cross-section and the confidence estimate is shown in figure 5.8c. The confidence value (C)N P of the non-parametric model is computed using eq.(5.8). The eigenvalue of the structure tensor are averaged along the tracked paths defined by the orientation and curvature estimate. The number of steps n of the tracking algorithm was set to 50, and the step size to one pixel. The result of the (C)N P estimation is shown in figure 5.8e.
94
Chapter 5 Structural analysis using a non-parametric description
(a) Channel image
(b) Manual detection
(c) (C)QGST
(d) Detection using (C)QGST
(e) (C)N P
(f) Detection using (C)N P
Figure 5.8: Channel detection based on the confidence estimate of the quadratic structure tensor (C)QGST (d) and the non-parametric model (C) N P (f). The visual comparison of these detections with the manual detection (b) shows that (C) N P yields a better detection.
5.3 Non-parametric confidence estimation
95
The confidence values do not give a segmentation of the image, but an estimate of the ‘channelness’ of each point in the image. For a segmentation or classification into the classes ‘channel’ and ‘non-channel’ based on the confidence estimation, a decision boundary or threshold value should be introduced. All the points with a confidence value above the threshold value, are classified as ‘channel’, and all other points as ‘non-channel’. The segmentation of the confidence values (C)QGST , (C)N P using a threshold value Ctr = 0.75 are shown in resp. figure 5.8d,f. The white regions form the channel class and the black region the non-channel class. The threshold value Ctr = 0.75 was chosen manually and optimized by visual comparison with result of the manual detection. This optimization is clearly very subjective. A more objective method to optimize the threshold value is to ‘learn’ it from examples. Many supervised statistical classifiers are developed in the pattern recognition community [Bis95]. classifiers do however require a set of examples of both classes. Since the ground truth in this detection problem is not known, these examples have to be selected manually by the human expert. Comparison of the channel detection in figures 5.8d and f, shows that the detection based on confidence value of the non-parametric model is less fragmented. The three channels of the manual detection are partly found in both cases, but in figure 5.8f they form three continuous regions. Furthermore, the number of false positives is reduced significantly due to the incorporation of the continuity of structures in the confidence estimation. We conclude that the hypothesis that (C)N P is more suitable for the detection of channels than the confidence value (C)QGST , is correct. For the determination of the false positive detections we have used the manual detection as the ground truth. The largest false positive in the top-left corner of figure 5.8f could however be a correct detection. This is not clear from this cross-section. Inspection of the 3D seismic data showed that in fact a part of another channel was found.
6. Coherency estimation Sedimentary rock consists of petrified layers of sediments that were deposited on top of each other in previous geological times. An individual layer forms a surface that approximately corresponds to the earth’s surface where the deposits were made upon. Therefore the layers usually form continuous1 surfaces. Discontinuities in the layered structure of sedimentary rock often reveal geologically interesting processes. Discontinuities can for example be created by faults1 or incisions made by channels. The continuity of the sedimental layers is often called the bedding continuity. The continuity that can be measured in a seismic image is the reflection continuity. Therefore we need to link the bedding continuity to the reflection continuity. The rock property that is measured by the seismic reflection method is the acoustic impedance. If we could measure this property directly we would obtain a 3D data volume with in each point the acoustic impedance. During the seismic acquisition an acoustic wave package is generated at the surface. This package travels down the subsurface that acts like an acoustic filter and partly reflects where the acoustic impedance changes. The reflected wave packages are recorded at the surface. The amplitude of this wave package can be used for the determination reflectivity and the travel time for the determination of the depth where the package is reflected. The conversion of an acoustic impedance volume to seismic amplitude volume consist of two steps. First the reflectivity is determined by computing the vertical derivative. Next, the filtering effect of the subsurface is taken into account by vertically convolving the reflectivity with the seismic wavelet that was measured at the surface. The two vertical operations of the conversion change the intensity values, but the geometry of the layered structure remains unchanged. The reflection continuity therefore corresponds directly to the bedding continuity. Coherency is a measure for the reflection continuity in a seismic image. It was first introduced in literature as a seismic attribute by Bahorich and Farmer [BF95]. Inspired by horizon tracking algorithms they used the cross-correlation between vertical lines of data for the estimation of coherency. Bahorich et al. showed that by making the continuity explicit, the delineation of both faults and channels becomes easier, especially in time slices. Time slices of seismic data are difficult to interpret because they often cut through different reflectors. Traditionally this is solved by tracking of the reflector of interest and creating 2D map view of the intensity values on created horizon. The drawback of horizon tracking 1
Here we refer to continuity in the mathematical sense. Consider a surface in 3D, z = f (x, y) = f (x). The surface is said to be continuous if limx→x0 f (x) = f (x0 ) for all x0 .
98
Chapter 6 Coherency estimation
is that this process is not yet fully automated, and therefore it remains time consuming. Due to the ever increasing computing power and size of computer memory, the computation of 3D seismic attributes at each point of 3D seismic data volume became feasible. It appeared that the interpretation of a time slice of a coherency attribute volume is much less hampered by the fact that it cuts through different reflectors. This can be explained by the fact that each point in the attribute volume is computed using a neighborhood in the seismic image around the point. The computation of a 3D attribute volume can be done ‘off-line’, i.e. without the intervention of the interpreter, and is therefore not time consuming for the interpreter. As mentioned before, reflectors are generally not parallel to the xy-planes or time slices of the seismic data. The estimation of coherency can therefore not be performed accurately without the estimation of the orientation of the reflector. This orientation estimation was incorporated into the cross-correlation algorithm by searching for the maximum crosscorrelation of a user-defined set of lagged or skewed vertical data line pairs. The drawback of this algorithm based on the cross-correlation between individual lines of data, is that its estimates for coherency and orientation are not robust with respect to noise. The estimation can be stabilized by increasing the size of the analysis window. Coherency can be defined as the resemblance of the local data to a planar reflector. Examples of mathematical resemblance measures used for the estimation of coherency are semblance [MKF98] and correlation [GM99]. A seismic image of a single constant planar reflector is just a stack of isophote planes, and it therefore has a plane-like linear structure. Since the confidence value Cplane of the GST is a measure for the resemblance of an image structure to a plane-like linear structure, it can be used as an estimate of coherency as well. In this chapter we will compare the confidence estimate of the GST with the coherency estimate based on the eigenstructure of the covariance matrix described in [GM99]. We start with the mathematical definition of both methods. Next, we will show the effect of structural dip on the coherency estimation. Two methods for the compensation of structural dip are discussed and compared. We finish this chapter with the experimental comparison of the coherency estimates applied to the detection of faults.
6.1
Coherency based on the eigenstructure of the covariance matrix
The first coherency measure was based on the cross-correlation between neighboring windowed traces, the one-dimensional time signals of a seismic image I(x = x0 , y = y0 , t). If we compute this coherency measure using a time window of N samples, we need to correlate two data vectors d1 and d2 , with N samples each. This yields the off-diagonal
6.1 Coherency based on the eigenstructure of the covariance matrix element C12 of a 2 by 2 covariance matrix PN 2 PN 1 d d C11 C12 d n1 n2 n1 n=1 PN n=1 P , C= = N 2 C21 C22 N n=1 dn1 dn2 n=1 dn2
99
(6.1)
where C11 and C22 are the auto-covariances of the data vectors d1 and d2 . The word trace can either refer to the sum of the diagonal elements of a matrix or a seismic time signal. Since the only data vector we use in this section are those corresponding to windowed seismic traces, we use ‘data vector’ as synonym for ‘windowed seismic trace’. For the generalization to the covariance matrix of J data vectors, we define the data matrix D. This JxN matrix is given by d11 d12 . . . d1J d21 d22 . . . d2J D = .. (6.2) .. .. , . . . . . . dN 1 dN 2 . . . dN J where each single column represents a single data vector with N samples. If we assume that each data vector has a zero mean, then the JxJ data covariance matrix is given by
DT · D , (6.3) N where the dot denotes the matrix product. Since the definition is based on the correlation between data vectors, which are one-dimensional signals, it can be applied to both 2D and 3D seismic images without change. Applied in 3D using a rectangular analysis window with sizes (wx , wy , wt ), we have N = wt and J = wx · wy . C=
If all J data vectors in matrix D are identical, which is the case for a flat planar reflector without noise, then the rank of C is one. This means that C has only one non-zero eigenvalue. An increase in the variability of the traces, e.g. due to a discontinuity, causes an increase in the number of non-zero eigenvalues. The resemblance of the local data to a flat continuous planar reflector can therefore be measured by the fraction of the total signal energy that is captured by the largest eigenvalue of C. It can easily be checked that the total signal energy E, is equal to the sum of the diagonal elements of C, known as the trace T r(C). From linear algebra [GvL96], we know that the trace of a matrix is invariant under orthonormal transformations. Therefore we have N J J J X X X X 2 dnj = cjj = T r(C) = λj , (6.4) E= j=1 n=1
j=1
j=1
with dnj , cjj the elements of D, C and λj the eigenvalues of C. If λ1 is the largest eigenvalue then the coherency estimate based on the eigenstructure of the covariance matrix is given by λ1 . (6.5) cˆcov = T r(C) This measure was presented as an estimate of seismic coherency in [GM99].
100
6.2
Chapter 6 Coherency estimation
Coherency estimation using the GST
An image of a single planar reflector can be described as a stack of isophote planes. The confidence value Cplane of the gradient structure tensor (GST) can therefore be used as a measure of the resemblance of the local data to a planar reflector. The GST was presented and discussed in detail in chapter 2. In this section we will discuss the GST only in relation with coherency estimation. The gradient structure tensor T is defined as the averaged dyadic product of the gradients g T = ggT .
(6.6)
For most practical purposes, the tensor can be treated as a N xN matrix, with N the dimensionality of the data. The eigenvalues of this tensor indicate the gradient energy in the orientations defined by the corresponding eigenvectors. In the case of a planar reflector the tensor has only one non-zero eigenvalue, and the corresponding eigenvector is the normal vector of the reflector. Any deviation of the data from a constant planar reflector leads to an increase of the gradient energy in the lateral direction. The coherency of the GST could therefore be estimated by cˆgst =
λ1 T r(T)
(6.7)
Although the definitions of cˆcov and cˆgst seem identical, they are not. The eigenvalues of the covariance matrix represent the correlation between seismic traces and the eigenvalues the GST represent the gradient energies of a geometrically ordered set of traces. This means that the reflector continuity is measured with the correlation between traces by the first method, while the second method uses the gradient energy in the lateral direction as a measurement for continuity. For simplicity, we will now consider only two traces (t1 , t2 ), both with a zero mean value X X t1i = 0 , t2i = 0 (6.8) i
i
The normalized correlation C between the traces is defined as C=
t1 · t2 . kt1 kkt2 k
(6.9)
The gradient in the lateral direction can be defined as the difference between the traces. The corresponding lateral gradient energy El is in that case given by El = kt1 − t2 k2
(6.10)
If the gradient energy El is zero then the correlation C is one. However, if the normalized correlation between the traces is one, then the gradient energy is not necessarily zero. In practice this means that for a perfectly flat reflector with an increasing amplitude cˆcov = 1, but cˆgst < 1.
6.3 The presence of structural dip
6.3
101
The presence of structural dip
The coherency estimate cˆcov is based on the correlation of windowed traces. The optimal correlation measurement is obtained when the correlation is computed between 1D subsets perpendicular to the reflector surface. The estimate cˆcov is therefore optimal if the reflector being analyzed is parallel to the xy-plane, and becomes less accurate in the presence of structural dip. Here structural means, consistent at a larger scale not determined by random fluctuations. The compensation of structural dip requires the estimation of this dip and the coherency estimate cˆcov should be computed after transforming the data in such a way that the reflector becomes horizontal. In figure 6.1 a tilted reflector is drawn. Window A in this figure is the xt-aligned window of cˆcov , and B is the effective window after dip compensation.
t
B A
x
Figure 6.1: Coherency estimation in the presence of structural dip in an xt-aligned analysis window (A) yields an sub-optimal result. The window can be adapted (B), using an estimation of the structural dip.
To demonstrate the negative effect of structural dip on the estimate cˆcov , we applied it to a 2D seismic image (128*128 pixel) with a lot of faults, see figure 6.2. We next estimate the orientation of the reflectors using the GST at scale (σg = 1, σt = 5). We then use this orientation estimate to locally rotate the data in such a way that the reflectors become horizontal before estimating the coherency. The result of this adaptive version of cˆcov is shown in figure 6.2e. A 5 by 11 pixel coherency window was used in both cases. For the non-adaptive windows this means that the coherency was computed over 5 traces using a time window of 11 pixels, and the widths of the adaptive windows were 5 pixels laterally and 11 pixels in the perpendicular direction. The results show a clear improvement due to the dip compensation. To show the difference between cˆcov and the semblance based coherency estimate presented in [MKF98], we have repeated the experiment using the semblance based coherency. The semblance based coherency estimate is given by 1 1 1 u·C·u (6.11) , u = √ .. , cˆsem = T r(C) J . 1
102
Chapter 6 Coherency estimation
(a) 2D seismic image
(b) Covariance
(c) Semblance
(d) GST Orientation
(e) Adaptive Covariance
(f) Adaptive Semblance
Figure 6.2: Coherency estimation of a 2D seismic image using the eigenstructure of the covariance matrix (b) and semblance (c) in an xt-aligned window. The orientation estimate (d) of the GST is used to compensate for the structural dip (e,f).
6.3 The presence of structural dip
103
with C the J by J covariance matrix defined in (6.3). The results of cˆsem and its orientation adaptive version are shown in figures 6.2c and f. The results show that the faults highlighted by cˆcov are highlighted by cˆsem as well. However, the eigenstructure based coherency estimate has a narrower fault response. A more detailed experimental comparison is given in [MSG+ 99] and a theoretical comparison is given in [GM99]. The choice of the window size for the coherency and dip estimation, is a trade-off due to the uncertainty principle in filtering [WG84, Dau85]. A small window yields an accurate, high resolution measurement, but is sensitive to noise. A larger window yields a more robust measurement with respect to noise, but with a lower resolution. The optimal window size for the coherency measurement is not necessarily the optimal window size for the dip measurement.
6.3.1
Dip estimation by sampling the dip dependency
We have demonstrated that compensating for structural dip significantly improves the coherency estimation. We used the GST to estimate the orientation of the reflector, and rotated the data before estimating coherency. The dip compensation for the semblance based method [MKF98] was done by a straight forward dip search. The idea of this dip search is to compute the coherency for all possible dips, and the trial dip that yields the highest coherency value corresponds to the actual dip of the reflector. The coherency as a function of both the position (x, y, t) and the apparent dips (p, q) is defined as c(x, y, t, p, q) = c(x, y, t − px − qy). The estimate of the reflector dip d is given by p dˆ = pˆ2 + qˆ2 , (ˆ p, qˆ) = arg max c(x, y, t, p, q). p,q
(6.12)
(6.13)
The sample spacing of the apparent dip parameters ∆p and ∆q, depends on the sizes (wx , wy ) of the coherency window and the highest temporal frequency component fmax contained in the seismic data. The Nyquist sampling criterion gives the following expressions 1 1 ∆p = , ∆q = . (6.14) 2wx fmax 2wy fmax In [MKF98] they assume that the interpreter is able to estimate the maximum true dip, dmax , to reduce the range of the dip search. Using this maximum dip, c(x, y, t, p, q) has to be computed over np · nq discrete dip pairs (p, q), where np =
2dmax 2dmax + 1 , nq = + 1. ∆p ∆q
(6.15)
For example, if the image is critically sampled 1/fmax ≈ 2 pixels, the maximum dip dmax = 1 pixel (45 degrees), and the window sizes are wx = wy = 5, then c(x, y, t, p, q) needs to computed over 121 dips.
104
Chapter 6 Coherency estimation
In [MSG+ 99], Marfurt et al. applied this dip compensation to the coherency estimate based on the eigenstructure of the covariance matrix. They found that the high fault discrimination and resolution of this algorithm reverted to that of the semblance based method. For high resolution fault highlighting a small lateral window size is needed, e.g. wx = wy = 5. Estimating the dip near a fault in such a small window does not yield the desired result. The apparent dips that give the maximum coherency value appeared to compensate the offset of the fault rather then the structural dip. The structural dip needs to be estimated at a larger scale.
t
2
1
x
Figure 6.3: Dip estimation near a fault at two scales (1,2), indicated by the gray discs. The white rectangles indicate the effective windows that yield the maximum coherency estimate c(x, y, t, p, q).
A schematic representation of dip estimation at two scales is depicted in figure 6.3. The gray discs indicate the desired scales for the coherency estimation (1), and the dip estimation (2). The white rectangles indicates the effective windows that yield the maximum coherency estimate c(x, y, t, p, q).
6.3.2
Comparison between the dip estimates of the GST and dip search
In this section we will compare the dip estimates of the GST dˆgst and of the dip search method dˆds . The GST does not yield an estimation for dip directly. The first eigenvector, belonging to the largest eigenvalue, of the GST gives an estimate of the normal of the reflector plane. The dip can be derived from this vector by φ = arccos
π e1 · ez , dˆgst = tan( − φ), ke1 kkez k 2
(6.16)
where e1 is the first eigenvector of the GST, ez is a vector along the z-axis, and φ is the angle between e1 and ez . The differences in the key features of the GST and dip search analysis are summarized in table 6.1, and will be discussed in descending order.
6.3 The presence of structural dip
105
GST uni-modal rotation invariant gradient based
Dip search multi-modal skew invariant covariance based
Table 6.1: Key features of the dip analysis of the GST and dip search.
For the determination of the modality of the methods, we can interpret both methods as the measurement of the resemblance R(d) between the local data and a planar reflector as a function of the dip2 . The dip search method samples the dip parameter di , and measures R for all instances i in the search range. In this way the distribution of R as a function of d is measured allowing multiple peaks or maxima as a function of the parameter at a single position in the image. This is per definition a multi-modal parameter estimation. The orientation analysis of the GST on the other hand, is uni-modal. It implicitly assumes that there is only one orientation at a single position in the image, see chapter 2. In figure 6.4 we have drawn three reflector configurations: a single reflector (a), two touching reflectors (b), and two intersecting reflectors (c). A multi-modal analysis is able to represent all three configurations, whereas an uni-modal analysis is limited to the case of the single reflector. However, in chapter 2 we introduced the generalized Kuwahara filter to avoid filtering across borders. Combined with the GST, this filter yields a correct orientation estimate of touching reflectors as well. Since intersecting reflectors have no geological meaning, it is not necessary to deploy a multi-modal analysis for the estimation of reflector orientation or dip.
(a) single
(b) touching
(c) intersecting
Figure 6.4: Three different reflector configurations.
The methods also differ in their invariance. The GST is rotation invariant and the dip search is skew invariant. The parameter estimate of the GST, orientation, indicates the rotation needed to level the reflector. The dip estimate is the amount of skewing necessary to level the reflector. The difference between a rotated and a skewed analysis window increases with the dip. In figure 6.5, we have drawn a skewed and a rotated window on top of each other for two different angles: 20 and 45 degrees. The advantage of skewing is that if wx = wy < wt the spatial and temporal window sizes do not dependent on the amount of skewing. Furthermore, transformation of the data based on skewing only requires a one-dimensional interpolation along the t-axis. Rotation on the other hand 2
For simplicity we assume that the azimuth is known. In practice the estimation of the reflector dip requires the simultaneous estimation of two parameters.
106
Chapter 6 Coherency estimation
treats all dimensions on equal footing. The advantage of rotation is that the shape of the window is unaffected. If an isotropic window is used, such as the Gaussian window used for the computation of the GST, then the part of the image data inside the analysis window is identical for all orientations.
Figure 6.5: The difference between rotated and skewed analysis windows, for 20 (left) and 45 (right) degree angles.
The third difference between the two dip estimators is the resemblance measure. The difference of gradient versus covariance is already discussed in section 6.2. The experimental comparison between the GST and the covariance based coherency measures will be performed in the next section. To see the practical implications of the differences between the dip-estimators, we have applied both the GST and the dip search using cˆcov to the dip estimation on a 2D seismic image. We then used the two dip estimates for the computation of the dip-compensated coherency based on the covariance eigenstructure and the results are shown in figure 6.6. The dip search and all coherency measurements are performed in a wx = 5 by wt = 11 window. The GST is applied at scale (σg = 1, σT = 5). The dip search estimate dˆds and its corresponding coherency measurement are shown in figures 6.6a and d. To increase the effective scale of the dip estimate dˆds , we filtered it using a isotropic Gaussian filter (σ = 5), see figures 6.6b. Contrary to orientation, dip is not periodic and may therefore be smoothed directly without the need of a mapping. If we analyze the results in figure 6.6, we see that the dip estimate dˆds computed at a small scale (a), shows rapid changes near the faults. As depicted in figure 6.3, the throw of the fault influences the dip estimate. Smoothing of the dip estimate dˆds (b) practically removes this influence. A comparison of the corresponding coherency estimates (d,e) shows that increasing the scale of the dip estimation improves both the discrimination and the resolution of the fault-highlighting. The coherency estimate using GST dip dˆgst (f), does not visibly differ from the result using the smoothed dip (e). In (h) we have shown the result of rotating the data to compensate the reflector dip using GST orientation (i). If we compare (f) with (h), we see that for the two faults in the right half of the image, rotation yields a more continuous fault-highlighting. The discrimination of the two faults in the bottom-left corner is less in (h).
6.3 The presence of structural dip
107
(a) Dip search
(b) Smoothed
(c) GST dip
(d)
(e)
(f)
(g) 2D seismic image
(h) OA Coherency
(i) GST Orientation
Figure 6.6: Dip estimation on a 2D seismic image (g) using both the dip search (a) as the GST (c). A Gaussian smoothed (σ = 5) version of (a) is shown in (b). Dip adaptive coherency estimates using the dip estimates (a,b,c) are shown in (d,e,f) respectively. The GST orientation and corresponding adaptive coherency are shown in (h,i) for comparison.
108
6.4
Chapter 6 Coherency estimation
An experimental comparison for fault detection
In the sections above we discussed the coherency estimation cˆcov based on the eigenstructure of the covariance matrix [GM99] and we showed that the eigenvalues of the GST can be used for the estimation of coherency as well. In chapter 2 we described the application of the GST to the detection of faults. In section 2.4, we showed that a combination of the eigenvalues of the GST, Cf ault , is able to highlight faults. We used this highlighting as the basis for the detection of faults in 3D seismic data. In this section we will compare Cf ault with the coherency estimate cˆcov . Furthermore, we will compare the automatic fault detection based on these to measures, by applying the same post-processing steps to cˆcov as we applied in the fault detection application. In the previous section we have addressed the effects of structural dip on the coherency estimation. We have shown that the orientation estimate of the GST can be used to compensate for this dip, and we will use this dip compensation for all coherency measurements in this section. We perform the fault analysis on the same 3D seismic image (1283 voxels) as we used in section 2.4. We start by highlighting the faults using both Cf ault and cˆcov . The results on the 3D image are visualized in figure 6.7 by three perpendicular 2D cross-sections, namely a xt, a yt, and a xy-slice or time-slice. To allow easy visual comparison between cˆcov and Cf ault , we have displayed (1 − Cf ault ). The coherency estimate cˆcov was computed in a (21 ∗ 5 ∗ 5) window, and Cf ault was computed at scale (σ1 = 6, σ2 = σ3 = 2). The elongated windows are steered perpendicular to the reflector planes using the orientation of the first eigenvector of the GST at scale (σg = 1, σT = 9). Visual comparison of the results in figure 6.7, shows that the GST based highlighting is smoother than the covariance based highlighting. This could mean that the resolution of cˆcov is higher. With the term resolution we refer to the ability to resolve faults that are in close proximity of each other. A lower resolution does not necessarily result in a worse localization of the faults. The two closest faults visible in figure 6.7 in the bottom-left quadrant of the time slice (g,h,i), are still resolved by the GST based highlighting. Therefore we are not able to confirm the lower resolution. With respect to the discrimination between faults and continuous reflectors, both methods perform equally well. In some cases Cf ault performs better, in some cases cˆcov . We next use the fault highlighting shown in figure 6.7 for the segmentation faults in the seismic image as described in 2.4. The ‘faultiness’ of a position in the image is given by Cf ault and (1−ˆ ccov ). First we apply a non-maximum suppresion. The candidate fault-points are reduced to the local maxima in the faultiness images. These maxima are computed in a one-dimensional window (5 ∗ 1 ∗ 1) along the gradient direction. We next take the value of the faultiness of the candidate points into account using two threshold values, Tl = 0.1, Th = 0.4. Points with a faultiness below the low threshold Tl value are rejected. The final segmentation is obtained by accepting only those connected segments whose maximum value is above the high threshold value Th . The resulting fault detections based on both the GST and the covariance analysis are shown in figure 6.8. The covariance based
6.4 An experimental comparison for fault detection
109
detection is more sensitive to small variations in the data than the GST based detection. This is consistent with the smoother fault highlighting of the GST.
110
Chapter 6 Coherency estimation
(a) Slice y = 34
(b) GST
(c) Covariance
(d) Slice x = 113
(e) GST
(f) Covariance
(g) Slice t = 64
(h) GST
(i) Covariance
Figure 6.7: Coherency estimation on a 3D seismic image containing many faults, visualized in three 2D cross-sections, a xt (a),a yt (d), and a xy-slice or time-slice (g). The GST based coherency was computed at scale (σ 1 = 6, σ2 = σ3 = 2), and the covariance based coherency in a window with sizes (w 1 = 21, w2 = w3 = 5). The elongated windows are steered perpendicular to the reflector planes using the first eigenvector of the GST applied at scale (σg = 1, σT = 9).
6.4 An experimental comparison for fault detection
111
(a) Slice y = 34
(b) GST
(c) Covariance
(d) Slice x = 113
(e) GST
(f) Covariance
(g) Slice t = 64
(h) GST
(i) Covariance
Figure 6.8: Fault detection based on the fault highlighting shown in figure 6.7.
7. Conclusions Is it possible to improve the efficiency of the seismic interpretation process by automating certain interpretation tasks? In this thesis an image processing approach to this question is presented. We limited the broad range of seismic interpretation tasks to horizon tracking, fault detection, and channel detection in 2D and 3D migrated seismic reflection images. The global structure of a seismic image of sedimentary rock is determined by reflector surfaces stacked on top of each other. The interpretation tasks mentioned above all require an estimation of the global structure. We found that the faults and the channels manifest themselves as deviations from the global structure. Futhermore, automatic horizon tracking could be improved by enhancing the global structure before tracking. The first step towards the automation of the interpretation tasks therefore is to model the global structure in the data and to find or create image processing algorithms for the estimation of the modelparameters. In addition to the inherent difficulty of finding the most effective structure model, we face a practical problem regarding the vast size of seismic datasets1 . This demands for a minimization of the computational complexity of the image processing algorithms used.
7.1
Image processing approach
The image processing approach chosen in this thesis is the local analysis of the image structures. We have presented several geometric structure models in both two and threedimensions. We started the structural analysis of images by assuming local shift invariance and we called an image structure with shift invariance in one or more orientations, a linear structure. The gradient structure tensor (GST) locally models an image as a linear structure. This is equivalent to describing the image structure as locally straight (line-like structures) or planar (plane-like structures). The GST yields estimates of both the parameters, orientation, and the confidence of the linear structure model. The confidence estimate decreases at image locations where the neighborhood deviates from the linear structure model. An 1
A seismic dataset can take up to a couple of Terabytes of computer storage space. During the time this thesis was written, an average workstation contains a storage capacity smaller than 0.1 Terabyte and approximately 2 Gigabytes of memory.
114
Chapter 7 Conclusions
example of such a neighborhood is the area around a border between two oriented textures. The confidence estimate can therefore be used to these highlight borders. The generalized Kuwahara filter avoids filtering across borders by allowing decentralized windows and selecting the decentralized window with the highest confidence value. The uni-modal orientation estimate of the GST gives near orientation borders a weighed average of the two orientations on both sides of the border. We have shown that this orientation estimate can be corrected by combining the GST with the generalized Kuwahara filter. Compared to a multi-orientation approach, this method has a low computational complexity. The corrected orientation estimate makes it possible to correctly steer an orientation adaptive filter over images composed of several touching oriented textures. An orientation adaptive filter can be made edge or border preserving by combining it with a steered version of the generalized Kuwahara filter. This edge preserving orientation adaptive filter gives results that are visually comparable to the results of anisotropic diffusion, but is at least one order of magnitude faster. As mentioned before, the confidence estimate of the GST decreases at image locations where the local structure deviates from the linear structure model. A frequently occurring deviation in linear structures is curvature. For the description of line-like and plane-like curvilinear structures, we introduced the curvature corrected GST based on a quadratic curve model and the quadratic surface model respectively. Both models yield a curvature corrected confidence estimate and an estimate of the curvature as a valuable by-product. We have used the orientation and curvature estimate of the corrected GST to control an adaptive filter for the reduction of noise in curved oriented textures. The curvature and orientation adaptive filter (COAF) yields both a visibly better result and a higher signal-to-noise ratio improvement than the orientation adaptive filter (OAF). Describing image structures, such as the meandering structure of a channel, at a ‘larger’ scale with a simple geometric model, is often not possible. Instead of introducing a more complex parametric model with more parameters, we turned to a non-parametric model. The local parameter estimates obtained from the simple geometric models can be used for a piece-wise description of curvilinear image structures. The shape of the structure can be found by the tracking of this piece-wise description. We have applied the corresponding non-parametric adaptive filter (NPAF) to the reduction of noise in curved oriented textures. The NPAF yields a signal-to-noise ratio improvement that is higher than that of the COAF, and approaches the improvement of the theoretically predicted ideal filter.
7.2
Application to seismic interpretation
The parameter estimates of the local geometric models can control the window shape of a filter in such a way that it approximates the shape of the local structure. A low-pass filter controlled in this way can enhance the structure of a layered texture, and thus the structural
7.2 Application to seismic interpretation
115
information in a seismic image. Automatic horizon trackers have been around for years and are available in almost every seismic interpretation software package. The problem however is that these trackers fail at many points in the image due to small deviations from the layered structure. Accurate horizon tracking still requires a lot of user interaction in the form of providing many starting points for the auto-tracker. The structural information in seismic images can be enhanced by removing the deviations using an adaptive smoothing filter controlled by the gradient structure tensor. We have shown that by combining the adaptive smoothing filter with the generalized Kuwahara filter, smoothing across faults can be avoided. The confidence estimate of the local straight model can be used for the detection of faults. The automatic detection of faults is presented as a two-step procedure. First an attribute is computed that measures the ‘faultiness’ at each point in the image. In the second step this faultiness is used for the segmentation of the faults. We have shown that a faultiness measure Cf ault can be constructed from the eigenvalues of the GST. We compared this measure with the coherency estimate cˆcov based on the eigenstructure of the covariance matrix [MSG+ 99]. We found that Cf ault yields a smoother fault-highlighting than cˆcov . For the manual interpretation of the fault-highlighting the more crisp fault responses of cˆcov are preferred, furthermore the erratic false responses of cˆcov are easily suppressed by the human visual system. The automatic detection of faults, on the other hand, is not hampered by a smoother fault-highlighting as long as all the faults are spatially resolved. Moreover, the automatic discrimination between small erratic responses and true fault responses requires an additional image processing step. This increases the computational demands of the detection based on the already more computational complex coherency estimate cˆcov . The automatic detection of channels appeared to be a more difficult task than the automatic detection of faults. Due to the large natural variability in the shape of channels, it is not possible to capture their characteristics in a simple geometric model. We have presented a channel detection method based on the confidence estimate of the structure tensor for 3D line-like curvilinear structures. This method was able to detect most parts of the target channels. However, the detection was fragmented and many other sedimentary structures were falsely detected as channels as well. The visually most distinguishing feature between the channels and the other sedimentary structures is the larger spatial continuity of the channels. we have shown that the exploitation of this continuity by elongating the analysis window along the channel using a piece-wise description of the channel shape, results in a more continuous detection and a significant reduction of the number of false-positive detections.
116
Chapter 7 Conclusions
Bibliography [BF95]
M. Bahorich and S. Farmer. The coherency cube. Leading edge, pages 1053– 1058, October 1995.
[Bis95]
C.M. Bishop. Neural Networks for Pattern Recognition. Oxford University Press, 1995.
[BSMM00]
I.N. Bronstein, K.A. Semendjajew, G. Musiol, and H. M¨ uhlig. Taschenbuch der Mathematik. Verlag Harri Deutsch, 2000.
[BVV99]
P. Bakker, P.W. Verbeek, and L.J. van Vliet. Edge preserving orientation adaptive filtering. In CVPR’99 (Colorado, U.S.A), volume 2, pages 535–540, June 1999.
[BVV01]
P. Bakker, P.W. Verbeek, and L.J. van Vliet. Confidence and curvature estimation of curvilinear structures in 3-d. In ICCV’01 (Vancouver, Canada), volume II, pages 139–144. IEEE, July 9-12 2001.
[Can86]
J. Canny. A computational approach to edge detection. IEEE transactions on Pattern Analysis and Machine Intelligence, 8(6):679–697, November 1986.
[CH89]
Y.S. Chen and W.H. Hsu. An interpretive model of line continuation in human visual perception. Pattern Recognition, 22(5):619–639, 1989.
[CST00]
J. Chen, Y. Sato, and S. Tamura. Orientation space filtering for multiple orientation line segmentation. Pattern Analysis and Machine Intelligence, 22(5):417–429, May 2000.
[Dau85]
J.G. Daugman. Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters. Journal of the Optical Society of America A, 2(7):1160–1169, July 1985.
[Des97]
F.J. Dessing. A wavelet transform approach to seismic processing. PhD thesis, Delft University of Technology, 1997.
[DGS+ 89]
R.M. Dalley, E.C.A. Gevers, G.M. Stampfli, D.J. Davies, C.N. Gastaldi, P.A. Ruitjenberg, and G.J.O. Vermeer. Dip and azimuth displays for 3d seismic interpretation. First Break, 7(3):86–95, 1989.
118
Bibliography
[FA91]
W.T. Freeman and E.H. Adelson. The design and use of steerable filters. IEEE transactions on Pattern Analysis and Machine Intelligence, 13(9):891– 906, September 1991.
[FIS89]
S.H. Friedberg, A.J. Insel, and L.E. Spence. Linear Algebra. Prentice-Hall, 1989.
[FJ97]
M. Frigo and S.G. Johnson. The fastest fourier transform in the west. Technical report, Massachusetts Institute of Technology, USA, September 1997. MIT-LCS-TR-728.
[Fre92]
W.T. Freeman. Steerable Filters and Local Analysis of Image Structure. PhD thesis, Massachusetts Institute of Technology, 1992.
[FS99]
B. Fischl and E.L. Schwartz. Adaptive nonlocal filtering: A fast alternative to anisotropic diffusion for image enhancement. IEEE Trans. on PAMI, VOL. 21(NO. 1), 1999.
[GK95]
G.H. Granlund and H. Knutsson. Signal processing for computer vision. Kluwer Academic Publishers, Boston/Dordrecht/London, 1995.
[GM99]
A. Gersztenkorn and J.K. Marfurt. Eigenstructure-based coherency computations as an aid to 3-d structural and stratigraphic mapping. Geophysics, 64(5):1468–1479, October 1999.
[Gol80]
H. Goldstein. Classical Mechanics. Addison Wesley, 1980.
[Gra78]
G.H. Granlund. In search of a general picture processing operator. Computer Graphics and Image Processing, 8:155–173, 1978.
[GvL96]
G.H. Golub and C.F. van Loan. Matrix Computations. John Hopkins University Press, Baltimore, 3 edition, 1996.
[GVV97]
M. van Ginkel, P.W. Verbeek, and L.J. van Vliet. Orientation selectivity for orientation estimation. In Proceedings of the 10th Scandinavian Conference on Image Analysis (Lappeenranta, Finland), volume I, pages 533–537, June 9-11 1997.
[GWVV99]
M. van Ginkel, J. van de Weijer, P.W. Verbeek, and L.J. van Vliet. Curvature estimation from orientation fields. In B.K. Ersboll and P. Johansen, editors, SCIA’99, Proc. 11th Scandinavian Conference on Image Analysis (Kangerlussuaq, Greenland), pages 545–551. Pattern Recognition Society of Denmark, Lyngby, June 7-11 1999.
[Hag92]
L. Haglund. Adaptive Multidimensional Filtering. PhD thesis, Link¨oping University, Link¨oping, Sweden, 1992.
[KBV+ 99]
G.M.P. van Kempen, N. van den Brink, L.J. van Vliet, M. van Ginkel, and P.W. Verbeek. The application of a local dimensionality estimator to
Bibliography
119 the analysis of 3d microscopic network structures. In SCIA’99, Proc. 11th Scandinavian Conference on Image Analysis (Kangerlussuaq, Greenland), pages 447–455, June 7-11 1999.
[Kea76]
M. Kuwahara and et al. Processing of ri-angiocardiographic images. Digital Processing of Biomedical Images, pages p. 187–203, 1976.
[KHRV99]
S.N. Kalitzin, B. ter Haar Romeny, and M.A. Viergever. Invertible apertured orientation filters in image analysis. International Journal of Computer Vision, 31(2/3):145–158, 1999.
[Knu85]
H. Knutsson. Producing a continuous and distance preserving 5-d vector representation of 3-d orientation. In IEEE Computer Society Workshop on Computer Architecture for Pattern Analysis and Image Database Management, pages 175–182. Miami Beach, Florida, November 18-20 1985.
[Knu89]
H. Knutsson. Representing local structure using tensors. In The 6th Scandinavian Conference in Image Analysis, pages 244–251, June 19-22 1989.
[KW87]
M. Kass and A. Witkin. Analyzing oriented patterns. Computer Vision, Graphics and Image Processing, 37:362–385, 1987.
[Mar97]
J.B. Martens. Local orientation analysis in images by means of the hermite transform. IEEE transactions on Image Processing, 6(8):1103–1116, August 1997.
[MK00]
J.K. Marfurt and R.L Kirlin. 3-d broad-band estimates of reflector dip and amplitude. Geophysics, 65(1):304–320, January 2000.
[MKF98]
J.K. Marfurt, R.L Kirlin, and S. Farmer. 3-d attributes using a semblancebased coherency algorithm. Geophysics, 63(4):1150–1165, July 1998.
[MSG+ 99]
J.K. Marfurt, V. Sudhaker, A. Gersztenkorn, K.D. Crawford, and S.E. Nissen. Coherency calculations in the presence of structural dip. Geophysics, 64(1):104–111, January 1999.
[NM79]
M. Nagao and T. Matsuyama. Edge preserving smoothing. Computer Graphics and Image Processing, 9:394–407, 1979.
[Per95]
P. Perona. Deformable models for early vision. IEEE transactions on Pattern Analysis and Machine Intelligence, 17(5):488–499, May 1995.
[PM90]
P. Perona and J. Malik. Scale-space and edge detection using anisotropic diffusion. IEEE-PAMI, 12(7):629–639, July 1990.
[Pra72]
W.K. Pratt. Generalized Wiener filtering computation techniques. IEEE Trans. Comput., C-21:297–303, 1972.
[RS91]
A.R. Rao and B.G. Schunck. Computing oriented texture fields. CVGIP: Graphical Models and Image Processing, 53(2):157–185, March 1991.
120
Bibliography
[SF96]
E.P. Simoncelli and H. Farid. Steerable wedge filters for local orientation analysis. IEEE transactions on Image Processing, 5(9):1377–1382, 1996.
[Sha49]
C.E. Shannon. Communication in the presence of noise. Proceedings of the IRE, 29(4):10–21, 1949.
[Spi79]
M. Spivak. A Comprehensive Introduction to Differential Geometry, volume 2. Publish or Perish, Berkeley, 1979.
[Ste97]
P. Steeghs. Local power spectra and seismic interpretation. PhD thesis, Delft University of Technology, 1997.
[TKS79]
T.H. Tanner, F. Koehler, and R.E. Sherrif. Complex seismic trace analysis. Geophysics, 44, 1979.
[VvVvdW98] P.W. Verbeek, L.J. van Vliet, and J. van de Weijer. Improved curvature and anisotropy estimation for curved line bundles. In A.K. Jain, S. Venkatesh, and B.C. Lovell, editors, ICPR’98, Proc. 14th Int. Conference on Pattern Recognition, pages 528–533. IEEE Computer Society Press, Los Alamitos, 17-20 August 1998. [Wal87]
D. Walters. Selection of image primitives for general-purpose visual processing. Computer Vision, Graphics, and Image Processing, 37:261–298, 1987.
[Wei98]
J. Weickert. Anisotropic diffusion in image processing. Stuttgart, 1998.
[WG84]
R. Wolson and G.H. Granlund. The uncertainty principle in image processing. IEEE transactions on Pattern Analysis and Machine Intelligence, 6(6), November 1984.
[Wie49]
N. Wiener. Extrapolation, interpolation and smoothing of stationary time series. Cambridge, MA : M.I.T. Press, 1949.
[WVVG01]
J. van de Weijer, L.J. van Vliet, P.W. Verbeek, and M. van Ginkel. Curvature estimation in oriented patterns using curvlinear models applied to gradient vector fields. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2001. in press.
[YGV98]
I.T. Young, J.J. Gerbrands, and L.J. van Vliet. Image processing fundamentals. In V.K. Madisetti and D.B. Williams, editors, The Digital Signal Processing handbook, chapter 51. IEEE Press and CRC Press, 1998.
B.G. Teubner
Summary In this thesis we have examined the possibility of automating 3-D seismic interpretation, by making use of image processing techniques. We limited the broad range of seismic interpretation tasks to horizon tracking, fault detection, and channel detection in 2D and 3D migrated seismic reflection images. All these tasks require a geometrical description of the earth layers, i.e. a quantitative structural interpretation. Earth strata manifest themselves in seismic images as layered structures. For the quantitative description of these structures we have decided to use a local differential geometric description by means of the structure tensor. The structure tensor models the data as locally straight and yields both an estimate of the orientation of the structure as a confidence measure for the model. This confidence measure can be used to highlight and detect faults that manifest themselves as discontinuities. We have compared our discontinuity based fault measure both theoretically and experimentally with a seismic coherency measure from recent literature. The orientation estimate can be used to steer an adaptive filter. This steered filter can enhance the structural information in the image by reducing the noise in seismic image without degrading the layered structure. By combining this filter with the generalized Kuwahara filter, it can be made fault preserving as well. Furthermore, the fault detection can be improved by elongating the tensor window along the faults. Since curved structures occur often in practice, we have introduced the curvature corrected structure tensor. This corrected tensor yields an estimate of the curvature as a valuable by-product. We have described the curvature correction to both line-like and plane-like curvilinear structures. The confidence measure of the line-like curvilinear model can be used to detect channels. To make better use of the continuity of channels, we have introduced a non-parametric structure description method. This method is able to describe complex structures such as channels on a large scale without the need for additional geometric parameters.
Samenvatting In dit proefschrift is onderzocht of het mogelijk is om, met behulp van beeldbewerkingstechnieken, driedimensionale seismische interpretatie te automatiseren. De interpretatietaken waar we ons in eerste instantie tot beperkt hebben, zijn het vinden van horizon- en breukvlakken en het vinden van geulen. Voor deze taken is een geometrische beschrijving van de aardlagen nodig, i.e. een structurele interpretatie. De aardlagen manifesteren zich in seismische beelden als gelaagde structuren. Voor een kwantitatieve beschrijving van deze gelaagde structuren, hebben we gekozen voor een lokaal differentieel geometrische beschrijving met behulp van de structuurtensor. De structuurtensor modelleert de data als lokaal recht en levert zowel een schatting van de ori¨entatie van de structuur als een betrouwbaarheidsmaat van dit model. Deze betrouwbaarheidsmaat kan gebruikt worden om breuken die zich manifesteren als discontinuiteiten te accentueren en te detecteren. We hebben onze op discontinuiteit gebaseerde breukmaat theoretisch en experimenteel vergeleken met een coherentiemaat uit de recente literatuur. De ori¨entatieschatting kan gebruikt worden om een adaptief filter te sturen. Dit gestuurde filter is in staat de structurele informatie in het beeld te versterken, door de ruis in seismische beelden te reduceren zonder de gelaagde structuur aan te tasten. Door dit filter te combineren met het gegeneraliseerde Kuwahara-filter, blijven ook de breuken intact. Verder kan de breukdetectie verbetert worden door het tensorvenster langs de breuk te verlengen. Omdat in de praktijk gekromde structuren vaak voorkomen, hebben we de curvatuurgecorrigeerde structuurtensor geintroduceerd. Deze gecorrigeerde structuurtensor levert een schatting van de curvatuur als waardevol bijproduct. We hebben de curvatuur correctie zowel voor vlak-achtige als lijn-achtige structuren beschreven. De betrouwbaarheidsmaat van het lijn-achtige curvilineaire model kan gebruikt worden voor het detecteren van geulen. Om beter gebruik te maken van de continuiteit van geulen, hebben we een nietparametrische structuurbeschrijvingsmethode geintroduceerd. Deze methode kan zonder toevoeging van parameters, complexe structuren zoals geulen op een grote schaal beschrijven.
Dankwoord Het leven als AIO, in Delft omgedoopt tot promovendus, is voor mij een heel interessante ”reis” geweest. De eerste drie¨eneenhalf jaar zijn voor mijn gevoel echt voorbij gevlogen. Toen kwam het onvermijdelijke moment in het leven van een promovendus: het schrijven van het proefschrift. De metamorfose van experimentator tot schrijver. In de praktijk betekent dit vaak een half jaar (of langer...) de laatste loodjes. Ik wil graag iedereen bedanken die direct en indirect heeft bijgedragen aan het voltooien van deze reis. Het werk wat ten grondslag ligt aan dit proefschrift is uitgevoerd in een samenwerkingsverband tussen TU Delft, TNO-TPD en Shell. Ik wil Gert van Antwerpen (TNO-TPD), Martin Kraaijveld (Shell), Gijs Fehmers (Shell) en de andere betrokkenen bedanken voor alle discussies. Speciaal wil ik bedanken mijn begeleiders Piet Verbeek en Lucas van Vliet. Jullie deur staat altijd open en jullie hebben mij tijdens deze hele reis enorm geholpen met al jullie goede idee¨en en adviezen. Verder wil ik mijn kamergenoten Michael, Judith, en Bernd bedanken voor alle wetenschappelijke en vooral ook de minder wetenschappelijke discussies. Jullie hebben samen met alle ph-aio’s ervoor gezorgd dat het niet alleen een productieve maar ook een gezellige reis was. De vele TPKV en de ASCI avonden, de ‘vlak na de lunch’ recreatieve aktiviteiten en de squashpartijen hebben hieraan zeker bijgedragen. Mijn ouders wil ik bedanken voor hun onophoudelijke steun en vertrouwen. Tenslotte, ‘last but definitely not least’, Jeanine, bedankt voor alles.
Curriculum vitae Peter Bakker was born in Linz on December 8, 1974. In 1993 he obtained the VWO diploma at the Bonaventura college in Leiden and started a study in physics at the University of Utrecht. After finishing the propedeuse he continued his study at the Rijksuniversiteit Leiden, and obtained his MSc. in the quantum optics group in 1997. In the same year he started a Ph.D at the Pattern Recognition Group of the Technical University Delft, in an IOP Beeldbewerking project called ‘Advanced texture segmentation in 2D and 3D subsurface images’. This work was carried out under supervision of Dr. Piet Verbeek and Prof.Dr.Ir. Lucas J. van Vliet. In 2002 he joined Shell International Exploration and Production in Rijswijk as a Research Geophysicist.
