Chapter 4
Building Simple Applications with OTB

Well, that’s it, you’ve just downloaded and installed OTB, lured by the promise that you will be able to do everything with it. That’s true, you will be able to do everything but - there is always a but - some effort is required.

OTB uses the very powerful systems of generic programming, many classes are already available, some powerful tools are defined to help you with recurrent tasks, but it is not an easy world to enter.

These tutorials are designed to help you enter this world and grasp the logic behind OTB. Each of these tutorials should not take more than 10 minutes (typing included) and each is designed to highlight a specific point. You may not be concerned by the latest tutorials but it is strongly advised to go through the first few which cover the basics you’ll use almost everywhere.

4.1 Hello world

4.1.1 Linux and Mac OS X

Let’s start by the typical Hello world program. We are going to compile this C++ program linking to your new OTB.

First, create a new folder to put your new programs (all the examples from this tutorial) in and go into this folder.

Since all programs using OTB are handled using the CMake system, we need to create a CMakeLists.txt that will be used by CMake to compile our program. An example of this file can be found in the OTB/Examples/Tutorials directory. The CMakeLists.txt will be very similar between your projects.

Open the CMakeLists.txt file and write in the few lines:

  PROJECT(Tutorials)
  
  cmake_minimum_required(VERSION 2.6)
  
  FIND_PACKAGE(OTB)
  IF(OTB_FOUND)
    INCLUDE(${OTB_USE_FILE})
  ELSE(OTB_FOUND)
    MESSAGE(FATAL_ERROR
        "Cannot build OTB project without OTB.  Please set OTB_DIR.")
  ENDIF(OTB_FOUND)
  
  ADD_EXECUTABLE(HelloWorldOTB HelloWorldOTB.cxx )
  TARGET_LINK_LIBRARIES(HelloWorldOTB ${OTB_LIBRARIES})

The first line defines the name of your project as it appears in Visual Studio (it will have no effect under UNIX or Linux). The second line loads a CMake file with a predefined strategy for finding OTB 1. If the strategy for finding OTB fails, CMake will prompt you for the directory where OTB is installed in your system. In that case you will write this information in the OTB_DIR variable. The line INCLUDE(${USE_OTB_FILE}) loads the UseOTB.cmake file to set all the configuration information from OTB.

The line ADD_EXECUTABLE defines as its first argument the name of the executable that will be produced as result of this project. The remaining arguments of ADD_EXECUTABLE are the names of the source files to be compiled and linked. Finally, the TARGET_LINK_LIBRARIES line specifies which OTB libraries will be linked against this project.

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

The following code is an implementation of a small OTB program. It tests including header files and linking with OTB libraries.

  #include "otbImage.h"
  #include <iostream>
  
  int main(int itkNotUsed(argc), char  itkNotUsed(argv)[])
  {
    typedef otb::Image<unsigned short, 2> ImageType;
  
    ImageType::Pointer image = ImageType::New();
  
    std::cout << "OTB Hello World !" << std::endl;
  
    return EXIT_SUCCESS;
  }

This code instantiates an image whose pixels are represented with type unsigned short. The image is then created and assigned to a itk::SmartPointer . Later in the text we will discuss SmartPointers in detail, for now think of it as a handle on an instance of an object (see section 3.2.4 for more information).

Once the file is written, run ccmake on the current directory (that is ccmake ./ under Linux/Unix). If OTB is on a non standard place, you will have to tell CMake where it is. Once your done with CMake (you shouldn’t have to do it anymore) run make.

You finally have your program. When you run it, you will have the OTB Hello World ! printed.

Ok, well done! You’ve just compiled and executed your first OTB program. Actually, using OTB for that is not very useful, and we doubt that you downloaded OTB only to do that. It’s time to move on to a more advanced level.

4.1.2 Windows

Create a directory (with write access) where to store your work (for example at C:\path\to\MyFirstCode). Organize your repository as it :

