PROJECTREPORT AutomaticSeismicImageAnalysisTool

(1)

Automatic Seismic Image Analysis Tool

Niklas Borwell Bo Lee

Klas Eurenius

Project in Computational Science: Report Januari 2012

PROJECT REPORT

(2)

(3)

Abstract

We have created a tool in C# that looks for geological structures in 2D seismic images, with the main motivation in its large potential for time and cost savings in the process of finding oil. Today, still, image screening is done manually at the same time as the amounts of aquired geological data grows faster than ever.

Using computer assisted image analysis to screen seismic images is a largely unexplored field of research, so we apply well known image processing and analysis theory. For image preprocessing this includes intensity level thresholding, frequency filtering, mathematical morphology etc. For analysis and classification we use skeletonization, distance transformation, property measuring etc.

Considering this works as an early stage prototype, the results are very promising. The produced software tool successes to find and classify objects correctly, except for some problems which occur when different image objects are tightly connected. Processing time for processing a

”standard sized” single image ranges from 0.2 s to 10 s depending on how much information it contains, noting that this is for an AMD Athlon X2

@ 2 GHz and that the code is not optimized.

(4)

Acknowledgements

We wish to thank supervisor Sergio Courtade from Schlumberger Ltd. for providing guidance and expertise in geology and, specifically, in seismic image interpretation. We also thank course coordinator Maya Neytcheva for making this project possible at all.

(5)

1 Introduction

As it becomes increasingly difficult to find profitable deposits of oil to explore, the accuracy of the information regarding promising geological sites becomes more and more important. Huge amounts of seismic data are collected every week but the manpower and time needed to analyze all the data manually imposes a huge cost. This means that the need for a computerized automated analysis tool is in great demand.

1.1 Project description

The aim of this project is to create a prototype program that can automatically find and classify interesting objects in seismic images. The prototype program should work as a first step towards a fully functional automated analysis tool.

Our prototype program is developed to work on a special type of images, namely colour coded amplitude spectrum images. It can easily be modified by changing a few parameters to handle different images. Automatic parameter adjustment is very complex, however, and lies outside our project scope.

We were given nine time slice images as our development dataset. They are all from the same place, only separated by some milliseconds in depth(see Section 2.1).

In the amplitude spectrum images the reflections from softer sediments are clearly differentiable from the reflections of solid rock, making interesting areas easily detectable for human eyes. However, present in all these images are virtual artefacts (see Section 2.3) that make the task of classification and dif- ferentiation much harder.

(a) Original image (b) Interestin geological features detected and marked

Figure 1: Example of seismic image before and after processing

1.2 Confusing terminology

In this report the words channel and object are used extensively. Both words carry two meanings.

(7)

Object refers to an image object, meaning a feature in the image that has been segmented (see Section 4.1). Geological objects are the actual physical objects represented in the unprocessed images such as channels, deposits and similar features that are interested to locate.

Channel will mainly refer to long elongated geological objects, like old riverbeds or valleys. Sometimes channel will refer to an information channel from the measuring equipment that was used for a particular picture. For example, the data we have used is the acoustic amplitude reflection channel in the seismological measuring equipment, or it might refer to the red color channel in a RGB image. In these cases we clarify what kind of channel we refer to, and if we don’t have a better name for it then it will be called “information channel”.

2 Background

2.1 Seismic data format

The subsurface is mapped to digital 3D volumes by sound waves - if onshore by special vibrator trucks attached with an array of receiving geophones, if offshore by boats with an air gun generating vibrations. The transmitted waves are reflected back when hitting density boundaries of, e.g., changes in rock types or fluid contents. Using this method the signal’s time of travel is closely related to subsurface depth (Satterfield/Whiteley).

Upon interpreting, the recorded volume is often sliced vertically, horizontally or along a ”horizon”. The latter represents a line or surface corresponding to a density border. Though, for this project we will solely work with horizontal slices which are more commonly known as ”time slices”. The data is colored red and blue where high amplitude reflections occur (Figure 2).

2.2 Computer aided search

Due to modern techniques and depleting oil supplies the oil industry encounters a very rapid growth in the amount of seismic data to analyze. But the increse in the number of employed doesn’t by far reflect this increase in collected data.

As for now, the images are screened manually. In order to streamline the process it does seem natural to consider computer based image analysis.

2.3 Image interpretation

