Chapter 14
Feature Extraction

Under the term Feature Extraction we include several techniques aiming to detect or extract information of low level of abstraction from images. These features can be objects : points, lines, etc. They can also be measures : moments, textures, etc.

14.1 Textures

14.1.1 Haralick Descriptors

This example illustrates the use of the otb::ScalarImageToTexturesFilter , which compute the standard Haralick’s textural features [56] presented in table 14.1.1, where μ_t and σ_t are the mean and standard deviation of the row (or column, due to symmetry) sums, μ = (weighted pixel average) = ∑_i,ji⋅g(i,j) = ∑_i,jj⋅g(i,j) due to matrix summetry, and σ = (weighted pixel variance) = ∑_i,j(i-μ)² ⋅g(i,j) = ∑_i,j(j-μ)² ⋅g(i,j) due to matrix symmetry.

---

Energy	f1 = ∑i,jg(i,j)2
Entropy	f2 = -∑i,jg(i,j)log2g(i,j), or 0 if g(i,j) = 0
Correlation	f3 = ∑i,j
Difference Moment	f4 = ∑i,jg(i,j)
Inertia (a.k.a. Contrast)	f5 = ∑i,j(i-j)2g(i,j)
Cluster Shade	f6 = ∑i,j((i-μ)+(j-μ))3g(i,j)
Cluster Prominence	f7 = ∑i,j((i-μ)+(j-μ))4g(i,j)
Haralick’s Correlation	f8 =

Table 14.1: Haralick features [56] available in otb::ScalarImageToTexturesFilter

---

More features are available in otb::ScalarImageToAdvancedTexturesFilter . The following classes provide similar functionality:

• otb::ScalarImageToAdvancedTexturesFilter
• otb::ScalarImageToPanTexTextureFilter
• otb::GreyLevelCooccurrenceIndexedList

The source code for this example can be found in the file
Examples/FeatureExtraction/TextureExample.cxx.

The first step required to use the filter is to include the header file.

After defining the types for the pixels and the images used in the example, we define the types for the textures filter. It is templated by the input and output image types.

We can now instantiate the filters.

The texture filters takes at least 2 parameters: the radius of the neighborhood on which the texture will be computed and the offset used. Texture features are bivariate statistics, that is, they are computed using pair of pixels. Each texture feature is defined for an offset defining the pixel pair.

The radius parameter can be passed to the filter as a scalar parameter if the neighborhood is square, or as SizeType in any case.

The offset is always an array of N values, where N is the number of dimensions of the image.

The textures filter will automatically derive the optimal bin size for co-occurences histogram, but they need to know the input image minimum and maximum. These values can be set like this :

To tune co-occurence histogram resolution, you can use the SetNumberOfBinsPerAxis() method.

We can now plug the pipeline.

Figure 14.1 shows the result of applying the contrast texture computation.

14.1.2 PanTex

The source code for this example can be found in the file
Examples/FeatureExtraction/PanTexExample.cxx.

This example illustrates the use of the otb::ScalarImageToPanTexTextureFilter . This texture parameter was first introduced in [107] and is very useful for urban area detection. The following classes provide similar functionality:

• otb::ScalarImageToTexturesFilter
• otb::ScalarImageToAdvancedTexturesFilter

The first step required to use this filter is to include its header file.

After defining the types for the pixels and the images used in the example, we define the type for the PanTex filter. It is templated by the input and output image types.

We can now instatiate the filter.

Then, we set the parameters of the filter.The radius of the neighborhood to compute the texture. The number of bins per axis for histogram generation (it is the size of the co-occurrence matrix). Moreover, we have to specify the Min/Max in the input image. In the example, image Min/Max is set by the user to 0 and 255. Alternatively you can use the class itk::MinimumMaximumImageCalculator to calculate these values.

We can now plug the pipeline and trigger the execution by calling the Update method of the writer.

Figure 14.2 shows the result of applying the PanTex computation.

14.1.3 Structural Feature Set

The source code for this example can be found in the file
Examples/FeatureExtraction/SFSExample.cxx.

This example illustrates the use of the otb::SFSTexturesImageFilter . This filter computes the Structural Feature Set as described in [61]. These features are textural parameters which give information about the structure of lines passing through each pixel of the image.

The first step required to use this filter is to include its header file.

As with every OTB program, we start by defining the types for the images, the readers and the writers.

The we can instantiate the type for the SFS filter, which is templated over the input and output pixel types.