Follow the following steps:

  1. Create a CMakeLists.txt into the src repository with the following lines:

      project(MyFirstProcessing)
      
      cmake_minimum_required(VERSION 2.8)
      
      find_package(OTB REQUIRED)
      include(${OTB_USE_FILE})
      
      add_executable(MyFirstProcessing MyFirstProcessing.cxx )
      
      target_link_libraries(MyFirstProcessing ${OTB_LIBRARIES} )
  2. Create a MyFirstProcessing.cxx into the src repository with the following lines:

      #include "otbImage.h"
      #include "otbVectorImage.h"
      #include "otbImageFileReader.h"
      #include "otbImageFileWriter.h"
      #include "otbMultiToMonoChannelExtractROI.h"
      
      int main(int argc, char argv[])
      {
        if (argc < 3)
        {
          std::cerr << "Usage: " << std::endl;
          std::cerr << argv[0] << "  inputImageFile  outputImageFile" << std::endl;
          return EXIT_FAILURE;
        }
      
        typedef unsigned short PixelType;
        typedef otb::Image <PixelType, 2> ImageType;
        typedef otb::VectorImage <PixelType, 2> VectorImageType;
        typedef otb::MultiToMonoChannelExtractROI <PixelType, PixelType> FilterType;
        typedef otb::ImageFileReader<VectorImageType> ReaderType;
        typedef otb::ImageFileWriter<ImageType> WriterType;
      
        FilterType::Pointer filter = FilterType::New();
        ReaderType::Pointer reader = ReaderType::New();
        WriterType::Pointer writer = WriterType::New();
      
        reader->SetFileName(argv[1]);
        filter->SetInput(reader->GetOutput());
        writer->SetFileName(argv[2]);
        writer->SetInput(filter->GetOutput());
      
        return EXIT_SUCCESS;
      }
  3. create a file named BuildMyFirstProcessing.bat into the MyFirstCode directory with the following lines:
      @echo off
      
      set /A ARGS_COUNT=0
      for %%A in (%⋆) do set /A ARGS_COUNT+=1
      if %ARGS_COUNT% NEQ 3 (goto :Usage)
      
      if NOT DEFINED OSGEO4W_ROOT (goto :NoOSGEO4W)
      
      set src_dir=%1
      set build_dir=%2
      set otb_install_dir=%3
      set current_dir=%CD%
      
      cd %build_dir%
      
      cmake %src_dir% ^
            -DCMAKE_INCLUDE_PATH:PATH="%OSGEO4W_ROOT%\include" ^
            -DCMAKE_LIBRARY_PATH:PATH="%OSGEO4W_ROOT%\lib" ^
            -DOTB_DIR:PATH=%otb_install_dir% ^
            -DCMAKE_CONFIGURATION_TYPES:STRING=Release
      
      cmake --build . --target INSTALL --config Release
      
      cd %current_dir%
      
      goto :END
      
      :Usage
      echo You need to provide 3 arguments to the script:
      echo   1. path to the source directory
      echo   2. path to the build directory
      echo   3. path to the installation directory
      GOTO :END
      
      :NoOSGEO4W
      echo You need to run this script from an OSGeo4W shell
      GOTO :END
      
      :END
  4. into a OSGEo4W shell, run the configure.bat with the right arguments: full path to your src directory, full path to your build directory, full path to the place where find OTBConfig.cmake file (should be C:\path\to\MyOTBDir\install\lib\otb).
  5. into the OSGeo4W shell, open the MyFirstProcessing.sln
  6. build the solution
  7. into the OSGeo4W shell, go to the bin\Release directory and run MyFirstProcessing.exe. You can try for example with the otb_logo.tif file which can be found into the OTB source.

4.2 Pipeline basics: read and write

OTB is designed to read images, process them and write them to disk or view the result. In this tutorial, we are going to see how to read and write images and the basics of the pipeline system.

First, let’s add the following lines at the end of the CMakeLists.txt file:

  ADD_EXECUTABLE(Pipeline Pipeline.cxx )
  TARGET_LINK_LIBRARIES(Pipeline ${OTB_LIBRARIES})

Now, create a Pipeline.cxx file.

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

Start by including some necessary headers and with the usual main declaration:

  #include "otbImage.h"
  #include "otbImageFileReader.h"
  #include "otbImageFileWriter.h"
  
  int main(int argc, char  argv[])
  {
    if (argc != 3)
      {
      std::cerr << "Usage: "
          << argv[0]
          << " <input_filename> <output_filename>"
          << std::endl;
      }

Declare the image as an otb::Image , the pixel type is declared as an unsigned char (one byte) and the image is specified as having two dimensions.

    typedef otb::Image<unsigned char, 2> ImageType;

To read the image, we need an otb::ImageFileReader which is templated with the image type.

    typedef otb::ImageFileReader<ImageType> ReaderType;
    ReaderType::Pointer reader = ReaderType::New();

Then, we need an otb::ImageFileWriter also templated with the image type.

    typedef otb::ImageFileWriter<ImageType> WriterType;
    WriterType::Pointer writer = WriterType::New();

The filenames are passed as arguments to the program. We keep it simple for now and we don’t check their validity.

    reader->SetFileName(argv[1]);
    writer->SetFileName(argv[2]);

Now that we have all the elements, we connect the pipeline, pluging the output of the reader to the input of the writer.

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

And finally, we trigger the pipeline execution calling the Update() method on the last element of the pipeline. The last element will make sure to update all previous elements in the pipeline.

    writer->Update();
  
    return EXIT_SUCCESS;
  }