As mentioned in Section 2.1, red and blue parts of an image are strong reflections. In general, the more reflections an image shows, the more interesting it is. The reflections build up shapes that represent different geological objects, and we learned that our image set has three different kinds of objects (Figure 3):

• channels,

• deposits,

• structural reflections.

(8)

(a) 3D volume illustration with stacked time slices. The time resolution of the original data set is 2 ms.

(b) Time slice where a channel can be seen.

We work with this kind of images.

(c) Vertical slice of the 3D volume. Horizons are the red and blue lines running through the image horizontally.

Figure 2: 3D volume and 2D slices.

Of these, the last one requires additional explanation. Structural reflections somtimes occur when slicing the dataset horizontally where there is a density region border that is at an angle compared to the horizontal plane of the 3D volume, for example a bedrock slope. This means that structural reflections are only artificial and not real geological objects and are only present in the 2D slices. They are very difficult to tell apart from other objects by looking at their shape. The outlining done in Figure 3(b) is made by an expert. An ideal seismic image analysis tool can manage to ignore structural reflections.

(9)

(a) Image with three different object types from our dataset.

(b) Objects outlined. Red: deposit, green:

structural reflection, blue: channcel.

Figure 3: Object types.

3 Methodology

This project is not based on or built on any previous work in this area, so we build our tool using well-known image analysis algorithms programmed in high level languages. In general the following is our method of development:

1. Use MATLAB and its Image Processing Toolbox to try different image analysis methods. The Image Processing Toolbox offers a comprehensive set of image analysis functions, and allows for a very quick startup phase in the development process.

2. When we feel confident that we have found a working method, the MAT- LAB code is migrated to C# code. For image processing and analysis we use EmguCV which is a .NET wrapper to the Intel OpenCV image processing library (for C/C++). While searching for suitable image processing libraries we found that OpenCV is, arguably, a de facto standard in image processing and computer vision implementations. EmguCV/OpenCV does actually not do everything MATLAB does, which forces us to implement many algorithms ourselves.

The choice of using C# (instead of e.g. C) is solely because Schlumberger’s software platform Petrel nicely allows for .NET plugins or apps. While our tool is not a Petrel application, it had been considered in an early stage.

Furthermore, making it C# will facilitate any attempts to integrate it within Petrel in the future.

3. Final evaluation of results and computation times. Even though speed is not a part of our goals, the problem formulation has an inherent speed concern - the tool needs to outperform human screening.

4 Theory

This section aims to explain certain image analysis concepts and methods used in the implementation.

(10)

4.1 Objects

We often describe a binary image (i.e. an image with only two colors: black and white) in terms of objects and background. An object in this context is defined as a connected white area, while all black pixels are considered part of the background. An example of this concept is shown in Figure 4.

Figure 4: A binary image with 3 objects.

4.2 Threshold

Threshold is the most common way to select parts of an image. In a normal 8-bit grayscale image, each pixel has a brightness value that ranges between 0 and 255. When we use a threshold, we simply select those pixels that have a brightness value larger than a certain threshold. Then we create a new binary image where the selected pixels are set to white and the rest set to black. An example of this technique with different threshold values is shown in Figure 5.

By using color space transformations this can be extended to a more general concept. By converting the image to a HLS (Hue-Lightness-Saturation) color space, the saturation and hue will also be represented by a grayscale image.

So by applying threshold to for instance the saturation channel, we can select saturated areas.

4.3 Closing operation

A closing operation is defined by two consecutive operations; a dilation followed by an erosion. The purpose of this operation is to connect neighbouring white areas in binary images while keeping as much of the original shape as possible.

The dilation step is used to extend the borders of white areas to connect neighbouring areas. The erosion step peels off the outer border of the white areas in the same way they were dilated to restore the original shape, while keeping the new connections between areas. The process is demonstrated in Figure 6.

(11)

(a) Original image (b) Result of a threshold at 167

(c) Result of a threshold at 200

Figure 5: Example of threshold technique

4.4 Bounding box

The bounding box of an object is the smallest rectangular box that can be drawn (with a vertical-horizontal alignment) that encapsulates every pixel of the object. An example of this is shown in Figure 7

4.5 Skeleton representation

A skeleton representation is a way of representing objects with thin lines, where the lines should represent the inner structure of the object. The method to achieve this used in the implementation works by iteravely finding and removing pixels which are found to be non-skeleton points, i.e. a pixel which is not important for the inner structure. An example of a skeleton representation is shown in Figure 8.

