Chapter 8
Basic Filtering

This chapter introduces the most commonly used filters found in OTB. Most of these filters are intended to process images. They will accept one or more images as input and will produce one or more images as output. OTB is based ITK’s data pipeline architecture in which the output of one filter is passed as input to another filter. (See Section 3.5 on page 57 for more information.)

8.1 Thresholding

The thresholding operation is used to change or identify pixel values based on specifying one or more values (called the threshold value). The following sections describe how to perform thresholding operations using OTB.

8.1.1 Binary Thresholding

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

PIC

This example illustrates the use of the binary threshold image filter. This filter is used to transform an image into a binary image by changing the pixel values according to the rule illustrated in Figure ??. The user defines two thresholds—Upper and Lower—and two intensity values—Inside and Outside. For each pixel in the input image, the value of the pixel is compared with the lower and upper thresholds. If the pixel value is inside the range defined by [Lower,Upper] the output pixel is assigned the InsideValue. Otherwise the output pixels are assigned to the OutsideValue. Thresholding is commonly applied as the last operation of a segmentation pipeline.

The first step required to use the itk::BinaryThresholdImageFilter is to include its header file.

#include "itkBinaryThresholdImageFilter.h"

The next step is to decide which pixel types to use for the input and output images.

  typedef  unsigned char InputPixelType; 
  typedef  unsigned char OutputPixelType;

The input and output image types are now defined using their respective pixel types and dimensions.

  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The filter type can be instantiated using the input and output image types defined above.

  typedef itk::BinaryThresholdImageFilter< 
      InputImageType, OutputImageType>  FilterType;

An otb::ImageFileReader class is also instantiated in order to read image data from a file. (See Section 6 on page 151 for more information about reading and writing data.)

  typedef otb::ImageFileReader<InputImageType> ReaderType;

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

  typedef otb::ImageFileWriter<InputImageType> WriterType;

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

  ReaderType::Pointer reader = ReaderType::New(); 
  FilterType::Pointer filter = FilterType::New();

The image obtained with the reader is passed as input to the BinaryThresholdImageFilter.

  filter->SetInput(reader->GetOutput());

The method SetOutsideValue() defines the intensity value to be assigned to those pixels whose intensities are outside the range defined by the lower and upper thresholds. The method SetInsideValue() defines the intensity value to be assigned to pixels with intensities falling inside the threshold range.

  filter->SetOutsideValue(outsideValue); 
  filter->SetInsideValue(insideValue);

The methods SetLowerThreshold() and SetUpperThreshold() define the range of the input image intensities that will be transformed into the InsideValue. Note that the lower and upper thresholds are values of the type of the input image pixels, while the inside and outside values are of the type of the output image pixels.

  filter->SetLowerThreshold(lowerThreshold); 
  filter->SetUpperThreshold(upperThreshold);

The execution of the filter is triggered by invoking the Update() method. If the filter’s output has been passed as input to subsequent filters, the Update() call on any posterior filters in the pipeline will indirectly trigger the update of this filter.

  filter->Update();


PIC PIC

Figure 8.1: Effect of the BinaryThresholdImageFilter on a ROI of a Spot 5 image.


Figure 8.1 illustrates the effect of this filter on a ROI of a Spot 5 image of an agricultural area. This figure shows the limitations of this filter for performing segmentation by itself. These limitations are particularly noticeable in noisy images and in images lacking spatial uniformity.

The following classes provide similar functionality:

8.1.2 General Thresholding

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


PIC PIC

Figure 8.2: ThresholdImageFilter using the threshold-below mode.



PIC PIC

Figure 8.3: ThresholdImageFilter using the threshold-above mode.



PIC PIC

Figure 8.4: ThresholdImageFilter using the threshold-outside mode.


This example illustrates the use of the itk::ThresholdImageFilter. This filter can be used to transform the intensity levels of an image in three different ways.

The following methods choose among the three operating modes of the filter.

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

#include "itkThresholdImageFilter.h"

Then we must decide what pixel type to use for the image. This filter is templated over a single image type because the algorithm only modifies pixel values outside the specified range, passing the rest through unchanged.

  typedef  unsigned char PixelType;

The image is defined using the pixel type and the dimension.

  typedef otb::Image<PixelType,  2> ImageType;

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

  typedef itk::ThresholdImageFilter<ImageType> FilterType;

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

  typedef otb::ImageFileReader<ImageType> ReaderType;

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

  typedef otb::ImageFileWriter<ImageType> WriterType;

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

  ReaderType::Pointer reader = ReaderType::New(); 
  FilterType::Pointer filter = FilterType::New();

The image obtained with the reader is passed as input to the itk::ThresholdImageFilter.

  filter->SetInput(reader->GetOutput());

The method SetOutsideValue() defines the intensity value to be assigned to those pixels whose intensities are outside the range defined by the lower and upper thresholds.

  filter->SetOutsideValue(0);

The method ThresholdBelow() defines the intensity value below which pixels of the input image will be changed to the OutsideValue.

  filter->ThresholdBelow(40);

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.

  filter->Update();

The output of this example is shown in Figure 8.2. The second operating mode of the filter is now enabled by calling the method ThresholdAbove().

  filter->ThresholdAbove(100); 
  filter->SetOutsideValue(255); 
  filter->Update();

Updating the filter with this new setting produces the output shown in Figure 8.3. The third operating mode of the filter is enabled by calling ThresholdOutside().

  filter->ThresholdOutside(40, 100); 
  filter->SetOutsideValue(0); 
  filter->Update();

The output of this third, “band-pass” thresholding mode is shown in Figure 8.4.

The examples in this section also illustrate the limitations of the thresholding filter for performing segmentation by itself. These limitations are particularly noticeable in noisy images and in images lacking spatial uniformity.

The following classes provide similar functionality:

8.1.3 Threshold to Point Set

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

Sometimes, it may be more valuable not to get an image from the threshold step but rather a list of coordinates. This can be done with the otb::ThresholdImageToPointSetFilter.

The following example illustrates the use of the otb::ThresholdImageToPointSetFilter which provide a list of points within given thresholds. Points set are described in section 5.2 on page 126.

The first step required to use this filter is to include the header

#include "otbThresholdImageToPointSetFilter.h" 
#include "itkPointSet.h"

The next step is to decide which pixel types to use for the input image and the Point Set as well as their dimension.

  typedef unsigned char PixelType; 
  const unsigned int Dimension = 2; 
 
  typedef otb::Image<PixelType, Dimension>    ImageType; 
  typedef itk::PointSet<PixelType, Dimension> PointSetType;

A reader is instanciated to read the input image

  typedef otb::ImageFileReader<ImageType> ReaderType; 
  ReaderType::Pointer reader = ReaderType::New(); 
 
  const char  filenamereader = argv[1]; 
  reader->SetFileName(filenamereader);

We get the parameters from the command line for the threshold filter. The lower and upper thresholds parameters are similar to those of the itk::BinaryThresholdImageFilter (see Section 8.1.1 on page 217 for more information).

  int lowerThreshold = atoi(argv[2]); 
  int upperThreshold = atoi(argv[3]);

Then we create the ThresholdImageToPointSetFilter and we pass the parameters.

  typedef otb::ThresholdImageToPointSetFilter 
  <ImageType, PointSetType> FilterThresholdType; 
  FilterThresholdType::Pointer filterThreshold = FilterThresholdType::New(); 
  filterThreshold->SetLowerThreshold(lowerThreshold); 
  filterThreshold->SetUpperThreshold(upperThreshold); 
  filterThreshold->SetInput(0, reader->GetOutput());

To manipulate and display the result of this filter, we manually instanciate a point set and we call the Update() method on the threshold filter to trigger the pipeline execution.

After this step, the pointSet variable contains the point set.

  PointSetType::Pointer pointSet     = PointSetType::New(); 
  pointSet = filterThreshold->GetOutput(); 
 
  filterThreshold->Update();

To display each point, we create an iterator on the list of points, which is accessible through the method GetPoints() of the PointSet.

  typedef PointSetType::PointsContainer ContainerType; 
  ContainerType pointsContainer = pointSet->GetPoints(); 
  typedef ContainerType::Iterator IteratorType; 
  IteratorType itList = pointsContainer->Begin();

A while loop enable us to through the list a display the coordinate of each point.

  while (itList != pointsContainer->End()) 
    { 
    std::cout << itList.Value() << std::endl; 
    ++itList; 
    }

8.2 Mathematical operations on images