Once this file is written you just have to run make. The ccmake call is not required anymore.

Get one image from the OTB-Data/Examples directory from the OTB-Data repository. You can get it either by cloning the OTB data repository (git clone https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb-data.git), but that might be quite long as this also gets the data to run the tests. Alternatively, you can get it from http://www.orfeo-_toolbox.org/packages/OTB-_Data-_Examples.tgz. Take for example get QB_Suburb.png.

Now, run your new program as Pipeline QB_Suburb.png output.png. You obtain the file output.png which is the same image as QB_Suburb.png. When you triggered the Update() method, OTB opened the original image and wrote it back under another name.

Well…that’s nice but a bit complicated for a copy program!

Wait a minute! We didn’t specify the file format anywhere! Let’s try Pipeline QB_Suburb.png output.jpg. And voila! The output image is a jpeg file.

That’s starting to be a bit more interesting: this is not just a program to copy image files, but also to convert between image formats.

You have just experienced the pipeline structure which executes the filters only when needed and the automatic image format detection.

Now it’s time to do some processing in between.

4.3 Filtering pipeline

We are now going to insert a simple filter to do some processing between the reader and the writer.

Let’s first add the 2 following lines to the CMakeLists.txt file:

  ADD_EXECUTABLE(FilteringPipeline FilteringPipeline.cxx )
  TARGET_LINK_LIBRARIES(FilteringPipeline ${OTB_LIBRARIES})

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

We are going to use the itk::GradientMagnitudeImageFilter to compute the gradient of the image. The beginning of the file is similar to the Pipeline.cxx.

We include the required headers, without forgetting to add the header for the itk::GradientMagnitudeImageFilter .

  #include "otbImage.h"
  #include "otbImageFileReader.h"
  #include "otbImageFileWriter.h"
  #include "itkGradientMagnitudeImageFilter.h"
  
  int main(int argc, char  argv[])
  {
    if (argc != 3)
      {
      std::cerr << "Usage: "
          << argv[0]
          << " <input_filename> <output_filename>"
          << std::endl;
      }

We declare the image type, the reader and the writer as before:

    typedef otb::Image<unsigned char, 2> ImageType;
  
    typedef otb::ImageFileReader<ImageType> ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
  
    typedef otb::ImageFileWriter<ImageType> WriterType;
    WriterType::Pointer writer = WriterType::New();
  
    reader->SetFileName(argv[1]);
    writer->SetFileName(argv[2]);

Now we have to declare the filter. It is templated with the input image type and the output image type like many filters in OTB. Here we are using the same type for the input and the output images:

    typedef itk::GradientMagnitudeImageFilter
    <ImageType, ImageType> FilterType;
    FilterType::Pointer filter = FilterType::New();

Let’s plug the pipeline:

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

And finally, we trigger the pipeline execution calling the Update() method on the writer

    writer->Update();
  
    return EXIT_SUCCESS;
  }

Compile with make and execute as FilteringPipeline QB_Suburb.png output.png.

You have the filtered version of your image in the output.png file.

Now, you can practice a bit and try to replace the filter by one of the 150+ filters which inherit from the itk::ImageToImageFilter class. You will definitely find some useful filters here!

4.4 Handling types: scaling output

If you tried some other filter in the previous example, you may have noticed that in some cases, it does not make sense to save the output directly as an integer. This is the case if you tried the itk::CannyEdgeDetectionImageFilter . If you tried to use it directly in the previous example, you will have some warning about converting to unsigned char from double.

The output of the Canny edge detection is a floating point number. A simple solution would be to used double as the pixel type. Unfortunately, most image formats use integer typed and you should convert the result to an integer image if you still want to visualize your images with your usual viewer (we will see in a tutorial later how you can avoid that using the built-in viewer).