After that, we can instantiate the filter. We will also instantiate the reader and one writer for each output image, since the SFS filter generates 6 different features.

The SFS filter has several parameters which have to be selected. They are:

1. a spectral threshold to decide if 2 neighboring pixels are connected;
2. a spatial threshold defining the maximum length for an extracted line;
3. the number of directions which will be analyzed (the first one is to the right and they are equally distributed between 0 and 2π);
4. the α parameter fort the ω-mean feature;
5. the RatioMax parameter fort the ω-mean feature.

In order to disable the computation of a feature, the SetFeatureStatus parameter can be used. The true value enables the feature (default behavior) and the false value disables the computation. Therefore, the following line is useless, but is given here as an example.

Now, we plug the pipeline using all the writers.

Figure 14.3 shows the result of applying the SFS computation to an image

14.2 Interest Points

14.2.1 Harris detector

The source code for this example can be found in the file
Examples/FeatureExtraction/HarrisExample.cxx.

This example illustrates the use of the otb::HarrisImageFilter .

The first step required to use this filter is to include its header file.

The otb::HarrisImageFilter is templated over the input and output image types, so we start by defining:

The otb::HarrisImageFilter needs some parameters to operate. The derivative computation is performed by a convolution with the derivative of a Gaussian kernel of variance σ_D (derivation scale) and the smoothing of the image is performed by convolving with a Gaussian kernel of variance σ_I (integration scale). This allows the computation of the following matrix:

(14.1)

The output of the detector is

Figure 14.4 shows the result of applying the interest point detector to a small patch extracted from a Spot 5 image.

The output of the otb::HarrisImageFilter is an image where, for each pixel, we obtain the intensity of the detection. Often, the user may want to get access to the set of points for which the output of the detector is higher than a given threshold. This can be obtained by using the otb::HarrisImageToPointSetFilter . This filter is only templated over the input image type, the output being a itk::PointSet with pixel type equal to the image pixel type.

We declare now the filter and a pointer to the output point set.

The otb::HarrisImageToPointSetFilter takes the same parameters as the otb::HarrisImageFilter and an additional parameter : the threshold for the point selection.

We can now iterate through the obtained pointset and access the coordinates of the points. We start by accessing the container of the points which is encapsulated into the point set (see section 5.2 for more information on using itk::PointSet s) and declaring an iterator to it.

And we get the points coordinates

14.2.2 SIFT detector

14.2.3 SURF detector

The source code for this example can be found in the file
Examples/FeatureExtraction/SURFExample.cxx.

This example illustrates the use of the otb::ImageToSURFKeyPointSetFilter . The Speed-Up Robust Features (or SURF) is an algorithm in computer vision to detect and describe local features in images. The algorithm is detailed in [10]. The applications of SURF are the same as those for SIFT.

The first step required to use this filter is to include its header file.

We will start by defining the required types. We will work with a scalar image of float pixels. We also define the corresponding image reader.

The SURF descriptors will be stored in a point set containing the vector of features.

The SURF filter itself is templated over the input image and the generated point set.

We instantiate the reader.

We instantiate the filter.

We plug the filter and set the number of scales for the SURF computation. We can afterwards run the processing with the Update() method.

Once the SURF are computed, we may want to draw them on top of the input image. In order to do this, we will create the following RGB image and the corresponding writer:

We set the regions of the image by copying the information from the input image and we allocate the memory for the output image.

We can now proceed to copy the input image into the output one using region iterators. The input image is a grey level one. The output image will be made of color crosses for each SURF on top of the grey level input image. So we start by copying the grey level values on each of the 3 channels of the color image.

We are now going to plot color crosses on the output image. We will need to define offsets (top, bottom, left and right) with respect to the SURF position in order to draw the cross segments.

Now, we are going to access the point set generated by the SURF filter. The points are stored into a points container that we are going to walk through using an iterator. These are the types needed for this task:

We set the iterator to the beginning of the point set.

We get the information about image size and spacing before drawing the crosses.

And we iterate through the SURF set:

We get the pixel coordinates for each SURF by using the Value() method on the point set iterator. We use the information about size and spacing in order to convert the physical coordinates of the point into pixel coordinates.

We create a green pixel.

We draw the crosses using the offsets and checking that we are inside the image, since SURFs on the image borders would cause an out of bounds pixel access.

Finally, we write the image.

Figure 14.5 shows the result of applying the SURF point detector to a small patch extracted from a Spot 5 image.