4.6 Constrained distance transform

The distance transform is an operation which outputs an image where each object pixel value is equal to the distance between the pixel and the closest background pixel in the original image. When using constrained distance transform, the distance is only allowed to be calculated over certain areas.

A special case of this is used in the implementation, where we only allow distance to be calculated on a line from an endpoint of the line. The result of this is a numbered line where each value represents the number of pixels between the pixel and the starting point. An example of this is shown in Figure 9.

(12)

(a) Input image with three disconnected white areas

(b) Image after dilation. The white areas have been connected

(c) Image after erosion. Notice that the white areas are connected while the major shape of the original areas are preserved

Figure 6: Steps for the closing operation

4.7 Elongatedness measure

The elongatedness measure is a way of differentiating between elongated shapes and more compact, circular shapes. It is defined by Sonka, Hlaval och Boyle (2008, 355) as

elongatedness = area (2d)²,

(13)

(a) Original image (b) Image with bounding box in red

Figure 7: Example of a bounding box

(a) Original image (b) Skeleton representation

Figure 8: Example of a skeleton representation

where d is the number of erosions that can be applied to an object before it disappears (see Section 4.3). It is a measure of the maximum thickness of a shape and the elongatedness measure therefore gives a measure of the ratio of the area to the squared maximum thickness of a shape. It is very useful since it works for any general shape.

5 Implementation

The input image used to demonstrate the steps of the program is shown in Figure 10. The main algorithm consists of the following steps.

1. Mark interesting areas.

2. Merge areas into objects.

3. Remove uninteresting objects.

4. Classifiy objects.

(14)

(a) Input image with a branching line and blue pixel as a starting point

(b) Result of constrained distance transform

Figure 9: Example of constrained distance transform

5. Vector representation.

The algorithm uses four user defined search parameters, which makes it possible to adapt the algorithm for different kinds of images:

• A search criterion which defines what to look for

• A sensitivity value for the search criterion

• A merging parameter which defines how dense the interesting areas are.

• A size limit that defines how small objects we are interested in.

The main steps of the algorithm are described below.

5.1 Mark interesting areas

In the first step interesting areas are marked based on a search parameter. This parameter consists of two parts; one or more search criterion combined with a sensitivity value for each. The search criterion in the given example is areas of blue, red and yellow color. This is achieved by converting the image to an appropriate color space and selecting areas in different information channels.

This method gives the algorithm a lot of possibilities due to the wide range of color space transformations available. Examples of possible search options could be areas with a certain color, saturated or unsaturated areas, dark or bright areas or any combination of the above. The sensitivity parameter for each channel is implemented through a simple binary threshold value (see Section 4.2).

The output of this step is a black and white image where white areas are considered interesting (see Figure 11).

5.2 Merge areas into objects

In this step interesting areas are combined into objects, where object in this case means a connected white area (see Section 4.1). The goal is for these objects to represent real life objects so that no separate objects are connected to each other and no single object is divided into two separate objects. This is however in some cases impossible, since different objects often overlap and share the

(15)

Figure 10: Original image. We want to find and mark the channel and the deposit.

same properties (see Section 4.3). The result is therefore often a compromise between the two. Possible solutions to this problem are discussed in Section 7.

The technique used to merge areas is a closing operation (see Section 4.3).

The size of the closing operation is controlled by a search parameter which defines the minimum distance between two areas to be considered part of the same object. The output of this step is a binary image consisting of objects and background (see Section 4.1), see Figure 12.

5.3 Classifying objects

In this step we perform some measurements on each object which are used to remove uninsteresting objects and classify the remaining. Currently we have two measurements; area and elongatedness (see Section 4.7).

The area measurement is used for removing objects considered too small and is controlled by a search parameter which defines a minimum area for interesting objects. One can also consider using more measurements to define interesting objects such as length, positon or base it on the classification.

The classification is based on the elongatedness measurement which differ- entiates channels from more round objects such as deposits. The classification uses a threshold value; objects with an elongatedness over 10 is considered a channel and below 10 a non-channel. The value 10 was chosen by comparing the value for several different shapes and choosing the value that best agreed

(16)

Figure 11: Output of step 1. Interesting areas are marked in white.