OTB and ITK provide a lot of filters allowing to perform basic operations on image layers (thresholding, ratio, layers combinations...). It allows to create a processing chain defining at each step operations and to combine them in the data pipeline. But the library offers also the possibility to perform more generic complex mathematical operation on images in a single filter: the otb::BandMathImageFilter.

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

This filter is based on the mathematical parser library muParser. The built in functions and operators list is available at: http://muparser.sourceforge.net/mup_features.html.

In order to use this filter, at least one input image is to be set. An associated variable name can be specified or not by using the corresponding SetNthInput method. For the nth input image, if no associated variable name has been specified, a default variable name is given by concatenating the letter ”b” (for band) and the corresponding input index.

The next step is to set the expression according to the variable names. For example, in the default case with three input images the following expression is valid : ”(b1+b2)*b3”.

We start by including the needed header file. The aim of this example is to compute the Normalized Difference Vegetation Index (NDVI) from a multispecral image and perform a threshold on this indice to extract area containing a dense vegetation canopy.

#include "otbBandMathImageFilter.h"

We start by the classical typedefs needed for reading and writing the images. The otb::BandMathImageFilter class works with otb::Image as input so we need to define additional filters to extract each layer of the multispectral image

  typedef double                                                          PixelType; 
  typedef otb::VectorImage<PixelType, 2>                                  InputImageType; 
  typedef otb::Image<PixelType, 2>                                        OutputImageType; 
  typedef otb::ImageList<OutputImageType>                                 ImageListType; 
  typedef OutputImageType::PixelType                                      VPixelType; 
  typedef otb::VectorImageToImageListFilter<InputImageType, ImageListType> 
  VectorImageToImageListType; 
  typedef otb::ImageFileReader<InputImageType>                            ReaderType; 
  typedef otb::ImageFileWriter<OutputImageType>                           WriterType;

We can now define the type for the filter:

  typedef otb::BandMathImageFilter<OutputImageType>   FilterType;

We instantiate the filter, the reader, and the writer:

  ReaderType::Pointer reader = ReaderType::New(); 
  WriterType::Pointer writer = WriterType::New(); 
 
  FilterType::Pointer filter = FilterType::New(); 
 
  writer->SetInput(filter->GetOutput()); 
  reader->SetFileName(argv[1]); 
  writer->SetFileName(argv[2]);

We need now to extract now each band from the input otb::VectorImage, it illustrates the use of the otb::VectorImageToImageList. Each extracted layer are inputs of the otb::BandMathImageFilter:

  VectorImageToImageListType::Pointer imageList = VectorImageToImageListType::New(); 
  imageList->SetInput(reader->GetOutput()); 
 
  imageList->UpdateOutputInformation(); 
 
  const unsigned int nbBands = reader->GetOutput()->GetNumberOfComponentsPerPixel(); 
 
  for(unsigned int j = 0; j < nbBands; ++j) 
      { 
      filter->SetNthInput(j, imageList->GetOutput()->GetNthElement(j)); 
      }

Now we can define the mathematical expression to perform on the layers (b1, b2, b3, b4). The filter takes advantage of the parsing capabilities of the muParser library and allows to set the expression as on a digital calculator.

The expression below returns 255 if the ratio (NIRRED)(NIR+RED) is greater than 0.4 and 0 if not.

  filter->SetExpression("if((b4-b3)/(b4+b3) > 0.4, 255, 0)");

We can now plug the pipeline and run it.

  writer->Update();

The muParser library offers also the possibility to extended existing built-in functions. For example, you can use the OTB expression ”ndvi(b3, b4)” with the filter. The mathematical expression would be in this case if(ndvi(b3,b4) > 0.4, 255, 0). It will return the same result.

Figure 8.5 shows the result of the threshold over the NDVI indice to a Quickbird image.


PIC PIC

Figure 8.5: From left to right: Original image, thresholded NDVI indice.


8.3 Gradients

Computation of gradients is a fairly common operation in image processing. The term “gradient” may refer in some contexts to the gradient vectors and in others to the magnitude of the gradient vectors. ITK filters attempt to reduce this ambiguity by including the magnitude term when appropriate. ITK provides filters for computing both the image of gradient vectors and the image of magnitudes.

8.3.1 Gradient Magnitude

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

The magnitude of the image gradient is extensively used in image analysis, mainly to help in the determination of object contours and the separation of homogeneous regions. The itk::GradientMagnitudeImageFilter computes the magnitude of the image gradient at each pixel location using a simple finite differences approach. For example, in the case of 2D the computation is equivalent to convolving the image with masks of type

PICT

then adding the sum of their squares and computing the square root of the sum.

This filter will work on images of any dimension thanks to the internal use of itk::NeighborhoodIterator and itk::NeighborhoodOperator.

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

#include "itkGradientMagnitudeImageFilter.h"

Types should be chosen for the pixels of the input and output images.

  typedef    float InputPixelType; 
  typedef    float OutputPixelType;

The input and output image types can be defined using the pixel types.

  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The type of the gradient magnitude filter is defined by the input image and the output image types.

  typedef itk::GradientMagnitudeImageFilter< 
      InputImageType, OutputImageType>  FilterType;

A filter object is created by invoking the New() method and assigning the result to a itk::SmartPointer.

  FilterType::Pointer filter = FilterType::New();

The input image can be obtained from the output of another filter. Here, the source is an image reader.

  filter->SetInput(reader->GetOutput());

Finally, the filter is executed by invoking the Update() method.

  filter->Update();

If the output of this filter has been connected to other filters in a pipeline, updating any of the downstream filters will also trigger an update of this filter. For example, the gradient magnitude filter may be connected to an image writer.

  rescaler->SetInput(filter->GetOutput()); 
  writer->SetInput(rescaler->GetOutput()); 
  writer->Update();


PIC PIC

Figure 8.6: Effect of the GradientMagnitudeImageFilter.


Figure 8.6 illustrates the effect of the gradient magnitude. The figure shows the sensitivity of this filter to noisy data.

Attention should be paid to the image type chosen to represent the output image since the dynamic range of the gradient magnitude image is usually smaller than the dynamic range of the input image. As always, there are exceptions to this rule, for example, images of man-made objects that contain high contrast objects.

This filter does not apply any smoothing to the image before computing the gradients. The results can therefore be very sensitive to noise and may not be best choice for scale space analysis.

8.3.2 Gradient Magnitude With Smoothing

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

Differentiation is an ill-defined operation over digital data. In practice it is convenient to define a scale in which the differentiation should be performed. This is usually done by preprocessing the data with a smoothing filter. It has been shown that a Gaussian kernel is the most convenient choice for performing such smoothing. By choosing a particular value for the standard deviation (σ) of the Gaussian, an associated scale is selected that ignores high frequency content, commonly considered image noise.

The itk::GradientMagnitudeRecursiveGaussianImageFilter computes the magnitude of the image gradient at each pixel location. The computational process is equivalent to first smoothing the image by convolving it with a Gaussian kernel and then applying a differential operator. The user selects the value of σ.

Internally this is done by applying an IIR 1 filter that approximates a convolution with the derivative of the Gaussian kernel. Traditional convolution will produce a more accurate result, but the IIR approach is much faster, especially using large σs [3637].

GradientMagnitudeRecursiveGaussianImageFilter will work on images of any dimension by taking advantage of the natural separability of the Gaussian kernel and its derivatives.

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

#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"

Types should be instantiated based on the pixels of the input and output images.

  typedef    float InputPixelType; 
  typedef    float OutputPixelType;

With them, the input and output image types can be instantiated.

  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The filter type is now instantiated using both the input image and the output image types.

  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< 
      InputImageType, OutputImageType>  FilterType;

A filter object is created by invoking the New() method and assigning the result to a itk::SmartPointer.

  FilterType::Pointer filter = FilterType::New();

The input image can be obtained from the output of another filter. Here, an image reader is used as source.

  filter->SetInput(reader->GetOutput());

The standard deviation of the Gaussian smoothing kernel is now set.

  filter->SetSigma(sigma);

Finally the filter is executed by invoking the Update() method.

  filter->Update();

If connected to other filters in a pipeline, this filter will automatically update when any downstream filters are updated. For example, we may connect this gradient magnitude filter to an image file writer and then update the writer.

  rescaler->SetInput(filter->GetOutput()); 
  writer->SetInput(rescaler->GetOutput()); 
  writer->Update();


PIC PIC