14.3 Alignments

The source code for this example can be found in the file
Examples/FeatureExtraction/AlignmentsExample.cxx.

This example illustrates the use of the otb::ImageToPathListAlignFilter . This filter allows extracting meaningful alignments. Alignments (that is edges and lines) are detected using the Gestalt approach proposed by Desolneux et al. [39]. In this context, an event is considered meaningful if the expectation of its occurrence would be very small in a random image. One can thus consider that in a random image the direction of the gradient of a given point is uniformly distributed, and that neighbouring pixels have a very low probability of having the same gradient direction. This algorithm gives a set of straight line segments defined by the two extremity coordinates under the form of a std::list of itk::PolyLineParametricPath.

The first step required to use this filter is to include its header.

In order to visualize the detected alignments, we will use the facility class otb::DrawPathFilter which draws a itk::PolyLineParametricPath on top of a given image.

The otb::ImageToPathListAlignFilter is templated over the input image type and the output path type, so we start by defining:

Next, we build the pipeline.

We can choose the number of accepted false alarms in the detection with the method SetEps() for which the parameter is of the form -log10(max. number of false alarms).

As stated, above, the otb::DrawPathFilter , is useful for drawint the detected alignments. This class is templated over the input image and path types and also on the output image type.

We will now go through the list of detected paths and feed them to the otb::DrawPathFilter inside a loop. We will use a list iterator inside a while statement.

We define a dummy image will be iteratively fed to the otb::DrawPathFilter after the drawing of each alignment.

We iterate through the list and write the result to a file.

Figure 14.6 shows the result of applying the alignment detection to a small patch extracted from a VHR image.

14.4 Lines

14.4.1 Line Detection

The source code for this example can be found in the file
Examples/FeatureExtraction/RatioLineDetectorExample.cxx.

This example illustrates the use of the otb::RatioLineDetectorImageFilter . This filter is used for line detection in SAR images. Its principle is described in [131]: a line is detected if two parallel edges are present in the images. These edges are detected with the ratio of means detector.

The first step required to use this filter is to include its header file.

Then we must decide what pixel type to use for the image. We choose to make all computations with floating point precision and rescale the results between 0 and 255 in order to export PNG images.

The images are defined using the pixel type and the dimension.

The filter can be instantiated using the image types defined above.

An otb::ImageFileReader class is also instantiated in order to read image data from a file.

An otb::ImageFileWriter is instantiated in order to write the output image to a file.

The intensity rescaling of the results will be carried out by the itk::RescaleIntensityImageFilter which is templated by the input and output image types.

Both the filter and the reader are created by invoking their New() methods and assigning the result to SmartPointers.

The same is done for the rescaler and the writer.

The itk::RescaleIntensityImageFilter needs to know which is the minimu and maximum values of the output generated image. Those can be chosen in a generic way by using the NumericTraits functions, since they are templated over the pixel type.

The image obtained with the reader is passed as input to the otb::LineRatioDetectorImageFilter . The pipeline is built as follows.

The methods SetLengthLine() and SetWidthLine() allow setting the minimum length and the typical witdh of the lines which are to be detected.

The filter is executed by invoking the Update() method. If the filter is part of a larger image processing pipeline, calling Update() on a downstream filter will also trigger update of this filter.

We can also obtain the direction of the lines by invoking the GetOutputDirection() method.

shows the result of applying the LineRatio edge detector filter to a SAR image.

The following classes provide similar functionality:

• otb::LineCorrelationDetectorImageFilter

The source code for this example can be found in the file
Examples/FeatureExtraction/CorrelationLineDetectorExample.cxx.

This example illustrates the use of the otb::CorrelationLineDetectorImageFilter . This filter is used for line detection in SAR images. Its principle is described in [131]: a line is detected if two parallel edges are present in the images. These edges are detected with the correlation of means detector.

The first step required to use this filter is to include its header file.

Then we must decide what pixel type to use for the image. We choose to make all computations with floating point precision and rescale the results between 0 and 255 in order to export PNG images.

The images are defined using the pixel type and the dimension.

The filter can be instantiated using the image types defined above.

An otb::ImageFileReader class is also instantiated in order to read image data from a file.

An otb::ImageFileWriter is instantiated in order to write the output image to a file.

The intensity rescaling of the results will be carried out by the itk::RescaleIntensityImageFilter which is templated by the input and output image types.

Both the filter and the reader are created by invoking their New() methods and assigning the result to SmartPointers.