To realize this conversion, we will use the itk::RescaleIntensityImageFilter .

Add the two lines to the CMakeLists.txt file:

  ADD_EXECUTABLE(ScalingPipeline ScalingPipeline.cxx )
  TARGET_LINK_LIBRARIES(ScalingPipeline ${OTB_LIBRARIES})

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

This example illustrates the use of the itk::RescaleIntensityImageFilter to convert the result for proper display.

We include the required header including the header for the itk::CannyEdgeDetectionImageFilter and the itk::RescaleIntensityImageFilter .

  #include "otbImage.h"
  #include "otbImageFileReader.h"
  #include "otbImageFileWriter.h"
  #include "itkUnaryFunctorImageFilter.h"
  #include "itkCannyEdgeDetectionImageFilter.h"
  #include "itkRescaleIntensityImageFilter.h"
  
  int main(int argc, char  argv[])
  {
    if (argc != 3)
      {
      std::cerr << "Usage: "
          << argv[0]
          << " <input_filename> <output_filename>"
          << std::endl;
      }

We need to declare two different image types, one for the internal processing and one to output the results:

    typedef double                   PixelType;
    typedef otb::Image<PixelType, 2> ImageType;
  
    typedef unsigned char                  OutputPixelType;
    typedef otb::Image<OutputPixelType, 2> OutputImageType;

We declare the reader with the image template using the pixel type double. It is worth noticing that this instantiation does not imply anything about the type of the input image. The original image can be anything, the reader will just convert the result to double.

The writer is templated with the unsigned char image to be able to save the result on one byte images (like png for example).

    typedef otb::ImageFileReader<ImageType> ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
  
    typedef otb::ImageFileWriter<OutputImageType> WriterType;
    WriterType::Pointer writer = WriterType::New();
  
    reader->SetFileName(argv[1]);
    writer->SetFileName(argv[2]);

Now we are declaring the edge detection filter which is going to work with double input and output.

    typedef itk::CannyEdgeDetectionImageFilter
    <ImageType, ImageType> FilterType;
    FilterType::Pointer filter = FilterType::New();

Here comes the interesting part: we declare the itk::RescaleIntensityImageFilter . The input image type is the output type of the edge detection filter. The output type is the same as the input type of the writer.

Desired minimum and maximum values for the output are specified by the methods SetOutputMinimum() and SetOutputMaximum().

This filter will actually rescale all the pixels of the image but also cast the type of these pixels.

    typedef itk::RescaleIntensityImageFilter
    <ImageType, OutputImageType> RescalerType;
    RescalerType::Pointer rescaler = RescalerType::New();
  
    rescaler->SetOutputMinimum(0);
    rescaler->SetOutputMaximum(255);

Let’s plug the pipeline:

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

And finally, we trigger the pipeline execution calling the Update() method on the writer

    writer->Update();
  
    return EXIT_SUCCESS;
  }

As you should be getting used to it by now, compile with make and execute as ScalingPipeline QB_Suburb.png output.png.

You have the filtered version of your image in the output.png file.

4.5 Working with multispectral or color images

So far, as you may have noticed, we have been working with grey level images, i.e. with only one spectral band. If you tried to process a color image with some of the previous examples you have probably obtained a deceiving grey result.

Often, satellite images combine several spectral band to help the identification of materials: this is called multispectral imagery. In this tutorial, we are going to explore some of the mechanisms used by OTB to process multispectral images.

Add the following lines in the CMakeLists.txt file:

  ADD_EXECUTABLE(Multispectral Multispectral.cxx )
  TARGET_LINK_LIBRARIES(Multispectral ${OTB_LIBRARIES})

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

First, we are going to use otb::VectorImage instead of the now traditional otb::Image . So we include the required header:

  