Figure 8.7: Effect of the GradientMagnitudeRecursiveGaussianImageFilter.


Figure 8.7 illustrates the effect of this filter using σ values of 3 (left) and 5 (right). The figure shows how the sensitivity to noise can be regulated by selecting an appropriate σ. This type of scale-tunable filter is suitable for performing scale-space analysis.

Attention should be paid to the image type chosen to represent the output image since the dynamic range of the gradient magnitude image is usually smaller than the dynamic range of the input image.

8.3.3 Derivative Without Smoothing

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

The itk::DerivativeImageFilter is used for computing the partial derivative of an image, the derivative of an image along a particular axial direction.

The header file corresponding to this filter should be included first.

#include "itkDerivativeImageFilter.h"

Next, the pixel types for the input and output images must be defined and, with them, the image types can be instantiated. Note that it is important to select a signed type for the image, since the values of the derivatives will be positive as well as negative.

  typedef   float InputPixelType; 
  typedef   float OutputPixelType; 
 
  const unsigned int Dimension = 2; 
 
  typedef otb::Image<InputPixelType,  Dimension> InputImageType; 
  typedef otb::Image<OutputPixelType, Dimension> OutputImageType;

Using the image types, it is now possible to define the filter type and create the filter object.

  typedef itk::DerivativeImageFilter< 
      InputImageType, OutputImageType>  FilterType; 
 
  FilterType::Pointer filter = FilterType::New();

The order of the derivative is selected with the SetOrder() method. The direction along which the derivative will be computed is selected with the SetDirection() method.

  filter->SetOrder(atoi(argv[4])); 
  filter->SetDirection(atoi(argv[5]));

The input to the filter can be taken from any other filter, for example a reader. The output can be passed down the pipeline to other filters, for example, a writer. An update call on any downstream filter will trigger the execution of the derivative filter.

  filter->SetInput(reader->GetOutput()); 
  writer->SetInput(filter->GetOutput()); 
  writer->Update();


PIC PIC

Figure 8.8: Effect of the Derivative filter.


Figure 8.8 illustrates the effect of the DerivativeImageFilter. The derivative is taken along the x direction. The sensitivity to noise in the image is evident from this result.

8.4 Second Order Derivatives

8.4.1 Laplacian Filters

Laplacian Filter Recursive Gaussian

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

This example illustrates how to use the itk::RecursiveGaussianImageFilter for computing the Laplacian of an image.

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

#include "itkRecursiveGaussianImageFilter.h"

Types should be selected on the desired input and output pixel types.

  typedef    float InputPixelType; 
  typedef    float OutputPixelType;

The input and output image types are instantiated using the pixel types.

  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The filter type is now instantiated using both the input image and the output image types.

  typedef itk::RecursiveGaussianImageFilter< 
      InputImageType, OutputImageType>  FilterType;

This filter applies the approximation of the convolution along a single dimension. It is therefore necessary to concatenate several of these filters to produce smoothing in all directions. In this example, we create a pair of filters since we are processing a 2D image. The filters are created by invoking the New() method and assigning the result to a itk::SmartPointer.

We need two filters for computing the X component of the Laplacian and two other filters for computing the Y component.

  FilterType::Pointer filterX1 = FilterType::New(); 
  FilterType::Pointer filterY1 = FilterType::New(); 
 
  FilterType::Pointer filterX2 = FilterType::New(); 
  FilterType::Pointer filterY2 = FilterType::New();

Since each one of the newly created filters has the potential to perform filtering along any dimension, we have to restrict each one to a particular direction. This is done with the SetDirection() method.

  filterX1->SetDirection(0);     // 0 --> X direction 
  filterY1->SetDirection(1);     // 1 --> Y direction 
 
  filterX2->SetDirection(0);     // 0 --> X direction 
  filterY2->SetDirection(1);     // 1 --> Y direction

The itk::RecursiveGaussianImageFilter can approximate the convolution with the Gaussian or with its first and second derivatives. We select one of these options by using the SetOrder() method. Note that the argument is an enum whose values can be ZeroOrder, FirstOrder and SecondOrder. For example, to compute the x partial derivative we should select FirstOrder for x and ZeroOrder for y. Here we want only to smooth in x and y, so we select ZeroOrder in both directions.

  filterX1->SetOrder(FilterType::ZeroOrder); 
  filterY1->SetOrder(FilterType::SecondOrder); 
 
  filterX2->SetOrder(FilterType::SecondOrder); 
  filterY2->SetOrder(FilterType::ZeroOrder);

There are two typical ways of normalizing Gaussians depending on their application. For scale-space analysis it is desirable to use a normalization that will preserve the maximum value of the input. This normalization is represented by the following equation.

  1
-√---
σ  2π
(8.1)

In applications that use the Gaussian as a solution of the diffusion equation it is desirable to use a normalization that preserve the integral of the signal. This last approach can be seen as a conservation of mass principle. This is represented by the following equation.

--1---
σ2√2π
(8.2)

The itk::RecursiveGaussianImageFilter has a boolean flag that allows users to select between these two normalization options. Selection is done with the method SetNormalizeAcrossScale(). Enable this flag to analyzing an image across scale-space. In the current example, this setting has no impact because we are actually renormalizing the output to the dynamic range of the reader, so we simply disable the flag.

  const bool normalizeAcrossScale = false; 
  filterX1->SetNormalizeAcrossScale(normalizeAcrossScale); 
  filterY1->SetNormalizeAcrossScale(normalizeAcrossScale); 
  filterX2->SetNormalizeAcrossScale(normalizeAcrossScale); 
  filterY2->SetNormalizeAcrossScale(normalizeAcrossScale);

The input image can be obtained from the output of another filter. Here, an image reader is used as the source. The image is passed to the x filter and then to the y filter. The reason for keeping these two filters separate is that it is usual in scale-space applications to compute not only the smoothing but also combinations of derivatives at different orders and smoothing. Some factorization is possible when separate filters are used to generate the intermediate results. Here this capability is less interesting, though, since we only want to smooth the image in all directions.

  filterX1->SetInput(reader->GetOutput()); 
  filterY1->SetInput(filterX1->GetOutput()); 
 
  filterY2->SetInput(reader->GetOutput()); 
  filterX2->SetInput(filterY2->GetOutput());

It is now time to select the σ of the Gaussian used to smooth the data. Note that σ must be passed to both filters and that sigma is considered to be in the units of the image spacing. That is, at the moment of applying the smoothing process, the filter will take into account the spacing values defined in the image.

  filterX1->SetSigma(sigma); 
  filterY1->SetSigma(sigma); 
  filterX2->SetSigma(sigma); 
  filterY2->SetSigma(sigma);

Finally the two components of the Laplacian should be added together. The itk::AddImageFilter is used for this purpose.

  typedef itk::AddImageFilter< 
      OutputImageType, 
      OutputImageType, 
      OutputImageType> AddFilterType; 
 
  AddFilterType::Pointer addFilter = AddFilterType::New(); 
 
  addFilter->SetInput1(filterY1->GetOutput()); 
  addFilter->SetInput2(filterX2->GetOutput());

The filters are triggered by invoking Update() on the Add filter at the end of the pipeline.

  try 
    { 
    addFilter->Update(); 
    } 
  catch (itk::ExceptionObject& err) 
    { 
    std::cout << "ExceptionObject caught !" << std::endl; 
    std::cout << err << std::endl; 
    return EXIT_FAILURE; 
    }

The resulting image could be saved to a file using the otb::ImageFileWriter class.

  typedef  float WritePixelType; 
 
  typedef otb::Image<WritePixelType, 2> WriteImageType; 
 
  typedef otb::ImageFileWriter<WriteImageType> WriterType; 
 
  WriterType::Pointer writer = WriterType::New(); 
 
  writer->SetInput(addFilter->GetOutput()); 
 
  writer->SetFileName(argv[2]); 
 
  writer->Update();


PIC PIC

Figure 8.9: Effect of the RecursiveGaussianImageFilter.


Figure 8.9 illustrates the effect of this filter using σ values of 3 (left) and 5 (right). The figure shows how the attenuation of noise can be regulated by selecting the appropriate standard deviation. This type of scale-tunable filter is suitable for performing scale-space analysis.

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

The previous exampled showed how to use the itk::RecursiveGaussianImageFilter for computing the equivalent of a Laplacian of an image after smoothing with a Gaussian. The elements used in this previous example have been packaged together in the itk::LaplacianRecursiveGaussianImageFilter in order to simplify its usage. This current example shows how to use this convenience filter for achieving the same results as the previous example.

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