The same is done for the rescaler and the writer.

The itk::RescaleIntensityImageFilter needs to know which is the minimu and maximum values of the output generated image. Those can be chosen in a generic way by using the NumericTraits functions, since they are templated over the pixel type.

The image obtained with the reader is passed as input to the otb::LineCorrelationDetectorImageFilter . The pipeline is built as follows.

The methods SetLengthLine() and SetWidthLine() allow setting the minimum length and the typical witdh of the lines which are to be detected.

The filter is executed by invoking the Update() method. If the filter is part of a larger image processing pipeline, calling Update() on a downstream filter will also trigger update of this filter.

We can also obtain the direction of the lines by invoking the GetOutputDirections() method.

shows the result of applying the LineCorrelation edge detector filter to a SAR image.

The following classes provide similar functionality:

• otb::LineCorrelationDetectorImageFilter

The source code for this example can be found in the file
Examples/FeatureExtraction/AsymmetricFusionOfLineDetectorExample.cxx.

This example illustrates the use of the otb::AssymmetricFusionOfLineDetectorImageFilter .

The first step required to use this filter is to include its header file.

Then we must decide what pixel type to use for the image. We choose to make all computations with floating point precision and rescale the results between 0 and 255 in order to export PNG images.

The images are defined using the pixel type and the dimension.

The filter can be instantiated using the image types defined above.

An otb::ImageFileReader class is also instantiated in order to read image data from a file.

An otb::ImageFileWriter is instantiated in order to write the output image to a file.

The intensity rescaling of the results will be carried out by the itk::RescaleIntensityImageFilter which is templated by the input and output image types.

Both the filter and the reader are created by invoking their New() methods and assigning the result to SmartPointers.

The same is done for the rescaler and the writer.

The itk::RescaleIntensityImageFilter needs to know which is the minimu and maximum values of the output generated image. Those can be chosen in a generic way by using the NumericTraits functions, since they are templated over the pixel type.

The image obtained with the reader is passed as input to the otb::AsymmetricFusionOfDetectorImageFilter . The pipeline is built as follows.

The methods SetLengthLine() and SetWidthLine() allow setting the minimum length and the typical witdh of the lines which are to be detected.

The filter is executed by invoking the Update() method. If the filter is part of a larger image processing pipeline, calling Update() on a downstream filter will also trigger update of this filter.

Figure 14.9 shows the result of applying the AsymmetricFusionOf edge detector filter to a SAR image.

The source code for this example can be found in the file
Examples/FeatureExtraction/ParallelLineDetectionExample.cxx.

This example illustrates the details of the otb::ParallelLinePathListFilter .

14.4.2 Segment Extraction

14.4.2.1 Local Hough Transform

The source code for this example can be found in the file
Examples/FeatureExtraction/LocalHoughExample.cxx.

This example illustrates the use of the otb::ExtractSegmentsImageFilter .

The first step required to use this filter is to include its header file.

Then we must decide what pixel type to use for the image. We choose to make all computations with floating point precision and rescale the results between 0 and 255 in order to export PNG images.

The images are defined using the pixel type and the dimension.

The filter can be instantiated using the image types defined above.

An otb::ImageFileReader class is also instantiated in order to read image data from a file.

An otb::ImageFileWriter is instantiated in order to write the output image to a file.

Both the filter and the reader are created by invoking their New() methods and assigning the result to SmartPointers.

The same is done for the writer.

The image obtained with the reader is passed as input to the otb::ExtractSegmentsImageFilter . The pipeline is built as follows.

Figure 14.10 shows the result of applying the otb::LocalHoughImageFilter .

14.5 Density Features

An interesting approach to feature extraction consists in computing the density of previously detected features as simple edges or interest points.

14.5.1 Edge Density

The source code for this example can be found in the file
Examples/FeatureExtraction/EdgeDensityExample.cxx.

This example illustrates the use of the otb::EdgeDensityImageFilter . This filter computes a local density of edges on an image and can be useful to detect man made objects or urban areas, for instance. The filter has been implemented in a generic way, so that the way the edges are detected and the way their density is computed can be chosen by the user.

The first step required to use this filter is to include its header file.

We will also include the header files for the edge detector (a Canny filter) and the density estimation (a simple count on a binary image).

The first step required to use this filter is to include its header file.

As usual, we start by defining the types for the images, the reader and the writer.