We also include some other header which will be useful later. Note that we are still using the otb::Image in this example for some of the output.

  #include "otbImageFileReader.h"
  #include "otbImageFileWriter.h"
  #include "otbMultiToMonoChannelExtractROI.h"
  #include "itkShiftScaleImageFilter.h"
  #include "otbPerBandVectorImageFilter.h"
  
  int main(int argc, char  argv[])
  {
    if (argc != 4)
      {
      std::cerr << "Usage: "
          << argv[0]
          << " <input_filename> <output_extract> <output_shifted_scaled>"
          << std::endl;
      }

We want to read a multispectral image so we declare the image type and the reader. As we have done in the previous example we get the filename from the command line.

    typedef unsigned short int             PixelType;
    typedef otb::VectorImage<PixelType, 2> VectorImageType;
  
    typedef otb::ImageFileReader<VectorImageType> ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
  
    reader->SetFileName(argv[1]);

Sometime, you need to process only one spectral band of the image. To get only one of the spectral band we use the otb::MultiToMonoChannelExtractROI . The declaration is as usual:

    typedef otb::MultiToMonoChannelExtractROI<PixelType, PixelType>
    ExtractChannelType;
    ExtractChannelType::Pointer extractChannel = ExtractChannelType::New();

We need to pass the parameters to the filter for the extraction. This filter also allow extracting only a spatial subset of the image. However, we will extract the whole channel in this case.

To do that, we need to pass the desired region using the SetExtractionRegion() (method such as SetStartX, SetSizeX are also available). We get the region from the reader with the GetLargestPossibleRegion() method. Before doing that we need to read the metadata from the file: this is done by calling the UpdateOutputInformation() on the reader’s output. The difference with the Update() is that the pixel array is not allocated (yet !) and reduce the memory usage.

    reader->UpdateOutputInformation();
    extractChannel->SetExtractionRegion(
      reader->GetOutput()->GetLargestPossibleRegion());

We chose the channel number to extract (starting from 1) and we plug the pipeline.

    extractChannel->SetChannel(3);
    extractChannel->SetInput(reader->GetOutput());

To output this image, we need a writer. As the output of the otb::MultiToMonoChannelExtractROI is a otb::Image , we need to template the writer with this type.

    typedef otb::Image<PixelType, 2>                 ImageType;
    typedef otb::ImageFileWriter<ImageType> WriterType;
    WriterType::Pointer writer = WriterType::New();
  
    writer->SetFileName(argv[2]);
    writer->SetInput(extractChannel->GetOutput());
  
    writer->Update();

After this, we have a one band image that we can process with most OTB filters.

In some situation, you may want to apply the same process to all bands of the image. You don’t have to extract each band and process them separately. There is several situations:

Let’s see how this filter is working. We chose to apply the itk::ShiftScaleImageFilter to each of the spectral band. We start by declaring the filter on a normal otb::Image . Note that we don’t need to specify any input for this filter.

    typedef itk::ShiftScaleImageFilter<ImageType, ImageType> ShiftScaleType;
    ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
    shiftScale->SetScale(0.5);
    shiftScale->SetShift(10);

We declare the otb::PerBandVectorImageFilter which has three template: the input image type, the output image type and the filter type to apply to each band.

The filter is selected using the SetFilter() method and the input by the usual SetInput() method.

    typedef otb::PerBandVectorImageFilter
    <VectorImageType, VectorImageType, ShiftScaleType> VectorFilterType;
    VectorFilterType::Pointer vectorFilter = VectorFilterType::New();
    vectorFilter->SetFilter(shiftScale);
  
    vectorFilter->SetInput(reader->GetOutput());

Now, we just have to save the image using a writer templated over an otb::VectorImage :

    typedef otb::ImageFileWriter<VectorImageType> VectorWriterType;
    VectorWriterType::Pointer writerVector = VectorWriterType::New();
  
    writerVector->SetFileName(argv[3]);
    writerVector->SetInput(vectorFilter->GetOutput());
  
    writerVector->Update();
  
    return EXIT_SUCCESS;
  }

Compile with make and execute as ./Multispectral qb_RoadExtract.tif qb_blue.tif qb_shiftscale.tif.

4.6 Parsing command line arguments

Well, if you play with some other filters in the previous example, you probably noticed that in many cases, you need to set some parameters to the filters. Ideally, you want to set some of these parameters from the command line.

In OTB, there is a mechanism to help you parse the command line parameters. Let try it!

Add the following lines in the CMakeLists.txt file:

  ADD_EXECUTABLE(SmarterFilteringPipeline SmarterFilteringPipeline.cxx )
  TARGET_LINK_LIBRARIES(SmarterFilteringPipeline ${OTB_LIBRARIES})

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

We are going to use the otb::HarrisImageFilter to find the points of interest in one image.

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:

                   [                      ]
μ(x,σI,σD)= σ2g(σI)⋆   L2x(x,σD )  LxL2y(x,σD)
            D        LxL2y(x,σD)   L2y(x,σD )

The output of the detector is det(μ)-αtrace2(μ).