#include "itkLaplacianRecursiveGaussianImageFilter.h"

Types should be selected on the desired input and output pixel types.

  typedef    float InputPixelType; 
  typedef    float OutputPixelType;

The input and output image types are instantiated using the pixel types.

  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The filter type is now instantiated using both the input image and the output image types.

  typedef itk::LaplacianRecursiveGaussianImageFilter< 
      InputImageType, OutputImageType>  FilterType;

This filter packages all the components illustrated in the previous example. The filter is created by invoking the New() method and assigning the result to a itk::SmartPointer.

  FilterType::Pointer laplacian = FilterType::New();

The option for normalizing across scale space can also be selected in this filter.

  laplacian->SetNormalizeAcrossScale(false);

The input image can be obtained from the output of another filter. Here, an image reader is used as the source.

  laplacian->SetInput(reader->GetOutput());

It is now time to select the σ of the Gaussian used to smooth the data. Note that σ must be passed to both filters and that sigma is considered to be in the units of the image spacing. That is, at the moment of applying the smoothing process, the filter will take into account the spacing values defined in the image.

  laplacian->SetSigma(sigma);

Finally the pipeline is executed by invoking the Update() method.

  try 
    { 
    laplacian->Update(); 
    } 
  catch (itk::ExceptionObject& err) 
    { 
    std::cout << "ExceptionObject caught !" << std::endl; 
    std::cout << err << std::endl; 
    return EXIT_FAILURE; 
    }


PIC PIC

Figure 8.10: Effect of the LaplacianRecursiveGaussianImageFilter.


Figure 8.10 illustrates the effect of this filter using σ values of 3 (left) and 5 (right). The figure shows how the attenuation of noise can be regulated by selecting the appropriate standard deviation. This type of scale-tunable filter is suitable for performing scale-space analysis.

8.5 Edge Detection

8.5.1 Canny Edge Detection

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

This example introduces the use of the itk::CannyEdgeDetectionImageFilter. This filter is widely used for edge detection since it is the optimal solution satisfying the constraints of good sensitivity, localization and noise robustness.

The first step required for using this filter is to include its header file

#include "itkCannyEdgeDetectionImageFilter.h"

As the Canny filter works with real values, we can instanciated the reader using an image with pixels as double. This does not imply anything on the real image coding format which will be cast into double.

  typedef otb::ImageFileReader<RealImageType> ReaderType;

The itk::CannyEdgeDetectionImageFilter is instantiated using the float image type.


PIC PIC

Figure 8.11: Effect of the CannyEdgeDetectorImageFilter on a ROI of a Spot 5 image.


Figure 8.11 illustrates the effect of this filter on a ROI of a Spot 5 image of an agricultural area.

8.5.2 Ratio of Means Detector

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

This example illustrates the use of the otb::TouziEdgeDetectorImageFilter. This filter belongs to the family of the fixed false alarm rate edge detectors but it is apropriate for SAR images, where the speckle noise is considered as multiplicative. By analogy with the classical gradient-based edge detectors which are suited to the additive noise case, this filter computes a ratio of local means in both sides of the edge [129]. In order to have a normalized response, the following computation is performed :

r= 1− min{μA,μB},
          μB μA
(8.3)

where μA and μB are the local means computed at both sides of the edge. In order to detect edges with any orientation, r is computed for the 4 principal directions and the maximum response is kept.

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

#include "otbTouziEdgeDetectorImageFilter.h"

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.

  typedef  float         InternalPixelType; 
  typedef  unsigned char OutputPixelType;

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

  typedef otb::Image<InternalPixelType,  2> InternalImageType; 
  typedef otb::Image<OutputPixelType,  2>   OutputImageType;

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

  typedef otb::TouziEdgeDetectorImageFilter<InternalImageType, 
      InternalImageType> FilterType;

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

  typedef otb::ImageFileReader<InternalImageType> ReaderType;

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

  typedef otb::ImageFileWriter<OutputImageType> WriterType;

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

  typedef itk::RescaleIntensityImageFilter<InternalImageType, 
      OutputImageType> RescalerType;

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

  ReaderType::Pointer reader = ReaderType::New(); 
  FilterType::Pointer filter = FilterType::New();

The same is done for the rescaler and the writer.

  RescalerType::Pointer rescaler = RescalerType::New(); 
  WriterType::Pointer   writer = WriterType::New();

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.

  rescaler->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min()); 
  rescaler->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());

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

  filter->SetInput(reader->GetOutput()); 
  rescaler->SetInput(filter->GetOutput()); 
  writer->SetInput(rescaler->GetOutput());

The method SetRadius() defines the size of the window to be used for the computation of the local means.

  FilterType::SizeType Radius; 
  Radius[0] = atoi(argv[4]); 
  Radius[1] = atoi(argv[4]); 
 
  filter->SetRadius(Radius);

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.

  filter->Update();

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

  rescaler->SetInput(filter->GetOutputDirection()); 
  writer->SetInput(rescaler->GetOutput()); 
  writer->Update();

Figure 8.12 shows the result of applying the Touzi edge detector filter to a SAR image.


PIC PIC PIC

Figure 8.12: Result of applying the otb::TouziEdgeDetectorImageFilter to a SAR image. From left to right : original image, edge intensity and edge orientation.


8.6 Neighborhood Filters

The concept of locality is frequently encountered in image processing in the form of filters that compute every output pixel using information from a small region in the neighborhood of the input pixel. The classical form of these filters are the 3×3 filters in 2D images. Convolution masks based on these neighborhoods can perform diverse tasks ranging from noise reduction, to differential operations, to mathematical morphology.

The Insight toolkit implements an elegant approach to neighborhood-based image filtering. The input image is processed using a special iterator called the itk::NeighborhoodIterator. This iterator is capable of moving over all the pixels in an image and, for each position, it can address the pixels in a local neighborhood. Operators are defined that apply an algorithmic operation in the neighborhood of the input pixel to produce a value for the output pixel. The following section describes some of the more commonly used filters that take advantage of this construction. (See Chapter 25 on page 1072 for more information about iterators.)

8.6.1 Mean Filter

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

The itk::MeanImageFilter is commonly used for noise reduction. The filter computes the value of each output pixel by finding the statistical mean of the neighborhood of the corresponding input pixel. The following figure illustrates the local effect of the MeanImageFilter. The statistical mean of the neighborhood on the left is passed as the output value associated with the pixel at the center of the neighborhood.

PICT

Note that this algorithm is sensitive to the presence of outliers in the neighborhood. This filter will work on images of any dimension thanks to the internal use of itk::SmartNeighborhoodIterator and itk::NeighborhoodOperator. The size of the neighborhood over which the mean is computed can be set by the user.

The header file corresponding to this filter should be included first.

#include "itkMeanImageFilter.h"

Then the pixel types for input and output image must be defined and, with them, the image types can be instantiated.

  typedef   unsigned char InputPixelType; 
  typedef   unsigned char OutputPixelType; 
 
  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

Using the image types it is now possible to instantiate the filter type and create the filter object.

  typedef itk::MeanImageFilter< 
      InputImageType, OutputImageType>  FilterType; 
 
  FilterType::Pointer filter = FilterType::New();

The size of the neighborhood is defined along every dimension by passing a SizeType object with the corresponding values. The value on each dimension is used as the semi-size of a rectangular box. For example, in 2D a size of 1,2  will result in a 3×5 neighborhood.

  InputImageType::SizeType indexRadius; 
 
  indexRadius[0] = 1; // radius along x 
  indexRadius[1] = 1; // radius along y 
 
  filter->SetRadius(indexRadius);

The input to the filter can be taken from any other filter, for example a reader. The output can be passed down the pipeline to other filters, for example, a writer. An update call on any downstream filter will trigger the execution of the mean filter.

  filter->SetInput(reader->GetOutput()); 
  writer->SetInput(filter->GetOutput()); 
  writer->Update();


PIC PIC

Figure 8.13: Effect of the MeanImageFilter.


Figure 8.13 illustrates the effect of this filter using neighborhood radii of 1,1  which corresponds to a 3×3 classical neighborhood. It can be seen from this picture that edges are rapidly degraded by the diffusion of intensity values among neighbors.

8.6.2 Median Filter

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