with our own observations. A lot of improvement in this area is possible and discussed further in Section 7. The output of this step is a list of interesting objects with measurements and classification result (see Figure 13).

5.4 Vector representation

In this step each object is processed individually based on its classification. Cur- rently only objects classified as channels will go through individual processing.

For these channels a vector representation is created, which has many potential benefits. Currently we use this to calculate the length of a channel, but since we have a mathematical representation of the channel’s path, many more advanced measurements are possible. We could for instance calculate whether it is straight or curved and how it bends.

Another important possibility is related to the segmentation problem of di- viding objects described in Section 5.5. With the vector representation we can see the direction of a channels path and by using that we could automatically connect accidentally disconnected channel segments.

The vector representation should have the following properties

• Runs along the entire length of the channel.

• Is centered in the channel in each point.

• Has a variable detail level, controlled by a parameter.

(17)

Figure 12: Output of step 2. Interesting areas are merged into objects.

It is accopmlished through the following steps 5.4.1 Skeletonization

The first step is to create a skeleton representation of the object (see Section 5.4.1). The resulting image has a centerline running through the object along with a number of branches (see Figure 14).

5.4.2 Path selection

The goal in this step is to select only the centerline of the skeleton, which is accomplished by making two asumptions

1. One of the endpoints of the centerline touches the bounding box of the object (see Section 4.4).

2. The centerline is the the longest possible line in the skeleton that agrees with 1.

Although an implementation without the need for touching the bounding box can be accomplished, it would make the algorithm considerably slower and was therefore not implemented. These asumptions are clearly not accurate for all types of possible elongated shapes, but are accurate enough when considering real-life channels. By using these asumptions we can find the centerline from the skeleton in the following way

(18)

Figure 13: Output of step 3. Uninteresting objects are removed and remaining objects are classified.

(19)

Figure 14: Skeleton representation of a channel.

1. Find all points connected to the bounding box of the object

2. Calculate a constrained distance transform image for each point found in 1 (see Section 4.6)

3. Select the point with the highest value found in the distance transform images.

4. Trace the path from the highest value back to zero.

The traced path will be the assumed centerline of the object (see Figure 15).

5.4.3 Path sampling

The final step is to sample the centerline of the object to create a simplified vector representation. This is done by selecting points from the centerline and measuring the error between the curves. The error is defined as the relative difference between the length of the original line and the length of the vector representation. More points are then added from the centerline until the error is below a threshold.

To improve the result, two accuracy limits were used; one for normal and one for short curve segments. This gives a total of three parameters to control the detail level of the vector representation: accuracy for normal segements, accuracy for short segments and length threshold for what is considered a short

(20)

Figure 15: Asumed centerline.

segment. Although it would be possible to control all these manually, we decided to simplify the usage by introducing a single parameter which controls the other three which ranges between 0.8 and 1.

An example of a final vector representation is shown in Figure 16.

5.5 Separation problem

The most difficult problem of the algorithm appears in Section 5.2: Merging of areas into objects. The different objects we wish to distinguish between usually share the same appearance and properties which makes a direct classification impossible. This is further complicated by the fact that they often overlap, see Figure17. Different objects are sometimes more connected to each other than objects are in themselves, which makes a perfect separation based on connect- edness also impossible.

The issue comes from the fact that we are looking at 2D slices from a complete 3D volume which sometimes even makes it hard for a person to distinguish between different objects. Possible solutions for this problem are discussed in Section 7.

(21)

Figure 16: Channel object with corresponding vector representation.

6 Results

In this section we consider the output of two different types of images. The first we call ”main images” because the program development is almost exclusively based on this type. The other image type is discussed in ”Other images” (Section 6.2) and we use it to demonstrate some flexibility of the tool. Last, a discussion about speed is included.

6.1 Main images

To evaluate results of the program we start with one of the nine images from Schlumberger. Because of the separation problem (Section 5.5) we use two versions of this image. The first version is the unmodified image (Figure 18(b)).

In the other version we have by hand disconnected connected objects (Figure 18(a)).

We see that in the modified image, objects are classfied correctly. There are a channel object and a deposit object.

The channel and the deposit are connected to each other in the original image, where the connection region is equally wide as the channel. This is a segmentation problem that requires a much more complex solution than we have implemented. We discuss possible improvements in Section 7.

Even though the deposit and the channel are segmented as one object, the vectorization step finds the right channel path. This can be considered as a