We want to set 3 parameters of this filter through the command line: σD (SigmaD), σI (SigmaI) and α (Alpha).

We are also going to do the things properly and catch the exceptions.

Let first add the two following headers:

  #include "itkMacro.h"
  #include "otbCommandLineArgumentParser.h"

The first one is to handle the exceptions, the second one to help us parse the command line.

We include the other required headers, without forgetting to add the header for the otb::HarrisImageFilter . Then we start the usual main function.

  #include "otbImage.h"
  #include "otbImageFileReader.h"
  #include "otbImageFileWriter.h"
  #include "itkUnaryFunctorImageFilter.h"
  #include "itkRescaleIntensityImageFilter.h"
  #include "otbHarrisImageFilter.h"
  
  int main(int argc, char  argv[])
  {

To handle the exceptions properly, we need to put all the instructions inside a try.

    try
      {

Now, we can declare the otb::CommandLineArgumentParser which is going to parse the command line, select the proper variables, handle the missing compulsory arguments and print an error message if necessary.

Let’s declare the parser:

      typedef otb::CommandLineArgumentParser ParserType;
      ParserType::Pointer parser = ParserType::New();

It’s now time to tell the parser what are the options we want. Special options are available for input and output images with the AddInputImage() and AddOutputImage() methods.

For the other options, we need to use the AddOption() method. This method allows us to specify

      parser->SetProgramDescription(
        "This program applies a Harris detector on the input image");
      parser->AddInputImage();
      parser->AddOutputImage();
      parser->AddOption("--SigmaD",
                        "Set the sigmaD parameter. Default is 1.0.",
                        "-d",
                        1,
                        false);
      parser->AddOption("--SigmaI",
                        "Set the sigmaI parameter. Default is 1.0.",
                        "-i",
                        1,
                        false);
      parser->AddOption("--Alpha",
                        "Set the alpha parameter. Default is 1.0.",
                        "-a",
                        1,
                        false);

Now that the parser has all this information, it can actually look at the command line to parse it. We have to do this within a try - catch loop to handle exceptions nicely.

      typedef otb::CommandLineArgumentParseResult ParserResultType;
      ParserResultType::Pointer parseResult = ParserResultType::New();
  
      try
        {
        parser->ParseCommandLine(argc, argv, parseResult);
        }
  
      catch (itk::ExceptionObject& err)
        {
        std::string descriptionException = err.GetDescription();
        if (descriptionException.find("ParseCommandLine(): Help Parser")
            != std::string::npos)
          {
          return EXIT_SUCCESS;
          }
        if (descriptionException.find("ParseCommandLine(): Version Parser")
            != std::string::npos)
          {
          return EXIT_SUCCESS;
          }
        return EXIT_FAILURE;
        }

Now, we can declare the image type, the reader and the writer as before:

      typedef double                   PixelType;
      typedef otb::Image<PixelType, 2> ImageType;
  
      typedef unsigned char                  OutputPixelType;
      typedef otb::Image<OutputPixelType, 2> OutputImageType;
  
      typedef otb::ImageFileReader<ImageType> ReaderType;
      ReaderType::Pointer reader = ReaderType::New();
  
      typedef otb::ImageFileWriter<OutputImageType> WriterType;
      WriterType::Pointer writer = WriterType::New();

We are getting the filenames for the input and the output images directly from the parser:

      reader->SetFileName(parseResult->GetInputImage().c_str());
      writer->SetFileName(parseResult->GetOutputImage().c_str());

Now we have to declare the filter. It is templated with the input image type and the output image type like many filters in OTB. Here we are using the same type for the input and the output images:

      typedef otb::HarrisImageFilter
      <ImageType, ImageType> FilterType;
      FilterType::Pointer filter = FilterType::New();

We set the filter parameters from the parser. The method IsOptionPresent() let us know if an optional option was provided in the command line.

      if (parseResult->IsOptionPresent("--SigmaD"))
        filter->SetSigmaD(
          parseResult->GetParameterDouble("--SigmaD"));
  
      if (parseResult->IsOptionPresent("--SigmaI"))
        filter->SetSigmaI(
          parseResult->GetParameterDouble("--SigmaI"));
  
      if (parseResult->IsOptionPresent("--Alpha"))
        filter->SetAlpha(
          parseResult->GetParameterDouble("--Alpha"));

We add the rescaler filter as before

      typedef itk::RescaleIntensityImageFilter
      <ImageType, OutputImageType> RescalerType;
      RescalerType::Pointer rescaler = RescalerType::New();
  
      rescaler->SetOutputMinimum(0);
      rescaler->SetOutputMaximum(255);

Let’s plug the pipeline:

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

We trigger the pipeline execution calling the Update() method on the writer

      writer->Update();
  
      }

Finally, we have to handle exceptions we may have raised before

    catch (itk::ExceptionObject& err)
      {
      std::cout << "Following otbException catch :" << std::endl;
      std::cout << err << std::endl;
      return EXIT_FAILURE;
      }
    catch (std::bad_alloc& err)
      {
      std::cout << "Exception bad_alloc : " << (char) err.what() << std::endl;
      return EXIT_FAILURE;
      }
    catch (...)
      {
      std::cout << "Unknown Exception found !" << std::endl;
      return EXIT_FAILURE;
      }
    return EXIT_SUCCESS;
  }