The itk::MedianImageFilter is commonly used as a robust approach for noise reduction. This filter is particularly efficient against salt-and-pepper noise. In other words, it is robust to the presence of gray-level outliers. MedianImageFilter computes the value of each output pixel as the statistical median of the neighborhood of values around the corresponding input pixel. The following figure illustrates the local effect of this filter. The statistical median of the neighborhood on the left is passed as the output value associated with the pixel at the center of the neighborhood.

PICT

This filter will work on images of any dimension thanks to the internal use of itk::NeighborhoodIterator and itk::NeighborhoodOperator. The size of the neighborhood over which the median is computed can be set by the user.

The header file corresponding to this filter should be included first.

#include "itkMedianImageFilter.h"

Then the pixel and image types of the input and output must be defined.

  typedef   unsigned char InputPixelType; 
  typedef   unsigned char OutputPixelType; 
 
  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

Using the image types, it is now possible to define the filter type and create the filter object.

  typedef itk::MedianImageFilter< 
      InputImageType, OutputImageType>  FilterType; 
 
  FilterType::Pointer filter = FilterType::New();

The size of the neighborhood is defined along every dimension by passing a SizeType object with the corresponding values. The value on each dimension is used as the semi-size of a rectangular box. For example, in 2D a size of 1,2  will result in a 3×5 neighborhood.

  InputImageType::SizeType indexRadius; 
 
  indexRadius[0] = 1; // radius along x 
  indexRadius[1] = 1; // radius along y 
 
  filter->SetRadius(indexRadius);

The input to the filter can be taken from any other filter, for example a reader. The output can be passed down the pipeline to other filters, for example, a writer. An update call on any downstream filter will trigger the execution of the median filter.

  filter->SetInput(reader->GetOutput()); 
  writer->SetInput(filter->GetOutput()); 
  writer->Update();


PIC PIC

Figure 8.14: Effect of the MedianImageFilter.


Figure 8.14 illustrates the effect of the MedianImageFilter filter a neighborhood radius of 1,1  , which corresponds to a 3×3 classical neighborhood. The filtered image demonstrates the moderate tendency of the median filter to preserve edges.

8.6.3 Mathematical Morphology

Mathematical morphology has proved to be a powerful resource for image processing and analysis [123]. ITK implements mathematical morphology filters using NeighborhoodIterators and itk::NeighborhoodOperators. The toolkit contains two types of image morphology algorithms, filters that operate on binary images and filters that operate on grayscale images.

Binary Filters

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

The following section illustrates the use of filters that perform basic mathematical morphology operations on binary images. The itk::BinaryErodeImageFilter and itk::BinaryDilateImageFilter are described here. The filter names clearly specify the type of image on which they operate. The header files required to construct a simple example of the use of the mathematical morphology filters are included below.

#include "itkBinaryErodeImageFilter.h" 
#include "itkBinaryDilateImageFilter.h" 
#include "itkBinaryBallStructuringElement.h"

The following code defines the input and output pixel types and their associated image types.

  const unsigned int Dimension = 2; 
 
  typedef unsigned char InputPixelType; 
  typedef unsigned char OutputPixelType; 
 
  typedef otb::Image<InputPixelType,  Dimension> InputImageType; 
  typedef otb::Image<OutputPixelType, Dimension> OutputImageType;

Mathematical morphology operations are implemented by applying an operator over the neighborhood of each input pixel. The combination of the rule and the neighborhood is known as structuring element. Although some rules have become de facto standards for image processing, there is a good deal of freedom as to what kind of algorithmic rule should be applied to the neighborhood. The implementation in ITK follows the typical rule of minimum for erosion and maximum for dilation.

The structuring element is implemented as a NeighborhoodOperator. In particular, the default structuring element is the itk::BinaryBallStructuringElement class. This class is instantiated using the pixel type and dimension of the input image.

  typedef itk::BinaryBallStructuringElement< 
      InputPixelType, 
      Dimension>             StructuringElementType;

The structuring element type is then used along with the input and output image types for instantiating the type of the filters.

  typedef itk::BinaryErodeImageFilter< 
      InputImageType, 
      OutputImageType, 
      StructuringElementType>  ErodeFilterType; 
 
  typedef itk::BinaryDilateImageFilter< 
      InputImageType, 
      OutputImageType, 
      StructuringElementType>  DilateFilterType;

The filters can now be created by invoking the New() method and assigning the result to itk::SmartPointers.

  ErodeFilterType::Pointer  binaryErode  = ErodeFilterType::New(); 
  DilateFilterType::Pointer binaryDilate = DilateFilterType::New();

The structuring element is not a reference counted class. Thus it is created as a C++ stack object instead of using New() and SmartPointers. The radius of the neighborhood associated with the structuring element is defined with the SetRadius() method and the CreateStructuringElement() method is invoked in order to initialize the operator. The resulting structuring element is passed to the mathematical morphology filter through the SetKernel() method, as illustrated below.

  StructuringElementType structuringElement; 
 
  structuringElement.SetRadius(1);    // 3x3 structuring element 
 
  structuringElement.CreateStructuringElement(); 
 
  binaryErode->SetKernel(structuringElement); 
  binaryDilate->SetKernel(structuringElement);

A binary image is provided as input to the filters. This image might be, for example, the output of a binary threshold image filter.

  thresholder->SetInput(reader->GetOutput()); 
 
  InputPixelType background =   0; 
  InputPixelType foreground = 255; 
 
  thresholder->SetOutsideValue(background); 
  thresholder->SetInsideValue(foreground); 
 
  thresholder->SetLowerThreshold(lowerThreshold); 
  thresholder->SetUpperThreshold(upperThreshold);
  binaryErode->SetInput(thresholder->GetOutput()); 
  binaryDilate->SetInput(thresholder->GetOutput());

The values that correspond to “objects” in the binary image are specified with the methods SetErodeValue() and SetDilateValue(). The value passed to these methods will be considered the value over which the dilation and erosion rules will apply.

  binaryErode->SetErodeValue(foreground); 
  binaryDilate->SetDilateValue(foreground);

The filter is executed by invoking its Update() method, or by updating any downstream filter, like, for example, an image writer.

  writerDilation->SetInput(binaryDilate->GetOutput()); 
  writerDilation->Update();


PIC PIC PIC

Figure 8.15: Effect of erosion and dilation in a binary image.


Figure 8.15 illustrates the effect of the erosion and dilation filters. The figure shows how these operations can be used to remove spurious details from segmented images.

Grayscale Filters

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

The following section illustrates the use of filters for performing basic mathematical morphology operations on grayscale images. The itk::GrayscaleErodeImageFilter and itk::GrayscaleDilateImageFilter are covered in this example. The filter names clearly specify the type of image on which they operate. The header files required for a simple example of the use of grayscale mathematical morphology filters are presented below.

#include "itkGrayscaleErodeImageFilter.h" 
#include "itkGrayscaleDilateImageFilter.h" 
#include "itkBinaryBallStructuringElement.h"

The following code defines the input and output pixel types and their associated image types.

  const unsigned int Dimension = 2; 
 
  typedef unsigned char InputPixelType; 
  typedef unsigned char OutputPixelType; 
 
  typedef otb::Image<InputPixelType,  Dimension> InputImageType; 
  typedef otb::Image<OutputPixelType, Dimension> OutputImageType;

Mathematical morphology operations are based on the application of an operator over a neighborhood of each input pixel. The combination of the rule and the neighborhood is known as structuring element. Although some rules have become the de facto standard in image processing there is a good deal of freedom as to what kind of algorithmic rule should be applied on the neighborhood. The implementation in ITK follows the typical rule of minimum for erosion and maximum for dilation.

The structuring element is implemented as a itk::NeighborhoodOperator. In particular, the default structuring element is the itk::BinaryBallStructuringElement class. This class is instantiated using the pixel type and dimension of the input image.

  typedef itk::BinaryBallStructuringElement< 
      InputPixelType, 
      Dimension>             StructuringElementType;

The structuring element type is then used along with the input and output image types for instantiating the type of the filters.

  typedef itk::GrayscaleErodeImageFilter< 
      InputImageType, 
      OutputImageType, 
      StructuringElementType>  ErodeFilterType; 
 
  typedef itk::GrayscaleDilateImageFilter< 
      InputImageType, 
      OutputImageType, 
      StructuringElementType>  DilateFilterType;

The filters can now be created by invoking the New() method and assigning the result to SmartPointers.

  ErodeFilterType::Pointer  grayscaleErode  = ErodeFilterType::New(); 
  DilateFilterType::Pointer grayscaleDilate = DilateFilterType::New();