(22)

Figure 17: Illustration of separation problem. The channel (orange arrow) is strongly connected to the deposit (blue arrow) which is further connected to the structural reflections (green arrows).

coincidental case because the correct channel path just happens to be slightly longer than the second longest endpoint-to-endpoint path. Unfortunately this can not be guaranteed to hold in the general case.

6.2 Other images

Here we try another image (Figure 19) with a visible channel. The common feature between this image and the main image type is that objects of interest have distinctly different color intensity compared to background objects. On the other hand notable differences are color space, contrast and size. By changing the parameters of the algorithm, this image can also be nicely segmented and classified for objects.

6.3 Speed

In the long run, beyond our prototype program, speed is of utmost importance.

Due to time constraints we haven’t optimized our code with respect to speed, but we will here present some figures to show that doing automatic segmentation and classification can become a huge benefit in future.

The computations are done on an Athlon X2 2.0 Ghz (2005) CPU, and the code is not in any way optimized for this CPU type. A single original size (557x475 pixels) image takes up to 10 s to process if there are objects present.

(23)

(a) Modified image where the different object types have been separated.

(b) Unmodified image where all objects are connected.

(c) Object 1 Class: channel, Elongated- ness: 18.4, Size: 43 · 10³, Length: 813. Ob- ject 2 Class: deposit, Elongatedness: 5.2, Size: 14 · 10³, Length: 813.

(d) Channel and deposit are segmented as one image object. Note that the vector representation is the same as in (a) because of longest path beteween endpoints.

Figure 18: Result images showing segmentation, vector representation, classification and measurements.

On the other hand an ”empty” image takes only a fraction of a second to process (see Figure 20).

We note that the images we are working with are high resolution compared to the amount of information (useful for us) they contain. As the computation time is linear to the image’s pixel count (see Table 1), we can make use of downscaling in order to gain speed. As an example, our main image (Figure 18) segments the same even when downscaled to 25% of its original size. This would have to be a rather dynamic algorithm because the width of channels vary greatly.

Finally we state that these speed results show an obvious advantage of using machines instead of humans for image screening. Also there are many options to optimize the performance of our implementation even further, e.g. optimizing for CPU, optimizing the image processing algorithms, and even changing to a more modern computer.

(24)

(a) Input image. Alternative image type from Schlum- berger.

(b) Output image. Class: channel, Elongatedness: 28.4, Size: 35 · 10³, Length: 1040

Figure 19: Results for the alternative image type.

(a) 4.5 s to process this image with objects present.

(b) 0.3 s to process image with nothing to classify.

Figure 20: Computation times for images with and without objects.

7 Further development

This chapter discusses some ideas we have for further development, some of the issues we have not been able to solve yet and some ideas we have for solving them.

(25)

Resolution (pixels) Relative size (%) Time (s)

278x237 25 ≈ 1

417x356 50 ≈ 2

557x475 100 ≈ 4 − 5

1114x950 200 ≈ 15 − 18

Table 1: Computation time is linear to pixel count.

7.1 3D

Our program is developed to work on 2D data slices from a large 3D data set, handling one image at a time segmenting and classifying the objects we think are interesting and marking them before moving on to the next one.

This process is repeated for every slice in the dataset until the whole volume is analysed.

The same methods that are used to analyse the 2D slices could be generalised to work on the 3D volume. The reason for doing the analysis in 2D in the first place is inherited from manual analysis. It is difficult for humans to see through a 3D volume so the volume is sliced to allow it to be properly examined. This limitation does not apply to computers.

We think that analysing the data set in 3D is preferable to 2D for several reasons. The main reason is that having the whole 3D object gives us much more information and would render otherwise quite difficult operations unnecessary.

When handling 2D slices we have to take into account that the seismological objects we are looking for are not flat but stretched over several time-slices, making it necessary to track the objects through the dataset in order to determine how many different objects that are actually present in the volume.

Using 3D, this would become unnecessary as each object would be segmented as a whole from the start. Having access to full 3D objects could also be helpful when trying to differentiate between different geological objects that are connected to each other as image objects for some reason.