We define now the type for the function which will be used by the edge density filter to estimate this density. Here we choose a function which counts the number of non null pixels per area. The function takes as template the type of the image to be processed.

These non null pixels will be the result of an edge detector. We use here the classical Canny edge detector, which is templated over the input and output image types.

Finally, we can define the type for the edge density filter which takes as template the input and output image types, the edge detector type, and the count function type..

We can now instantiate the different processing objects of the pipeline using the New() method.

The edge detection filter needs to be instantiated because we need to set its parameters. This is what we do here for the Canny filter.

After that, we can pass the edge detector to the filter which will be used it internally.

Finally, we set the file names for the input and the output images and we plug the pipeline. The Update() method of the writer will trigger the processing.

Figure 14.11 shows the result of applying the edge density filter to an image.

14.5.2 SIFT Density

14.6 Geometric Moments

14.6.1 Complex Moments

The complex geometric moments are defined as:

(14.2)

where x and y are the coordinates of the image f(x,y), i is the imaginary unit and p+q is the order of c_pq. The geometric moments are particularly useful in the case of scale changes.

14.6.1.1 Complex Moments for Images

The source code for this example can be found in the file
Examples/FeatureExtraction/ComplexMomentsImageFunctionExample.cxx.

This example illustrates the use of the otb::ComplexMomentsImageFunction .

The first step required to use this filter is to include its header file.

The otb::ComplexMomentImageFunction is templated over the input image type and the output complex type value, so we start by defining:

Next, we plug the input image into the complex moment function and we set its parameters.

We can chose the pixel of the image which will used as center for the moment computation

We can also choose the size of the neighborhood around the center pixel for the moment computation.

In order to get the value of the moment, we call the EvaluateAtIndex method.

14.6.1.2 Complex Moments for Paths

The source code for this example can be found in the file
Examples/FeatureExtraction/ComplexMomentPathExample.cxx.

The complex moments can be computed on images, but sometimes we are interested in computing them on shapes extracted from images by segmentation algorithms. These shapes can be represented by itk::Path s. This example illustrates the use of the otb::ComplexMomentPathFunction for the computation of complex geometric moments on ITK paths.

The first step required to use this filter is to include its header file.

The otb::ComplexMomentPathFunction is templated over the input path type and the output complex type value, so we start by defining:

Next, we set the parameters of the plug the input path into the complex moment function and we set its parameters.

Since the paths are defined in physical coordinates, we do not need to set the center for the moment computation as we did with the otb::ComplexMomentImageFunction . The same applies for the size of the neighborhood around the center pixel for the moment computation. The moment computation is triggered by calling the Evaluate method.

14.6.2 Hu Moments

Using the algebraic moment theory, H. Ming-Kuel obtained a family of 7 invariants with respect to planar transformations called Hu invariants, [60]. Those invariants can be seen as nonlinear combinations of the complex moments. Hu invariants have been very much used in object recognition during the last 30 years, since they are invariant to rotation, scaling and translation. [46] gives their expressions :

(14.3)

[42] have used these invariants for the recognition of aircraft silhouettes. Flusser and Suk have used them for image registration, [73].

14.6.2.1 Hu Moments for Images

The source code for this example can be found in the file
Examples/FeatureExtraction/HuMomentsImageFunctionExample.cxx.

This example illustrates the use of the otb::HuMomentsImageFunction .

The first step required to use this filter is to include its header file.

The otb::HuImageFunction is templated over the input image type and the output (real) type value, so we start by defining:

We can choose the region and the pixel of the image which will used as coordinate origin for the moment computation

Next, we plug the input image into the complex moment function and we set its parameters.

In order to get the value of the moment, we call the EvaluateAtIndex method.

The following classes provide similar functionality:

• otb::HuPathFunction

14.6.3 Flusser Moments

The Hu invariants have been modified and improved by several authors. Flusser used these moments in order to produce a new family of descriptors of order higher than 3, [46]. These descriptors are invariant to scale and rotation. They have the following expressions:

(14.4)

Examples

14.6.3.1 Flusser Moments for Images

The source code for this example can be found in the file
Examples/FeatureExtraction/FlusserMomentsImageFunctionExample.cxx.

This example illustrates the use of the otb::FlusserMomentsImageFunction .

The first step required to use this filter is to include its header file.

The otb::FlusserMomentsImageFunction is templated over the input image type and the output (real) type value, so we start by defining:

We can choose the region and the pixel of the image which will used as coordinate origin for the moment computation