The structuring element is not a reference counted class. Thus it is created as a C++ stack object instead of using New() and SmartPointers. The radius of the neighborhood associated with the structuring element is defined with the SetRadius() method and the CreateStructuringElement() method is invoked in order to initialize the operator. The resulting structuring element is passed to the mathematical morphology filter through the SetKernel() method, as illustrated below.

  StructuringElementType structuringElement; 
 
  structuringElement.SetRadius(1);    // 3x3 structuring element 
 
  structuringElement.CreateStructuringElement(); 
 
  grayscaleErode->SetKernel(structuringElement); 
  grayscaleDilate->SetKernel(structuringElement);

A grayscale image is provided as input to the filters. This image might be, for example, the output of a reader.

  grayscaleErode->SetInput(reader->GetOutput()); 
  grayscaleDilate->SetInput(reader->GetOutput());

The filter is executed by invoking its Update() method, or by updating any downstream filter, like, for example, an image writer.

  writerDilation->SetInput(grayscaleDilate->GetOutput()); 
  writerDilation->Update();


PIC PIC PIC

Figure 8.16: Effect of erosion and dilation in a grayscale image.


Figure 8.16 illustrates the effect of the erosion and dilation filters. The figure shows how these operations can be used to remove spurious details from segmented images.

8.7 Smoothing Filters

Real image data has a level of uncertainty that is manifested in the variability of measures assigned to pixels. This uncertainty is usually interpreted as noise and considered an undesirable component of the image data. This section describes several methods that can be applied to reduce noise on images.

8.7.1 Blurring

Blurring is the traditional approach for removing noise from images. It is usually implemented in the form of a convolution with a kernel. The effect of blurring on the image spectrum is to attenuate high spatial frequencies. Different kernels attenuate frequencies in different ways. One of the most commonly used kernels is the Gaussian. Two implementations of Gaussian smoothing are available in the toolkit. The first one is based on a traditional convolution while the other is based on the application of IIR filters that approximate the convolution with a Gaussian [3637].

Discrete Gaussian

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

PIC

The itk::DiscreteGaussianImageFilter computes the convolution of the input image with a Gaussian kernel. This is done in ND by taking advantage of the separability of the Gaussian kernel. A one-dimensional Gaussian function is discretized on a convolution kernel. The size of the kernel is extended until there are enough discrete points in the Gaussian to ensure that a user-provided maximum error is not exceeded. Since the size of the kernel is unknown a priori, it is necessary to impose a limit to its growth. The user can thus provide a value to be the maximum admissible size of the kernel. Discretization error is defined as the difference between the area under the discrete Gaussian curve (which has finite support) and the area under the continuous Gaussian.

Gaussian kernels in ITK are constructed according to the theory of Tony Lindeberg [87] so that smoothing and derivative operations commute before and after discretization. In other words, finite difference derivatives on an image I that has been smoothed by convolution with the Gaussian are equivalent to finite differences computed on I by convolving with a derivative of the Gaussian.

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

#include "itkDiscreteGaussianImageFilter.h"

Types should be chosen for the pixels of the input and output images. Image types can be instantiated using the pixel type and dimension.

  typedef    float InputPixelType; 
  typedef    float OutputPixelType; 
 
  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The discrete Gaussian filter type is instantiated using the input and output image types. A corresponding filter object is created.

  typedef itk::DiscreteGaussianImageFilter< 
      InputImageType, OutputImageType>  FilterType; 
 
  FilterType::Pointer filter = FilterType::New();

The input image can be obtained from the output of another filter. Here, an image reader is used as its input.

  filter->SetInput(reader->GetOutput());

The filter requires the user to provide a value for the variance associated with the Gaussian kernel. The method SetVariance() is used for this purpose. The discrete Gaussian is constructed as a convolution kernel. The maximum kernel size can be set by the user. Note that the combination of variance and kernel-size values may result in a truncated Gaussian kernel.

  filter->SetVariance(gaussianVariance); 
  filter->SetMaximumKernelWidth(maxKernelWidth);

Finally, the filter is executed by invoking the Update() method.

  filter->Update();

If the output of this filter has been connected to other filters down the pipeline, updating any of the downstream filters would have triggered the execution of this one. For example, a writer could have been used after the filter.

  rescaler->SetInput(filter->GetOutput()); 
  writer->SetInput(rescaler->GetOutput()); 
  writer->Update();


PIC PIC

Figure 8.17: Effect of the DiscreteGaussianImageFilter.


Figure 8.17 illustrates the effect of this filter.

Note that large Gaussian variances will produce large convolution kernels and correspondingly slower computation times. Unless a high degree of accuracy is required, it may be more desirable to use the approximating itk::RecursiveGaussianImageFilter with large variances.

8.7.2 Edge Preserving Smoothing

Introduction to Anisotropic Diffusion

The drawback of image denoising (smoothing) is that it tends to blur away the sharp boundaries in the image that help to distinguish between the larger-scale anatomical structures that one is trying to characterize (which also limits the size of the smoothing kernels in most applications). Even in cases where smoothing does not obliterate boundaries, it tends to distort the fine structure of the image and thereby changes subtle aspects of the anatomical shapes in question.

Perona and Malik [106] introduced an alternative to linear-filtering that they called anisotropic diffusion. Anisotropic diffusion is closely related to the earlier work of Grossberg [53], who used similar nonlinear diffusion processes to model human vision. The motivation for anisotropic diffusion (also called nonuniform or variable conductance diffusion) is that a Gaussian smoothed image is a single time slice of the solution to the heat equation, that has the original image as its initial conditions. Thus, the solution to

∂g(x,y,t)-
  ∂t   = ∇ ⋅∇g(x,y,t),
(8.4)

where g(x,y,0) = f(x,y) is the input image, is g(x,y,t) = G(√ --
  2t)f(x,y), where G(σ) is a Gaussian with standard deviation σ.

Anisotropic diffusion includes a variable conductance term that, in turn, depends on the differential structure of the image. Thus, the variable conductance can be formulated to limit the smoothing at “edges” in images, as measured by high gradient magnitude, for example.

gt = ∇ ⋅c(|∇g |)∇g,
(8.5)

where, for notational convenience, we leave off the independent parameters of g and use the subscripts with respect to those parameters to indicate partial derivatives. The function c(|∇g|) is a fuzzy cutoff that reduces the conductance at areas of large |∇g|, and can be any one of a number of functions. The literature has shown

c(|∇g|)= e− |∇g2k|22
(8.6)

to be quite effective. Notice that conductance term introduces a free parameter k, the conductance parameter, that controls the sensitivity of the process to edge contrast. Thus, anisotropic diffusion entails two free parameters: the conductance parameter, k, and the time parameter, t, that is analogous to σ, the effective width of the filter when using Gaussian kernels.

Equation 8.5 is a nonlinear partial differential equation that can be solved on a discrete grid using finite forward differences. Thus, the smoothed image is obtained only by an iterative process, not a convolution or non-stationary, linear filter. Typically, the number of iterations required for practical results are small, and large 2D images can be processed in several tens of seconds using carefully written code running on modern, general purpose, single-processor computers. The technique applies readily and effectively to 3D images, but requires more processing time.

In the early 1990’s several research groups [49137] demonstrated the effectiveness of anisotropic diffusion on medical images. In a series of papers on the subject [142139141137138140], Whitaker described a detailed analytical and empirical analysis, introduced a smoothing term in the conductance that made the process more robust, invented a numerical scheme that virtually eliminated directional artifacts in the original algorithm, and generalized anisotropic diffusion to vector-valued images, an image processing technique that can be used on vector-valued medical data (such as the color cryosection data of the Visible Human Project).

For a vector-valued input ⃗F : U↦→m the process takes the form

⃗Ft = ∇⋅c(D⃗F )F⃗,
(8.7)

where D⃗F is a dissimilarity measure of ⃗F, a generalization of the gradient magnitude to vector-valued images, that can incorporate linear and nonlinear coordinate transformations on the range of ⃗F. In this way, the smoothing of the multiple images associated with vector-valued data is coupled through the conductance term, that fuses the information in the different images. Thus vector-valued, nonlinear diffusion can combine low-level image features (e.g. edges) across all “channels” of a vector-valued image in order to preserve or enhance those features in all of image “channels”.

