HillShadingExample.cxx

Visualization of digital elevation models (DEM) is often more intuitive by simulating a lighting source and generating the corresponding shadows. This principle is called hill shading.

Using otb::HillShadingFilter and the DEM image generated using the otb::DEMToImageGenerator, you can easily obtain a representation of the DEM. Better yet, using the itk::ScalarToRGBColormapImageFilter combined with the ReliefColormapFunctor you can easily generate the classic elevation maps.

This example will focus on the shading itself.

Hill shading obtained from SRTM data (left) and combined with the color representation (right)

Example usage:

./HillShadingExample Output/HillShadingExample.png Output/HillShadingColorExample.png 6.5 45.5 500 500 0.002 -0.002 Input/DEM/srtm_directory

Example source code (HillShadingExample.cxx):

#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"

#include "otbDEMToImageGenerator.h"
#include "otbHillShadingFilter.h"

#include "itkScalarToRGBColormapImageFilter.h"
#include "otbReliefColormapFunctor.h"
#include "itkMultiplyImageFilter.h"
#include "itkShiftScaleImageFilter.h"
#include "otbWorldFile.h"

int main(int argc, char* argv[])
{
  if (argc < 10)
  {
    std::cout << argv[0] << " <output_filename> <output_color_filename> "
              << " <Longitude Output Origin point> <Latitude Output Origin point>"
              << " <X Output Size> <Y Output size>"
              << " <X Spacing> <Y Spacing> <DEM folder path>"
              << " [Projection Ref]" << std::endl;
    return EXIT_FAILURE;
  }

  using PixelType        = double;
  using UCharPixelType   = unsigned char;
  using RGBPixelType     = itk::RGBPixel<UCharPixelType>;
  using ImageType        = otb::Image<PixelType, 2>;
  using RGBImageType     = otb::Image<RGBPixelType, 2>;
  using ScalarImageType  = otb::Image<UCharPixelType, 2>;
  using WriterType       = otb::ImageFileWriter<RGBImageType>;
  using ScalarWriterType = otb::ImageFileWriter<ScalarImageType>;

  ScalarWriterType::Pointer writer = ScalarWriterType::New();
  writer->SetFileName(argv[1]);

  WriterType::Pointer writer2 = WriterType::New();
  writer2->SetFileName(argv[2]);

  using DEMToImageGeneratorType = otb::DEMToImageGenerator<ImageType>;

  DEMToImageGeneratorType::Pointer demToImage = DEMToImageGeneratorType::New();

  using SizeType    = DEMToImageGeneratorType::SizeType;
  using SpacingType = DEMToImageGeneratorType::SpacingType;
  using PointType   = DEMToImageGeneratorType::PointType;

  otb::DEMHandler::GetInstance().OpenDEMDirectory(argv[9]);

  PointType origin;
  origin[0] = ::atof(argv[3]);
  origin[1] = ::atof(argv[4]);

  demToImage->SetOutputOrigin(origin);

  SizeType size;
  size[0] = ::atoi(argv[5]);
  size[1] = ::atoi(argv[6]);

  demToImage->SetOutputSize(size);

  SpacingType spacing;
  spacing[0] = ::atof(argv[7]);
  spacing[1] = ::atof(argv[8]);

  demToImage->SetOutputSpacing(spacing);
  double res = 0;

  if (argc > 10)
  {
    demToImage->SetOutputProjectionRef(argv[10]);
    res = spacing[0];
  }
  else
  {
    // Compute the resolution (Vincenty formula)
    double lon1 = origin[0];
    double lon2 = origin[0] + size[0] * spacing[0];
    double lat1 = origin[1];
    double lat2 = origin[1] + size[1] * spacing[1];
    double R    = 6371; // km
    double d    = std::acos(std::sin(lat1) * std::sin(lat2) + std::cos(lat1) * std::cos(lat2) * std::cos(lon2 - lon1)) * R;
    res         = d / std::sqrt(2.0);
  }
  // After generating the DEM image as in the DEMToImageGenerator example, you can declare
  // the hill shading mechanism. The hill shading is implemented as a functor doing some
  // operations in its neighborhood. A convenient filter called \doxygen{otb}{HillShadingFilter}
  // is defined around this mechanism.

  using HillShadingFilterType                = otb::HillShadingFilter<ImageType, ImageType>;
  HillShadingFilterType::Pointer hillShading = HillShadingFilterType::New();
  hillShading->SetRadius(1);
  hillShading->SetInput(demToImage->GetOutput());

  hillShading->GetFunctor().SetXRes(res);
  hillShading->GetFunctor().SetYRes(res);

  using RescalerType             = itk::ShiftScaleImageFilter<ImageType, ScalarImageType>;
  RescalerType::Pointer rescaler = RescalerType::New();
  rescaler->SetScale(255.0);
  rescaler->SetInput(hillShading->GetOutput());

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

  using ColorMapFilterType                = itk::ScalarToRGBColormapImageFilter<ImageType, RGBImageType>;
  ColorMapFilterType::Pointer colormapper = ColorMapFilterType::New();
  colormapper->UseInputImageExtremaForScalingOff();

  using ColorMapFunctorType             = otb::Functor::ReliefColormapFunctor<PixelType, RGBPixelType>;
  ColorMapFunctorType::Pointer colormap = ColorMapFunctorType::New();
  colormap->SetMinimumInputValue(0);
  colormap->SetMaximumInputValue(4000);
  colormapper->SetColormap(colormap);

  colormapper->SetInput(demToImage->GetOutput());

  using MultiplyFilterType =
      itk::BinaryFunctorImageFilter<RGBImageType, ImageType, RGBImageType, otb::Functor::HillShadeModulationFunctor<RGBPixelType, PixelType, RGBPixelType>>;

  MultiplyFilterType::Pointer multiply = MultiplyFilterType::New();
  multiply->SetInput1(colormapper->GetOutput());
  multiply->SetInput2(hillShading->GetOutput());

  writer2->SetInput(multiply->GetOutput());

  writer->Update();
  writer2->Update();

  otb::WorldFile::Pointer worldFile = otb::WorldFile::New();
  worldFile->SetLonOrigin(origin[0]);
  worldFile->SetLatOrigin(origin[1]);
  worldFile->SetLonSpacing(spacing[0]);
  worldFile->SetLatSpacing(spacing[1]);

  worldFile->SetImageFilename(argv[1]);
  worldFile->Update();
  worldFile->SetImageFilename(argv[2]);
  worldFile->Update();
}