Compile with make as usual. The execution is a bit different now as we have an automatic parsing of the command line. First, try to execute as SmarterFilteringPipeline without any argument.

The usage message (automatically generated) appears:

'--InputImage' option is obligatory !!!  
 
 Usage : ./SmarterFilteringPipeline  
      [--help|-h]           :  Help  
      [--version|-v]        :  Version  
       --InputImage|-in     :  input image file name   (1 parameter)  
       --OutputImage|-out   :  output image file name   (1 parameter)  
      [--SigmaD|-d]         :  Set the sigmaD parameter of the Harris points of  
interest  algorithm. Default is 1.0.  (1 parameter)  
      [--SigmaI|-i]         :  Set the SigmaI parameter of the Harris points of  
interest  algorithm. Default is 1.0.  (1 parameter)  
      [--Alpha|-a]          :  Set the alpha parameter of the Harris points of  
interest  algorithm. Default is 1.0.  (1 parameter)

That looks a bit more professional: another user should be able to play with your program. As this is automatic, that’s a good way not to forget to document your programs.

So now you have a better idea of the command line options that are possible. Try SmarterFilteringPipeline -in QB_Suburb.png -out output.png for a basic version with the default values.

If you want a result that looks a bit better, you have to adjust the parameter with SmarterFilteringPipeline -in QB_Suburb.png -out output.png -d 1.5 -i 2 -a 0.1 for example.

4.7 Going from raw satellite images to useful products

Quite often, when you buy satellite images, you end up with several images. In the case of optical satellite, you often have a panchromatic spectral band with the highest spatial resolution and a multispectral product of the same area with a lower resolution. The resolution ratio is likely to be around 4.

To get the best of the image processing algorithms, you want to combine these data to produce a new image with the highest spatial resolution and several spectral band. This step is called fusion and you can find more details about it in 13. However, the fusion suppose that your two images represents exactly the same area. There are different solutions to process your data to reach this situation. Here we are going to use the metadata available with the images to produce an orthorectification as detailed in 11.

First you need to add the following lines in the CMakeLists.txt file:

  ADD_EXECUTABLE(OrthoFusion  OrthoFusion.cxx)
  TARGET_LINK_LIBRARIES(OrthoFusion ${OTB_LIBRARIES})

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

