HarrisExample.cxx 6.31 KB
Newer Older
1
/*
Julien Michel's avatar
Julien Michel committed
2
 * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
Jordi Inglada's avatar
Jordi Inglada committed
20

Patrick Imbo's avatar
Patrick Imbo committed
21

22
#include "otbMacro.h"
23
#include "otbImage.h"
Jordi Inglada's avatar
Jordi Inglada committed
24 25 26 27

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

28 29 30 31
/* Example usage:
./HarrisExample Input/ROISpot5.png Output/ROISpot5Harris.png 1.5 2 0.1
*/

Jordi Inglada's avatar
Jordi Inglada committed
32

33
// This example illustrates the use of the \doxygen{otb}{HarrisImageFilter}.
Jordi Inglada's avatar
Jordi Inglada committed
34
//
35
// The first step required to use this filter is to include its header file.
Jordi Inglada's avatar
Jordi Inglada committed
36

37
#include "otbHarrisImageToPointSetFilter.h"
Jordi Inglada's avatar
Jordi Inglada committed
38 39
#include "itkRescaleIntensityImageFilter.h"

40
int main(int argc, char* argv[])
Jordi Inglada's avatar
Jordi Inglada committed
41
{
OTB Bot's avatar
OTB Bot committed
42
  if (argc != 6)
43
  {
Jordi Inglada's avatar
Jordi Inglada committed
44
    std::cerr << "Usage: " << argv[0] << " inputImageFile ";
45
    std::cerr << " outputHarrisImageFile sigmaD sigmaI alpha" << std::endl;
Jordi Inglada's avatar
Jordi Inglada committed
46
    return EXIT_FAILURE;
47
  }
48

49 50
  const char* inputFilename  = argv[1];
  const char* outputFilename = argv[2];
51

52 53 54
  double SigmaD((double)::atof(argv[3]));
  double SigmaI((double)::atof(argv[4]));
  double Alpha((double)::atof(argv[5]));
55

56 57 58
  using InputPixelType         = float;
  const unsigned int Dimension = 2;
  using OutputPixelType        = unsigned char;
59

60 61
  using InputImageType  = otb::Image<InputPixelType, Dimension>;
  using OutputImageType = otb::Image<OutputPixelType, Dimension>;
62

63
  using ReaderType = otb::ImageFileReader<InputImageType>;
64 65 66 67 68

  //  The \doxygen{otb}{HarrisImageFilter} is templated over the
  //  input and output image types, so we start by
  //  defining:

69 70
  using HarrisFilterType = otb::HarrisImageFilter<InputImageType, InputImageType>;
  using RescalerType     = itk::RescaleIntensityImageFilter<InputImageType, OutputImageType>;
71

72
  using WriterType = otb::ImageFileWriter<OutputImageType>;
73

OTB Bot's avatar
OTB Bot committed
74 75
  ReaderType::Pointer       reader   = ReaderType::New();
  WriterType::Pointer       writer   = WriterType::New();
76
  HarrisFilterType::Pointer harris   = HarrisFilterType::New();
OTB Bot's avatar
OTB Bot committed
77
  RescalerType::Pointer     rescaler = RescalerType::New();
78

OTB Bot's avatar
OTB Bot committed
79 80
  reader->SetFileName(inputFilename);
  writer->SetFileName(outputFilename);
81

OTB Bot's avatar
OTB Bot committed
82
  harris->SetInput(reader->GetOutput());
83 84 85 86 87 88 89 90 91 92 93 94

  // The \doxygen{otb}{HarrisImageFilter} needs some parameters to
  // operate. The derivative computation is performed by a
  // convolution with the derivative of a Gaussian kernel of
  // variance $\sigma_D$ (derivation scale) and
  // the smoothing of the image is performed by convolving with a
  // Gaussian kernel of variance $\sigma_I$ (integration
  // scale). This allows the computation of the following matrix:
  // \begin{equation}
  // \mu(\mathbf{x},\sigma_I,\sigma_D) = \sigma_D^2 g(\sigma_I)\star
  // \left[\begin{array}{cc} L_x^2(\mathbf{x},\sigma_D) &
  // L_xL_y(\mathbf{x},\sigma_D)\\ L_xL_y(\mathbf{x},\sigma_D)&
95 96
  // L_y^2(\mathbf{x},\sigma_D) \end{array}\right]
  // \end{equation}
97 98
  // The output of the detector is $$det(\mu) - \alpha trace^2(\mu).$$

OTB Bot's avatar
OTB Bot committed
99 100 101
  harris->SetSigmaD(SigmaD);
  harris->SetSigmaI(SigmaI);
  harris->SetAlpha(Alpha);
102

OTB Bot's avatar
OTB Bot committed
103 104
  rescaler->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  rescaler->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
105

OTB Bot's avatar
OTB Bot committed
106 107
  rescaler->SetInput(harris->GetOutput());
  writer->SetInput(rescaler->GetOutput());
108 109
  writer->Update();

Jordi Inglada's avatar
Jordi Inglada committed
110 111 112 113 114 115
  // Figure~\ref{fig:Harris} shows the result of applying the interest
  // point detector to a small patch extracted from a Spot 5 image.
  // \begin{figure}
  // \center
  // \includegraphics[width=0.25\textwidth]{ROISpot5.eps}
  // \includegraphics[width=0.25\textwidth]{ROISpot5Harris.eps}
Jordi Inglada's avatar
Jordi Inglada committed
116
  // \itkcaption[Harris Filter Application]{Result of applying the
117
  // \doxygen{otb}{HarrisImageFilter} to a Spot 5 image.}
Jordi Inglada's avatar
Jordi Inglada committed
118 119 120
  // \label{fig:Harris}
  // \end{figure}
  //
121
  // The output of the \doxygen{otb}{HarrisImageFilter} is an image
OTB Bot's avatar
OTB Bot committed
122
  // where, for each pixel, we obtain the intensity of the
123 124 125
  // detection. Often, the user may want to get access to the set of
  // points for which the output of the detector is higher than a
  // given threshold. This can be obtained by using the
126
  // \doxygen{otb}{HarrisImageToPointSetFilter}. This filter is only
127
  // templated over the input image type, the output being a
128
  // \doxygen{itk}{PointSet} with pixel type equal to the image pixel type.
Jordi Inglada's avatar
Jordi Inglada committed
129

130
  using FunctionType = otb::HarrisImageToPointSetFilter<InputImageType>;
131 132

  //  We declare now the filter and a pointer to the output point set.
133
  using OutputPointSetType = FunctionType::OutputPointSetType;
134

135 136
  FunctionType::Pointer       harrisPoints = FunctionType::New();
  OutputPointSetType::Pointer pointSet     = OutputPointSetType::New();
137 138 139 140 141

  //  The \doxygen{otb}{HarrisImageToPointSetFilter} takes the same
  // parameters as the \doxygen{otb}{HarrisImageFilter} and an
  // additional parameter : the threshold for the point selection.

OTB Bot's avatar
OTB Bot committed
142 143 144 145 146
  harrisPoints->SetInput(0, reader->GetOutput());
  harrisPoints->SetSigmaD(SigmaD);
  harrisPoints->SetSigmaI(SigmaI);
  harrisPoints->SetAlpha(Alpha);
  harrisPoints->SetLowerThreshold(10);
147 148 149 150 151 152 153 154 155 156 157
  pointSet = harrisPoints->GetOutput();

  harrisPoints->Update();

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

158 159 160 161
  using ContainerType            = OutputPointSetType::PointsContainer;
  ContainerType* pointsContainer = pointSet->GetPoints();
  using IteratorType             = ContainerType::Iterator;
  IteratorType itList            = pointsContainer->Begin();
162 163 164

  //  And we get the points coordinates

OTB Bot's avatar
OTB Bot committed
165
  while (itList != pointsContainer->End())
166
  {
167 168
    using OutputPointType       = OutputPointSetType::PointType;
    OutputPointType pCoordinate = (itList.Value());
169
    otbLogMacro(Debug, << pCoordinate);
170
    ++itList;
171
  }
172

Jordi Inglada's avatar
Jordi Inglada committed
173 174
  return EXIT_SUCCESS;
}