Vector-valued anisotropic diffusion is useful for denoising data from devices that produce multiple values such as MRI or color photography. When performing nonlinear diffusion on a color image, the color channels are diffused separately, but linked through the conductance term. Vector-valued diffusion it is also useful for processing registered data from different devices or for denoising higher-order geometric or statistical features from scalar-valued images [140149].

The output of anisotropic diffusion is an image or set of images that demonstrates reduced noise and texture but preserves, and can also enhance, edges. Such images are useful for a variety of processes including statistical classification, visualization, and geometric feature extraction. Previous work has shown [138] that anisotropic diffusion, over a wide range of conductance parameters, offers quantifiable advantages over linear filtering for edge detection in medical images.

Since the effectiveness of nonlinear diffusion was first demonstrated, numerous variations of this approach have surfaced in the literature [128]. These include alternatives for constructing dissimilarity measures [121], directional (i.e., tensor-valued) conductance terms [1356] and level set interpretations [143].

Gradient Anisotropic Diffusion

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

The itk::GradientAnisotropicDiffusionImageFilter implements an N-dimensional version of the classic Perona-Malik anisotropic diffusion equation for scalar-valued images [106].

The conductance term for this implementation is chosen as a function of the gradient magnitude of the image at each point, reducing the strength of diffusion at edge pixels.

       −(∥∇U(x)∥)2
C (x)= e    K
(8.8)

The numerical implementation of this equation is similar to that described in the Perona-Malik paper [106], but uses a more robust technique for gradient magnitude estimation and has been generalized to N-dimensions.

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

#include "itkGradientAnisotropicDiffusionImageFilter.h"

Types should be selected based on the pixel types required for the input and output images. The image types are defined using the pixel type and the dimension.

  typedef    float InputPixelType; 
  typedef    float OutputPixelType; 
 
  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The filter type is now instantiated using both the input image and the output image types. The filter object is created by the New() method.

  typedef itk::GradientAnisotropicDiffusionImageFilter< 
      InputImageType, OutputImageType>  FilterType; 
  FilterType::Pointer filter = FilterType::New();

The input image can be obtained from the output of another filter. Here, an image reader is used as source.

  filter->SetInput(reader->GetOutput());

This filter requires three parameters, the number of iterations to be performed, the time step and the conductance parameter used in the computation of the level set evolution. These parameters are set using the methods SetNumberOfIterations(), SetTimeStep() and SetConductanceParameter() respectively. The filter can be executed by invoking Update().

  filter->SetNumberOfIterations(numberOfIterations); 
  filter->SetTimeStep(timeStep); 
  filter->SetConductanceParameter(conductance); 
 
  filter->Update();

A typical value for the time step is 0.125. The number of iterations is typically set to 5; more iterations result in further smoothing and will increase the computing time linearly.


PIC PIC

Figure 8.18: Effect of the GradientAnisotropicDiffusionImageFilter.


Figure 8.18 illustrates the effect of this filter. In this example the filter was run with a time step of 0.125, and 5 iterations. The figure shows how homogeneous regions are smoothed and edges are preserved.

The following classes provide similar functionality:

Mean Shift filtering and clustering

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

This example demonstrates the use of the otb::MeanShiftSegmentationFilter class which implements filtering and clustering using the mean shift algorithm [29]. For a given pixel, the mean shift will build a set of neighboring pixels within a given spatial radius and a color range. The spatial and color center of this set is then computed and the algorithm iterates with this new spatial and color center. The Mean Shift can be used for edge-preserving smoothing, or for clustering.

We start by including the needed header file.

#include "otbMeanShiftSegmentationFilter.h"

We start by the classical typedefs needed for reading and writing the images.

  const unsigned int Dimension = 2; 
 
  typedef float                        PixelType; 
  typedef unsigned int                 LabelPixelType; 
  typedef itk::RGBPixel<unsigned char> ColorPixelType; 
 
  typedef otb::VectorImage<PixelType, Dimension> ImageType; 
  typedef otb::Image<LabelPixelType, Dimension>  LabelImageType; 
  typedef otb::Image<ColorPixelType, Dimension>  RGBImageType; 
 
  typedef otb::ImageFileReader<ImageType>               ReaderType; 
  typedef otb::ImageFileWriter<ImageType>      WriterType; 
  typedef otb::ImageFileWriter<LabelImageType> LabelWriterType; 
 
  typedef otb::MeanShiftSegmentationFilter<ImageType, LabelImageType, ImageType> FilterType;

We instantiate the filter, the reader, and 2 writers (for the labeled and clustered images).

  FilterType::Pointer filter = FilterType::New(); 
  ReaderType::Pointer reader = ReaderType::New(); 
  WriterType::Pointer writer1 = WriterType::New(); 
  LabelWriterType::Pointer writer2 = LabelWriterType::New();

We set the file names for the reader and the writers:

  reader->SetFileName(infname); 
  writer1->SetFileName(clusteredfname); 
  writer2->SetFileName(labeledfname);

We can now set the parameters for the filter. There are 3 main parameters: the spatial radius used for defining the neighborhood, the range radius used for defining the interval in the color space and the minimum size for the regions to be kept after clustering.

  filter->SetSpatialBandwidth(spatialRadius); 
  filter->SetRangeBandwidth(rangeRadius); 
  filter->SetMinRegionSize(minRegionSize);

Two another parameters can be set : the maximum iteration number, which defines maximum number of iteration until convergence. Algorithm iterative scheme will stop if convergence hasn’t been reached after the maximum number of iterations. Threshold parameter defines mean-shift vector convergence value. Algorithm iterative scheme will stop if mean-shift vector is below this threshold or if iteration number reached maximum number of iterations.

  filter->SetMaxIterationNumber(maxiter); 
  filter->SetThreshold(thres);

We can now plug the pipeline and run it.

  filter->SetInput(reader->GetOutput()); 
  writer1->SetInput(filter->GetClusteredOutput()); 
  writer2->SetInput(filter->GetLabelOutput()); 
 
  writer1->Update(); 
  writer2->Update();

Figure 8.19 shows the result of applying the mean shift to a Quickbird image.


PIC PIC PIC

Figure 8.19: From top to bottom and left to right: Original image, image filtered by mean shift after clustering , and labeled image.


8.7.3 Edge Preserving Speckle Reduction Filters

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

This example illustrates the use of the otb::LeeImageFilter. This filter belongs to the family of the edge-preserving smoothing filters which are usually used for speckle reduction in radar images. The Lee filter [84] aplies a linear regression which minimizes the mean-square error in the frame of a multiplicative speckle model.

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

#include "otbLeeImageFilter.h"

Then we must decide what pixel type to use for the image.

  typedef  unsigned char PixelType;

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

  typedef otb::Image<PixelType,  2> InputImageType; 
  typedef otb::Image<PixelType,  2> OutputImageType;

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

  typedef otb::LeeImageFilter<InputImageType, OutputImageType> FilterType;

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

  typedef otb::ImageFileReader<InputImageType> ReaderType;

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

  typedef otb::ImageFileWriter<OutputImageType> WriterType;

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

  ReaderType::Pointer reader = ReaderType::New(); 
  FilterType::Pointer filter = FilterType::New();

The image obtained with the reader is passed as input to the otb::LeeImageFilter.

  filter->SetInput(reader->GetOutput());

The method SetRadius() defines the size of the window to be used for the computation of the local statistics. The method SetNbLooks() sets the number of looks of the input image.

  FilterType::SizeType Radius; 
  Radius[0] = atoi(argv[3]); 
  Radius[1] = atoi(argv[3]); 
 
  filter->SetRadius(Radius); 
  filter->SetNbLooks(atoi(argv[4]));

Figure 8.20 shows the result of applying the Lee filter to a SAR image.


PIC PIC

Figure 8.20: Result of applying the otb::LeeImageFilter to a SAR image.


The following classes provide similar functionality:

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

This example illustrates the use of the otb::FrostImageFilter. This filter belongs to the family of the edge-preserving smoothing filters which are usually used for speckle reduction in radar images.

This filter uses a negative exponential convolution kernel. The output of the filter for pixel p is: s = pηpmpIp

where : mp = --KC2sexp2(−KC2sds,p)2--
∑p∈ηpKCsexp(−KCsds,p) and ds,p = ∘ ----------------
  (i− ip)2+(j− jp)2

Most of this example is similar to the previous one and only the differences will be highlighted.

First, we need to include the header:

#include "otbFrostImageFilter.h"

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

  typedef otb::FrostImageFilter<InputImageType, OutputImageType> FilterType;