Start by including some necessary headers and with the usual main declaration. Apart from the classical header related to image input and output. We need the headers related to the fusion and the orthorectification. One header is also required to be able to process vector images (the XS one) with the orthorectification.

  #include "otbImageFileReader.h"
  #include "otbImageFileWriter.h"
  
  #include "otbOrthoRectificationFilter.h"
  #include "otbGenericMapProjection.h"
  
  #include "otbSimpleRcsPanSharpeningFusionImageFilter.h"
  #include "otbStandardFilterWatcher.h"
  
  int main(int argc, char argv[])
  {

We initialize ossim which is required for the orthorectification and we check that all parameters are provided. Basically, we need:

We check that all those parameters are provided.

    if (argc != 12)
      {
      std::cout << argv[0] << " <input_pan_filename> <input_xs_filename> ";
      std::cout << "<output_filename> <utm zone> <hemisphere N/S>  ";
      std::cout << "<x_ground_upper_left_corner> <y_ground_upper_left_corner> ";
      std::cout << "<x_Size> <y_Size> ";
      std::cout << "<x_groundSamplingDistance> ";
      std::cout << "<y_groundSamplingDistance "
                << "(negative since origin is upper left)>"
                << std::endl;
  
      return EXIT_FAILURE;
      }

We declare the different images, readers and writer:

    typedef otb::Image<unsigned int, 2>                    ImageType;
    typedef otb::VectorImage<unsigned int, 2>              VectorImageType;
    typedef otb::Image<double, 2>                          DoubleImageType;
    typedef otb::VectorImage<double, 2>                    DoubleVectorImageType;
    typedef otb::ImageFileReader<ImageType>                ReaderType;
    typedef otb::ImageFileReader<VectorImageType>          VectorReaderType;
    typedef otb::ImageFileWriter<VectorImageType> WriterType;
  
    ReaderType::Pointer       readerPAN = ReaderType::New();
    VectorReaderType::Pointer readerXS = VectorReaderType::New();
    WriterType::Pointer       writer = WriterType::New();
  
    readerPAN->SetFileName(argv[1]);
    readerXS->SetFileName(argv[2]);
    writer->SetFileName(argv[3]);

We declare the projection (here we chose the UTM projection, other choices are possible) and retrieve the parameters from the command line:

    typedef otb::GenericMapProjection<otb::TransformDirection::INVERSE> InverseProjectionType;
    InverseProjectionType::Pointer utmMapProjection = InverseProjectionType::New();
    utmMapProjection->SetWkt("Utm");
    utmMapProjection->SetParameter("Zone", argv[4]);
    utmMapProjection->SetParameter("Hemisphere", argv[5]);

We will need to pass several parameters to the orthorectification concerning the desired output region:

    ImageType::IndexType start;
    start[0] = 0;
    start[1] = 0;
  
    ImageType::SizeType size;
    size[0] = atoi(argv[8]);
    size[1] = atoi(argv[9]);
  
    ImageType::SpacingType spacing;
    spacing[0] = atof(argv[10]);
    spacing[1] = atof(argv[11]);
  
    ImageType::PointType origin;
    origin[0] = strtod(argv[6], ITK_NULLPTR);
    origin[1] = strtod(argv[7], ITK_NULLPTR);

We declare the orthorectification filter. And provide the different parameters:

    typedef otb::OrthoRectificationFilter<ImageType, DoubleImageType,
        InverseProjectionType>
    OrthoRectifFilterType;
  
    OrthoRectifFilterType::Pointer orthoRectifPAN =
      OrthoRectifFilterType::New();
    orthoRectifPAN->SetMapProjection(utmMapProjection);
  
    orthoRectifPAN->SetInput(readerPAN->GetOutput());
  
    orthoRectifPAN->SetOutputStartIndex(start);
    orthoRectifPAN->SetOutputSize(size);
    orthoRectifPAN->SetOutputSpacing(spacing);
    orthoRectifPAN->SetOutputOrigin(origin);

Now we are able to have the orthorectified area from the PAN image. We just have to follow a similar process for the XS image.

    typedef otb::OrthoRectificationFilter<VectorImageType,
        DoubleVectorImageType, InverseProjectionType>
    VectorOrthoRectifFilterType;
  
  
    VectorOrthoRectifFilterType::Pointer orthoRectifXS =
      VectorOrthoRectifFilterType::New();
  
    orthoRectifXS->SetMapProjection(utmMapProjection);
  
    orthoRectifXS->SetInput(readerXS->GetOutput());
  
    orthoRectifXS->SetOutputStartIndex(start);
    orthoRectifXS->SetOutputSize(size);
    orthoRectifXS->SetOutputSpacing(spacing);
    orthoRectifXS->SetOutputOrigin(origin);

It’s time to declare the fusion filter and set its inputs:

    typedef otb::SimpleRcsPanSharpeningFusionImageFilter
    <DoubleImageType, DoubleVectorImageType, VectorImageType> FusionFilterType;
    FusionFilterType::Pointer fusion = FusionFilterType::New();
    fusion->SetPanInput(orthoRectifPAN->GetOutput());
    fusion->SetXsInput(orthoRectifXS->GetOutput());

And we can plug it to the writer. To be able to process the images by tiles, we use the SetAutomaticTiledStreaming() method of the writer. We trigger the pipeline execution with the Update() method.

    writer->SetInput(fusion->GetOutput());
  
    writer->SetAutomaticTiledStreaming();
  
    otb::StandardFilterWatcher watcher(writer, "OrthoFusion");
  
    writer->Update();
  
    return EXIT_SUCCESS;
  
  }