We discussed in Section 6.1 the problem of having two different kinds of objects closely connected, with no image information that clearly separates them, (Figure 17). This connection however is not present in all of the images in our dataset meaning that a human analyst can with relative certainty say that those are two separate objects. If we examine the 3D mock up (Figure 21), created by interpolating our nine time slice dataset to give it a bit more volume, we can quite easily see that the deposit, even though it is connected to the channel is probably not a part of the channel. More importantly, it will be easier for the computer to make this discovery. This kind of interpolation would of course be unnecessary if the whole 3D volume is used.

Structural reflections are a major problem when processing and analysing 2D slices (see Section 2.3), handling the data volume in 3D will ease this problem.

7.2 Use more data

One major difficulty we have faced during the development of this program is lack of data. Due to restrictive security classification of the data volumes we where only allowed to work with 9 time slices of somewhat preprocessed

(26)

Figure 21: 3D mock up of the channel and the deposit. In this image it is possible to see that the deposit (blue area on the right) is probably not a part of the channel.

images (the acoustic amplitude reflection information-channel was chosen from the dataset).

If the program could use all attributes and information channels from the dataset we think that both segmentation and classification could be made much more accurate. Comparisons between as well as combination of different information channels will be possible. This information can then be used to discover and remove noise and to help separate the kind of objects that we are not able to separate in this version of the program.

7.3 Complementing segmentation methods

In this version of the program, we employ thresholding in the acoustic amplitude information channel as the main segmentation method, but thresholding is not the only segmentation method available that we have considered. It is by far the most effective method for the images in our dataset though, which is why it was chosen. There might however be characteristics in other information channels and images where it is desirable to employ other segmentation method to complement thresholding. There are two additional methods for segmentation that are often used, edge detection and texture analysis (Gonzales, Woods 2008, 447-450, 675-679). Both of these methods should be considered again when implementing further developments.

7.4 Employ an advanced classification tool

Our prototype program uses a very simple classifier. First all objects that are too small to be of interest are removed. The remaining objects are then classified based on how elongated they are. If the object is long enough it is classified as a channel. This simple classifier is the result of our very limited dataset.

Implementing an advanced classifier based solely on nine images was never considered. We have also held back on the amount of geometric measurements we do in the segmented image as we don’t really know what measurements are actually useful and it will be very easy to implement them later.

(27)

Principal component analysis (PCA) can be used in conjunction with as suitable classifier to create a powerful classification tool. PCA transforms the parameter space into an eigenvector space based om maximum variance along each eigenvector. This transformation reduces the number of dimensions to consider when classifying a dataset, making it easier and faster to do so (Jolliffe 2002, 1-6).

We mentioned “suitable” classifier above instead of naming one. This is because different methods are suitable for different types of data. Before a classifier is chosen we need to know very well how the measured data is distributed.

A supervised classifier is probably the best choice though as it can be tailored to the data. This means that the classifier needs to be taught how to behave, and for that end a large quantity of data is needed.

7.5 Parallelization and optimization

At the moment it takes about 0.2 - 10 seconds to process an image on our test machine (see Section 6.3). This is a long time if you intend to process thousands of pictures, so like in all applications parallelization is necessary. If the program is to continue processing 2D time slices, parallelization could be as simple as to divide the images between the threads as each image can be processed individually.

The code we have written is not optimized for speed meaning that there is much room for improvement that can reduce the processing time considerably.

8 Conclusion

We have met the goals we set at the beginning of this project. We have created a prototype program that can segment, detect and, to some degree, classify interesting objects in seismic images. The program is tuned to work on color coded acoustic amplitude reflection images but can be adapted to different images by changing a few parameters.

Two major concerns are still present that we have not been able to solve in the scope of this project. These are separation of geological objects that should not be connected and an advanced classification tool. The program is by no means a complete analysis tool and there is a lot of room for further development. We have given a few suggestions to this end in this report.

9 References

Dorothy Satterfield, Martin Whiteley. Seismic Interpretation and Subsurface Mapping. http://www.derby.ac.uk/files/seismic interpretation.ppt (January 12, 2012)

Gonzales, Rafael C, Woods, Richard E. 2008. Digital image processing. 3rd Efition. New Jersey: Pearson Education Inc.

Jolliffe, I.T. 2002 Principal component analysis. 2nd Edition. New York:

Springer-verlag New York Inc.

(28)

Sonka, Milan, Hlavac, Vaclav och Boyle, Roger. 2008. Image processing, Anal- ysis, and Machine Vision. 3rd Edition. Stamford: CENGAGE Learning