The image obtained with the reader is passed as input to the otb::FrostImageFilter.

  filter->SetInput(reader->GetOutput());

The method SetRadius() defines the size of the window to be used for the computation of the local statistics. The method SetDeramp() sets the K coefficient.

  FilterType::SizeType Radius; 
  Radius[0] = atoi(argv[3]); 
  Radius[1] = atoi(argv[3]); 
 
  filter->SetRadius(Radius); 
  filter->SetDeramp(atof(argv[4]));

Figure 8.21 shows the result of applying the Frost filter to a SAR image.


PIC PIC

Figure 8.21: Result of applying the otb::FrostImageFilter to a SAR image.


The following classes provide similar functionality:

8.7.4 Edge preserving Markov Random Field

The Markov Random Field framework for OTB is more detailled in 19.2.6 (p. 895).

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

The Markov Random Field framework can be used to apply an edge preserving filtering, thus playing a role of restauration.

This example applies the otb::MarkovRandomFieldFilter for image restauration. The structure of the example is similar to the other MRF example. The original image is assumed to be coded in one byte, thus 256 states are possible for each pixel. The only other modifications reside in the energy function chosen for the fidelity and for the regularization.

For the regularization energy function, we choose an edge preserving function:

         2
Φ(u)= --u--
      1 +u2
(8.9)

and for the fidelity function, we choose a gaussian model.

The starting state of the Markov Random Field is given by the image itself as the final state should not be too far from it.

The first step toward the use of this filter is the inclusion of the proper header files:

#include "otbMRFEnergyEdgeFidelity.h" 
#include "otbMRFEnergyGaussian.h" 
#include "otbMRFOptimizerMetropolis.h" 
#include "otbMRFSamplerRandom.h"

We declare the usual types:

  const unsigned int Dimension = 2; 
 
  typedef double                                   InternalPixelType; 
  typedef unsigned char                            LabelledPixelType; 
  typedef otb::Image<InternalPixelType, Dimension> InputImageType; 
  typedef otb::Image<LabelledPixelType, Dimension> LabelledImageType;

We need to declare an additional reader for the initial state of the MRF. This reader has to be instantiated on the LabelledImageType.

  typedef otb::ImageFileReader<InputImageType>    ReaderType; 
  typedef otb::ImageFileReader<LabelledImageType> ReaderLabelledType; 
  typedef otb::ImageFileWriter<LabelledImageType> WriterType; 
 
  ReaderType::Pointer         reader = ReaderType::New(); 
  ReaderLabelledType::Pointer reader2 = ReaderLabelledType::New(); 
  WriterType::Pointer         writer = WriterType::New(); 
 
  const char  inputFilename  = argv[1]; 
  const char  labelledFilename  = argv[2]; 
  const char  outputFilename = argv[3]; 
 
  reader->SetFileName(inputFilename); 
  reader2->SetFileName(labelledFilename); 
  writer->SetFileName(outputFilename);

We declare all the necessary types for the MRF:

  typedef otb::MarkovRandomFieldFilter 
  <InputImageType, LabelledImageType> MarkovRandomFieldFilterType; 
 
  typedef otb::MRFSamplerRandom<InputImageType, LabelledImageType> SamplerType; 
 
  typedef otb::MRFOptimizerMetropolis OptimizerType;

The regularization and the fidelity energy are declared and instanciated:

  typedef otb::MRFEnergyEdgeFidelity 
  <LabelledImageType, LabelledImageType>  EnergyRegularizationType; 
  typedef otb::MRFEnergyGaussian 
  <InputImageType, LabelledImageType>  EnergyFidelityType;
  MarkovRandomFieldFilterType::Pointer markovFilter = 
    MarkovRandomFieldFilterType::New(); 
 
  EnergyRegularizationType::Pointer energyRegularization = 
    EnergyRegularizationType::New(); 
  EnergyFidelityType::Pointer energyFidelity = EnergyFidelityType::New(); 
 
  OptimizerType::Pointer optimizer = OptimizerType::New(); 
  SamplerType::Pointer   sampler = SamplerType::New();

The number of possible states for each pixel is 256 as the image is assumed to be coded on one byte and we pass the parameters to the markovFilter.

  unsigned int nClass = 256; 
 
  optimizer->SetSingleParameter(atof(argv[6])); 
  markovFilter->SetNumberOfClasses(nClass); 
  markovFilter->SetMaximumNumberOfIterations(atoi(argv[5])); 
  markovFilter->SetErrorTolerance(0.0); 
  markovFilter->SetLambda(atof(argv[4])); 
  markovFilter->SetNeighborhoodRadius(1); 
 
  markovFilter->SetEnergyRegularization(energyRegularization); 
  markovFilter->SetEnergyFidelity(energyFidelity); 
  markovFilter->SetOptimizer(optimizer); 
  markovFilter->SetSampler(sampler);

The original state of the MRF filter is passed through the SetTrainingInput() method:

  markovFilter->SetTrainingInput(reader2->GetOutput());

And we plug the pipeline:

  markovFilter->SetInput(reader->GetOutput()); 
 
  typedef itk::RescaleIntensityImageFilter 
  <LabelledImageType, LabelledImageType> RescaleType; 
  RescaleType::Pointer rescaleFilter = RescaleType::New(); 
  rescaleFilter->SetOutputMinimum(0); 
  rescaleFilter->SetOutputMaximum(255); 
 
  rescaleFilter->SetInput(markovFilter->GetOutput()); 
 
  writer->SetInput(rescaleFilter->GetOutput()); 
 
  writer->Update();

Figure 8.22 shows the output of the Markov Random Field restauration.


PIC PIC

Figure 8.22: Result of applying the otb::MarkovRandomFieldFilter to an extract from a PAN Quickbird image for restauration. From left to right : original image, restaured image with edge preservation.


8.8 Distance Map

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

This example illustrates the use of the itk::DanielssonDistanceMapImageFilter. This filter generates a distance map from the input image using the algorithm developed by Danielsson [32]. As secondary outputs, a Voronoi partition of the input elements is produced, as well as a vector image with the components of the distance vector to the closest point. The input to the map is assumed to be a set of points on the input image. Each point/pixel is considered to be a separate entity even if they share the same gray level value.

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

#include "itkConnectedComponentImageFilter.h" 
#include "itkDanielssonDistanceMapImageFilter.h"

Then we must decide what pixel types to use for the input and output images. Since the output will contain distances measured in pixels, the pixel type should be able to represent at least the width of the image, or said in N D terms, the maximum extension along all the dimensions. The input and output image types are now defined using their respective pixel type and dimension.

  typedef  unsigned char                 InputPixelType; 
  typedef  unsigned short                OutputPixelType; 
  typedef otb::Image<InputPixelType,  2> InputImageType; 
  typedef otb::Image<OutputPixelType, 2> OutputImageType;

The filter type can be instantiated using the input and output image types defined above. A filter object is created with the New() method.

  typedef itk::ConnectedComponentImageFilter< 
    InputImageType, InputImageType> ConnectedType; 
  ConnectedType::Pointer connectedComponents = ConnectedType::New(); 
 
  typedef itk::DanielssonDistanceMapImageFilter< 
    InputImageType, OutputImageType, OutputImageType> FilterType; 
  FilterType::Pointer filter = FilterType::New();

The input to the filter is taken from a reader and its output is passed to a itk::RescaleIntensityImageFilter and then to a writer.

  connectedComponents->SetInput(reader->GetOutput()); 
  filter->SetInput(connectedComponents->GetOutput()); 
  scaler->SetInput(filter->GetOutput()); 
  writer->SetInput(scaler->GetOutput());

The type of input image has to be specified. In this case, a binary image is selected.

  filter->InputIsBinaryOff();


PIC PIC PIC

Figure 8.23: DanielssonDistanceMapImageFilter output. Set of pixels, distance map and Voronoi partition.


Figure 8.23 illustrates the effect of this filter on a binary image with a set of points. The input image is shown at left, the distance map at the center and the Voronoi partition at right. This filter computes distance maps in N-dimensions and is therefore capable of producing N D Voronoi partitions.

The Voronoi map is obtained with the GetVoronoiMap() method. In the lines below we connect this output to the intensity rescaler and save the result in a file.

  scaler->SetInput(filter->GetVoronoiMap()); 
  writer->SetFileName(voronoiMapFileName); 
  writer->Update();

Execution of the writer is triggered by the invocation of the Update() method. Since this method can potentially throw exceptions it must be placed in a try/catch block.